some corrections in the extensions

extensions/base_problems/BaseProblem.hh

@@ -18,7 +18,7 @@ BaseProblem<ProblemType>::BaseProblem(const std::string &name_, bool createProbl
   dow = Global::getGeo(WORLD);
   
   Initfile::get(name + "->secure iteration", secureIteration);
-};
+}
 
 
 template<typename ProblemType>
@@ -163,7 +163,7 @@ Flag BaseProblem<ProblemType>::initDataFromFile(AdaptInfo *adaptInfo)
   initFlag.setFlag(MESH_ADOPTED);
 
   return initFlag;
-};
+}
 
 
 template<typename ProblemType>
@@ -187,7 +187,7 @@ void BaseProblem<ProblemType>::beginIteration(AdaptInfo *adaptInfo)
       fixedMatrixTimestep[i].second,
       fixedMatrixCell);
   }
-};
+}
 
 
 template<typename ProblemType>
@@ -247,7 +247,7 @@ Flag BaseProblem<ProblemType>::oneIteration(AdaptInfo *adaptInfo,
   oldMeshChangeIdx = getMesh()->getChangeIndex();
 
   return flag;
-};
+}
 
 
 template<typename ProblemType>
@@ -256,7 +256,7 @@ void BaseProblem<ProblemType>::endIteration(AdaptInfo *adaptInfo)
 	
   MSG("\n");
   MSG(("[[ end of <"+name+"> iteration ]]\n").c_str());
-};
+}
 
 
 template<typename ProblemType>
@@ -268,8 +268,7 @@ void BaseProblem<ProblemType>::closeTimestep(AdaptInfo *adaptInfo)
     && adaptInfo->getStartTime() < adaptInfo->getEndTime()) {
     writeFiles(adaptInfo, false);
   }
-};
-
+}
 
 
 template<typename ProblemType>
@@ -282,11 +281,11 @@ void BaseProblem<ProblemType>::addTimeOperator(ProblemStat *prob, int i, int j)
 
   prob->addMatrixOperator(op, i, j, getInvTau(), getInvTau());
   prob->addVectorOperator(opRhs, i, getInvTau(), getInvTau());
-};
+}
 
 
 template<typename ProblemType>
 void BaseProblem<ProblemType>::addTimeOperator(RosenbrockStationary *prob, int i, int j)
 {
   prob->addTimeOperator(i,j);
-};
+}

extensions/base_problems/NavierStokes_TH_MultiPhase.cc

@@ -50,10 +50,12 @@ void NavierStokes_TH_MultiPhase::initTimestep(AdaptInfo *adaptInfo)
 { FUNCNAME("NavierStokes_TH_MultiPhase::initTimestep()");
 
   super::initTimestep(adaptInfo);
-  LinearInterpolation1 dLI(density1, density2);
-  LinearInterpolation1 vLI(viscosity1, viscosity2);
-  transformDOF(multiPhase, densityPhase, &dLI);
-  transformDOF(multiPhase, viscosityPhase, &vLI);
+  if (multiPhase) {
+    LinearInterpolation1 dLI(density1, density2);
+    LinearInterpolation1 vLI(viscosity1, viscosity2);
+    transformDOF(multiPhase, densityPhase, &dLI);
+    transformDOF(multiPhase, viscosityPhase, &vLI);
+  }
 };
 
 

extensions/base_problems/NavierStokes_TH_MultiPhase.h

@@ -50,10 +50,36 @@ public: // methods
    **/
   void setMultiPhase(DOFVector<double>* p) { multiPhase=p; }
 
-protected: // methods
-
   virtual void fillOperators();
   virtual void addLaplaceTerm(int i);
+  
+  // get-methods
+  
+  DOFVector<double>* getViscosityPhase()
+  {
+    return viscosityPhase;
+  }
+  
+  DOFVector<double>* getDensityPhase()
+  {
+    return densityPhase;
+  }
+  
+  double getViscosity(int i=0)
+  {
+    if (i == 0)
+      return viscosity1;
+    else
+      return viscosity2;
+  }
+  
+  double getDensity(int i=0)
+  {
+    if (i == 0)
+      return density1;
+    else
+      return density2;
+  }
 
 protected: // variables
 

extensions/base_problems/PhaseFieldCrystal_Base.cc

@@ -8,13 +8,11 @@ PhaseFieldCrystal::PhaseFieldCrystal(const std::string &name_) :
   useMobility(false),
   tempParameter(-0.6),
   r(-0.4),		// temperature deviation
-  rho0(1.0),		// liquid density
   density(-0.3),	// mean density
   two(2.0),
   minus2(-2.0)
 {  
   Parameters::get(name + "->r",r);
-  Parameters::get(name + "->rho0", rho0);
   Parameters::get(name + "->density", density);
   Parameters::get(name + "->use mobility", useMobility);
   tempParameter= -(1.0+r);
@@ -34,57 +32,53 @@ void PhaseFieldCrystal::fillOperators()
   Operator *opL = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
   opL->addTerm(new Simple_SOT);
   
+  double M0 = 1.0;
+  Parameters::get(name + "->M0", M0);
   Operator *opLM = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
   if (useMobility) // non-constant mobility
     opLM->addTerm(new VecAtQP_SOT(
       prob->getSolution()->getDOFVector(0),
-      new MobilityPfc(density, rho0)));
+      new MobilityPfc(density, M0)));
   else // constant mobility
