diff --git a/extensions/base_problems/BaseProblem_RB.h b/extensions/base_problems/BaseProblem_RB.h
new file mode 100644
index 0000000000000000000000000000000000000000..b395524d4cc020c8e7a3dfe6fac47f31379b4dd5
--- /dev/null
+++ b/extensions/base_problems/BaseProblem_RB.h
@@ -0,0 +1,44 @@
+/** \file BaseProblem_RB.h */
+
+#ifndef BASE_PROBLEM_RB_H
+#define BASE_PROBLEM_RB_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "time/ExtendedRosenbrockStationary.h"
+
+using namespace AMDiS;
+
+/** BaseProblems for Rosenbrock time-discretization
+ */
+
+class BaseProblem_RB : public BaseProblem<ExtendedRosenbrockStationary>
+{
+public: // typedefs
+
+ typedef BaseProblem<ExtendedRosenbrockStationary> super;
+
+// Rosenbrock methods
+//_________________________________________________________________________
+
+public:
+
+ BaseProblem_RB(std::string name_) : super(name_), minus1(-1.0) {}
+
+ 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(); }
+
+protected:
+ double minus1;
+};
+
+#endif // BASE_PROBLEM_RB_H
\ No newline at end of file
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.cc b/extensions/base_problems/CahnHilliardNavierStokes.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6ef2b8ce261301bc96b758330866fcf4023c178d
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes.cc
@@ -0,0 +1,302 @@
+#include "CahnHilliardNavierStokes.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "POperators.h"
+
+using namespace AMDiS;
+
+CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
+ super(name_),
+ useMobility(false),
+ useConservationForm(false),
+ doubleWell(0),
+ gamma(1.0),
+ eps(0.1),
+ minusEps(-0.1),
+ epsInv(10.0),
+ minusEpsInv(-10.0),
+ epsSqr(0.01),
+ minusEpsSqr(-0.01),
+ laplaceType(0),
+ nonLinTerm(2),
+ oldTimestep(0.0),
+ viscosity(1.0),
+ density(1.0),
+ c(0.0),
+ sigma(1.0),
+ surfaceTension(1.0),
+ reinit(NULL),
+ velocity(NULL),
+ fileWriter(NULL)
+{
+ // 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 for navier-stokes
+ Initfile::get(name + "->viscosity", viscosity);
+ 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)
+ Initfile::get(name + "->non-linear term", nonLinTerm);
+
+ // Parameters for CH-Coupling
+ Initfile::get(name + "->use conservation form", useConservationForm);
+
+ // Parameters for NS-Coupling
+ Initfile::get(name + "->sigma", sigma); // surface tension
+ surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
+
+
+ force.set(0.0);
+ Initfile::get(name + "->force", force);
+
+ // transformation of the parameters
+ minusEps = -eps;
+ epsInv = 1.0/eps;
+ minusEpsInv = -epsInv;
+ epsSqr = sqr(eps);
+ minusEpsSqr = -epsSqr;
+}
+
+
+CahnHilliardNavierStokes::~CahnHilliardNavierStokes()
+{ FUNCNAME("CahnHilliardNavierStokes::~CahnHilliardNavierStokes()");
+
+ if (reinit != NULL) {
+ delete reinit;
+ reinit = NULL;
+ }
+ if (velocity != NULL) {
+ delete velocity;
+ velocity = NULL;
+ }
+ delete fileWriter;
+ fileWriter = NULL;
+}
+
+
+void CahnHilliardNavierStokes::initData()
+{ FUNCNAME("CahnHilliardNavierStokes::initData()");
+
+ dim = getMesh()->getDim();
+ // create instance redistancing class
+ reinit = new HL_SignedDistTraverse("reinit", dim);
+
+ if (velocity == NULL)
+ velocity = new DOFVector<WorldVector<double> >(getFeSpace(2), "velocity");
+ fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace(2)->getMesh(), velocity);
+
+ super::initData();
+}
+
+
+void CahnHilliardNavierStokes::fillOperators()
+{ FUNCNAME("CahnHilliardNavierStokes::fillOperators()");
+ MSG("CahnHilliardNavierStokes::fillOperators()\n");
+
+ // variable order:
+ // (c, mu, u0, u1 [, u2], p)
+
+ WorldVector<DOFVector<double>* > vel;
+ for (size_t i = 0; i < dow; i++)
+ vel[i] = prob->getSolution()->getDOFVector(2+i);
+
+
+ // dt(c) = laplace(mu) - u*grad(c)
+ // -----------------------------------
+ /// < phi*c/tau , psi >
+ Operator *opChMnew = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+ opChMnew->addZeroOrderTerm(new Simple_ZOT);
+ prob->addMatrixOperator(*opChMnew, 0, 0, getInvTau());
+
+ /// < phi*c^old/tau , psi >
+ Operator *opChMold = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+ opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0)));
+ prob->addVectorOperator(*opChMold,0, getInvTau());
+
+ /// < phi*grad(mu) , grad(psi) >
+ // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
+ Operator *opChLM = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ if (useMobility) {
+ if (doubleWell == 0)
+ opChLM->addSecondOrderTerm(new VecAtQP_SOT(
+ prob->getSolution()->getDOFVector(0),
+ new MobilityCH0(gamma)));
+ else
+ opChLM->addSecondOrderTerm(new VecAtQP_SOT(
+ prob->getSolution()->getDOFVector(0),
+ new MobilityCH1(gamma)));
+ } else
+ opChLM->addSecondOrderTerm(new Simple_SOT(gamma));
+ prob->addMatrixOperator(*opChLM, 0, 1);
+
+ /// < u_old * grad(c) - u_old * grad(c_old), psi >
+ Operator *uGradC = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ uGradC->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PSI);
+ else
+ uGradC->addTerm(new WorldVec_FOT(vel), GRD_PHI);
+ uGradC->setUhOld(prob->getSolution()->getDOFVector(0));
+ prob->addVectorOperator(*uGradC, 0);
+ prob->addMatrixOperator(*uGradC, 0, 0);
+
+ /// < u * grad(c_old) , psi >
+ for (size_t i = 0; i < dow; i++) {
+ Operator *uGradC2 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ uGradC2->addTerm(new VecAndPartialDerivative_FOT(prob->getSolution()->getDOFVector(0), i, -1.0), GRD_PSI);
+ else {
+ uGradC2->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(0), i));
+ }
+ prob->addMatrixOperator(*uGradC2, 0, 2+i);
+ }
+
+
+ // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
+ // ----------------------------------------------------------------------
+ /// < -3*phi*c*c_old^2 , psi >
+ // -3*c_old^2 * c
+ Operator *opChMPowImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new Pow2Functor(-3.0)));
+ if (doubleWell == 0) {
+ opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new AMDiS::Factor<double>(3.0)));
+ opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(-0.5));
+ } else {
+ opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
+ }
+ prob->addMatrixOperator(*opChMPowImpl, 1, 0);
+
+ /// < -2*phi*c_old^3 , psi >
+ // -2*c_old^3 + 3/2*c_old^2
+ Operator *opChMPowExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new Pow3Functor(-2.0)));
+ if (doubleWell == 0) {
+ opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new Pow2Functor(3.0/2.0)));
+ }
+ prob->addVectorOperator(*opChMPowExpl, 1);
+
+ /// < -eps^2*phi*grad(c) , grad(psi) >
+ Operator *opChL = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChL->addSecondOrderTerm(new Simple_SOT);
+ prob->addMatrixOperator(*opChL, 1, 0, &minusEpsSqr);
+
+ /// < phi*mu , psi >
+ Operator *opChM = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+ opChM->addZeroOrderTerm(new Simple_ZOT);
+ prob->addMatrixOperator(*opChM, 1, 1);
+
+ // -----------------------------------------------------------------------------------
+ for (size_t i = 0; i < dow; ++i) {
+ /// < (1/tau)*u'_i , psi >
+ Operator *opTime = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opTime->addTerm(new Simple_ZOT(density));
+ prob->addMatrixOperator(*opTime, 2+i, 2+i, getInvTau(), getInvTau());
+ /// < (1/tau)*u_i^old , psi >
+ Operator *opTimeOld = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opTime->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(2+i), density));
+ prob->addVectorOperator(*opTimeOld, 2+i, getInvTau(), getInvTau());
+
+ /// < u^old*grad(u_i^old) , psi >
+ Operator *opUGradU0 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(2+i));
+ if (nonLinTerm == 0) {
+ prob->addVectorOperator(*opUGradU0, 2+i);
+ }
+
+ if (nonLinTerm == 1) {
+ /// < u'*grad(u_i^old) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU1 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
+ opUGradU1->addTerm(new PartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(2+i), j, density));
+ prob->addMatrixOperator(*opUGradU1, 2+i, 2+j);
+ }
+ } else if (nonLinTerm == 2) {
+ /// < u^old*grad(u'_i) , psi >
+ for(size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU2 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
+ opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
+ prob->getSolution()->getDOFVector(2+j), j, density), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, 2+i, 2+i);
+ }
+ }
+
+// for (size_t j = 0; j < dow; ++j) {
+// Operator *opNull = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+// opNull->addTerm(new Simple_ZOT(0.0));
+// prob->addMatrixOperator(*opNull, 2+i, 2+j);
+// }
+
+ /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
+ addLaplaceTerm(i);
+
+ /// < p , d_i(psi) >
+ Operator *opGradP = new Operator(getFeSpace(2+i),getFeSpace(2+dow));
+ opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
+ prob->addMatrixOperator(*opGradP, 2+i, 2+dow);
+
+ /// external force, i.e. gravitational force
+ if (norm(force) > DBL_TOL) {
+ Operator *opForce = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opForce->addTerm(new Simple_ZOT(force[i]));
+ prob->addVectorOperator(*opForce, 2+i);
+ }
+
+ /// forces by Cahn-Hilliard terms
+ Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
+ opMuGradC->addTerm(new VecAndPartialDerivative_FOT(
+ prob->getSolution()->getDOFVector(1), i), GRD_PHI);
+ prob->addMatrixOperator(*opMuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
+
+ Operator *opMuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
+ opMuGradC2->addTerm(new PartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(0), i));
+ prob->addMatrixOperator(*opMuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
+ opMuGradC2->setUhOld(prob->getSolution()->getDOFVector(1));
+ prob->addVectorOperator(*opMuGradC2, 2+i, &surfaceTension, &surfaceTension);
+
+ }
+
+ /// div(u) = 0
+ for (size_t i = 0; i < dow; ++i) {
+ /// < d_i(u'_i) , psi >
+ Operator *opDivU = new Operator(getFeSpace(2+dow),getFeSpace(2+i));
+ opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+ prob->addMatrixOperator(*opDivU, 2+dow, 2+i);
+ }
+}
+
+
+void CahnHilliardNavierStokes::addLaplaceTerm(int i)
+{ FUNCNAME("CahnHilliardNavierStokes::addLaplaceTerm()");
+
+ /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+ if (laplaceType == 1) {
+ for (size_t j = 0; j < dow; ++j) {
+ 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);
+ }
+ }
+
+ /// < alpha*grad(u'_i) , grad(psi) >
+ Operator *opLaplaceUi = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+ prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
+}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.h b/extensions/base_problems/CahnHilliardNavierStokes.h
new file mode 100644
index 0000000000000000000000000000000000000000..56686fdd5789e67e57e7e3422ff8af9e83953842
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes.h
@@ -0,0 +1,175 @@
+/** \file CahnHilliardNavierStokes.h */
+
+#ifndef CAHN_HILLIARD_H
+#define CAHN_HILLIARD_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "ExtendedProblemStat.h"
+#include "HL_SignedDistTraverse.h"
+
+using namespace AMDiS;
+
+class CahnHilliardNavierStokes : public BaseProblem<ExtendedProblemStat>
+{
+public: // definition of types
+
+ typedef BaseProblem<ExtendedProblemStat> super;
+
+public: // public methods
+
+ CahnHilliardNavierStokes(const std::string &name_);
+ ~CahnHilliardNavierStokes();
+
+ virtual void initData();
+
+ double getEpsilon() { return eps; }
+ int getDoubleWellType() { return doubleWell; }
+
+ DOFVector<WorldVector<double> >* getVelocity()
+ {
+ return velocity;
+ }
+
+protected: // protected methods
+
+ virtual void fillOperators();
+ virtual void fillBoundaryConditions() {}
+ virtual void addLaplaceTerm(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 viscosity;
+ double density;
+ 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;
+};
+
+
+/** \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;
+};
+
+#endif // CAHN_HILLIARD_H
diff --git a/extensions/base_problems/CahnHilliard_RB.cc b/extensions/base_problems/CahnHilliard_RB.cc
new file mode 100644
index 0000000000000000000000000000000000000000..94cc4798c328fd827562ba340327ee696f99b30c
--- /dev/null
+++ b/extensions/base_problems/CahnHilliard_RB.cc
@@ -0,0 +1,226 @@
+#include "CahnHilliard_RB.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "Tools.h"
+
+using namespace AMDiS;
+
+CahnHilliard_RB::CahnHilliard_RB(const std::string &name_) :
+ super(name_),
+ 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),
+ reinit(NULL)
+{
+ // 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);
+
+ // transformation of the parameters
+ minusEps = -eps;
+ epsInv = 1.0/eps;
+ minusEpsInv = -epsInv;
+ epsSqr = sqr(eps);
+ minusEpsSqr = -epsSqr;
+}
+
+
+CahnHilliard_RB::~CahnHilliard_RB()
+{ FUNCNAME("CahnHilliard_RB::~CahnHilliard_RB()");
+
+ if (reinit)
+ delete reinit;
+}
+
+
+void CahnHilliard_RB::initData()
+{ FUNCNAME("CahnHilliard_RB::initData()");
+
+ // create instance redistancing class
+ reinit = new HL_SignedDistTraverse("reinit", getMesh()->getDim());
+
+ super::initData();
+}
+
+
+void CahnHilliard_RB::solveInitialProblem(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliard_RB::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 == -1) {
+ /// random values
+ prob->getSolution()->getDOFVector(0)->interpol(new tools::Random(0.5 - doubleWell*0.5, 1.0 + doubleWell));
+ }
+ else 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(0)->interpol(new Plane(a, dir));
+ }
+ else if (initialInterface == 1) {
+ /// schraege Linie
+ double theta = m_pi/4.0;
+ prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
+ transformDOFInterpolation(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(name + "->ellipse->a", a);
+ Initfile::get(name + "->ellipse->b", b);
+ 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(name + "->lines->pos1", a);
+ Initfile::get(name + "->lines->pos2", b);
+ prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
+ transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
+ }
+ else if (initialInterface == 4) {
+ /// Kreis
+ double radius= 1.0;
+ Initfile::get(name + "->kreis->radius", radius);
+ 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(name + "->rectangle->width", width);
+ Initfile::get(name + "->rectangle->height", height);
+ Initfile::get(name + "->rectangle->center", center);
+ prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center));
+ }
+
+ /// create phase-field from signed-dist-function
+ if (doubleWell == 0) {
+ transformDOF(prob->getSolution()->getDOFVector(0),
+ new SignedDistToPhaseField(initialEps));
+ } else {
+ transformDOF(prob->getSolution()->getDOFVector(0),
+ new SignedDistToCh(initialEps));
+ }
+ }
+}
+
+
+void CahnHilliard_RB::fillOperators()
+{ FUNCNAME("CahnHilliard_RB::fillOperators()");
+ MSG("CahnHilliard_RB::fillOperators()\n");
+
+ // dt(c) = laplace(mu) - u*grad(c)
+ // -----------------------------------
+ prob->addTimeOperator(0,0); /// < dt(c) , psi >
+
+ // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
+ Operator *opChLM = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ if (useMobility) {
+ if (doubleWell == 0) {
+ // jacobian operators
+ opChLM->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new MobilityCH0(gamma)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH0Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH0(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
+ } else {
+ // jacobian operators
+ opChLM->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new MobilityCH1(gamma)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH1Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH1(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
+ }
+ } else {
+ opChLM->addTerm(new Simple_SOT(gamma));
+ opChLM->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM, 0, &minus1, &minus1);
+ prob->addMatrixOperator(*opChLM, 0, 1); /// < phi*grad(mu) , grad(psi) >
+ }
+
+ // mu + eps^2*laplace(c) - M'(c)
+ // ----------------------------------------------------------------------
+ /// < phi*mu , psi >
+ Operator *opChMu = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+ opChMu->addZeroOrderTerm(new Simple_ZOT);
+ opChMu->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChMu, 1);
+ prob->addMatrixOperator(*opChMu, 1, 1, &minus1, &minus1);
+
+ /// < -eps^2*phi*grad(c) , grad(psi) >
+ Operator *opChLC = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChLC->addSecondOrderTerm(new Simple_SOT(minusEpsSqr));
+ opChLC->setUhOld(prob->getStageSolution(0));
+ prob->addVectorOperator(*opChLC, 1);
+ prob->addMatrixOperator(*opChLC, 1, 0, &minus1, &minus1);
+
+ /// < -M'(c) , psi >
+ if (doubleWell == 0) {
+ // jacobian operators
+ Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getUnVec(0),
+ new DoubleWell0Diff(-1.0))); // < (3c^2-3c+0.5)*c' , psi >
+ prob->addMatrixOperator(*opChMImpl, 1, 0, &minus1, &minus1);
+ // rhs operators
+ Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getStageSolution(0),
+ new DoubleWell0(-1.0))); // < (c^3-3/2*c^2+0.5) , psi >
+ prob->addVectorOperator(*opChMExpl, 1);
+ }
+ else if (doubleWell == 1) {
+ // jacobian operators
+ Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getUnVec(0),
+ new DoubleWell1Diff(-1.0))); // < (3c^2-1)*c' , psi >
+ prob->addMatrixOperator(*opChMImpl, 1, 0, &minus1, &minus1);
+ // rhs operators
+ Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getStageSolution(0),
+ new DoubleWell1(-1.0))); // < (c^3-c) , psi >
+ prob->addVectorOperator(*opChMExpl, 1);
+ }
+}
diff --git a/extensions/base_problems/CahnHilliard_RB.h b/extensions/base_problems/CahnHilliard_RB.h
new file mode 100644
index 0000000000000000000000000000000000000000..4614b33fcb2c7a3f3315ef3f180710ed044bb39a
--- /dev/null
+++ b/extensions/base_problems/CahnHilliard_RB.h
@@ -0,0 +1,207 @@
+/** \file CahnHilliard_RB.h */
+
+#ifndef CAHN_HILLIARD_RB_H
+#define CAHN_HILLIARD_RB_H
+
+#include "AMDiS.h"
+#include "BaseProblem_RB.h"
+#include "ExtendedProblemStat.h"
+#include "HL_SignedDistTraverse.h"
+
+using namespace AMDiS;
+
+class CahnHilliard_RB : public BaseProblem_RB
+{
+public: // definition of types
+
+ typedef BaseProblem_RB super;
+
+public: // public methods
+
+ CahnHilliard_RB(const std::string &name_);
+ ~CahnHilliard_RB();
+
+ virtual void initData();
+ virtual void solveInitialProblem(AdaptInfo *adaptInfo);
+
+ double getEpsilon() { return eps; }
+ int getDoubleWellType() { return doubleWell; }
+
+protected: // protected methods
+
+ virtual void fillOperators();
+ virtual void fillBoundaryConditions() {}
+
+protected: // protected variables
+
+ HL_SignedDistTraverse *reinit;
+
+ bool useMobility;
+
+ unsigned dim;
+
+ int doubleWell;
+
+ double gamma;
+ double eps;
+ double minusEps;
+ double epsInv;
+ double minusEpsInv;
+ double epsSqr;
+ double minusEpsSqr;
+};
+
+
+/** \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 MobilityCH0Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH0Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(0.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - 1.5*sqr(phase) + 0.5*phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class MobilityCH1Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH1Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(-1.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - 1.5*sqr(c) + 0.5*c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 3.0*c + 0.5);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 1.0);
+ }
+
+ protected:
+ double factor;
+};
+
+#endif // CAHN_HILLIARD_RB_H