a lot of changes

extensions/POperators_FOT.h

@@ -180,6 +180,45 @@ double fac;
 int numVecs;
 };
 
+
+
+/* -------------------------------------------------------------- */
+
+/// < f(vec0)*vecs*grad(u) , psi > or < f(vec0)*vecs*u , grad(psi) > where vecs in R^n
+struct VecAndWorldVec_FOT : FirstOrderTerm
+{
+  VecAndWorldVec_FOT(DOFVector<double>* vec0, 
+					 WorldVector<DOFVector<double>*> vecs, 
+					 AbstractFunction<double, double>* fct, 
+					 double fac = 1.0);
+
+  void initElement(const ElInfo*, 
+				   SubAssembler*,
+				   Quadrature *quad= NULL);
+  
+  void initElement(const ElInfo* largeElInfo, 
+				   const ElInfo* smallElInfo,
+				   SubAssembler*,
+				   Quadrature *quad = NULL);
+
+  void getLb(const ElInfo *elInfo,
+			 vector<mtl::dense_vector<double> >& result) const;
+
+  void eval(int nPoints,
+			const mtl::dense_vector<double>&,
+			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+			mtl::dense_vector<double>& result,
+			double opFactor);
+
+private:
+  DOFVector<double> *vDV, *vec0DV, *vec1DV, *vec2DV;
+  mtl::dense_vector<double> v, vec0, vec1, vec2;
+  AbstractFunction<double, double>* fct;
+  double fac;
+  int numVecs;
+};
+
 /* -------------------------------------------------------------- */
 
 /// < factor*d_i(u) , psi > or < factor*u , d_i(psi) >
@@ -242,6 +281,96 @@ mtl::dense_vector<double> vec1;
 double fac;
 int component;
 };
+
+/* -------------------------------------------------------------- */
+
+/// < d_i(u) * factor*v1*d_j(v2) , psi > or < u * factor*v1*d_j(v2) , d_i(psi) >
+struct VecAndPartialDerivativeIJ_FOT : FirstOrderTerm
+{
+  VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
+				DOFVectorBase<double> *dv2,
+				int i_,
+				int j_,
+				double fac_ = 1.0);
+
+  VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
+				DOFVectorBase<double> *dv2,
+				int i_,
+				int j_,
+				BinaryAbstractFunction<double, double, double>* fct, // = f(v1, d_j(v2))
+				double fac_ = 1.0);
+
+  void initElement(const ElInfo* elInfo,
+	SubAssembler* subAssembler,
+	Quadrature *quad = NULL);
+
+  void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+	SubAssembler* subAssembler,
+	Quadrature *quad = NULL);
+
+  void getLb(const ElInfo *elInfo, vector<mtl::dense_vector<double> >& result) const;
+
+  void eval(int nPoints,
+	const mtl::dense_vector<double>&,
+	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	mtl::dense_vector<double>& result,
+	double opFactor);
+
+private:
+  DOFVectorBase<double>* vec1DV;
+  DOFVectorBase<double>* vec2DV;
+  mtl::dense_vector<double> vec1;
+  mtl::dense_vector<WorldVector<double> > grad;
+  BinaryAbstractFunction<double, double, double>* fct;
+  double fac;
+  int xi, xj;
+};
+
+
+/* -------------------------------------------------------------- */
+
+/// < factor*v1*v2*d_i(u) , psi > or < factor*v1*v2*u , d_i(psi) >
+struct Vec2AndPartialDerivative_FOT : public FirstOrderTerm
+{
+  Vec2AndPartialDerivative_FOT(DOFVectorBase<double> *dv1, 
+			       DOFVectorBase<double> *dv2,
+			       int component_, 
+			       double fac_=1.0);
+  
+  Vec2AndPartialDerivative_FOT(DOFVectorBase<double> *dv1, 
+			       DOFVectorBase<double> *dv2,
+			       int component_,
+			       BinaryAbstractFunction<double, double, double>* fct, // = f(v1, v2)
+			       double fac_=1.0);
+
+  void initElement(const ElInfo* elInfo,
+	  SubAssembler* subAssembler,
+	  Quadrature *quad = NULL);
+  
+  void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+	  SubAssembler* subAssembler,
+	  Quadrature *quad = NULL);
+
+  void getLb(const ElInfo *elInfo, vector<mtl::dense_vector<double> >& result) const;
+
+  void eval(int nPoints,
+	  const mtl::dense_vector<double>&,
+	  const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	  const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	  mtl::dense_vector<double>& result,
+	  double opFactor);
+
+private:
+  DOFVectorBase<double>* vec1DV;
+  DOFVectorBase<double>* vec2DV;
+  mtl::dense_vector<double> vec1,vec2;
+  BinaryAbstractFunction<double, double, double>* fct;
+  double fac;
+  int component;
+};
+
+
 /* -------------------------------------------------------------- */
 
 

