codestyle and implicit meshes

AMDiS/CMakeLists.txt

-project("AMDiS")
+project(AMDiS)
 
 cmake_minimum_required(VERSION 2.6)
 
@@ -49,7 +49,8 @@ SET(AMDIS_SRC ${SOURCE_DIR}/DOFIndexed.cc
 	      ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/StandardProblemIteration.cc 
 	      ${SOURCE_DIR}/ProblemScal.cc 
-	      ${SOURCE_DIR}/ProblemVec.cc 
+	      ${SOURCE_DIR}/ProblemVec.cc
+	      ${SOURCE_DIR}/ProblemImplicit.cc 
 	      ${SOURCE_DIR}/ProblemVecDbg.cc 
 	      ${SOURCE_DIR}/DualTraverse.cc 
 	      ${SOURCE_DIR}/ElementData.cc 

AMDiS/src/MathFunctions.h

@@ -5,26 +5,26 @@
 
 namespace AMDiS {
   //converts signed distance to phasefield
-  inline double Phi1(double r, double eps) { return 0.5*(1-tanh(3*r/eps)); }
-  inline double Phi2(double r, double eps) { return 0.5*(1+tanh(3*r/eps)); }
+  inline double Phi1(double r, double eps) { return 0.5 * (1 - tanh(3 * r / eps)); }
+  inline double Phi2(double r, double eps) { return 0.5 * (1 + tanh(3 * r / eps)); }
   
   //levelset: positive (1) in the set, negative (-1) outside, zero on the boundary
   inline double LevelSet(double r)
   {
-    if(r<0)
+    if (r < 0)
       return 1;
-    if(r>0)
+    if (r > 0)
       return -1;
     return 0;
   }
 
   //convert Phi1 to r
   inline double Phi1ToR(double p1, double eps) {
-	return eps/3.0 * atanh( max(-1+1.e-14, min(1-1.e-14, 1-2*p1)) );
+	return eps / 3.0 * atanh( max(-1 + 1.0e-14, min(1 - 1.0e-14, 1 - 2 * p1)) );
   }
   //convert Phi2 to r
   inline double Phi2ToR(double p2, double eps) {
-	return eps/3.0 * atanh( max(-1+1.e-14, min(1-1.e-14, 1+2*p2)) );
+	return eps / 3.0 * atanh( max(-1 + 1.0e-14, min(1 - 1.0e-14, 1 + 2 * p2)) );
   }
 }
 #endif

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;
+      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());
+      }
+    }
+
+    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();
+
+      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("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]);
+      }
+    }
+
+    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();
+
+      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("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 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();
+
+      TEST_EXIT(vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size())("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]);
+      }
+    }
+
+
+  }
+
+  bool ProblemImplicitScal::createImplicitMesh() 
+  {
+    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)
+      return false;
+
+    std::string dofFilename("");
+    GET_PARAMETER(0, name + "->implicit mesh->dof file", &dofFilename); 
+    if (dofFilename.length() == 0)
+      return false;
+
+    double eps(-1);
+    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
+    if (eps<0)
+      return false;
+
+    int serType(-1);
+    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
+    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());
+    switch (serType) {
+    case 0:
+      calcspace::readR(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      break;
+    case 1:
+      calcspace::readPhi1(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      break;
+    case 2:
+      calcspace::readPhi2(inStream, eps, r, phi1, phi2, levelSet, mesh);
+      break;
+    default:
+      WARNING("unknown implicit mesh type\n");
+    }
+    inStream.close();
+    return true;
+  }
+
+  void ProblemImplicitScal::initialize(Flag initFlag, ProblemScal* adaptProblem, Flag adoptFlag)
+  {
+     FUNCNAME("ImplicitScal::initialize()");
+     ProblemScal::initialize(CREATE_MESH);
+     createImplicitMesh();
+     initFlag = initFlag & ~CREATE_MESH;
+     ProblemScal::initialize(initFlag);
+  }
+
+  bool ProblemImplicitVec::createImplicitMesh()
+  {
+    //check each mesh if it's an implicit one
+    r.resize(meshes.size());
+    phi1.resize(meshes.size());
+    phi2.resize(meshes.size());
+    levelSet.resize(meshes.size());
+    for(int i=0;i<meshes.size();++i) {
+      createImplicitMesh(i);
+    }
+  }
+
+  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)
+      return false;
+    std::string dofFilename("");
+    GET_PARAMETER(0,path + "dof file", &dofFilename); 
+    if (dofFilename.length() == 0)
+      return false;
+
+    double eps(-1);
+    GET_PARAMETER(0, path + "eps", "%d", &eps);
+    if(eps < 0)
+        return false;
+    int serType(-1);
+    GET_PARAMETER(0, path + "type", "%d", &serType);
+    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());
+    switch (serType) 
+    {
+      case 0:
+        calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
+        break;
+      case 1:
+	calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
+        break;
+      case 2:
+	calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p], levelSet[p], meshes[p]);
+        break;
+      default:
+        WARNING("unknown implicit mesh type\n");
+    }
+    inStream.close();
+    return true;
+  }
+
+  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, Flag adoptFlag)
+  {
+     FUNCNAME("ImplicitScal::initialize()");
+     ProblemVec::initialize(CREATE_MESH);
+     createImplicitMesh();
+     initFlag = initFlag & ~CREATE_MESH;
+     ProblemVec::initialize(initFlag);    
+  }
+}

