added a coupled NavierStokesCahnHilliard preconditioner and demo

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -509,7 +509,28 @@ namespace AMDiS {
     PCFieldSplitSetIS(pc, splitName, is);
     ISDestroy(&is);
   }
+  
+  
+  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
+  {
+    FUNCNAME("PetscSolverGlobalMatrix::extractVectorComponent()");
+    IS is;
+    interiorMap->createIndexSet(is, i, numberOfComponents);
+    VecGetSubVector(input, is, output);
+    ISDestroy(&is);
+  }
 
+  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
+  {
+    FUNCNAME("PetscSolverGlobalMatrix::extractMatrixComponent()");
+    IS isrow, iscol;
+    interiorMap->createIndexSet(isrow, startRow, numberOfRows);
+    interiorMap->createIndexSet(iscol, startCol, numberOfCols);
+    MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
+    ISDestroy(&iscol);
+    ISDestroy(&isrow);
+  }
+  
 
   void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
   {

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -62,6 +62,10 @@ namespace AMDiS {
     void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);    
 
     void solveGlobal(Vec &rhs, Vec &sol);
+    
+    void extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents=1);
+    
+    void extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output);
 
     void destroyMatrixData();
 
@@ -96,6 +100,8 @@ namespace AMDiS {
       components.push_back(component);
       createFieldSplit(pc, splitName, components);
     }
+    
+    
 
     virtual void initSolver(KSP &ksp);
 

extensions/base_problems/NavierStokesCahnHilliard.cc