extensions/POperators_ZOT.h

@@ -336,65 +336,87 @@ int component, component2;
 /* -------------------------------------------------------------- */
 
 /// < factor*v*d_i(w)*u , psi >
-class VecAndPartialDerivative_ZOT : public ZeroOrderTerm
+struct VecAndPartialDerivative_ZOT : ZeroOrderTerm
 {
-public:
-VecAndPartialDerivative_ZOT(DOFVectorBase<double>* vecDV_, DOFVectorBase<double>* gradDV_, int component_, double fac_=1.0);
-
-void initElement(const ElInfo*, 
-	SubAssembler*,
-	Quadrature *quad= NULL);
-void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
-	SubAssembler*,
-	Quadrature *quad = NULL);
-
-void getC(const ElInfo *, int nPoints, ElementVector &C);
-
-void eval(int nPoints,
-		const mtl::dense_vector<double>&,
-		const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
-		const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
-		mtl::dense_vector<double>& result,
-		double opFactor);
-
-protected:
-DOFVectorBase<double> *vecDV, *gradDV;
-mtl::dense_vector<double> vec;
-double fac;
-mtl::dense_vector<WorldVector<double> > grad;
-int component;
+  VecAndPartialDerivative_ZOT(DOFVectorBase<double>* vecDV_, 
+			      DOFVectorBase<double>* gradDV_, 
+			      int component_, 
+			      double fac_ = 1.0);
+  
+  VecAndPartialDerivative_ZOT(DOFVectorBase<double>* vecDV_, 
+			      DOFVectorBase<double>* gradDV_, 
+			      int component_, 
+			      AbstractFunction<double, double>* fct, // f(v)
+			      double fac_ = 1.0);
+
+  void initElement(const ElInfo*, 
+	  SubAssembler*,
+	  Quadrature *quad= NULL);
+  
+  void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+	  SubAssembler*,
+	  Quadrature *quad = NULL);
+
+  void getC(const ElInfo *, int nPoints, ElementVector &C);
+
+  void eval(int nPoints,
+		  const mtl::dense_vector<double>&,
+		  const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+		  const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+		  mtl::dense_vector<double>& result,
+		  double opFactor);
+
+private:
+  DOFVectorBase<double> *vecDV, *gradDV;
+  mtl::dense_vector<double> vec;
+  AbstractFunction<double, double>* fct;
+  double fac;
+  mtl::dense_vector<WorldVector<double> > grad;
+  int component;
 };
 
 /* -------------------------------------------------------------- */
 
 /// < factor*v1*v2*d_i(w)*u , psi >
