new demo for pfc added

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.")

demo/init/pfc.dat.2d

+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

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 

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);
 

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

demo/src/pfc.cc

+#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;
+};