AMDiS/src/ProblemImplicit.h

+#ifndef AMDIS_PROBLEMIMPLICIT_H
+#define AMDIS_PROBLEMIMPLICIT_H
+
+#include "AMDiS.h"
+#include "ProblemScal.h"
+#include "ProblemVec.h"
+
+namespace AMDiS {
+  class ProblemImplicitScal : public ProblemScal 
+  {
+
+  public:
+    ProblemImplicitScal(std::string name, ProblemIterationInterface* pis = NULL)
+      : ProblemScal(name, pis),
+        r(NULL),
+        phi1(NULL),
+        phi2(NULL),
+        levelSet(NULL)
+        {}
+
+    /// @return: true if implicit mesh type is ok (see below) and false in all other cases
+    /// initfile entries: 
+    /// name+"->implicit mesh->mesh file": name of serialized mesh
+    /// name+"->implicit mesh->dof file": name of serialized dofvector
+    /// name+"->implicit mesh->type": type of serialized dofvector. 
+    /// type can be: 
+    ///    0: dofvector is a signed distance r
+    ///    1: dofvector is a phasefield phi=0.5*(1-tanh(3*r/eps))
+    ///    2: dofvector is a phasefield phi=0.5*(1+tanh(3*r/eps))
+    /// assumes:
+    ///    a) already created, but not yet initialized mesh
+    ///    b) already created feSpace
+    virtual bool createImplicitMesh();
+    virtual void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL,
+			    Flag adoptFlag = INIT_NOTHING);
+
+  protected:
+    /// DOFVector for a signed distance
+    DOFVector< double > *r;
+
+    /// DOFVector for the phasefield function 0.5*(1-tanh(3*r/eps))
+    DOFVector< double > *phi1;
+
+    /// DOFVector for the phasefield function 0.5*(1+tanh(3*r/eps))
+    DOFVector< double > *phi2;
+
+    /// DOFVector for the levelset function (levelSet(x): x \in \Omega: 1, x \not \in Omega: -1, x \in \Gamma: 0)
+    DOFVector< double > *levelSet;
+
+  };
+
+  class ProblemImplicitVec : public ProblemVec
+  {
+  public:
+    ProblemImplicitVec(std::string name, ProblemIterationInterface* problem = NULL)
+      : ProblemVec(name, problem),
+	r(0),
+	phi1(0),
+	phi2(0),
+	levelSet(0)
+        {}
+
+    virtual bool createImplicitMesh();
+    virtual void initialize(Flag initFlag, ProblemScal *adoptProblem=NULL, 
+		            Flag adoptFlag = INIT_NOTHING);
+
+  protected:
+
+    bool createImplicitMesh(int p); 
+    /// DOFVector for a signed distance
+    std::vector< DOFVector< double >* > r;
+
+    /// DOFVector for the phasefield function 0.5*(1-tanh(3*r/eps))
+    std::vector< DOFVector< double >* > phi1;
+
+    /// DOFVector for the phasefield function 0.5*(1+tanh(3*r/eps))
+    std::vector< DOFVector< double >* > phi2;
+
+    /// DOFVector for the levelset function (levelSet(x): x \in \Omega: 1, x \not \in Omega: -1, x \in \Gamma: 0)
+    std::vector< DOFVector< double >* > levelSet;
+
+  };
+}
+#endif

AMDiS/src/ProblemScal.cc

@@ -151,110 +151,6 @@ namespace AMDiS {
       rhs->getBoundaryManager()->addBoundaryCondition(periodic);
   }
 