-class Vec2AndPartialDerivative_ZOT : public ZeroOrderTerm
+struct Vec2AndPartialDerivative_ZOT : ZeroOrderTerm
 {
-public:
-Vec2AndPartialDerivative_ZOT(DOFVectorBase<double>* vec1DV_, DOFVectorBase<double>* vec2DV_, DOFVectorBase<double>* gradDV_, int component_, double fac_=1.0);
-
-void initElement(const ElInfo*,
-	SubAssembler*,
-	Quadrature *quad= NULL);
-void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
-	SubAssembler*,
-	Quadrature *quad = NULL);
-
-void getC(const ElInfo *, int nPoints, ElementVector &C);
-
-void eval(int nPoints,
-		const mtl::dense_vector<double>&,
-		const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
-		const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
-		mtl::dense_vector<double>& result,
-		double opFactor);
-
-protected:
-DOFVectorBase<double> *vec1DV, *vec2DV, *gradDV;
-mtl::dense_vector<double> vec1, vec2;
-double fac;
-mtl::dense_vector<WorldVector<double> > grad;
-int component;
+  Vec2AndPartialDerivative_ZOT(DOFVectorBase<double>* vec1DV_, 
+							  DOFVectorBase<double>* vec2DV_, 
+							  DOFVectorBase<double>* gradDV_, 
+							  int component_, 
+							  double fac_ = 1.0);
+  
+  Vec2AndPartialDerivative_ZOT(DOFVectorBase<double>* vec1DV_, 
+							  DOFVectorBase<double>* vec2DV_, 
+							  DOFVectorBase<double>* gradDV_, 
+							  int component_, 
+							  BinaryAbstractFunction<double, double, double>* fct,
+							  double fac_ = 1.0);
+
+  void initElement(const ElInfo*,
+	  SubAssembler*,
+	  Quadrature *quad= NULL);
+
+  void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+	  SubAssembler*,
+	  Quadrature *quad = NULL);
+
+  void getC(const ElInfo *, int nPoints, ElementVector &C);
+
+  void eval(int nPoints,
+		  const mtl::dense_vector<double>&,
+		  const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+		  const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+		  mtl::dense_vector<double>& result,
+		  double opFactor);
+
+private:
+  DOFVectorBase<double> *vec1DV, *vec2DV, *gradDV;
+  mtl::dense_vector<double> vec1, vec2;
+  BinaryAbstractFunction<double, double, double>* fct;
+  double fac;
+  mtl::dense_vector<WorldVector<double> > grad;
+  int component;
 };
 
 /* -------------------------------------------------------------- */

extensions/SignedDistFunctors.h

@@ -247,17 +247,17 @@ class Plane : public AbstractFunction<double, WorldVector<double> >
 {
 public:
 
-  Plane(double shift_, double factor_=1.0) :
-    shift(shift_), factor(factor_) {}
+  Plane(double shift_, double factor_=1.0, int component_=1) :
+    shift(shift_), factor(factor_), component(component_) {}
   double operator()(const WorldVector<double>& x) const
   {
-    return factor*(shift-x[1]);
+    return factor*(shift-x[component]);
   };
 
 private:
 
   double shift;
-  int comp;
+  int component;
   double factor;
 };
 

extensions/base_problems/BaseProblem_RB.h

@@ -30,7 +30,7 @@ public:
   
   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); }
+  void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3, double *ptr4) { prob->setTauGamma(ptr0,ptr1,ptr2,ptr3,ptr4); }
 
   double* getTauGamma() { return prob->getTauGamma(); }
   double* getMinusTauGamma() { return prob->getMinusTauGamma(); }

extensions/base_problems/CahnHilliardNavierStokes.cc

@@ -59,7 +59,7 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
   // Parameters for NS-Coupling
   // cahn-hiliard-force: sigma*mu*grad(c)
   Initfile::get(name + "->sigma", sigma); // surface tension
-  surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
+  surfaceTension = -sigma*3.0/(2.0*sqrt(2.0)) / eps;
 
     
   force.set(0.0);

extensions/base_problems/CahnHilliardNavierStokes_RB.cc

@@ -59,7 +59,7 @@ CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name
   // Parameters for NS-Coupling
   // cahn-hiliard-force: sigma*mu*grad(c)
   Initfile::get(name + "->sigma", sigma); // surface tension
