From ae8b0b7179b14055a47d88b6c3cb71433e2ee99b Mon Sep 17 00:00:00 2001
From: Andreas Naumann <andreas.naumann@tu-dresden.de>
Date: Tue, 25 May 2010 12:53:35 +0000
Subject: [PATCH] implicit meshes for scalar problems
---
AMDiS/src/AMDiS.h | 1 +
AMDiS/src/MathFunctions.h | 23 +++++++++
AMDiS/src/ProblemScal.cc | 106 ++++++++++++++++++++++++++++++++++++++
AMDiS/src/ProblemScal.h | 38 +++++++++++++-
4 files changed, 167 insertions(+), 1 deletion(-)
create mode 100644 AMDiS/src/MathFunctions.h
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 65dc56f9..a39c5a0a 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -1,5 +1,6 @@
#ifndef AMDIS_H
#define AMDIS_H
+#include "MathFunctions.h"
#include "AbstractFunction.h"
#include "AdaptInfo.h"
#include "AdaptInstationary.h"
diff --git a/AMDiS/src/MathFunctions.h b/AMDiS/src/MathFunctions.h
new file mode 100644
index 00000000..356bfa20
--- /dev/null
+++ b/AMDiS/src/MathFunctions.h
@@ -0,0 +1,23 @@
+#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
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 40d8ce4d..aaf4fefb 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/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");
diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h
index 55fb1ce0..eb694961 100644
--- a/AMDiS/src/ProblemScal.h
+++ b/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);
};
}
--
GitLab