implicit meshes tested and ready in gui

AMDiS/src/AMDiS.h

@@ -72,6 +72,7 @@
 #include "ProblemVec.h"
 #include "ProblemInstat.h"
 #include "ProblemTimeInterface.h"
+#include "ProblemImplicit.h"
 #include "ProblemInterpolScal.h"
 #include "ProblemNonLin.h"
 #include "ProblemStatBase.h"

AMDiS/src/ProblemImplicit.cc

 #include "ProblemImplicit.h"
 
 namespace AMDiS {
-  namespace calcspace {
-    void readDofvec(std::istream& in, DOFVector<double>* vec, Mesh* mesh) 
-    {
-      long size = mesh->getNumberOfVertices();
-      std::string buf;
+  void ProblemImplicitBase::readDofVec(std::istream& in, DOFVector<double>* vec, 
+                                       Mesh* mesh) 
+  {
+    long size = mesh->getNumberOfVertices();
+
+    TEST_EXIT(vec->getSize()>=size)("dofvector smaller than vertex data vector\n");
+    std::string buf;
+    getline(in, buf);
+    while (!in.eof() && buf != "vertex values: solution") 
       getline(in, buf);
-      while (buf != "vertex values: solution") 
-	getline(in, buf);
-      for (long i = 0; i < size; i++) {
-        in >> buf;
-        (*vec)[i] = atof(buf.c_str());
-      }
+
+    TEST_EXIT(!in.eof())("no vertex data\n");
+
+    for (long i = 0; i < size ; ++i) {
+      in >> buf;
+      (*vec)[i] = atof(buf.c_str());
     }
+  }
 
+  void ProblemImplicitBase::readR(std::istream& inStream, double eps, 
+	                          Mesh* mesh, int implMesh , int comp) 
+  {
+    FUNCNAME("ProblemImplicitBase::readR()");
 
-    void readR(std::istream& inStream, double eps, 
-	       DOFVector<double>* r, DOFVector<double>* phi1, 
-	       DOFVector<double>* phi2, DOFVector<double>* levelSet, 
-	       Mesh* mesh) 
-    {
-      readDofvec(inStream, r, mesh);
-      std::vector<double>& vecR = r->getVector();
-      std::vector<double>& vecPhi1 = phi1->getVector();
-      std::vector<double>& vecPhi2 = phi2->getVector();
-      std::vector<double>& vecLevelSet = levelSet->getVector();
-      bool checkSize = 
-	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
-      TEST_EXIT(checkSize)("something went wrong\n");
-      for (int i = vecR.size() - 1; i >= 0; --i) {
-        vecPhi1[i] = Phi1(vecR[i], eps);
-        vecPhi2[i] = Phi2(vecR[i], eps);
-        vecLevelSet[i] = LevelSet(vecR[i]);
-      }
+    DOFVector< double >* r = getSignedDistance(implMesh, comp);
+    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
+    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
+    DOFVector< double >* levelSet = getLevelset(implMesh, comp);
+
+    TEST_EXIT( r != NULL )("no signed distance vector\n");
+    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
+    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
+    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");
+
+    bool checkSize = r->getSize() == phi1->getSize() && 
+	             r->getSize() == phi2->getSize();
+
+    TEST_EXIT(checkSize)("something went wrong\n");
+    
+
+    readDofVec(inStream, r, mesh);
+
+    for (int i = r->getSize() - 1; i >= 0; --i) {
+      (*phi1)[i] = Phi1((*r)[i], eps);
+      (*phi2)[i] = Phi2((*r)[i], eps);
+      (*levelSet)[i] = LevelSet((*r)[i]);
     }
+  }
 
 
-    void readPhi1(istream& inStream, double eps, 
-		  DOFVector<double>* r, DOFVector<double>* phi1, 
-		  DOFVector<double>* phi2, DOFVector<double>* levelSet, 
-		  Mesh* mesh) 
-    {
-      readDofvec(inStream, phi1, mesh);
-      std::vector<double>& vecR = r->getVector();
-      std::vector<double>& vecPhi1 = phi1->getVector();
-      std::vector<double>& vecPhi2 = phi2->getVector();
-      std::vector<double>& vecLevelSet = levelSet->getVector();
-      
-      bool checkSize = 
-	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
-      TEST_EXIT(checkSize)("something went wrong\n");
-      for (int i = vecR.size() - 1; i>= 0; --i) {
-        vecR[i] = Phi1ToR(vecPhi1[i], eps);
-        //vecPhi2[i] = Phi2(vecR[i], eps);
-        vecPhi2[i] = 1 - vecPhi1[i];
-        vecLevelSet[i] = LevelSet(vecR[i]);
-      }
+  void ProblemImplicitBase::readPhi1(istream& inStream, double eps, 
+                                     Mesh* mesh, int implMesh, int comp )
+  {
+    DOFVector< double >* r = getSignedDistance(implMesh, comp);
+    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
+    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
+    DOFVector< double >* levelSet = getLevelset(implMesh, comp);
+
+    TEST_EXIT( r != NULL )("no signed distance vector\n");
+    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
+    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
+    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");
+
+    bool checkSize = r->getSize() == phi1->getSize() && 
+	             r->getSize() == phi2->getSize();
+
+    TEST_EXIT(checkSize)("something went wrong\n");
+
+    readDofVec(inStream, phi1, mesh);
+
+    for (int i = r->getSize() - 1; i>= 0; --i) {
+      (*r)[i] = Phi1ToR((*phi1)[i], eps);
+      (*phi2)[i] = 1 - (*phi1)[i];
+      (*levelSet)[i] = LevelSet((*r)[i]);
     }
+  }
 
 
-    void readPhi2(istream& inStream, double eps, 
-		    DOFVector<double>* r, DOFVector<double>* phi1, 
-		    DOFVector<double>* phi2, DOFVector<double>* levelSet, 
-		    Mesh* mesh) 
-    {
-      readDofvec(inStream, phi2, mesh);
-      std::vector<double>& vecR = r->getVector();
-      std::vector<double>& vecPhi1 = phi1->getVector();
-      std::vector<double>& vecPhi2 = phi2->getVector();
-      std::vector<double>& vecLevelSet = levelSet->getVector();
-
-      bool checkSize = vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
-      TEST_EXIT(checkSize)("something went wrong\n");
-      for (int i = vecR.size() - 1; i >= 0; --i) {
-        vecR[i] = Phi2ToR(vecPhi2[i], eps);
-        //vecPhi1[i] = Phi1(vecR[i], eps);
-        vecPhi1[i] = 1 - vecPhi2[i];
-        vecLevelSet[i] = LevelSet(vecR[i]);
-      }
+  void ProblemImplicitBase::readPhi2(istream& inStream, double eps, 
+	                             Mesh* mesh, int implMesh, int comp )
+  { 
+    DOFVector< double >* r = getSignedDistance(implMesh, comp);
+    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
+    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
+    DOFVector< double >* levelSet = getLevelset(implMesh, comp);
+
+    TEST_EXIT( r != NULL )("no signed distance vector\n");
+    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
+    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
+    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");
+
+    bool checkSize = r->getSize() == phi1->getSize() &&
+	             r->getSize() == phi2->getSize();
+    TEST_EXIT(checkSize)("something went wrong\n");
+
+    readDofVec(inStream, phi2, mesh);
+
+    for (int i = r->getSize() - 1; i >= 0; --i) {
+      (*r)[i] = Phi2ToR((*r)[i], eps);
+      (*phi1)[i] = 1 - (*phi2)[i];
+      (*levelSet)[i] = LevelSet((*r)[i]);
     }
   }
+  
+  std::string ProblemImplicitBase::getDofFilename(std::string path, int implMesh)
+  {
+    std::string suffix = "[" + 
+	                 boost::lexical_cast< std::string >(implMesh) +
+			 "]";
+
+    std::string dofFilename("");
+    GET_PARAMETER(0, path + "dof file" + suffix, &dofFilename);
+    if (implMesh == 0 && dofFilename.length() == 0)
+      GET_PARAMETER(0, path + "dof file", &dofFilename);
+    return dofFilename;    
+  }
 
+  double ProblemImplicitBase::getEpsilon(std::string path, int implMesh)
+  {
+    std::string suffix = "[" + 
+	                 boost::lexical_cast< std::string >(implMesh) +
+			 "]";
 
-  bool ProblemImplicitScal::createImplicitMesh() 
+    double eps(-1);
+    GET_PARAMETER(0, path + "eps" + suffix, "%g", &eps);
+    if (implMesh == 0 && eps < 0)
+      GET_PARAMETER(0, path + "eps", "%g", &eps);
+    return eps;
+  }
+
+  int ProblemImplicitBase::getType(std::string path, int implMesh)
   {
+    std::string suffix = "[" + 
+	                 boost::lexical_cast< std::string >(implMesh) +
+			 "]";
+    int serType(-1);
+    GET_PARAMETER(0, path + "type" + suffix, "%d", &serType);
+    if (implMesh == 0 && serType < 0)
+      GET_PARAMETER(0, path + "type", "%d", &serType);
+    return serType;
+  }
+
+  DOFVector< double >* ProblemImplicitScal::getSignedDistance(int im , int )
+  { 
+    if (readImplMesh)
+      return r[im]; 
+    return NULL;
+  }
+
+  DOFVector< double >* ProblemImplicitScal::getPhi1(int im , int )
+  { 
+    if (readImplMesh)
+      return phi1[im];
+    return NULL;
+  }
+
+  DOFVector< double >* ProblemImplicitScal::getPhi2(int im , int )
+  { 
+    if (readImplMesh)
+      return phi2[im]; 
+    return NULL;
+  }
+
+  DOFVector< double >* ProblemImplicitScal::getLevelset(int im , int ) 
+  { 
+    if (readImplMesh)
+      return levelSet[im]; 
+    return NULL;
+  }
+
+  void ProblemImplicitScal::createMesh()
+  {
+    FUNCNAME("ProblemImplicitScal::createMesh()");
+    ProblemScal::createMesh();
+    std::string path = name + "->implicit mesh";
     std::string meshFilename("");
-    GET_PARAMETER(0, name + "->implicit mesh->mesh file", &meshFilename);
-    //if no file name is given, there's no implicit mesh
-    if (meshFilename.length() == 0)
+
+    GET_PARAMETER(0, path + "->macro file name", &meshFilename);
+    std::string serFilename("");
+    GET_PARAMETER(0, path + "->serialization file name", &serFilename);
+    if ( meshFilename.length() > 0)
+      mesh->setName(path);
+    else 
+      if ( serFilename.length() > 0 ) {
+        //INFO(5,1)("loading implicit mesh %s \n", serFilename.c_str());
+	std::string oldName = mesh->getName();
+        std::ifstream inStream(serFilename.c_str());
+        mesh->deserialize(inStream);
+        inStream.close();
+	mesh->setName(oldName);
+        readImplMesh = true;
+      }
+
+  }
+
+  void ProblemImplicitScal::initialize(Flag initFlag, ProblemScal* adaptProblem, 
+		                       Flag adoptFlag)
+  {
+     FUNCNAME("ProblemImplicitScal::initialize()");
+     ProblemScal::initialize(initFlag);
+     createImplicitMesh();
+  }
+
+  bool ProblemImplicitScal::createImplicitMesh() 
+  {
+    if (!isImplicitMesh())
       return false;
+    std::string path = name + "->implicit mesh->";
+    int nImplMesh(1);
+    GET_PARAMETER(0, name + "nr meshes", "%d", &nImplMesh);
 
-    std::string dofFilename("");
-    GET_PARAMETER(0, name + "->implicit mesh->dof file", &dofFilename); 
+    r.resize(nImplMesh);
+    phi1.resize(nImplMesh);
+    phi2.resize(nImplMesh);
+    levelSet.resize(nImplMesh);
+    for ( int i = 0; i < nImplMesh ; ++i) {
+      r[i] = new DOFVector< double >(getFeSpace(), "r");
+      phi1[i] = new DOFVector< double >(getFeSpace(), "phi1");
+      phi2[i] = new DOFVector< double >(getFeSpace(), "phi2");
+      levelSet[i] = new DOFVector< double >(getFeSpace(), "levelSet");
+      createImplicitMesh(path, i);
+    }
+    return true;
+  }
+
+  bool ProblemImplicitScal::createImplicitMesh(string path, int p)
+  {
+    
+    std::string dofFilename = getDofFilename(path, p);
     if (dofFilename.length() == 0)
       return false;
 
-    double eps(-1.0);
-    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
-    if (eps < 0)
+    double eps = getEpsilon(path, p);
+    if (eps < 0.0)
       return false;
-
-    int serType(-1);
-    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
+    int serType = getType(path, p);
     if (serType < 0)
       return false;
 
     TEST_EXIT(mesh != NULL)("the mesh was not created\n");
-
-    std::ifstream inStream(meshFilename.c_str());
-    mesh->deserialize(inStream);
-    inStream.close();
-    //create the fespace with the correct admin
-
-    createFeSpace();
-
-    r = new DOFVector< double >(getFeSpace(), "r");
-    phi1 = new DOFVector<double>(getFeSpace(), "phi1");
-    phi2 = new DOFVector<double>(getFeSpace(), "phi2");
-    levelSet = new DOFVector< double >(getFeSpace(), "levelSet");
-    inStream.open(dofFilename.c_str());
+    
+    std::ifstream inStream(dofFilename.c_str());
     switch (serType) {
     case 0:
-      calcspace::readR(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      readR(inStream, eps, mesh, p);
       break;
     case 1:
-      calcspace::readPhi1(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      readPhi1(inStream, eps, mesh, p);
       break;
     case 2:
-      calcspace::readPhi2(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      readPhi2(inStream, eps, mesh, p);
       break;
     default:
       WARNING("unknown implicit mesh type\n");
@@ -138,16 +261,36 @@ namespace AMDiS {
   }
 
 
-  void ProblemImplicitScal::initialize(Flag initFlag,
-		  ProblemScal* adaptProblem, Flag adoptFlag)
+  DOFVector< double >* ProblemImplicitVec::getSignedDistance(int im , int m) 
+  { 
+    if ( m >= r.size() || im >= r[m].size())
+	return NULL;
+      return (r[m])[im]; 
+    }
+
+  DOFVector< double >* ProblemImplicitVec::getPhi1(int im, int m)
   {
-     FUNCNAME("ProblemImplicitScal::initialize()");
-     ProblemScal::initialize(CREATE_MESH);
-     createImplicitMesh();
-     initFlag = initFlag & ~CREATE_MESH;
-     ProblemScal::initialize(initFlag);
+    if (m >= phi1.size() || im >= phi1[m].size())
+      return NULL;
+
+    return (phi1[m])[im];
   }
 
+  DOFVector< double >* ProblemImplicitVec::getPhi2(int im, int m)
+  {
+    if (m >= phi2.size() || im >= phi2[m].size())
+      return NULL;
+
+    return (phi2[m])[im];
+  }
+
+  DOFVector< double >* ProblemImplicitVec::getLevelset(int im, int m)
+  {
+    if (m >= levelSet.size() || im >= levelSet[m].size())
+      return NULL;
+
+    return (levelSet[m])[im];
+  }
 
   bool ProblemImplicitVec::createImplicitMesh()
   {
@@ -157,55 +300,64 @@ namespace AMDiS {
     phi2.resize(meshes.size());
     levelSet.resize(meshes.size());
     for (int i = 0; i < meshes.size(); ++i)
-      createImplicitMesh(i);    
+      createImplicitMesh(i);
+    return true;
   }
 
 
   bool ProblemImplicitVec::createImplicitMesh(int p) 
   {
-    std::string path = 
-      name + "->implicit mesh[" + boost::lexical_cast< std::string >(p) + "]->";
-    std::string meshFilename("");
-    GET_PARAMETER(0, path + "mesh file", &meshFilename);
-    //if no file name is given, there's no implicit mesh
-    if (meshFilename.length() == 0)
+    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
+    //INFO(5,1)("initing mesh %d\n", p);
+    if (!isImplicitMesh(p))
       return false;
-    std::string dofFilename("");
-    GET_PARAMETER(0,path + "dof file", &dofFilename); 
+    std::string path = name + "->implicit mesh[" 
+	             + boost::lexical_cast< std::string >(p) + "]->";
+    int nImplMeshes(0);
+    GET_PARAMETER(0, path + "nr meshes", "%d", &nImplMeshes);
+    //INFO(5,1)("has %d meshe(s)\n", nImplMeshes);
+    r[p].resize(nImplMeshes);
+    phi1[p].resize(nImplMeshes);
+    phi2[p].resize(nImplMeshes);
+    levelSet[p].resize(nImplMeshes);
+
+    for ( int i = 0; i < nImplMeshes ; ++i ) {
+      (r[p])[i] = new DOFVector< double >(getFeSpace(p), "r");
+      (phi1[p])[i] = new DOFVector< double >(getFeSpace(p), "phi1");
+      (phi2[p])[i] = new DOFVector< double >(getFeSpace(p), "phi2");
+      (levelSet[p])[i] = new DOFVector< double >(getFeSpace(p), "levelSet");
+      createImplicitMesh(path, i, p);
+    }
+    return true;
+  }
+  
+  bool ProblemImplicitVec::createImplicitMesh(std::string path, int implMesh, 
+		                              int comp)
+  {
+    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
+    std::string dofFilename = getDofFilename(path, implMesh);
     if (dofFilename.length() == 0)
       return false;
 
-    double eps(-1.0);
-    GET_PARAMETER(0, path + "eps", "%d", &eps);
+    double eps = getEpsilon(path, implMesh);
     if (eps < 0.0)
       return false;
-    int serType(-1);
-    GET_PARAMETER(0, path + "type", "%d", &serType);
+    int serType = getType(path, implMesh);
     if (serType < 0)
       return false;
-    TEST_EXIT(meshes[p] != NULL)("the mesh was not created\n");
-    std::ifstream inStream(meshFilename.c_str());
-    meshes[p]->deserialize(inStream);
-    inStream.close();
-    //create the fespace with the correct admin
-    createFeSpace(NULL);
-    r[p] = new DOFVector< double >(getFeSpace(p), "r");
-    phi1[p] = new DOFVector<double>(getFeSpace(p), "phi1");
-    phi2[p] = new DOFVector<double>(getFeSpace(p), "phi2");
-    levelSet[p] = new DOFVector< double >(getFeSpace(p), "levelSet");
-    inStream.open(dofFilename.c_str());
+
+    TEST_EXIT(meshes[comp] != NULL)("the mesh was not created\n");
+    
+    std::ifstream inStream(dofFilename.c_str());
     switch (serType) {
     case 0:
-      calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p],
-		       levelSet[p], meshes[p]);
+      readR(inStream, eps, meshes[comp], implMesh, comp);
       break;
     case 1:
-      calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p],
-			  levelSet[p], meshes[p]);
+      readPhi1(inStream, eps, meshes[comp], implMesh, comp);
       break;
     case 2:
-      calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p],
-			  levelSet[p], meshes[p]);
+      readPhi2(inStream, eps, meshes[comp], implMesh, comp);
       break;
     default:
       WARNING("unknown implicit mesh type\n");
@@ -214,14 +366,39 @@ namespace AMDiS {
     return true;
   }
 
+  void ProblemImplicitVec::createMesh()
+  {
+    FUNCNAME("ProblemImplicitVec::createMesh()");
+    ProblemVec::createMesh();
+    implMesh.resize(meshes.size());
+    for (int i = 0 ; i < meshes.size() ; ++i )
+    {
+      std::string path = name + "->implicit mesh[" +
+                         boost::lexical_cast< std::string >(i) + "]";
+      std::string meshFilename("");
+      GET_PARAMETER(0, path + "->macro file name", &meshFilename);
+      std::string serFilename("");
+      GET_PARAMETER(0, path + "->serialization file name", &serFilename);
+      implMesh[i] = true;
+      if (meshFilename.length() > 0)
+        meshes[i]->setName(path);
+      else if (serFilename.length() > 0) 
+      {
+        //INFO(5,1)("reading serialization file name %s \n", serFilename.c_str());
+	std::ifstream inStream(serFilename.c_str());
+	meshes[i]->deserialize(inStream);
+	inStream.close();
+      } else
+        implMesh[i] = false;
+    }
+
+  }
 
-  void ProblemImplicitVec::initialize(Flag initFlag,
-		  ProblemScal* adaptProblem, Flag adoptFlag)
+  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, 
+		                      Flag adoptFlag)
   {
      FUNCNAME("ProblemImplicitVec::initialize()");
-     ProblemVec::initialize(CREATE_MESH);
-     createImplicitMesh();
-     initFlag = initFlag & ~CREATE_MESH;
      ProblemVec::initialize(initFlag);    
+     createImplicitMesh();
   }
 }

AMDiS/src/ProblemImplicit.h

@@ -27,15 +27,36 @@
 #include "ProblemVec.h"
 
 namespace AMDiS {
-  class ProblemImplicitScal : public ProblemScal 
+  class ProblemImplicitBase 
+  {
+  public:
+    virtual bool createImplicitMesh() = 0;
+    virtual DOFVector< double >* getSignedDistance(int im = 0, int m = 0) = 0;
+    virtual DOFVector< double >* getPhi1(int im = 0, int m = 0) = 0;
+    virtual DOFVector< double >* getPhi2(int im = 0, int m = 0) = 0;
+    virtual DOFVector< double >* getLevelset(int im = 0, int m = 0) = 0;
+
+  protected:
+    void readDofVec(std::istream& , DOFVector<double>* , Mesh* );
+    void readR(std::istream& , double , Mesh* , int implMesh = 0, int comp = 0);
+    void readPhi1(std::istream& , double , Mesh* , int implMesh = 0, int comp = 0);
+    void readPhi2(std::istream& , double , Mesh* , int implMesh = 0, int comp = 0);
+    std::string getDofFilename(std::string path, int implMesh = 0);
+    double getEpsilon(std::string path, int implMesh = 0);
+    int getType(std::string path, int implMesh = 0);
+    virtual bool isImplicitMesh(int comp = 0 ) = 0;
+  };
+ 
+  class ProblemImplicitScal : public ProblemScal , public ProblemImplicitBase 
   {
   public:
     ProblemImplicitScal(std::string name, ProblemIterationInterface* pis = NULL)
       : ProblemScal(name, pis),
-        r(NULL),
-        phi1(NULL),
-        phi2(NULL),
-        levelSet(NULL)
+        r(0),
+        phi1(0),
+        phi2(0),
+        levelSet(0),
+        readImplMesh(false)
         {}
 
     /// @return: true if implicit mesh type is ok (see below) and false in 
@@ -52,26 +73,37 @@ namespace AMDiS {