+#include "NavierStokesCahnHilliard.h"
+// #include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+using namespace AMDiS;
+
+NavierStokesCahnHilliard::NavierStokesCahnHilliard(const std::string &name_, bool createProblem) :
+  super(name_, createProblem),
+  useMobility(false),
+  doubleWell(0),
+  gamma(1.0),
+  eps(0.1),
+  minusEps(-0.1),
+  epsInv(10.0),
+  minusEpsInv(-10.0),
+  epsSqr(0.01),
+  minusEpsSqr(-0.01),
+  multiPhase(NULL),
+  densityPhase(NULL),
+  viscosityPhase(NULL),
+  viscosity1(1.0),
+  viscosity2(1.0),
+  density1(1.0),
+  density2(1.0),    
+  refFunction(NULL),
+  refinement(NULL),
+  sigma(0.0),
+  surfaceTension(0.0)
+  {
+    dow = Global::getGeo(WORLD);
+  Initfile::get(name + "->viscosity1", viscosity1); // viscosity of fluid 1
+  Initfile::get(name + "->viscosity2", viscosity2); // viscosity of fluid 2
+  Initfile::get(name + "->density1", density1); // density of fluid 1
+  Initfile::get(name + "->density2", density2); // density of fluid 2
+
+  Initfile::get(name + "->non-linear term", nonLinTerm); 
+  Initfile::get(name + "->theta", theta); 
+  Initfile::get("adapt->timestep", tau); 
+    
+    
+    Parameters::get(name + "->epsilon", eps); // interface width
+    Parameters::get(name + "->sigma", sigma); 
+    double a_fac;
+    Parameters::get(name + "->a_factor", a_fac);
+    a=1.0/sigma*sqr(eps)*a_fac* (std::max(density1,density2)/tau);
+    b=1.0;
+    ab=a*b;   
+    surfaceTension = -sigma/eps * 3.0/2.0/sqrt(2.0) * a ;
+
+  theta1 = 1.0-theta;
+  minusTheta1 = -theta1;
+
+  force.set(0.0);
+  Initfile::get(name + "->force", force);
+  
+  // parameters for CH
+  Parameters::get(name + "->use mobility", useMobility); // mobility
+  Parameters::get(name + "->gamma", gamma); // mobility
+
+  
+  // type of double well: 0= [0,1], 1= [-1,1]
+  Parameters::get(name + "->double-well type", doubleWell); 
+  
+  // transformation of the parameters
+  minusEps = -eps;
+  minus1 = -1.0;
+  epsInv = 1.0/eps;
+  minusEpsInv = -epsInv;
+  epsSqr = sqr(eps);
+  minusEpsSqr = -epsSqr*b;
+  
+  
+}
+
+
+NavierStokesCahnHilliard::~NavierStokesCahnHilliard() 
+{ FUNCNAME("NavierStokesCahnHilliard::~NavierStokesCahnHilliard()");   
+  delete densityPhase;
+  delete viscosityPhase;
+}
+
+
+void NavierStokesCahnHilliard::setTime(AdaptInfo* adaptInfo) 
+  {
+    super::setTime(adaptInfo);
+    delta = gamma*adaptInfo->getTimestep();
+  }
+  
+  
+
+
+void NavierStokesCahnHilliard::initData()
+{ FUNCNAME("NavierStokes_TH_MultiPhase::initTimeInterface()");
+
+  #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    string initFileStr = name + "->space->solver", solverType = "";
+    Parameters::get(initFileStr, solverType);
+    if (solverType == "petsc-nsch") {
+      PetscSolverNSCH* solver = dynamic_cast<PetscSolverNSCH*>(prob->getSolver()); 
+      if (solver)
+      {
+	solver->setChData(&eps, &delta);
+	solver->setStokesData( getInvTau(), prob->getSolution(), &viscosity1, &viscosity2, &density1, &density2);
+      }
+    }
+    #endif
+
+  densityPhase = new DOFVector<double>(prob->getFeSpace(0), "densityPhase");
+  viscosityPhase = new DOFVector<double>(prob->getFeSpace(0), "viscosityPhase");
+
+  densityPhase->set(density1);
+  viscosityPhase->set(viscosity1);
+  
+   if (velocity == NULL)
+    velocity = new DOFVector<WorldVector<double> >(getFeSpace(0), "velocity");
+ 
+  //fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
+
+  super::initData();
+};
+
+
+void NavierStokesCahnHilliard::closeTimestep(AdaptInfo *adaptInfo)
+{ super::closeTimestep(adaptInfo);
+
+}
+
+void NavierStokesCahnHilliard::initTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokes_TH_MultiPhase::initTimestep()");
+  
+  super::initTimestep(adaptInfo);
+  refinement->refine(2);
+  LinearInterpolation1 dLI(density1, density2);
+  LinearInterpolation1 vLI(viscosity1, viscosity2);
+  transformDOF(prob->getSolution()->getDOFVector(dow+2), densityPhase, &dLI);
+  transformDOF(prob->getSolution()->getDOFVector(dow+2), viscosityPhase, &vLI);
+};
+
+
+
+void NavierStokesCahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo)
+  {
+    // meshFunction for klocal refinement around the interface of the phasefield-DOFVector
+    refFunction = new PhaseFieldRefinement(prob->getMesh());
+
+    if (getDoubleWellType() == 0) {
+      refinement = new RefinementLevelDOF(
+	prob->getFeSpace(),
+	refFunction,
+	prob->getSolution()->getDOFVector(dow+2)); // phaseField-DOFVector in [0,1]
+    } else {
+      refinement = new RefinementLevelDOF(
+	prob->getFeSpace(),
+	refFunction,
+	new PhaseDOFView<double>(prob->getSolution()->getDOFVector(dow+2))); // phaseField-DOFVector in [-1,1]
+    }
+
+    // initial refinement
+    refinement->refine(0);
+
+    for (int i = 0; i < 3; i++) {
+      solveInitialProblem2(adaptInfo); 	// update phaseField-DOFVector
+      refinement->refine((i < 4 ? 4 : 10));	// do k steps of local refinement
+    }
+
+    // solve all initial problems of the problems added to the CouplingTimeInterface    
+  }
+
+
+
+void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo) 
+{ FUNCNAME("NavierStokesCahnHilliard::solveInitialProblem()");
+
+  Flag initFlag = initDataFromFile(adaptInfo);
+
+  if (!initFlag.isSet(DATA_ADOPTED)) {
+    int initialInterface = 0;
+    Initfile::get(name + "->initial interface", initialInterface);
+    double initialEps = eps;
+    Initfile::get(name + "->initial epsilon", initialEps);
+
+    if (initialInterface == 0) {
+      /// horizontale Linie
+      double a= 0.0, dir= -1.0;
+      Initfile::get(name + "->line->pos", a);
+      Initfile::get(name + "->line->direction", dir);
+      prob->getSolution()->getDOFVector(1+3)->interpol(new Plane(a, dir));
+    }
+    else if (initialInterface == 1) {
+      /// schraege Linie
+      double theta = m_pi/4.0;
+      prob->getSolution()->getDOFVector(1+3)->interpol(new PlaneRotation(0.0, theta, 1.0));
+      transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
+    }
+    else if (initialInterface == 2) {
+      /// Ellipse
+      double a= 1.0, b= 1.0;
+      Initfile::get(name + "->ellipse->a", a);
+      Initfile::get(name + "->ellipse->b", b);
+      prob->getSolution()->getDOFVector(1+3)->interpol(new Ellipse(a,b));
+    }
+    else if (initialInterface == 3) {
+      /// zwei horizontale Linien
+      double a= 0.75, b= 0.375;
+      Initfile::get(name + "->lines->pos1", a);
+      Initfile::get(name + "->lines->pos2", b);
+      prob->getSolution()->getDOFVector(1+3)->interpol(new Plane(a, -1.0));
+      transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),new Plane(b, 1.0), new AMDiS::Max<double>);
+    }
+    else if (initialInterface == 4) {
+      /// Kreis
+      double radius= 1.0, x=0, y=0;
+      Initfile::get(name + "->circle->radius", radius);
+      Initfile::get(name + "->circle->center_x", x);
+      Initfile::get(name + "->circle->center_y", y);
+      WorldVector<double> w;
+      w[0]=x; w[1]=y+adaptInfo->getTime();
+      prob->getSolution()->getDOFVector(1+3)->interpol(new Circle(radius, w));
+    } else if (initialInterface == 5) {
+      /// Rechteck
+      double width = 0.5;
+      double height = 0.3;
+      WorldVector<double> center; center.set(0.5);
+      Initfile::get(name + "->rectangle->width", width);
+      Initfile::get(name + "->rectangle->height", height);
+      Initfile::get(name + "->rectangle->center", center);
+      prob->getSolution()->getDOFVector(1+3)->interpol(new Rectangle(width, height, center));
+    }
+
+    /// create phase-field from signed-dist-function
+    if (doubleWell == 0) {
+      transformDOF(prob->getSolution()->getDOFVector(1+3),
+        new SignedDistToPhaseField(initialEps));
+    } else {
+      transformDOF(prob->getSolution()->getDOFVector(1+3),
+        new SignedDistToCh(initialEps));
+    }
+  }
+}
+
+
+void NavierStokesCahnHilliard::fillOperators()
+{ FUNCNAME("NavierStokesCahnHilliard::fillOperators()");
+  MSG("NavierStokesCahnHilliard::fillOperators()\n");
+
+  const FiniteElemSpace* feSpace = prob->getFeSpace(0);
+ 
+  // c
+  Operator *opChMnew = new Operator(feSpace,feSpace);
+  opChMnew->addZeroOrderTerm(new Simple_ZOT);
+  Operator *opChMold = new Operator(feSpace,feSpace);
+  opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1+3)));
+  // -nabla*(grad(c))
+  Operator *opChL = new Operator(feSpace,feSpace);
+  opChL->addSecondOrderTerm(new Simple_SOT);
+  
+  // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
+  Operator *opChLM = new Operator(feSpace,feSpace);
+  opChLM->addSecondOrderTerm(new Simple_SOT(gamma*ab));
+  
+  // -2*c_old^3 + 3/2*c_old^2
+  Operator *opChMPowExpl = new Operator(feSpace,feSpace);
+  opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+    prob->getSolution()->getDOFVector(1+3),
+    new Pow3Functor(-2.0)));
+  if (doubleWell == 0) {
+    opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+      prob->getSolution()->getDOFVector(1+3),
+      new Pow2Functor(3.0/2.0)));
+  }
+
+  // -3*c_old^2 * c
+  Operator *opChMPowImpl = new Operator(feSpace,feSpace);
+  opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+    prob->getSolution()->getDOFVector(1+3),
+    new Pow2Functor(-3.0)));
+  if (doubleWell == 0) {
+    opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+      prob->getSolution()->getDOFVector(1+3),
+      new AMDiS::Factor<double>(3.0)));
+    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(-0.5));
+  } else {
+    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
+  }
+
+  // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
+  // ----------------------------------------------------------------------
+  prob->addMatrixOperator(*opChMnew,0+3,0+3, &ab);              /// < phi*mu , psi >
+  
+  prob->addMatrixOperator(*opChMPowImpl,0+3,1+3, &b);          /// < -3*phi*c*c_old^2 , psi >
+  prob->addMatrixOperator(*opChL,0+3,1+3, &minusEpsSqr);   /// < -eps^2*phi*grad(c) , grad(psi) >
+  // . . . vectorOperators . . . . . . . . . . . . . . .
+  prob->addVectorOperator(*opChMPowExpl,0+3, &b);            /// < -2*phi*c_old^3 , psi >
+
+//   setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
+  
+  // dt(c) = laplace(mu) - u*grad(c)
+  // -----------------------------------
+  prob->addMatrixOperator(*opChMnew,1+3,1+3, &b); /// < phi*c/tau , psi >
+  prob->addMatrixOperator(*opChLM,1+3,0+3, getTau());                /// < phi*grad(mu) , grad(psi) >
+  // . . . vectorOperators . . . . . . . . . . . . . . .
+  prob->addVectorOperator(*opChMold,1+3, &b );   /// < phi*c^old/tau , psi >
+  
+  /**/
+  
+  
+  /// Navier-Stokes part
+    
+  WorldVector<DOFVector<double>* > vel;
+  for (unsigned i = 0; i < dow; i++){
+    vel[i]= prob->getSolution()->getDOFVector(i+2); 
+  }
+
+
+  // fill operators for prob
+  for (unsigned i = 0; i < dow; ++i) {
+    
+    /// < (1/tau)*rho*u'_i , psi >
+    Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
+    if (density1==density2)
+      opTime->addTerm(new Simple_ZOT(density1));  
+    else
+      opTime->addTerm(new VecAtQP_ZOT(densityPhase, NULL));
+    opTime->setUhOld(prob->getSolution()->getDOFVector(i));
+    prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());    
+    prob->addVectorOperator(*opTime, i, getInvTau(), getInvTau());    
+    
+    
+ 
+    /// < u^old*grad(u_i^old) , psi >
+/*    Operator *opUGradU0 = new Operator(getFeSpace(i),getFeSpace(i));
+    if (density1==density2)
+        opUGradU0->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PHI);
+    else
+        opUGradU0->addTerm(new WorldVecPhase_FOT(densityPhase, vel, -1.0), GRD_PHI);
+    opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
+    if (nonLinTerm == 0) {
+      prob->addVectorOperator(*opUGradU0, i);
+    } else {
+      prob->addVectorOperator(*opUGradU0, i, &theta1, &theta1);
+    }
+
+    if (nonLinTerm == 1) {
+      /// < u'*grad(u_i^old) , psi >
+      for (unsigned j = 2; j < 2+dow; ++j) {
+        Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(i));
+	if (density1==density2)
+	  opUGradU1->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(i), j-2));
+	else
+	  opUGradU1->addTerm(new VecAndPartialDerivative_ZOT(     densityPhase,    prob->getSolution()->getDOFVector(i), j-2));
+        prob->addMatrixOperator(*opUGradU1, i, j, &theta, &theta);
+      }
+    } else if (nonLinTerm == 2) {
+      /// < u^old*grad(u'_i) , psi >
+*/   for(unsigned j = 2; j < 2+dow; ++j) {
+        Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
+        opUGradU2->addTerm(new Vec2ProductPartial_FOT(   densityPhase,    prob->getSolution()->getDOFVector(j-2), j-2), GRD_PHI);
+        prob->addMatrixOperator(*opUGradU2, i, i);
+      }
+    
+  /**/
+    /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
+    Operator *opLaplaceUi = new Operator(getFeSpace(i), getFeSpace(i));
+    opLaplaceUi->addTerm(new VecAtQP_SOT(viscosityPhase, NULL));
+    prob->addMatrixOperator(*opLaplaceUi, i, i);  
+  
+    /// < p , d_i(psi) >
+    Operator *opGradP = new Operator(getFeSpace(i),getFeSpace(2));
+    opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
+    prob->addMatrixOperator(*opGradP, i, 2);
+    
+    /// external force, i.e. gravitational force
+    if (force[i]*force[i] > DBL_TOL) {
+      Operator *opForce = new Operator(getFeSpace(i), getFeSpace(i));
+      opForce->addZeroOrderTerm(new VecAtQP_ZOT(densityPhase, new Force(force[i])));
+      prob->addVectorOperator(*opForce, i);
+    }
+  
+    /// < d_i(u'_i) , psi >
+    Operator *opDivU = new Operator(getFeSpace(2),getFeSpace(i));
+    opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+    prob->addMatrixOperator(*opDivU, 2, i);
+    
+    
+    
+     ///coupling Operators 
+      Operator *opNuGradC = new Operator(getFeSpace(i), getFeSpace(dow+1));
+      opNuGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i));
+      prob->addMatrixOperator(opNuGradC, i, dow+1, &surfaceTension, &surfaceTension);
+    
+      Operator *opVGradC = new Operator(getFeSpace(dow+2), getFeSpace(i));
+      opVGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i, b));
+      prob->addMatrixOperator(opVGradC, dow+2, i,  getTau());
+      /**/
+     
+  }
+  
+  /**/
+};
+
+
+void NavierStokesCahnHilliard::addLaplaceTerm(int i)
+{ FUNCNAME("NavierStokes_TH_MultiPhase::addLaplaceTerm()");
+
+  /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+  if (viscosity1!=viscosity2) {
+    for (unsigned j = 0; j < dow; ++j) {
+      Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
+      opLaplaceUi1->addTerm(new VecAtQP_IJ_SOT(  viscosityPhase, NULL, j, i));      
+      prob->addMatrixOperator(*opLaplaceUi1, i, j);      
+    }
+  }/**/
+  
+  /// < alpha*grad(u'_i) , grad(psi) >
+  
+ 
+};