-  void ProblemScal::readR(std::istream& inStream,double eps) {
-    readDofvec(inStream,r);
-    std::vector< double >& vecR=r->getVector();
-    std::vector< double >& vecPhi1=phi1->getVector();
-    std::vector< double >& vecPhi2=phi2->getVector();
-    std::vector< double >& vecLevelSet=levelSet->getVector();
-
-    TEST_EXIT(vecR.size()!=vecPhi1.size() && vecPhi2.size()!=vecR.size())("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]);
-    }
-  }
-
-  void ProblemScal::readPhi1(istream& inStream,double eps) {
-    readDofvec(inStream,phi1);
-    std::vector< double >& vecR=r->getVector();
-    std::vector< double >& vecPhi1=phi1->getVector();
-    std::vector< double >& vecPhi2=phi2->getVector();
-    std::vector< double >& vecLevelSet=levelSet->getVector();
-
-    TEST_EXIT(vecR.size()!=vecPhi1.size() && vecPhi2.size()!=vecR.size())("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);
-      vecLevelSet[i]=LevelSet(vecR[i]);
-    }
-  }
-
-  void ProblemScal::readPhi2(istream& inStream,double eps) {
-    readDofvec(inStream,phi2);
-    std::vector< double >& vecR=r->getVector();
-    std::vector< double >& vecPhi1=phi1->getVector();
-    std::vector< double >& vecPhi2=phi2->getVector();
-    std::vector< double >& vecLevelSet=levelSet->getVector();
-
-    TEST_EXIT(vecR.size()!=vecPhi1.size() && vecPhi2.size()!=vecR.size())("something went wrong\n");
-    for(int i=vecR.size()-1;i>=0;--i) {
-      vecR[i]=Phi2ToR(vecPhi1[i],eps);
-      vecPhi1[i]=Phi1(vecR[i],eps);
-      vecLevelSet[i]=LevelSet(vecR[i]);
-    }
-  }
-
-  void ProblemScal::readDofvec(std::istream& in, DOFVector< double >* vec) 
-  {
-    long size = mesh->getNumberOfVertices();
-    std::string buf;
-    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());
-    }
-  }
-
-  bool ProblemScal::createImplicitMesh() 
-  {
-    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)
-	return false;
-    std::string dofFilename("");
-    GET_PARAMETER(0,name+"->implicit mesh->dof file",&dofFilename); 
-    if(dofFilename.length()==0)
-        return false;
-
-    double eps(-1);
-    GET_PARAMETER(0,name+"->implicit mesh->eps","%d",&eps);
-    if(eps<0)
-        return false;
-    int serType(-1);
-    GET_PARAMETER(0,name+"->implicit mesh->type","%d",&serType);
-    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());
-    switch(serType) {
-      case 0:
-        readR(inStream,eps);
-        break;
-      case 1:
-        readPhi1(inStream,eps);
-        break;
-      case 2:
-        readPhi2(inStream,eps);
-        break;
-      default:
-        WARNING("unknown implicit mesh type\n");
-    }
-    inStream.close();
-    return true;
-  }
 
   void ProblemScal::createMesh()
   {
@@ -364,7 +260,6 @@ namespace AMDiS {
 			       Flag adoptFlag) 
   {
     FUNCNAME("Problem::initialize()");
-
     // === create mesh ===
     if (mesh) { 
       WARNING("mesh already created\n");
@@ -389,7 +284,9 @@ namespace AMDiS {
     if (!mesh) 
       WARNING("no mesh created\n");
 
-    createImplicitMesh();
+    if (refinementManager == NULL)
+	WARNING("refinement manager not set\n");
+
     // === create fespace ===
     if (feSpace) {
       WARNING("feSpace already created\n");
@@ -529,6 +426,7 @@ namespace AMDiS {
 	}
       }
     }
+	
   }
 
   void ProblemScal::createFESpace()

AMDiS/src/ProblemScal.h

@@ -52,12 +52,7 @@ namespace AMDiS {
 	useGetBound(true),
 	refinementManager(NULL),
 	coarseningManager(NULL),
-	info(10),
-	r(NULL),
-	phi1(NULL),
-	phi2(NULL),
-	levelSet(NULL)
-    {}
+	info(10)    {}
 
     /// Destructor
     virtual ~ProblemScal();
@@ -70,22 +65,6 @@ namespace AMDiS {
     /// Used in \ref initialize().
     virtual void createMesh();
 
-    /// Used in \ref createMesh().
-    /// @return: true if implicit mesh type is ok (see below) and false in all other cases
-    /// initfile entries: 
-    /// name+"->implicit mesh->mesh file": name of serialized mesh
-    /// name+"->implicit mesh->dof file": name of serialized dofvector
-    /// name+"->implicit mesh->type": type of serialized dofvector. 
-    /// type can be: 
-    ///    0: dofvector is a signed distance r
-    ///    1: dofvector is a phasefield phi=0.5*(1-tanh(3*r/eps))
-    ///    2: dofvector is a phasefield phi=0.5*(1+tanh(3*r/eps))
-    /// assumes:
-    ///    a) already created, but not yet initialized mesh
-    ///    b) already created feSpace
-    virtual bool createImplicitMesh();
-    void readDofvec(std::istream&, DOFVector< double >*);
-
     /// Used in \ref initialize().
     virtual void createFESpace();
 
@@ -384,18 +363,6 @@ namespace AMDiS {
     /// DOFVector for the right hand side 
     DOFVector<double> *rhs;
 
-    /// DOFVector for a signed distance
-    DOFVector< double > *r;
-
-    /// DOFVector for the phasefield function 0.5*(1-tanh(3*r/eps))
-    DOFVector< double > *phi1;
-
-    /// DOFVector for the phasefield function 0.5*(1+tanh(3*r/eps))
-    DOFVector< double > *phi2;
-
-    /// DOFVector for the levelset function (levelSet(x): x \in \Omega: 1, x \not \in Omega: -1, x \in \Gamma: 0)
-    DOFVector< double > *levelSet;
-
     /// System matrix
     DOFMatrix *systemMatrix;
 
@@ -422,9 +389,7 @@ namespace AMDiS {
     /// Info level.
     int info;
 
-    void readR(std::istream&, double);
-    void readPhi1(std::istream&,double);
-    void readPhi2(std::istream&,double);
+  
   };
 
 }