some baseproblems extended

extensions/base_problems/CahnHilliard.cc

@@ -21,6 +21,7 @@
 
 #include "HL_SignedDistTraverse.h"
 #include "Recovery.h"
+#include "GenericZeroOrderTerm.h"
 
 using namespace AMDiS;
 
@@ -153,53 +154,67 @@ void CahnHilliard<P>::fillOperators()
 {  
   const FiniteElemSpace* feSpace = self::prob->self::getFeSpace();
   int degree = feSpace->getBasisFcts()->getDegree();
+  
+  DOFVector<double>* c = self::prob->getSolution()->getDOFVector(0);
+  DOFVector<double>* mu = self::prob->getSolution()->getDOFVector(1);
 
   // c
   Operator *opChMnew = new Operator(feSpace, feSpace);
-  opChMnew->addTerm(new Simple_ZOT);
+  addZOT(opChMnew, constant(1.0));
+//   opChMnew->addTerm(new Simple_ZOT);
   Operator *opChMold = new Operator(feSpace, feSpace);
-  opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
+  addZOT(opChMold, valueOf(c));
+//   opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
   // -nabla*(grad(c))
   Operator *opChL = new Operator(feSpace, feSpace);
-  opChL->addTerm(new Simple_SOT);
+  addSOT(opChL, constant(1.0) );
+//   opChL->addTerm(new Simple_SOT);
   
   // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
   Operator *opChLM = new Operator(feSpace, feSpace);
   if (useMobility) {
     if (doubleWell == 0)
-      opChLM->addTerm(new VecAtQP_SOT(
-	self::prob->getSolution()->getDOFVector(0),
-	new MobilityCH0(gamma, degree)));
+      addSOT(opChLM_, function_(MobilityCH0_(gamma), valueOf(c)) );
+//       opChLM->addTerm(new VecAtQP_SOT(
+// 	self::prob->getSolution()->getDOFVector(0),
+// 	new MobilityCH0(gamma, degree)));
     else
-      opChLM->addTerm(new VecAtQP_SOT(
-	self::prob->getSolution()->getDOFVector(0),
-	new MobilityCH1(gamma, degree)));
+      addSOT(opChLM, function_(MobilityCH1_(gamma), valueOf(c)) );
+//       opChLM->addTerm(new VecAtQP_SOT(
+// 	self::prob->getSolution()->getDOFVector(0),
+// 	new MobilityCH1(gamma, degree)));
   } else
-    opChLM->addTerm(new Simple_SOT(gamma));
+    addSOT(opChLM, constant(gamma) );
+//     opChLM->addTerm(new Simple_SOT(gamma));
   
   // -2*c_old^3 + 3/2*c_old^2
   Operator *opChMPowExpl = new Operator(feSpace, feSpace);
-  opChMPowExpl->addTerm(new VecAtQP_ZOT(
-    self::prob->getSolution()->getDOFVector(0),
-    new AMDiS::Pow<3>(-2.0, 3*degree)));
+  addZOT(opChMPowExpl, -2.0 * pow<3>(valueOf(c)));
+//   opChMPowExpl->addTerm(new VecAtQP_ZOT(
+//     self::prob->getSolution()->getDOFVector(0),
+//     new AMDiS::Pow<3>(-2.0, 3*degree)));
   if (doubleWell == 0) {
-    opChMPowExpl->addTerm(new VecAtQP_ZOT(
-      self::prob->getSolution()->getDOFVector(0),
-      new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
+    addZOT(opChMPowExpl, (3.0/2.0) * pow<2>(valueOf(c)));
+//     opChMPowExpl->addTerm(new VecAtQP_ZOT(
+//       self::prob->getSolution()->getDOFVector(0),
+//       new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
   }
 
   // -3*c_old^2 * c
   Operator *opChMPowImpl = new Operator(feSpace, feSpace);
-  opChMPowImpl->addTerm(new VecAtQP_ZOT(
-    self::prob->getSolution()->getDOFVector(0),
-    new AMDiS::Pow<2>(-3.0, 2*degree)));
+  addZOT(opChMPowImpl, -3.0 * pow<2>(valueOf(c)));
+//   opChMPowImpl->addTerm(new VecAtQP_ZOT(
+//     self::prob->getSolution()->getDOFVector(0),
+//     new AMDiS::Pow<2>(-3.0, 2*degree)));
   if (doubleWell == 0) {
-    opChMPowImpl->addTerm(new VecAtQP_ZOT(
-      self::prob->getSolution()->getDOFVector(0),
-      NULL, 3.0));
+    addZOT(opChMPowImpl, 3.0 * valueOf(c) - 0.5);
+//     opChMPowImpl->addTerm(new VecAtQP_ZOT(
+//       self::prob->getSolution()->getDOFVector(0),
+//       NULL, 3.0));
     opChMPowImpl->addTerm(new Simple_ZOT(-0.5));
   } else {
-    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
+    addZOT(opChMPowImpl, constant(1.0));
+//     opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
   }
 
   // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]

extensions/base_problems/CahnHilliard.h

@@ -21,6 +21,7 @@
 #include "AMDiS.h"
 #include "BaseProblem.h"
 #include "chns.h"
+#include "GenericOperatorTerm.h"
     
 namespace detail {
   
@@ -64,6 +65,51 @@ namespace detail {
     double epsSqr;
     double minusEpsSqr;
   };
+  
+  
+  class MobilityCH0_ : public FunctorBase //: public AbstractFunction<double,double>
+  {
+    public:
+      MobilityCH0_(double gamma_=1.0) : 
+	gamma(gamma_),
+	delta(1.e-6) { }
+	
+      typedef double value_type;
+      int getDegree(int d0) const { return d0 * 4; }
+  
+      double operator()(const double &ch) const 
+      {
+	double phase = std::max(0.0, std::min(1.0, ch));
+	double mobility = 0.25*sqr(phase)*sqr(phase-1.0);
+	return gamma * std::max(mobility, delta);
+      }
+
+    protected:
+      double gamma;
+      double delta;
+  };
+
+  class MobilityCH1_ : public FunctorBase // : public AbstractFunction<double,double>
+  {
+    public:
+      MobilityCH1_(double gamma_=1.0) :
+	gamma(gamma_),
+	delta(1.e-6) { }
+	
+      typedef double value_type;
+      int getDegree(int d0) const { return d0 * 4; }
+      
+      double operator()(const double &ch) const
+      {
+	double phase = std::max(-1.0, std::min(1.0, ch));
+	double mobility = 0.25*sqr(sqr(phase)-1.0);
+	return gamma * std::max(mobility, delta);
+      }
+
+    protected:
+      double gamma;
+      double delta;
+  };
 
 } // end namespace detail
 

extensions/base_problems/CahnHilliard.hh

@@ -153,55 +153,43 @@ void CahnHilliard<P>::fillOperators()
 {  
   const FiniteElemSpace* feSpace = self::prob->getFeSpace();
   int degree = feSpace->getBasisFcts()->getDegree();
+  
+  DOFVector<double>* c = self::prob->getSolution()->getDOFVector(0);
+  DOFVector<double>* mu = self::prob->getSolution()->getDOFVector(1);
 
   // c
   Operator *opChMnew = new Operator(feSpace, feSpace);
-  opChMnew->addTerm(new Simple_ZOT);
+  addZOT(opChMnew, constant(1.0));
   Operator *opChMold = new Operator(feSpace, feSpace);
-  opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
+  addZOT(opChMold, valueOf(c));
   // -nabla*(grad(c))
   Operator *opChL = new Operator(feSpace, feSpace);
-  opChL->addTerm(new Simple_SOT);
+  addSOT(opChL, constant(1.0) );
   
   // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
   Operator *opChLM = new Operator(feSpace, feSpace);
   if (useMobility) {
     if (doubleWell == 0)
-      opChLM->addTerm(new VecAtQP_SOT(
-	self::prob->getSolution()->getDOFVector(0),
-	new MobilityCH0(gamma, degree)));
+      addSOT(opChLM, function_(MobilityCH0_(gamma), valueOf(c)) );
     else
-      opChLM->addTerm(new VecAtQP_SOT(
-	self::prob->getSolution()->getDOFVector(0),
-	new MobilityCH1(gamma, degree)));
+      addSOT(opChLM, function_(MobilityCH1_(gamma), valueOf(c)) );
   } else
-    opChLM->addTerm(new Simple_SOT(gamma));
+    addSOT(opChLM, constant(gamma) );
   
   // -2*c_old^3 + 3/2*c_old^2
   Operator *opChMPowExpl = new Operator(feSpace, feSpace);
-  opChMPowExpl->addTerm(new VecAtQP_ZOT(
-    self::prob->getSolution()->getDOFVector(0),
-    new AMDiS::Pow<3>(-2.0, 3*degree)));
-  if (doubleWell == 0) {
-    opChMPowExpl->addTerm(new VecAtQP_ZOT(
-      self::prob->getSolution()->getDOFVector(0),
-      new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
-  }
+  addZOT(opChMPowExpl, -2.0 * pow<3>(valueOf(c)));
+  if (doubleWell == 0)
+    addZOT(opChMPowExpl, (3.0/2.0) * pow<2>(valueOf(c)));
 
   // -3*c_old^2 * c
   Operator *opChMPowImpl = new Operator(feSpace, feSpace);
-  opChMPowImpl->addTerm(new VecAtQP_ZOT(
-    self::prob->getSolution()->getDOFVector(0),
-    new AMDiS::Pow<2>(-3.0, 2*degree)));
-  if (doubleWell == 0) {
-    opChMPowImpl->addTerm(new VecAtQP_ZOT(
-      self::prob->getSolution()->getDOFVector(0),
-      NULL, 3.0));
-    opChMPowImpl->addTerm(new Simple_ZOT(-0.5));
-  } else {
-    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
-  }
-
+  addZOT(opChMPowImpl, -3.0 * pow<2>(valueOf(c)));
+  if (doubleWell == 0)
+    addZOT(opChMPowImpl, 3.0 * valueOf(c) - 0.5);
+  else
+    addZOT(opChMPowImpl, constant(1.0));
+  
   // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
   // ----------------------------------------------------------------------
   self::prob->addMatrixOperator(*opChMPowImpl,0,0);          /// < -3*phi*c*c_old^2 , psi >
@@ -217,7 +205,6 @@ void CahnHilliard<P>::fillOperators()
   self::prob->addMatrixOperator(*opChLM,1,1);                /// < phi*grad(mu) , grad(psi) >
   // . . . vectorOperators . . . . . . . . . . . . . . .
   self::prob->addVectorOperator(*opChMold,1, self::getInvTau());   /// < phi*c^old/tau , psi >
-  
 }
 
 

extensions/base_problems/CouplingBaseProblem.h

@@ -47,8 +47,8 @@ public:
   ~CouplingBaseProblem() { }  
   
   void initialize(Flag initFlag,
-			  ProblemStatSeq *adoptProblem = NULL,
-			  Flag adoptFlag = INIT_NOTHING) override
+		  ProblemStatSeq *adoptProblem = NULL,
+		  Flag adoptFlag = INIT_NOTHING) override
   {  
     for (size_t i = 0; i < baseProblems.size(); i++) {
       for (size_t j = 0; j < baseProblems[i]->getNumProblems(); j++)

extensions/base_problems/CouplingBaseProblem2.h

@@ -66,6 +66,30 @@ namespace detail {
   /// Functor for generic loops. Method fillBoundaryConditions() is called for each element in a sequence.
   struct FillBoundaryConditions { template<typename B> static void call(B& b) { b.fillBoundaryConditions(); } };
   
+  struct FindProblem {
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType& prob) {
+      if (baseProblem.getName() == name)
+	prob = baseProblem.getProblem();
+    }
+  };
+  
+  struct FindBaseProblem {
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType& prob) {
+      typedef typename boost::mpl::if_<typename boost::is_same<BaseProblemType, ProblemType>::type,
+				      boost::mpl::int_<true>,
+				      boost::mpl::bool_<false> >::type assign;
+      call(baseProblem, name, prob,  assign());
+    }
+    template<typename BaseProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, BaseProblemType& prob, boost::mpl::int_<true>) {
+      if (baseProblem.getName() == name)
+	prob = &baseProblem;
+    }
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, BaseProblemType& prob, boost::mpl::int_<false>) {}
+  };
 } // end namespace detail
 
 /**
@@ -156,7 +180,7 @@ public:
   template<int i>
   DOFVector<double> *getSolution(int j)
   { FUNCNAME("CouplingBaseProblem::getSolution<i>(j)");
-    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "BaseProblem-index out of range");
+    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "********** ERROR: BaseProblem-index out of range **********");
     
     TEST_EXIT(0 <= j && j <= _GET_<i>(baseProblems).getNumComponents())("Indices out of range!\n");
     return _GET_<i>(baseProblems).getSolution()->getDOFVector(j);
@@ -167,13 +191,35 @@ public:
   template<int i>
   inline const FiniteElemSpace* getFeSpace(int j=0) 
   { FUNCNAME("CouplingBaseProblem::getFeSpace<i>(j)");
-    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "BaseProblem index out of range");
+    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "********** ERROR: BaseProblem index out of range **********");
     
     TEST_EXIT(0 <= j && j <= _GET_<i>(baseProblems).getNumComponents())("Indices out of range!\n");
     return _GET_<i>(baseProblems).getFeSpace(j); 
   }
   
   std::string getName() { return name; }
+  
+
+  template<typename BaseProblemType>
+  BaseProblemType *getBaseProblem(std::string name_)
+  {
+    BaseProblemType *prob = NULL;
+    tools::FOR_EACH< detail::FindBaseProblem >::loop1(baseProblems, name_, prob);
+    if (prob)
+      return prob;
+    else
+      throw(std::runtime_error("problem with given name '" + name_ + "' does not exist"));
+  }
+
+  ProblemType *getProblem(std::string name_)
+  {
+    ProblemType *prob = NULL;
+    tools::FOR_EACH< detail::FindProblem >::loop1(baseProblems, name_, prob);
+    if (prob)
+      return prob;
+    else
+      throw(std::runtime_error("problem with given name '" + name_ + "' does not exist"));
+  }
 
 protected:
   BaseProblemsTupleType baseProblems;

extensions/base_problems/CouplingBaseProblem2_cxx11.h

@@ -68,6 +68,32 @@ namespace detail {
     static void call(BaseProblemType& baseProblem) { baseProblem.fillBoundaryConditions(); }
   };
   
+  
+  struct FindProblem {
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType& prob) {
+      if (baseProblem.getName() == name)
+	prob = baseProblem.getProblem();
+    }
+  };
+  
+  struct FindBaseProblem {
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType*& prob) {
+      typedef typename boost::mpl::if_<typename boost::is_same<BaseProblemType, ProblemType>::type,
+				      boost::mpl::bool_<true>,
+				      boost::mpl::bool_<false> >::type assign;
+      call(baseProblem, name, prob,  assign());
+    }
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType*& prob, boost::mpl::bool_<true>) {
+      if (baseProblem.getName() == name)
+	prob = &baseProblem;
+    }
+    template<typename BaseProblemType, typename ProblemType>
+    static void call(BaseProblemType& baseProblem, const std::string& name, ProblemType*& prob, boost::mpl::bool_<false>) {}
+  };
+  
 } // end namespace detail
 
 /**
@@ -151,7 +177,7 @@ public:
   template<int i>
   DOFVector<double> *getSolution(int j)
   { FUNCNAME("CouplingBaseProblem::getSolution<i>(j)");
-    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "BaseProblem-index out of range");
+    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "********** ERROR: BaseProblem-index out of range **********");
     
     TEST_EXIT(0 <= j && j <= _GET_<i>(baseProblems).getNumComponents())("Indices out of range!\n");
     return _GET_<i>(baseProblems).getSolution()->getDOFVector(j);
@@ -162,13 +188,35 @@ public:
   template<int i>
   inline const FiniteElemSpace* getFeSpace(int j=0) 
   { FUNCNAME("CouplingBaseProblem::getFeSpace<i>(j)");
-    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "BaseProblem index out of range");
+    BOOST_STATIC_ASSERT_MSG(0 <= i && i < _LENGTH_<BaseProblemsTupleType>::value , "********** ERROR: BaseProblem index out of range **********");
     
     TEST_EXIT(0 <= j && j <= _GET_<i>(baseProblems).getNumComponents())("Indices out of range!\n");
     return _GET_<i>(baseProblems).getFeSpace(j); 
   }
   
   std::string getName() { return name; }
+  
+
+  ProblemType *getProblem(std::string name_)
+  {
+    ProblemType *prob = NULL;
+    tools::FOR_EACH< detail::FindProblem >::loop1(baseProblems, name_, prob);
+    if (prob)
+      return prob;
+    else
+      throw(std::runtime_error("problem with given name '" + name_ + "' does not exist"));
+  }
+
+  template<typename BaseProblemType>
+  BaseProblemType *getBaseProblem(std::string name_)
+  {
+    BaseProblemType *prob = NULL;
+    tools::FOR_EACH< detail::FindBaseProblem >::loop1(baseProblems, name_, prob);
+    if (prob)
+      return prob;
+    else
+      throw(std::runtime_error("problem with given name '" + name_ + "' does not exist"));
+  }
 
 protected:
   BaseProblemsTupleType baseProblems;

extensions/base_problems/NavierStokesPhase_TaylorHood.h

@@ -19,6 +19,7 @@
 #define NAVIER_STOKES_PHASE_TAYLOR_HOOD_H
 
 #include "NavierStokes_TaylorHood.h"
+#include "ExtendedProblemStat.h"
 
 using namespace AMDiS;
 
@@ -29,11 +30,11 @@ using namespace AMDiS;
  * dt(phase*u_j) + phase*u*grad(u_j) = phase*d_j(p) - nu*div(phase*grad(u_j)) - beta/eps^alpha (1-phase)*(u_j - g_j); j=1...d
  * div(phase*u) = grad(phase)*g
  */
-class NavierStokesPhase_TaylorHood : public NavierStokes_TaylorHood
+class NavierStokesPhase_TaylorHood : public ::detail::NavierStokes_TaylorHood<ExtendedProblemStat>
 {
 public: // typedefs
 
-  typedef NavierStokes_TaylorHood super;
+  typedef ::detail::NavierStokes_TaylorHood<ExtendedProblemStat> super;
 
 public: // methods
 

extensions/base_problems/NavierStokes_TH_MultiPhase.h

@@ -19,147 +19,158 @@
 
 #include "NavierStokes_TaylorHood.h"
 
+
+namespace detail {
+  
 using namespace AMDiS;
 
-/** \ingroup NavierStokes_TaylorHood
-* \brief
-* Navier-Stokes multi-phase problem, using Taylor Hood elements
-* Standard implementation for two phases, but can be extended to
-* many more, by defining the multiPhase variable and overloading
-* the initTimestep, where densityPhase and viscosityPhase are defined
-* depending on the phases
-*/
-class NavierStokes_TH_MultiPhase : public NavierStokes_TaylorHood
-{
-public: // typedefs
+  /** \ingroup NavierStokes_TaylorHood
+  * \brief
+  * Navier-Stokes multi-phase problem, using Taylor Hood elements
+  * Standard implementation for two phases, but can be extended to
+  * many more, by defining the multiPhase variable and overloading
+  * the initTimestep, where densityPhase and viscosityPhase are defined
+  * depending on the phases
+  */
+  template<typename ProblemType>
+  class NavierStokes_TH_MultiPhase : public NavierStokes_TaylorHood<ProblemType>
+  {
+  public: // typedefs
 
-  typedef NavierStokes_TaylorHood super;
+    typedef NavierStokes_TaylorHood<ProblemType> super;
+    typedef NavierStokes_TH_MultiPhase<ProblemType> self;