BaseProblems updated and some Warnings removed

extensions/base_problems/PhaseFieldCrystal_Base.cc

@@ -30,7 +30,7 @@ void PhaseFieldCrystal::fillOperators()
 
   // -nabla*(phi*grad(rho))
   Operator *opL = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  opL->addTerm(new Simple_SOT);
+  opL->addTerm(new CoordsAtQP_SOT(new AMDiS::Const<double, WorldVector<double> >(2.0)));
   
   double M0 = 1.0;
   Parameters::get(name + "->M0", M0);

extensions/base_problems/PhaseFieldCrystal_RB.cc

@@ -35,28 +35,35 @@ void PhaseFieldCrystal_RB::fillOperators()
 
   // -nabla*(grad(rho))
   Operator *opL = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  opL->addTerm(new Simple_SOT);
+  opL->addTerm(new CoordsAtQP_SOT(new AMDiS::Const<double, WorldVector<double> >(4.0)));
   Operator *opL2 = new Operator(prob->getFeSpace(2), prob->getFeSpace(0));
-  opL2->addTerm(new Simple_SOT(2.0));
+  opL2->addTerm(new CoordsAtQP_SOT(new AMDiS::Const<double, WorldVector<double> >(8.0)));
   Operator *opLnu = new Operator(prob->getFeSpace(2), prob->getFeSpace(2));
-  opLnu->addTerm(new Simple_SOT);
+  opLnu->addTerm(new CoordsAtQP_SOT(new AMDiS::Const<double, WorldVector<double> >(4.0)));
   
-  Operator *opLM = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  if (useMobility) // non-constant mobility
+  Operator *opLM = new Operator(prob->getFeSpace(1), prob->getFeSpace(1));  
+  Operator *opLMdiff = new Operator(prob->getFeSpace(1), prob->getFeSpace(1));
+  if (useMobility) {
+   // non-constant mobility
     opLM->addTerm(new VecAtQP_SOT(
-      prob->getSolution()->getDOFVector(0),
+      prob->getStageSolution(0),
       new MobilityPfc(density, rho0)));
-  else // constant mobility
+    opLMdiff->addTerm(new VecAtQP_SOT(
+      prob->getUnVec(0),
+      new MobilityPfcDiff(density, rho0)));
+  } else {
     opLM->addTerm(new Simple_SOT);
+    opLMdiff->addTerm(new Simple_SOT);
+  }
   
   // -rho_old^3
   Operator *opMPowExpl = new Operator(prob->getFeSpace(2), prob->getFeSpace(0));
   opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
-    prob->getStageSolution(0), new Pow3Functor(-1.0)));
+    prob->getStageSolution(0), new AMDiS::Pow<3>(-1.0)));
   // 3*rho_old^2
   Operator *opMPowImpl = new Operator(prob->getFeSpace(2), prob->getFeSpace(0));
   opMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
-    prob->getUnVec(0), new Pow2Functor(3.0)));
+    prob->getUnVec(0), new AMDiS::Pow<2>(3.0)));
   
   // 0 = nu-laplace(rho)
   // -------------------
@@ -74,7 +81,7 @@ void PhaseFieldCrystal_RB::fillOperators()
 
   opLM->setUhOld(prob->getStageSolution(1));
   prob->addVectorOperator(opLM, 1, &minus1); 	// -laplace(mu)
-  prob->addMatrixOperator(opLM, 1, 1); 		// -laplace(mu)
+  prob->addMatrixOperator(opLMdiff, 1, 1); 		// -laplace(mu)
   
   // mu-2*nu-laplace(nu)-(1+r)*rho = (rho_old^3) + ExtPot - eps^2/(rho_old+0.9)
   // ----------------------------------------------------------------------

extensions/base_problems/PhaseFieldCrystal_RB.h

@@ -4,18 +4,17 @@
 #define PHASE_FIELD_CRYSTAL_H
 
 #include "AMDiS.h"
-#include "BaseProblem.h"
-#include "time/RosenbrockStationary.h"
+#include "BaseProblem_RB.h"
 
 using namespace AMDiS;
 
 /** Phase-field Crystal problem
  */
-class PhaseFieldCrystal_RB : public BaseProblem<RosenbrockStationary>
+class PhaseFieldCrystal_RB : public BaseProblem_RB
 {
 public: // typedefs
 
-  typedef BaseProblem<RosenbrockStationary> super;
+  typedef BaseProblem_RB super;
 
 public:
 
@@ -36,23 +35,6 @@ protected:
   double two;
   double minus2;
   double minus1;
-
-// Rosenbrock methods
-//_________________________________________________________________________
-
-public:
-  
-  void acceptTimestep() { prob->acceptTimestep(); }
-  void reset() { prob->reset(); }
-  
-  void setRosenbrockMethod(RosenbrockMethod *method) { prob->setRosenbrockMethod(method); }
-  void setTau(double *ptr) { prob->setTau(ptr); }
-  void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3) { prob->setTauGamma(ptr0,ptr1,ptr2,ptr3); }
-
-  double* getTauGamma() { return prob->getTauGamma(); }
-  double* getMinusTauGamma() { return prob->getMinusTauGamma(); }
-  double* getInvTauGamma() { return prob->getInvTauGamma(); }
-  double* getMinusInvTauGamma() { return prob->getMinusInvTauGamma(); }
 };
 
 
@@ -77,48 +59,48 @@ class MobilityPfc : public AbstractFunction<double,double>
       AbstractFunction<double,double>(1),
       density(density_),
       factor(factor_),
-      delta(1.e-6) { }
+      delta(1.e-2) {
+	rho0 = 3.0*density;
+      }
 
     double operator()(const double &rho) const 
     {
-      double mobility= abs(rho - 3.0*density)*factor;
-      return std::max(mobility, delta);
+      double mobility= (rho - rho0)*factor;
+      return std::max(delta, mobility);
     }
   
   protected:
     double density;
     double factor;
     double delta;
+    
+    double rho0;
 };
 
-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>
+class MobilityPfcDiff : public AbstractFunction<double,double>
 {
   public:
-    Pow3Functor(double factor_=1.0) :
-      AbstractFunction<double,double>(3),
-      factor(factor_) { }
-    double operator()(const double &ch) const
+    MobilityPfcDiff(double density_ = -0.3, double factor_ = 1.0) : 
+      AbstractFunction<double,double>(0),
+      density(density_),
+      factor(factor_),
+      delta(1.e-2) {
+	rho0 = 3.0*density;	
+      }
+
+    double operator()(const double &rho) const 
     {
-      return factor * sqr(ch) * ch;
+      double mobility= (rho > rho0+delta ? factor : delta);
+      return mobility;
     }
-
+  
   protected:
+    double density;
     double factor;
+    double delta;
+    
+    double rho0;
 };
 
 #endif // PHASE_FIELD_CRYSTAL_H