hh files for NavierStokes and CahnHilliard addded

extensions/base_problems/CahnHilliard.hh

+/******************************************************************************
+ *
+ * Extension of AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: Simon Praetorius et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+// #include "CahnHilliard.h"
+// #include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+#include "HL_SignedDistTraverse.h"
+#include "Recovery.h"
+
+using namespace AMDiS;
+
+namespace detail {
+
+template<typename P>
+CahnHilliard<P>::CahnHilliard(const std::string &name_) :
+  super(name_),
+  useMobility(false),
+  useReinit(false),
+  doubleWell(0),
+  gamma(1.0),
+  eps(0.1),
+  minusEps(-0.1),
+  epsInv(10.0),
+  minusEpsInv(-10.0),
+  epsSqr(0.01),
+  minusEpsSqr(-0.01)
+{
+  // parameters for CH
+  Parameters::get(name_ + "->use mobility", useMobility); // mobility
+  Parameters::get(name_ + "->gamma", gamma); // mobility
+  Parameters::get(name_ + "->epsilon", eps); // interface width
+  
+  // type of double well: 0= [0,1], 1= [-1,1]
+  Parameters::get(name_ + "->double-well type", doubleWell); 
+  
+  Parameters::get(name_ + "->use reinit", useReinit);
+
+  // transformation of the parameters
+  minusEps = -eps;
+  epsInv = 1.0/eps;
+  minusEpsInv = -epsInv;
+  epsSqr = sqr(eps);
+  minusEpsSqr = -epsSqr;
+}
+
+
+template<typename P>
+void CahnHilliard<P>::solveInitialProblem(AdaptInfo *adaptInfo) 
+{
+  Flag initFlag = self::initDataFromFile(adaptInfo);
+
+  if (!initFlag.isSet(DATA_ADOPTED)) {
+    int initialInterface = 0;
+    Initfile::get(self::name + "->initial interface", initialInterface);
+    double initialEps = eps;
+    Initfile::get(self::name + "->initial epsilon", initialEps);
+
+    if (initialInterface == 0) {
+      /// horizontale Linie
+      double a= 0.0, dir= -1.0;
+      Initfile::get(self::name + "->line->pos", a);
+      Initfile::get(self::name + "->line->direction", dir);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
+    }
+    else if (initialInterface == 1) {
+      /// schraege Linie
+      double theta = m_pi/4.0;
+      self::prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
+      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),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(self::name + "->ellipse->a", a);
+      Initfile::get(self::name + "->ellipse->b", b);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b));
+    }
+    else if (initialInterface == 3) {
+      /// zwei horizontale Linien
+      double a= 0.75, b= 0.375;
+      Initfile::get(self::name + "->lines->pos1", a);
+      Initfile::get(self::name + "->lines->pos2", b);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
+      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
+    }
+    else if (initialInterface == 4) {
+      /// Kreis
+      double radius= 1.0;
+      Initfile::get(self::name + "->kreis->radius", radius);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius));
+    } else if (initialInterface == 5) {
+      /// Rechteck
+      double width = 0.5;
+      double height = 0.3;
+      WorldVector<double> center; center.set(0.5);
+      Initfile::get(self::name + "->rectangle->width", width);
+      Initfile::get(self::name + "->rectangle->height", height);
+      Initfile::get(self::name + "->rectangle->center", center);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center));
+    }
+    
+    // TODO: Redistancing einfügen!
+    if (useReinit) {
+      FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace(
+	const_cast<DOFAdmin*>(self::getMesh()->getVertexAdmin()),
+	Lagrange::getLagrange(self::getMesh()->getDim(), 1),
+	self::getMesh(),
+	"P1");
+      DOFVector<double> tmp(feSpace, "tmp");
+      tmp.interpol(self::prob->getSolution()->getDOFVector(0));
+    
+      HL_SignedDistTraverse reinit("reinit", self::getMesh()->getDim());
+      reinit.calcSignedDistFct(adaptInfo, &tmp);
+    
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
+      Recovery recovery(L2_NORM, 1);
+      recovery.recoveryUh(&tmp, *self::prob->getSolution()->getDOFVector(0));
+#else
+      self::prob->getSolution()->getDOFVector(0)->interpol(&tmp);
+#endif
+    }
+    
+
+    /// create phase-field from signed-dist-function
+    if (doubleWell == 0) {
+      forEachDOF(self::prob->getSolution()->getDOFVector(0),
+        new SignedDistToPhaseField(initialEps));
+    } else {
+      forEachDOF(self::prob->getSolution()->getDOFVector(0),
+        new SignedDistToCh(initialEps));
+    }
+  }
+}
+
+
+template<typename P>
+void CahnHilliard<P>::fillOperators()
+{  
+  const FiniteElemSpace* feSpace = self::prob->getFeSpace();
+  int degree = feSpace->getBasisFcts()->getDegree();
+
+  // c
+  Operator *opChMnew = new Operator(feSpace, feSpace);
+  opChMnew->addTerm(new Simple_ZOT);
+  Operator *opChMold = new Operator(feSpace, feSpace);
+  opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
+  // -nabla*(grad(c))
+  Operator *opChL = new Operator(feSpace, feSpace);
+  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)));
+    else
+      opChLM->addTerm(new VecAtQP_SOT(
+	self::prob->getSolution()->getDOFVector(0),
+	new MobilityCH1(gamma, degree)));
+  } else
+    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)));
+  if (doubleWell == 0) {
+    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)));
+  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));
+  }
+
+  // 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 >
+  self::prob->addMatrixOperator(*opChL,0,0, &minusEpsSqr);   /// < -eps^2*phi*grad(c) , grad(psi) >
+  self::prob->addMatrixOperator(*opChMnew,0,1);              /// < phi*mu , psi >
+  // . . . vectorOperators . . . . . . . . . . . . . . .
+  self::prob->addVectorOperator(*opChMPowExpl,0);            /// < -2*phi*c_old^3 , psi >
+
+  
+  // dt(c) = laplace(mu) - u*grad(c)
+  // -----------------------------------
+  self::prob->addMatrixOperator(*opChMnew,1,0, self::getInvTau()); /// < phi*c/tau , psi >
+  self::prob->addMatrixOperator(*opChLM,1,1);                /// < phi*grad(mu) , grad(psi) >
+  // . . . vectorOperators . . . . . . . . . . . . . . .
+  self::prob->addVectorOperator(*opChMold,1, self::getInvTau());   /// < phi*c^old/tau , psi >
+  
+}
+
+
+template<typename P>
+void CahnHilliard<P>::finalizeData()
+{
+  self::setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
+  self::setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
+  if (!useMobility)
+    self::setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
+}
+
+} // end namespace detail

extensions/base_problems/NavierStokes_TaylorHood.hh

+/******************************************************************************
+ *
+ * Extension of AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: Simon Praetorius et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+#include "Helpers.h"
+#include "POperators.h"
+
+namespace detail {
+
+using namespace AMDiS;
+
+template<typename P> 
+NavierStokes_TaylorHood<P>::NavierStokes_TaylorHood(const std::string &name_, bool createProblem) :
+  super(name_, createProblem),
+  velocity(NULL),
+  laplaceType(0),
+  nonLinTerm(2),
+  oldTimestep(0.0),
+  viscosity(1.0),
+  density(1.0),
+  minus1(-1.0),
+  theta(0.5+1.e-2),
+  theta1(1-theta),
+  minusTheta1(-theta1),
+  fileWriter(NULL)
+{  
+  // fluid viscosity
+  Initfile::get(name_ + "->viscosity", viscosity);
+  // fluid density
+  Initfile::get(name_ + "->density", density);
+
+  // type of laplace operator: 
+  //   0... div(nu*grad(u)), 
+  //   1... div(0.5*nu*(grad(u)+grad(u)^T))
+  Initfile::get(name_ + "->laplace operator", laplaceType); 
+
+  // type of non-linear term: 
+  //   0... u^old*grad(u_i^old), 
+  //   1... u*grad(u_i^old), 
+  //   2... u^old*grad(u_i)
+  //   3... (1)+(2)-(0)
+  Initfile::get(name_ + "->non-linear term", nonLinTerm); 
+
+  // theta-scheme parameter
+  Initfile::get(name_ + "->theta", theta); 
+  theta1 = 1.0-theta;
+  minusTheta1 = -theta1;
+
+  // volume force
+  force.set(0.0);
+  Initfile::get(name_ + "->force", force);
+
+  oldSolution.resize(self::dow);
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i] = NULL;
+}
+
+template<typename P> 
+NavierStokes_TaylorHood<P>::~NavierStokes_TaylorHood() 
+{ FUNCNAME("NavierStokes_TaylorHood::~NavierStokes_TaylorHood()");
+  
+  if (velocity != NULL) {
+    delete velocity;
+    velocity = NULL;
+  }
+  
+  for (size_t i = 0; i < self::dow; i++) {
+    if (oldSolution[i] != NULL)
+      delete oldSolution[i];
+    oldSolution[i] = NULL;
+  }
+ 
+  if (fileWriter) { 
+    delete fileWriter;
+    fileWriter = NULL;
+  }
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::initData()
+{ FUNCNAME("NavierStokes_TaylorHood::initTimeInterface()");
+
+  if (velocity == NULL)
+    velocity = new DOFVector<WorldVector<double> >(self::getFeSpace(0), "velocity");
+  
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i] = new DOFVector<double>(self::getFeSpace(i), "old(v_"+Helpers::toString(i)+")");
+  
+  fileWriter = new FileVectorWriter(self::name + "->velocity->output", self::getFeSpace()->getMesh(), velocity);
+
+  super::initData();
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::transferInitialSolution(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokes_TaylorHood::transferInitialSolution()");
+
+  calcVelocity();
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i]->copy(*self::prob->getSolution()->getDOFVector(i));
+  
+  super::transferInitialSolution(adaptInfo);
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::fillOperators()
+{ FUNCNAME("NavierStokes_TaylorHood::fillOperators()");
+
+  WorldVector<DOFVector<double>* > vel;
+  for (size_t k = 0; k < self::dow; k++)
+    vel[k] = self::prob->getSolution()->getDOFVector(k);
+  
+  const FiniteElemSpace* feSpace_u = self::getFeSpace(0);
+  const FiniteElemSpace* feSpace_p = self::getFeSpace(self::dow);
+  
+  // fill operators for prob
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < (1/tau)*u'_i , psi >
+    Operator *opTime = new Operator(feSpace_u, feSpace_u);
+    opTime->addTerm(new Simple_ZOT(density));
+    self::prob->addMatrixOperator(*opTime, i, i, self::getInvTau(), self::getInvTau());
+    
+    /// < (1/tau)*u_i^old , psi >
+    Operator *opTimeOld = new Operator(feSpace_u, feSpace_u);
+    opTimeOld->addTerm(new Phase_ZOT(getOldSolution(i), density));
+    self::prob->addVectorOperator(*opTimeOld, i, self::getInvTau(), self::getInvTau());
+
+    /// < u^old*grad(u_i^old) , psi >
+    if (nonLinTerm == 0 || nonLinTerm == 3) {
+      Operator *opUGradU0 = new Operator(feSpace_u, feSpace_u);
+      opUGradU0->addTerm(new WorldVec_FOT(vel, density), GRD_PHI);
+      opUGradU0->setUhOld(self::prob->getSolution()->getDOFVector(i));
+      if (nonLinTerm == 0)
+        self::prob->addVectorOperator(*opUGradU0, i);
+      else
+        self::prob->addVectorOperator(*opUGradU0, i, &minus1, &minus1);
+    }
+    
+    /// < u*grad(u_i^old) , psi >
+    if (nonLinTerm == 1 || nonLinTerm == 3) {
+      for (size_t k = 0; k < self::dow; ++k) {
+        Operator *opUGradU1 = new Operator(feSpace_u, feSpace_u);
+        opUGradU1->addTerm(new PartialDerivative_ZOT(
+	  self::prob->getSolution()->getDOFVector(i), k, density));
+        self::prob->addMatrixOperator(*opUGradU1, i, k);
+      }
+    }
+
+    /// < u^old*grad(u_i) , psi >
+    if (nonLinTerm == 2 || nonLinTerm == 3) {
+      for (size_t k = 0; k < self::dow; ++k) {
+        Operator *opUGradU2 = new Operator(feSpace_u, feSpace_u);
+        opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
+	  vel[k], k, density), GRD_PHI);
+        self::prob->addMatrixOperator(*opUGradU2, i, i);
+      }
+    }
+
+    /// Diffusion-Operator
+    addLaplaceTerm(i);
+  
+    /// < p , d_i(psi) >
+    Operator *opGradP = new Operator(feSpace_u, feSpace_p);
+    opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
+    self::prob->addMatrixOperator(*opGradP, i, self::dow);
+  
+    /// external force, i.e. gravitational force
+    if (std::abs(force[i]) > DBL_TOL) {
+      Operator *opForce = new Operator(feSpace_u, feSpace_u);
+      opForce->addTerm(new Simple_ZOT(density*force[i]));
+      self::prob->addVectorOperator(*opForce, i);
+    }
+  }
+
+  /// div(u) = 0
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < d_i(u'_i) , psi >
+    Operator *opDivU = new Operator(feSpace_p, feSpace_u);
+    opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+    self::prob->addMatrixOperator(*opDivU, self::dow, i);
+  }
+
+  Operator *opNull = new Operator(feSpace_p, feSpace_p);
+  opNull->addTerm(new Simple_ZOT(0.0));
+  self::prob->addMatrixOperator(*opNull, self::dow, self::dow);
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::addLaplaceTerm(int i)
+{ FUNCNAME("NavierStokes_TaylorHood::addLaplaceTerm()");
+
+  const FiniteElemSpace* feSpace_u = self::getFeSpace(0);
+  
+  /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+  if (self::laplaceType == 1) {
+    for (size_t j = 0; j < self::dow; ++j) {
+      Operator *opLaplaceUi1 = new Operator(feSpace_u, feSpace_u);
+      opLaplaceUi1->addTerm(new MatrixIJ_SOT(j, i, viscosity));
+      self::prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta, &theta);
+
+      if (std::abs(minusTheta1) > DBL_TOL) {
+	opLaplaceUi1->setUhOld(self::prob->getSolution()->getDOFVector(j));
+	self::prob->addVectorOperator(*opLaplaceUi1, i, &minusTheta1, &minusTheta1);
+      }
+    }
+  }
+  
+  /// < alpha*grad(u'_i) , grad(psi) >
+  Operator *opLaplaceUi = new Operator(feSpace_u, feSpace_u);
+  opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+  self::prob->addMatrixOperator(*opLaplaceUi, i, i, &theta, &theta);
+
+  if (std::abs(minusTheta1) > DBL_TOL) {
+    opLaplaceUi->setUhOld(self::prob->getSolution()->getDOFVector(i));
+    self::prob->addVectorOperator(*opLaplaceUi, i, &minusTheta1, &minusTheta1);
+  }
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokes_TaylorHood::closeTimestep()");
+
+  calcVelocity();
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i]->copy(*self::prob->getSolution()->getDOFVector(i));
+  
+  super::closeTimestep(adaptInfo);
+}
+
+
+template<typename P> 
+void NavierStokes_TaylorHood<P>::writeFiles(AdaptInfo *adaptInfo, bool force)
+{ FUNCNAME("NavierStokesPhase_TaylorHood::closeTimestep()");
+
+  super::writeFiles(adaptInfo, force);
+  self::fileWriter->writeFiles(adaptInfo, false);
+}
+
+} // end namespace detail
\ No newline at end of file