extensions/base_problems/NavierStokesCahnHilliard.h

+/** \file NavierStokesCahnHilliard.h */
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "POperators.h"
+#include "ExtendedProblemStat.h"
+#include "chns.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+  #include "parallel/PetscSolverNSCH.h"
+#endif
+
+
+
+using namespace AMDiS;
+
+class NavierStokesCahnHilliard : public BaseProblem<ProblemStat>
+{
+public: // definition of types
+
+  typedef BaseProblem<ProblemStat> super;
+
+public: // public methods
+
+  NavierStokesCahnHilliard(const std::string &name_, bool createProblem = true);
+  ~NavierStokesCahnHilliard();
+
+  virtual void initData();
+  
+  /** initTimestep of super and
+   * initialization of densityPhase and viscosityPhase
+   **/
+  virtual void initTimestep(AdaptInfo *adaptInfo);
+  virtual void closeTimestep(AdaptInfo *adaptInfo);
+  
+  virtual void setTime(AdaptInfo* adaptInfo);
+  
+/*   * indicates the different phases. Will be initialized
+   * in \ref initTimeInterface with const 1
+   **/
+  void setMultiPhase(DOFVector<double>* p) { multiPhase=p; }
+
+  
+  virtual void solveInitialProblem(AdaptInfo *adaptInfo);
+  virtual void solveInitialProblem2(AdaptInfo *adaptInfo);
+
+  double getEpsilon() { return eps; }
+  int getDoubleWellType() { return doubleWell; }
+
+protected: // protected methods
+
+  virtual void fillOperators();
+  virtual void addLaplaceTerm(int i);
+
+
+protected: // protected variables
+  ///viscosity of phase 1
+  double viscosity1;
+  
+  ///viscosity of phase 2
+  double viscosity2;
+  
+  ///density of phase 1
+  double density1;
+  
+  ///density of phase 2
+  double density2;
+  
+  double a, b, ab;
+  
+  double tau;
+
+  /// phase dependent density
+  DOFVector<double> *densityPhase;
+  
+  /// phase dependent viscosity
+  DOFVector<double> *viscosityPhase;
+  
+  double delta;
+  
+  /// phase inticator
+  DOFVector<double> *multiPhase;
+  
+  DOFVector<WorldVector<double> >* velocity;
+
+  void calcVelocity()
+  {
+    if (dow == 1)
+      transformDOF(prob->getSolution()->getDOFVector(0), velocity, new AMDiS::Vec1WorldVec<double>);
+    else if (dow == 2)
+      transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), velocity, new AMDiS::Vec2WorldVec<double>);
+    else if (dow == 3)
+      transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec3WorldVec<double>);
+  }
+
+  FileVectorWriter* fileWriter;
+  
+  bool useMobility;
+
+  unsigned dow;		// dimension of the world
+  unsigned dim;
+  PhaseFieldRefinement* refFunction;
+  RefinementLevelDOF *refinement;
+
+  double sigma, minus1;		// coupling parameter to calculate the surface tension
+  double surfaceTension;// := sigma/epsilon
+
+  int doubleWell;
+
+  double gamma;
+  double eps;
+  double minusEps;
+  double epsInv;
+  double minusEpsInv;
+  double epsSqr;
+  double minusEpsSqr;
+  
+  int nonLinTerm;
+  double theta;
+  double theta1;
+  double minusTheta1;
+  WorldVector<double> force;
+};
+
+
+
+class Force : public AbstractFunction<double, double>
+{
+public:
+  Force(double c_) :
+    AbstractFunction<double, double>(2), c(c_) {};
+    
+  double operator()(const double &x) const {    
+          return x*c;    
+  }
+
+private:
+  double c;