implicit meshes for scalar problems

AMDiS/src/AMDiS.h

 #ifndef AMDIS_H
 #define AMDIS_H
+#include "MathFunctions.h"
 #include "AbstractFunction.h"
 #include "AdaptInfo.h"
 #include "AdaptInstationary.h"

AMDiS/src/MathFunctions.h

+#ifndef AMDIS_MATHFUNCTIONS_H
+#define AMDIS_MATHFUNCTIONS_H
+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)); }
+  
+  //levelset: positive (1) in the set, negative (-1) outside, zero on the boundary
+  inline double LevelSet(double r)
+  {
+    if(r<0)
+      return 1;
+    if(r>0)
+      return -1;
+    return 0;
+  }
+
+  //convert Phi1 to r
+  inline double Phi1ToR(double p1, double eps) { return atanh(1-2*p1)*eps/3; }
+  //convert Phi2 to r
+  inline double Phi2ToR(double p2, double eps) { return atanh(1+2*p2)*eps/3; }
+}
+#endif

AMDiS/src/ProblemScal.cc

@@ -151,6 +151,111 @@ 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()
   {
     TEST_EXIT(Parameters::initialized())("parameters not initialized\n");
@@ -284,6 +389,7 @@ namespace AMDiS {
     if (!mesh) 
       WARNING("no mesh created\n");
 
+    createImplicitMesh();
     // === create fespace ===
     if (feSpace) {
       WARNING("feSpace already created\n");

AMDiS/src/ProblemScal.h

@@ -52,7 +52,11 @@ namespace AMDiS {
 	useGetBound(true),
 	refinementManager(NULL),
 	coarseningManager(NULL),
-	info(10)
+	info(10),
+	r(NULL),
+	phi1(NULL),
+	phi2(NULL),
+	levelSet(NULL)
     {}
 
     /// Destructor
@@ -66,6 +70,22 @@ 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();
 
@@ -364,6 +384,18 @@ 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;
 
@@ -389,6 +421,10 @@ namespace AMDiS {
   
     /// Info level.
     int info;
+
+    void readR(std::istream&, double);
+    void readPhi1(std::istream&,double);
+    void readPhi2(std::istream&,double);
   };
 
 }