diff --git a/demo/CMakeLists.txt b/demo/CMakeLists.txt
index b9da147ef0c6c42d2b72de16eba207814eccbaaf..3dd1a835ffddd901d2e54bb4de0b576700a6e68c 100644
--- a/demo/CMakeLists.txt
+++ b/demo/CMakeLists.txt
@@ -13,8 +13,6 @@ endif(AMDIS_FOUND)
set(ball src/ball.cc)
set(bunny src/bunny.cc)
set(ellipt src/ellipt.cc)
-set(elliptBase src/elliptBaseProblem.cc)
-set(elliptImplicit src/elliptImplicit.cc)
set(heat src/heat.cc)
set(neumann src/neumann.cc)
set(nonlin src/nonlin.cc)
@@ -24,8 +22,6 @@ set(sphere src/sphere.cc)
set(torus src/torus.cc)
set(vecellipt src/vecellipt.cc)
set(vecheat src/vecheat.cc)
-set(rhsLaplace src/rhsLaplace.cc)
-set(laplace src/laplace.cc)
add_executable("ball" ${ball})
target_link_libraries("ball" ${BASIS_LIBS})
@@ -85,11 +81,22 @@ if (TOOLS_DIR)
include_directories(${TOOLS_DIR}/diffuseDomain)
include_directories(${TOOLS_DIR}/baseProblems)
+ set(elliptBase src/elliptBaseProblem.cc)
+ add_executable("elliptBase" ${elliptBase})
+ target_link_libraries("elliptBase" ${BASIS_LIBS})
+
+ set(elliptImplicit src/elliptImplicit.cc)
+ add_executable("elliptImplicit" ${elliptImplicit})
+ target_link_libraries("elliptImplicit" ${BASIS_LIBS})
set(cahn_hilliard src/cahn_hilliard.cc ${TOOLS_DIR}/baseProblems/CahnHilliard.cc)
add_executable("cahn_hilliard" ${cahn_hilliard})
target_link_libraries("cahn_hilliard" ${BASIS_LIBS})
+ set(pfc src/pfc.cc ${TOOLS_DIR}/baseProblems/PhaseFieldCrystal_Base.cc ${TOOLS_DIR}/diffuseDomain/Helpers.cc)
+ add_executable("pfc" ${pfc})
+ target_link_libraries("pfc" ${BASIS_LIBS})
+
set(drivenCavity ${TOOLS_DIR}/diffuseDomain/POperators.cc
${TOOLS_DIR}/baseProblems/CahnHilliard.cc
${TOOLS_DIR}/baseProblems/NavierStokes_TaylorHood.cc
@@ -97,11 +104,6 @@ if (TOOLS_DIR)
add_executable("drivenCavity" ${drivenCavity})
target_link_libraries("drivenCavity" ${BASIS_LIBS})
- add_executable("elliptBase" ${elliptBase})
- target_link_libraries("elliptBase" ${BASIS_LIBS})
-
- add_executable("elliptImplicit" ${elliptImplicit})
- target_link_libraries("elliptImplicit" ${BASIS_LIBS})
endif()
else()
message(WARNING "No tools directory specified! Some demos will not be build.")
diff --git a/demo/init/pfc.dat.2d b/demo/init/pfc.dat.2d
new file mode 100644
index 0000000000000000000000000000000000000000..9af20ab72952348e8a92944e9839d1e5367f4e30
--- /dev/null
+++ b/demo/init/pfc.dat.2d
@@ -0,0 +1,55 @@
+dimension of world: 2
+
+% =================== MESH ================================
+mesh->macro file name: ./macro/periodic.square.2d
+mesh->periodic file: ./macro/periodic.square.per
+mesh->global refinements: 10
+mesh->check: 0
+
+lattice: 4*M_PI/sqrt(3)
+mesh->scale mesh: 1
+mesh->dimension: [4*${lattice}, 2*sqrt(3)*${lattice}]
+
+% ============== USER-PARAMETER ==========================
+pfc->r: -0.4
+pfc->rho0: 1.0
+pfc->density: -0.3
+pfc->use mobility: 0
+
+pfc->density amplitude: 0.2
+
+% ============= PROBLEM-SPACE ==================================
+pfc->space->mesh: mesh
+pfc->space->components: 3
+pfc->space->polynomial degree[0]: 1
+pfc->space->polynomial degree[1]: 1
+pfc->space->polynomial degree[2]: 1
+pfc->space->dim: 2
+
+
+% ================== SOLVER ======================================
+pfc->space->solver: umfpack
+pfc->space->solver->symmetric strategy: 0
+pfc->space->solver->store symbolic: 0
+pfc->space->solver->info: 5
+
+
+% ==================== TIMESTEPS ===============================
+adapt->timestep: 1.e-1
+adapt->max timestep: 1e+10
+adapt->min timestep: 1e-10
+adapt->start time: 0.0
+adapt->end time: 50.0
+
+
+% =================== OUTPUT =========================================
+pfc->space->output->filename: ./output/pfc.2d
+pfc->space->output->ParaView animation: 1
+pfc->space->output->ParaView format: 1
+pfc->space->output->write every i-th timestep: 1
+pfc->space->output->append index: 1
+pfc->space->output->index length: 9
+pfc->space->output->index decimals: 7
+
+
+WAIT: 1
diff --git a/demo/macro/macro.stand.2d b/demo/macro/macro.stand.2d
index 8caa1a56795ca498a88cf5742929bd9f0729f1a5..a764b1ce6c7958c252d63a95464f789f1dcaa2ec 100644
--- a/demo/macro/macro.stand.2d
+++ b/demo/macro/macro.stand.2d
@@ -17,11 +17,11 @@ element boundaries:
0 0 2
vertex coordinates:
- 0.0 0.0
- 1.0 0.0
+ -1.0 -1.0
+ 1.0 -1.0
1.0 1.0
- 0.0 1.0
- 0.5 0.5
+ -1.0 1.0
+ 0.0 0.0
element neighbours:
1 3 -1
diff --git a/demo/src/cahn_hilliard.cc b/demo/src/cahn_hilliard.cc
index bb198d1c4e7fc096395554fba316f0b459575dd1..bb66191ce748c7295e52546554ce81d46644f30a 100644
--- a/demo/src/cahn_hilliard.cc
+++ b/demo/src/cahn_hilliard.cc
@@ -82,7 +82,7 @@ int main(int argc, char** argv)
RefinementTimeInterface timeInterface(&chProb);
- // Adapt-Infos
+ // Adapt-Infos
AdaptInfo adaptInfo("adapt", chProb.getNumComponents());
AdaptInstationary adaptInstat("adapt", chProb, adaptInfo, timeInterface, adaptInfo);
diff --git a/demo/src/chns/CahnHilliardNavierStokes.h b/demo/src/chns/CahnHilliardNavierStokes.h
index 39d7c4078243bda54126f71904e692bdfc5bb18e..0caa2e2b2a58297c8f1bf52ccd861031fcdbd7a1 100644
--- a/demo/src/chns/CahnHilliardNavierStokes.h
+++ b/demo/src/chns/CahnHilliardNavierStokes.h
@@ -13,7 +13,6 @@
#include "MeshFunction_Level.h"
#include "POperators.h"
-#include "CouplingOperators.h"
using namespace AMDiS;
@@ -22,7 +21,8 @@ using namespace AMDiS;
*
* \brief
*/
-template<typename CH_Type=CahnHilliard, typename NS_Type=NavierStokes_TaylorHood>
+template<typename CH_Type = CahnHilliard,
+ typename NS_Type = NavierStokes_TaylorHood>
class CahnHilliardNavierStokes : public CouplingIterationInterface,
public CouplingTimeInterface,
public CouplingProblemStat
@@ -53,6 +53,17 @@ public:
}
+ /**
+ * Add the problems to the iterationInterface, timeInterface and couplingProblemStat.
+ * As a consequence all problem can be initialized one after another and in the
+ * adaption loop they are solved in rotation.
+ * At the end the problems are filled with operators and coupling operators as well as
+ * boundary conditions are added.
+ *
+ * In the adaption loop the problems are solved the same order as they are added to the
+ * iterationInterface in this method. This order can be changed manually in the oneIteration
+ * method.
+ **/
void initialize(AdaptInfo *adaptInfo)
{
for (size_t i = 0; i < chProb->getNumProblems(); i++)
@@ -73,10 +84,15 @@ public:
// fillOperators and fillBoundaryConditions for chProb and nsProb
nsProb->initTimeInterface();
chProb->initTimeInterface();
- fillCouplingBoundaryConditions();
+// fillCouplingBoundaryConditions();
}
+ /**
+ * In this method the operators responsible for the coupling of all problems are defined
+ * and added to the problem. The access to the problems is through the baseBroblem->getProblem()
+ * method. All baseProblem provide FiniteElemSpaces and solution vectors.
+ **/
void fillCouplingOperators()
{ FUNCNAME("CahnHilliardNavierStokes::fillCouplingOperators()");
MSG("CahnHilliardNavierStokes::fillCouplingOperators()");
@@ -86,12 +102,6 @@ public:
Operator *opNuGradC = new Operator(nsProb->getFeSpace(i), chProb->getFeSpace(0));
opNuGradC->addTerm(new VecAndPartialDerivative_ZOT(chProb->getSolution()->getDOFVector(1), chProb->getSolution()->getDOFVector(0), i));
nsProb->getProblem(0)->addVectorOperator(opNuGradC, i, &surfaceTension, &surfaceTension);
-//
-// // stabilizing term
-// Operator *stabilBeltrami = new Operator(nsProb->getFeSpace(i), nsProb->getFeSpace(i));
-// stabilBeltrami->addTerm(new MatrixGradient_SOT(
-// chProb->getSolution()->getDOFVector(0), new ProjectionMatrix(surfaceTension), new DivFct())); // factor = sigma ?
-// nsProb->getProblem(0)->addMatrixOperator(stabilBeltrami, i, i, nsProb->getTau(), nsProb->getTau());
}
// < v * grad(c) , theta >
@@ -99,50 +109,48 @@ public:
opVGradC->addTerm(new WorldVector_FOT(nsProb->getVelocity(), -1.0), GRD_PSI);
chProb->getProblem()->addMatrixOperator(opVGradC, 0, 0);
}
-
-
- void fillCouplingBoundaryConditions()
- { FUNCNAME("CahnHilliardNavierStokes::fillCouplingBoundaryConditions()");
-
- }
- /// Solves the initial problem.
- virtual void solveInitialProblem(AdaptInfo *adaptInfo)
+ /**
+ * Solves the initial problems.Before this is done the mesh is refined by the Refinement-class
+ * This provides local refinement functions for phaseField-refinement or signed-distance-refinement.
+ * See MeshFunction_Level.h for some details.
+ **/
+ void solveInitialProblem(AdaptInfo *adaptInfo)
{
- refFunction= new PhaseFieldRefinement(chProb->getMesh());
- refinement= new RefinementLevelDOF(
- chProb->getFeSpace(),
- refFunction,
- new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(0))); // phaseField-DOFVector
+ // meshFunction for klocal refinement around the interface of the phasefield-DOFVector
+ refFunction = new PhaseFieldRefinement(chProb->getMesh());
+
+ if (chProb->getDoubleWellType() == 0) {
+ refinement = new RefinementLevelDOF(
+ chProb->getFeSpace(),
+ refFunction,
+ chProb->getSolution()->getDOFVector(0)); // phaseField-DOFVector in [0,1]
+ } else {
+ refinement = new RefinementLevelDOF(
+ chProb->getFeSpace(),
+ refFunction,
+ new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(0))); // phaseField-DOFVector in [-1,1]
+ }
- // set initial values
+ // initial refinement
refinement->refine(0);
- for (int i= 0; i < 3; ++i) {
- chProb->solveInitialProblem(adaptInfo); // initial phaseField
- refinement->refine((i < 4 ? 4 : 10));
+ for (int i = 0; i < 3; i++) {
+ chProb->solveInitialProblem(adaptInfo); // update phaseField-DOFVector
+ refinement->refine((i < 4 ? 4 : 10)); // do k steps of local refinement
}
-
- CouplingTimeInterface::solveInitialProblem(adaptInfo);
- }
-
- /// Called at the beginning of each timestep
- virtual void initTimestep(AdaptInfo *adaptInfo)
- {
- CouplingTimeInterface::initTimestep(adaptInfo);
- }
-
-
- Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
- {
- CouplingIterationInterface::oneIteration(adaptInfo, toDo);
+ // solve all initial problems of the problems added to the CouplingTimeInterface
+ CouplingTimeInterface::solveInitialProblem(adaptInfo);
}
- /// Called at the end of each timestep.
- virtual void closeTimestep(AdaptInfo *adaptInfo)
+ /**
+ * Called at the end of each timestep.
+ * If the phase-field has changed the mesh is updated by the refinement-method
+ **/
+ void closeTimestep(AdaptInfo *adaptInfo)
{
CouplingTimeInterface::closeTimestep(adaptInfo);
@@ -151,17 +159,17 @@ public:
protected:
- CH_Type *chProb;
- NS_Type *nsProb;
+ CH_Type *chProb; // CahnHilliard baseProblem
+ NS_Type *nsProb; // Navier-Stokes baseProblem
PhaseFieldRefinement* refFunction;
RefinementLevelDOF *refinement;
- unsigned dim;
- unsigned dow;
+ unsigned dim; // dimension of the mesh
+ unsigned dow; // dimension of the world
- double sigma;
- double surfaceTension;
+ double sigma; // coupling parameter to calculate the surface tension
+ double surfaceTension;// := sigma/epsilon
};
#endif // CAHN_HILLIARD_NAVIER_STOKES_H
diff --git a/demo/src/pfc.cc b/demo/src/pfc.cc
new file mode 100644
index 0000000000000000000000000000000000000000..76220ba8d6f9816fed352d928814f1b6575cdccd
--- /dev/null
+++ b/demo/src/pfc.cc
@@ -0,0 +1,93 @@
+#include "AMDiS.h"
+#include "PhaseFieldCrystal_Base.h"
+#include "Helpers.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+/// Functor that gives random values for each coordinate argument
+struct RandomValues : public AbstractFunction<double, WorldVector<double> >
+{
+ RandomValues(double mean_, double amplitude_) : mean(mean_), amplitude(amplitude_)
+ {
+ std::srand(time(NULL));
+ }
+
+ double operator()(const WorldVector<double> &x) const
+ {
+ return mean + amplitude * (std::rand() / static_cast<double>(RAND_MAX-0.5));
+ }
+
+private:
+ double mean;
+ double amplitude;
+};
+
+/// derived class from BaseProblem->PhaseFieldCrystal
+struct PFC_Demo : public PhaseFieldCrystal
+{
+ PFC_Demo(const std::string &name_) : PhaseFieldCrystal(name_) {}
+
+ // generate initial solution for evolution equation
+ void solveInitialProblem(AdaptInfo *adaptInfo)
+ { FUNCNAME("PFC_Demo::solveInitialProblem()");
+
+ Flag initFlag = initDataFromFile(adaptInfo);
+ if (initFlag.isSet(DATA_ADOPTED))
+ return;
+
+ double amplitude = 0.1;
+ Parameters::get(name + "->density amplitude",amplitude);
+
+ // random data
+ prob->getSolution()->getDOFVector(0)->interpol(new RandomValues(density, amplitude));
+ }
+
+ // add periodic boundary condition on all boundaries
+ void fillBoundaryConditions()
+ {
+ for (size_t i = 0; i < prob->getNumComponents(); i++) {
+ for (size_t j = 0; j < prob->getNumComponents(); j++) {
+ prob->addPeriodicBC(-1, i, j);
+ prob->addPeriodicBC(-2, i, j);
+ }
+ }
+ }
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+ AMDiS::init(argc, argv);
+
+ PFC_Demo pfcProb("pfc");
+ pfcProb.initialize(INIT_ALL);
+
+ // Adapt-Infos
+ AdaptInfo adaptInfo("adapt", pfcProb.getNumComponents());
+ // adaption loop - solve problems up to stationary solution
+ AdaptInstationary adaptInstat("adapt", pfcProb, adaptInfo, pfcProb, adaptInfo);
+
+ // Scale Mesh
+ bool scaleMesh = false;
+ Initfile::get("mesh->scale mesh",scaleMesh);
+ if (scaleMesh) {
+ WorldVector<double> scale; scale.set(1.0);
+ Initfile::get("mesh->dimension",scale);
+ Helpers::scaleMesh(pfcProb.getMesh(), scale);
+ }
+
+ ptime start_time = microsec_clock::local_time();
+ pfcProb.initTimeInterface(); // fill operators and BC
+ int error_code = adaptInstat.adapt();
+ time_duration td = microsec_clock::local_time()-start_time;
+
+ MSG("elapsed time= %d sec\n", td.total_seconds());
+
+ AMDiS::finalize();
+
+ return error_code;
+};