diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 109edee7b5832bc948e2ba27755fb440e8f131d6..af080f62dbe0b50b98c6d87add33238223c00883 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -1,4 +1,4 @@
-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
diff --git a/AMDiS/src/MathFunctions.h b/AMDiS/src/MathFunctions.h
index 0b440198f243888afb4bb801fe1911a364371228..335d144c3e8784efb8a356e6b72c7f468e16bff3 100644
--- a/AMDiS/src/MathFunctions.h
+++ b/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
diff --git a/AMDiS/src/ProblemImplicit.cc b/AMDiS/src/ProblemImplicit.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fb11fd9e6311b0bbb83907cb1835aa6f34f41182
--- /dev/null
+++ b/AMDiS/src/ProblemImplicit.cc
@@ -0,0 +1,216 @@
+#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);
+ }
+}
diff --git a/AMDiS/src/ProblemImplicit.h b/AMDiS/src/ProblemImplicit.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ee6e7a12a0fda014ec4812f202c138c0c13d921
--- /dev/null
+++ b/AMDiS/src/ProblemImplicit.h
@@ -0,0 +1,84 @@
+#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
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index aaf4fefb5ff5fb5c4bc055cc9fd6d35abe6ca198..f9714bd7d7edccf1d59050c3c366b9fd54cdcd92 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/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()
diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h
index eb6949614324490d89e5532bf831eb83c0a2c461..e1eab402210c00047ce13a708b3f6d23af427221 100644
--- a/AMDiS/src/ProblemScal.h
+++ b/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);
+
};
}