-  surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
+  surfaceTension = -sigma*3.0/(2.0*sqrt(2.0)) / eps;
 
     
   force.set(0.0);
@@ -148,35 +148,43 @@ void CahnHilliardNavierStokes_RB::fillOperators()
       opChLM->addTerm(new VecAtQP_SOT(
 	prob->getUnVec(0),
 	new MobilityCH0(gamma)));
-      opChLM->addTerm(new VecGrad_FOT(
+      prob->addMatrixOperator(*opChLM, 0, 1);
+      
+      Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+      opChLM2->addTerm(new VecGrad_FOT(
 	prob->getUnVec(0),
 	prob->getUnVec(1),
 	new MobilityCH0Diff(gamma)), GRD_PSI);
-      prob->addMatrixOperator(*opChLM, 0, 1);
+      prob->addMatrixOperator(*opChLM2, 0, 0);
+      
       // rhs operators
-      Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
-      opChLM2->addTerm(new VecAtQP_SOT(
+      Operator *opChLMe = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+      opChLMe->addTerm(new VecAtQP_SOT(
 	prob->getStageSolution(0),
 	new MobilityCH0(gamma)));
-      opChLM2->setUhOld(prob->getStageSolution(1));
-      prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
-    } else {
+      opChLMe->setUhOld(prob->getStageSolution(1));
+      prob->addVectorOperator(*opChLMe, 0, &minus1, &minus1);
+    } else if (doubleWell == 1) {
       // jacobian operators
       opChLM->addTerm(new VecAtQP_SOT(
 	prob->getUnVec(0),
 	new MobilityCH1(gamma)));
-      opChLM->addTerm(new VecGrad_FOT(
+      prob->addMatrixOperator(*opChLM, 0, 1);
+      
+      Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+      opChLM2->addTerm(new VecGrad_FOT(
 	prob->getUnVec(0),
 	prob->getUnVec(1),
 	new MobilityCH1Diff(gamma)), GRD_PSI);
-      prob->addMatrixOperator(*opChLM, 0, 1);
+      prob->addMatrixOperator(*opChLM2, 0, 0);
+      
       // rhs operators
-      Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
-      opChLM2->addTerm(new VecAtQP_SOT(
+      Operator *opChLMe = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+      opChLMe->addTerm(new VecAtQP_SOT(
 	prob->getStageSolution(0),
 	new MobilityCH1(gamma)));
-      opChLM2->setUhOld(prob->getStageSolution(1));
-      prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
+      opChLMe->setUhOld(prob->getStageSolution(1));
+      prob->addVectorOperator(*opChLMe, 0, &minus1, &minus1);
     }
   } else {
     opChLM->addTerm(new Simple_SOT(gamma));
@@ -319,7 +327,7 @@ void CahnHilliardNavierStokes_RB::fillOperators()
 	/// < mu_s * grad(c_s) , psi >
     Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(1));
     opMuGradC->addTerm(new VecAndPartialDerivative_ZOT(
-	  prob->getStageSolution(1),
+      prob->getStageSolution(1),
       prob->getStageSolution(0), i));
     prob->addVectorOperator(*opMuGradC, 2+i, &surfaceTension, &surfaceTension);
   }
@@ -347,8 +355,8 @@ void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
         Operator *opLaplaceUi1 = new Operator(getFeSpace(2+i), getFeSpace(2+j));
         opLaplaceUi1->addTerm(new MatrixIJ_SOT(1-i, 1-j, viscosity));
         prob->addMatrixOperator(*opLaplaceUi1, 2+i, 2+j);
-		opLaplaceUi1->setUhOld(prob->getStageSolution(2+j));
-		prob->addVectorOperator(*opLaplaceUi1, 2+i, &minus1, &minus1);
+	opLaplaceUi1->setUhOld(prob->getStageSolution(2+j));
+	prob->addVectorOperator(*opLaplaceUi1, 2+i, &minus1, &minus1);
       }
     }
     

extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h

+/** \file CahnHilliardNavierStokes_TwoPhase.h */
+
+#ifndef CAHN_HILLIARD_H
+#define CAHN_HILLIARD_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "ExtendedProblemStat.h"
+#include "HL_SignedDistTraverse.h"
+#include "chns.h"
+
+using namespace AMDiS;
+
+class CahnHilliardNavierStokes_TwoPhase : public BaseProblem<ExtendedProblemStat>
+{
+public: // definition of types
+
+  typedef BaseProblem<ExtendedProblemStat> super;
+
+public: // public methods
+
+  CahnHilliardNavierStokes_TwoPhase(const std::string &name_);
+  ~CahnHilliardNavierStokes_TwoPhase();
+
+  virtual void initData();
+
+  virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
+  virtual void closeTimestep(AdaptInfo *adaptInfo);
+  
+  double getEpsilon() { return eps; }
+  int getDoubleWellType() { return doubleWell; }
+
+  DOFVector<WorldVector<double> >* getVelocity()
+  {
+    return velocity;
+  }
+  
+  virtual void fillOperators();
+  virtual void fillBoundaryConditions() {}
+  
+  virtual void addMaterialDerivativeConst(int i);
+  virtual void addMaterialDerivativeNonConst(int i);
+  virtual void addLaplaceTermConst(int i);
+  virtual void addLaplaceTermNonConst(int i);
+
+protected: // protected variables
+
+  HL_SignedDistTraverse *reinit;
+
+  DOFVector<WorldVector<double> >* velocity;
+
+  void calcVelocity()
+  {
+    if (dow == 1)
+      transformDOF(prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec1WorldVec<double>);
+    else if (dow == 2)
+      transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), velocity, new AMDiS::Vec2WorldVec<double>);
+    else if (dow == 3)
+      transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
+  }
+  
+  bool useMobility;
+  bool useConservationForm;
+
+  unsigned dim;
+  
+  int laplaceType;
+  int nonLinTerm;
+  int doubleWell;
+
+  // navierStokes parameters
+  double oldTimestep;
+  double viscosity1;
+  double viscosity2;
+  double density1;
+  double density2;
+  double c;
+  
+  // Cahn-Hilliard parameters
+  double gamma;
+  double eps;
+  double minusEps;
+  double epsInv;
+  double minusEpsInv;
+  double epsSqr;
+  double minusEpsSqr;
+  
+  double sigma;		// coupling parameter to calculate the surface tension
+  double surfaceTension;// := sigma/epsilon
+  
+  WorldVector<double> force;
+
+  FileVectorWriter *fileWriter;
+};
+
+#endif // CAHN_HILLIARD_H

extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h

@@ -6,6 +6,7 @@
 #include "AMDiS.h"
 #include "BaseProblem_RB.h"
 #include "HL_SignedDistTraverse.h"
+#include "chns.h"
 
 using namespace AMDiS;
 
@@ -89,238 +90,4 @@ protected: // protected variables
 };
 
 
-/** \brief
- * Abstract function for Cahn-Hilliard mobility
- */
-class MobilityCH0 : public AbstractFunction<double,double>
-{
-  public:
-    MobilityCH0(double gamma_=1.0) : 
-      AbstractFunction<double,double>(4), 
-      gamma(gamma_),
-      delta(1.e-6) { }
-    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 AbstractFunction<double,double>
-{
-  public:
-    MobilityCH1(double gamma_=1.0) :
-      AbstractFunction<double,double>(4),
-      gamma(gamma_),
-      delta(1.e-6) { }
-    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;
-};
-
-class ScaleFunctor : public AbstractFunction<double,double>
-{
-  public:
-    ScaleFunctor(double factor_=1.0) :
-      AbstractFunction<double,double>(1),
-      factor(factor_) { }
-    double operator()(const double &ch) const
-    {
-      return factor * ch;
-    }
-
-  protected:
-    double factor;
-};
-
-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;
-};