-    opLM->addTerm(new Simple_SOT);
+    opLM->addTerm(new Simple_SOT(M0));
   
   // -(1+r)*phi*rho
   Operator *opMTemp = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  opMTemp->addZeroOrderTerm(new Simple_ZOT());
+  opMTemp->addTerm(new Simple_ZOT());
   // -2*rho_old^3
   Operator *opMPowExpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0),new Pow3Functor(-2.0)));
+  opMPowExpl->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0),new AMDiS::Pow<3>(-2.0)));
   // -3*rho_old^2
   Operator *opMPowImpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  opMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0),new Pow2Functor(-3.0)));
+  opMPowImpl->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0),new AMDiS::Pow<2>(-3.0)));
   
   // 0 = nu-laplace(rho)
   // -------------------
-  prob->addMatrixOperator(opL, 0, 0); // -laplace(rho)
-  prob->addMatrixOperator(opMnew, 0, 2); // nu
+  prob->addMatrixOperator(opL, 0, 0); 						// -laplace(rho)
+  prob->addMatrixOperator(opMnew, 0, 2); 					// nu
   setAssembleMatrixOnlyOnce_butTimestepChange(0,0);
   setAssembleMatrixOnlyOnce_butTimestepChange(0,2);
   
   // dt(rho) = laplace(mu) - u*grad(rho)
   // -----------------------------------
   prob->addMatrixOperator(opMnew, 1, 0, getInvTau(), getInvTau());
-  prob->addMatrixOperator(opLM, 1, 1); // -laplace(mu)
+  prob->addMatrixOperator(opLM, 1, 1); 					// -laplace(mu)
   // . . . vectorOperators . . . . . . . . . . . . . . .
   prob->addVectorOperator(opMold, 1, getInvTau(), getInvTau());
-  setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
+//   setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
   if (!useMobility)
     setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
   
   // mu-2*nu-laplace(nu)-(1+r)*rho = (rho_old^3) + ExtPot - eps^2/(rho_old+0.9)
   // ----------------------------------------------------------------------
-  prob->addMatrixOperator(opMTemp, 2, 0, getTempParameter(), getTempParameter()); // -phi*(1+r)*rho
-  prob->addMatrixOperator(opMPowImpl, 2, 0); // -3*rho*rho_old^2
-  prob->addMatrixOperator(opL, 2, 0, &two, &two); // -2*phi*laplace(rho) * psi
-  prob->addMatrixOperator(opMnew, 2, 1); // phi*mu * psi
-  prob->addMatrixOperator(opL, 2, 2); // phi*grad(nu) * grad(psi)
-  // prob.addMatrixOperator(opMnew,2,2,&minus2);
+  prob->addMatrixOperator(opMTemp, 2, 0, &tempParameter, &tempParameter); 	// -phi*(1+r)*rho
+  prob->addMatrixOperator(opMPowImpl, 2, 0); 					// -3*rho*rho_old^2
+  prob->addMatrixOperator(opL, 2, 0, &two, &two); 				// -2*phi*laplace(rho) * psi
+  prob->addMatrixOperator(opMnew, 2, 1); 					// phi*mu * psi
+  prob->addMatrixOperator(opL, 2, 2); 						// phi*grad(nu) * grad(psi)
   // . . . vectorOperators . . . . . . . . . . . . . . .
-  prob->addVectorOperator(opMPowExpl, 2); // -2*phi^old*rho_old^3
+  prob->addVectorOperator(opMPowExpl, 2); 					// -2*phi^old*rho_old^3
   
   setAssembleMatrixOnlyOnce_butTimestepChange(2,1);
   setAssembleMatrixOnlyOnce_butTimestepChange(2,2);
 };
-
-
-void PhaseFieldCrystal::fillBoundaryConditions()
-{ FUNCNAME("PhaseFieldCrystal::fillBoundaryConditions()");
-};

extensions/base_problems/PhaseFieldCrystal_Base.h

@@ -21,10 +21,7 @@ public:
   PhaseFieldCrystal(const std::string &name_);
   ~PhaseFieldCrystal() {};
 
-  double *getTempParameter() { return &tempParameter; }
-
   virtual void fillOperators();
-  virtual void fillBoundaryConditions();
 
 protected:
 
@@ -74,34 +71,4 @@ class MobilityPfc : public AbstractFunction<double,double>
     double delta;
 };
 
-class Pow2Functor : public AbstractFunction<double,double>
-{
-  public:
-    Pow2Functor(double factor_=1.0) :
-      AbstractFunction<double,double>(2),
-      factor(factor_) { }
-    double operator()(const double &ch) const
-    {
-      return factor * sqr(ch);
-    }
-
-  protected:
-    double factor;
-};
-
-class Pow3Functor : public AbstractFunction<double,double>
-{
-  public:
-    Pow3Functor(double factor_=1.0) :
-      AbstractFunction<double,double>(3),
-      factor(factor_) { }
-    double operator()(const double &ch) const
-    {
-      return factor * sqr(ch) * ch;
-    }
-
-  protected:
-    double factor;
-};
-
 #endif // PHASE_FIELD_CRYSTAL_H

extensions/base_problems/PhaseFieldCrystal_Phase.h

@@ -15,7 +15,7 @@ class PhaseFieldCrystal_Phase : public PhaseFieldCrystal
 {
 public: // typedefs
 
-  typedef NavierStokes_Chorin super;
+  typedef PhaseFieldCrystal super;
 
 public: // methods