diff --git a/extensions/POperators.cc b/extensions/POperators.cc
index 0023fa8bf55d503fee12cab258316f5cd630eaec..91bd2b32f739af40d0fbcc87dc2f1cc9188ef814 100644
--- a/extensions/POperators.cc
+++ b/extensions/POperators.cc
@@ -827,6 +827,86 @@ void WorldVecPhase_FOT::eval(int nPoints,
}
};
+
+/* ----------------------------------------------------------- */
+
+VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_,
+ WorldVector<DOFVector<double>*> vecs_,
+ AbstractFunction<double, double>* fct_, double fac_)
+ : FirstOrderTerm(2), vDV(vec0_), fct(fct_), fac(fac_)
+{
+ numVecs = vecs_.size();
+ TEST_EXIT(numVecs == 2 || numVecs == 3)
+ ("Only Dim=2 or Dim=3 possible\n");
+
+ for (int i = 0; i < numVecs; i++) {
+ TEST_EXIT(vecs_[i])("One vector is NULL!\n");
+ auxFeSpaces.insert(vecs_[i]->getFeSpace());
+ }
+
+ vec0DV = vecs_[0];
+ vec1DV = vecs_[1];
+ if (numVecs >= 3)
+ vec2DV = vecs_[2];
+};
+void VecAndWorldVec_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vDV, elInfo, subAssembler, quad, v);
+ getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
+ getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
+ if (numVecs >= 3)
+ getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
+};
+void VecAndWorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vDV, smallElInfo, largeElInfo, subAssembler, quad, v);
+ getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
+ getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
+ if (numVecs >= 3)
+ getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
+};
+void VecAndWorldVec_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& result) const
+{
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(result.size());
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ double f = (*fct)(v[iq]);
+ WorldVector<double> vec;
+ vec[0] = vec0[iq];
+ vec[1] = vec1[iq];
+ if (numVecs >= 3)
+ vec[2] = vec2[iq];
+ lb(Lambda, f*vec, result[iq], fac);
+ }
+};
+void VecAndWorldVec_FOT::eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double opFactor)
+{
+ double factor = fac * opFactor;
+ for (int iq = 0; iq < nPoints; iq++) {
+ double resultQP = 0.0;
+
+ resultQP += grdUhAtQP[iq][0] * vec0[iq];
+ resultQP += grdUhAtQP[iq][1] * vec1[iq];
+ if (numVecs >= 3)
+ resultQP += grdUhAtQP[iq][2] * vec2[iq];
+
+ result[iq] += factor * resultQP * (*fct)(v[iq]);
+ }
+};
+
+
+
/* ----------------------------------------------------------- */
// vec1*vec2*d/dxi(...)
@@ -872,6 +952,164 @@ void VecAndPartialDerivative_FOT::eval(int nPoints,
result[iq] += opFactor * vec1[iq] * grdUhAtQP[iq][component] * fac;
};
+/* ----------------------------------------------------------- */
+
+// vec1*vec2*d/dxi(...)
+VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ int i_,
+ int j_,
+ double fac_)
+: FirstOrderTerm(4), vec1DV(dv1), vec2DV(dv2), fct(NULL), fac(fac_), xi(i_), xj(j_)
+{
+ auxFeSpaces.insert(vec1DV->getFeSpace());
+ auxFeSpaces.insert(vec2DV->getFeSpace());
+
+ setB(xi);
+}
+
+VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ int i_,
+ int j_,
+ BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, d_j(v2))
+ double fac_)
+: FirstOrderTerm(4), vec1DV(dv1), vec2DV(dv2), fct(fct_), fac(fac_), xi(i_), xj(j_)
+{
+ auxFeSpaces.insert(vec1DV->getFeSpace());
+ auxFeSpaces.insert(vec2DV->getFeSpace());
+
+ setB(xi);
+}
+
+void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
+ getGradientsAtQPs(vec2DV, elInfo, subAssembler, quad, grad);
+}
+
+void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
+ getGradientsAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, grad);
+}
+
+void VecAndPartialDerivativeIJ_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& result) const
+{
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(result.size());
+
+ if (fct == NULL) {
+ for(int iq = 0; iq < nPoints; iq++) {
+ lb_one(Lambda, result[iq], vec1[iq] * grad[iq][xj] * fac);
+ }
+ } else {
+ for(int iq = 0; iq < nPoints; iq++) {
+ lb_one(Lambda, result[iq], (*fct)(vec1[iq], grad[iq][xj]) * fac);
+ }
+ }
+}
+
+void VecAndPartialDerivativeIJ_FOT::eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double opFactor)
+{
+ if (fct == NULL) {
+ for(int iq = 0; iq < nPoints; iq++)
+ result[iq] += opFactor * vec1[iq] * grad[iq][xj] * grdUhAtQP[iq][xi] * fac;
+ } else {
+ for(int iq = 0; iq < nPoints; iq++)
+ result[iq] += opFactor * (*fct)(vec1[iq], grad[iq][xj]) * grdUhAtQP[iq][xi] * fac;
+ }
+}
+
+
+/* ----------------------------------------------------------- */
+
+// vec1*vec2*d/dxi(...)
+Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
+ DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ int component_,
+ double fac_)
+: FirstOrderTerm(3), vec1DV(dv1), vec2DV(dv2), fct(NULL), fac(fac_), component(component_)
+{
+ auxFeSpaces.insert(vec1DV->getFeSpace());
+ auxFeSpaces.insert(vec2DV->getFeSpace());
+
+ setB(component);
+}
+
+
+Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
+ DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ int component_,
+ BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, v2)
+ double fac_)
+: FirstOrderTerm(3), vec1DV(dv1), vec2DV(dv2), fct(fct_), fac(fac_), component(component_)
+{
+ auxFeSpaces.insert(vec1DV->getFeSpace());
+ auxFeSpaces.insert(vec2DV->getFeSpace());
+
+ setB(component);
+}
+
+void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
+ getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
+}
+
+void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+{
+ getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
+ getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
+}
+
+void Vec2AndPartialDerivative_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& result) const
+{
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(result.size());
+
+ if (fct == NULL) {
+ for(int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, result[iq], vec1[iq] * vec2[iq] * fac);
+ } else {
+ for(int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, result[iq], (*fct)(vec1[iq], vec2[iq]) * fac);
+ }
+};
+void Vec2AndPartialDerivative_FOT::eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double opFactor)
+{
+ if (fct == NULL) {
+ for(int iq = 0; iq < nPoints; iq++)
+ result[iq] += opFactor * vec1[iq] * vec2[iq] * grdUhAtQP[iq][component] * fac;
+ } else {
+ for(int iq = 0; iq < nPoints; iq++)
+ result[iq] += opFactor * (*fct)(vec1[iq], vec2[iq]) * grdUhAtQP[iq][component] * fac;
+ }
+};
+
+
/* ----------------------------------------------------------- */
// vec1*vec2*d/dxi(...)
@@ -1161,6 +1399,27 @@ VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *
vecDV(vecDV_),
gradDV(gradDV_),
component(component_),
+ fct(NULL),
+ fac(fac_)
+{
+ TEST_EXIT(vecDV_)("No value vector!\n");
+ TEST_EXIT(gradDV_)("No gradient vector!\n");
+
+ auxFeSpaces.insert(vecDV_->getFeSpace());
+ auxFeSpaces.insert(gradDV_->getFeSpace());
+}
+
+
+VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_,
+ DOFVectorBase<double> *gradDV_,
+ int component_,
+ AbstractFunction<double, double>* fct_,
+ double fac_)
+: ZeroOrderTerm(2),
+ vecDV(vecDV_),
+ gradDV(gradDV_),
+ component(component_),
+ fct(fct_),
fac(fac_)
{
TEST_EXIT(vecDV_)("No value vector!\n");
@@ -1192,8 +1451,14 @@ void VecAndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo,
void VecAndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
- for (int iq = 0; iq < nPoints; iq++)
- C[iq] += vec[iq] * grad[iq][component] * fac;
+ if (fct == NULL) {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += vec[iq] * grad[iq][component] * fac;
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*fct)(vec[iq]) * grad[iq][component] * fac;
+ }
+
}
@@ -1205,42 +1470,73 @@ void VecAndPartialDerivative_ZOT::eval(int nPoints,
double opFactor)
{
double fac_ = fac * opFactor;
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] += vec[iq] * grad[iq][component] * uhAtQP[iq] * fac_;
+ if (fct == NULL) {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += vec[iq] * grad[iq][component] * uhAtQP[iq] * fac_;
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += (*fct)(vec[iq]) * grad[iq][component] * uhAtQP[iq] * fac_;
+ }
}
/* ----------------------------------------------------------- */
-Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, DOFVectorBase<double> *vec2DV_, DOFVectorBase<double> *gradDV_,int component_, double fac_)
- : ZeroOrderTerm(3), vec1DV(vec1DV_), vec2DV(vec2DV_), gradDV(gradDV_), fac(fac_), component(component_)
+Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_,
+ DOFVectorBase<double> *vec2DV_,
+ DOFVectorBase<double> *gradDV_,
+ int component_,
+ double fac_)
+ : ZeroOrderTerm(3), vec1DV(vec1DV_), vec2DV(vec2DV_), gradDV(gradDV_), fct(NULL), fac(fac_), component(component_)
{
- TEST_EXIT(vec1DV_)("No first vector!\n");
- TEST_EXIT(vec2DV_)("No second vector!\n");
- TEST_EXIT(gradDV_)("No gradient vector!\n");
- auxFeSpaces.insert(vec1DV_->getFeSpace());
- auxFeSpaces.insert(vec2DV_->getFeSpace());
- auxFeSpaces.insert(gradDV_->getFeSpace());
+ TEST_EXIT(vec1DV_)("No first vector!\n");
+ TEST_EXIT(vec2DV_)("No second vector!\n");
+ TEST_EXIT(gradDV_)("No gradient vector!\n");
+ auxFeSpaces.insert(vec1DV_->getFeSpace());
+ auxFeSpaces.insert(vec2DV_->getFeSpace());
+ auxFeSpaces.insert(gradDV_->getFeSpace());
+};
+
+Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_,
+ DOFVectorBase<double> *vec2DV_,
+ DOFVectorBase<double> *gradDV_,
+ int component_,
+ BinaryAbstractFunction<double, double, double>* fct_,
+ double fac_)
+ : ZeroOrderTerm(3), vec1DV(vec1DV_), vec2DV(vec2DV_), gradDV(gradDV_), fct(fct_), fac(fac_), component(component_)
+{
+ TEST_EXIT(vec1DV_)("No first vector!\n");
+ TEST_EXIT(vec2DV_)("No second vector!\n");
+ TEST_EXIT(gradDV_)("No gradient vector!\n");
+ auxFeSpaces.insert(vec1DV_->getFeSpace());
+ auxFeSpaces.insert(vec2DV_->getFeSpace());
+ auxFeSpaces.insert(gradDV_->getFeSpace());
};
+
void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
- getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
- getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
- getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
+ getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
+ getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
+ getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
};
void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
- getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
- getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
- getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
+ getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
+ getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
+ getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
};
void Vec2AndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
- for (int iq = 0; iq < nPoints; iq++)
- C[iq] += vec1[iq]*vec2[iq]*grad[iq][component]*fac;
+ if (fct == NULL) {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += vec1[iq]*vec2[iq]*grad[iq][component]*fac;
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*fct)(vec1[iq],vec2[iq])*grad[iq][component]*fac;
+ }
};
void Vec2AndPartialDerivative_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
@@ -1249,9 +1545,15 @@ void Vec2AndPartialDerivative_ZOT::eval(int nPoints,
mtl::dense_vector<double>& result,
double opFactor)
{
- for (int iq = 0; iq < nPoints; iq++) {
- result[iq] += (vec1[iq] * vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
- }
+ if (fct == NULL) {
+ for (int iq = 0; iq < nPoints; iq++) {
+ result[iq] += (vec1[iq] * vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
+ }
+ } else {
+ for (int iq = 0; iq < nPoints; iq++) {
+ result[iq] += (*fct)(vec1[iq], vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
+ }
+ }
};
/* ----------------------------------------------------------- */
diff --git a/extensions/POperators_FOT.h b/extensions/POperators_FOT.h
index af232352f38032ed3e175b2fdc75a649eb9a7659..c496b92837f494f28d53821cb9faec4ad0bf3e16 100644
--- a/extensions/POperators_FOT.h
+++ b/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;
+};
+
+
/* -------------------------------------------------------------- */
diff --git a/extensions/POperators_ZOT.h b/extensions/POperators_ZOT.h
index 6671c8133f2004da8d44219755efa33f5e040b2d..8a950b2af8d8eb898b2e1ba1c9000f4cc30e3029 100644
--- a/extensions/POperators_ZOT.h
+++ b/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;
};
/* -------------------------------------------------------------- */
diff --git a/extensions/SignedDistFunctors.h b/extensions/SignedDistFunctors.h
index 0d35e9c715adbae345fb8a1b589999685e6720eb..821e83f3571c66a8dee716f58f655a8ded84775d 100644
--- a/extensions/SignedDistFunctors.h
+++ b/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;
};
diff --git a/extensions/base_problems/BaseProblem_RB.h b/extensions/base_problems/BaseProblem_RB.h
index b395524d4cc020c8e7a3dfe6fac47f31379b4dd5..2ab0c40e9441c167ee637297a92cb9043fa62353 100644
--- a/extensions/base_problems/BaseProblem_RB.h
+++ b/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(); }
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.cc b/extensions/base_problems/CahnHilliardNavierStokes.cc
index c194ad5abbf7664c03653a67556b894d4e719d84..aabe82027c4aa1663af89794d93fd75861cd659a 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes.cc
+++ b/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);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
index 70a1b7be0896a9a7c4076892e003ae72d2535506..9a23aead05689a10c600ed2798ba7d5f7cff45c6 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
+++ b/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);
}
}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e77b6244f6505337765cdc557e96f13a5b878736
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
@@ -0,0 +1,436 @@
+#include "CahnHilliardNavierStokes_TwoPhase.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "POperators.h"
+
+using namespace AMDiS;
+
+CahnHilliardNavierStokes_TwoPhase::CahnHilliardNavierStokes_TwoPhase(const std::string &name_) :
+ super(name_),
+ useMobility(false),
+ useConservationForm(false),
+ doubleWell(1),
+ gamma(1.0),
+ eps(0.1),
+ minusEps(-0.1),
+ epsInv(10.0),
+ minusEpsInv(-10.0),
+ epsSqr(0.01),
+ minusEpsSqr(-0.01),
+ laplaceType(0),
+ nonLinTerm(3),
+ oldTimestep(0.0),
+ viscosity1(1.0),
+ viscosity2(1.0),
+ density1(1.0),
+ density2(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); // dt(c) = div(M(c)*grad(c)) + ...
+ Parameters::get(name + "->gamma", gamma); // M(c) = gamma * Double-Well
+ 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 + "->viscosity1", viscosity1);
+ Initfile::get(name + "->viscosity2", viscosity2);
+ Initfile::get(name + "->density1", density1);
+ Initfile::get(name + "->density2", density2);
+
+ // 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... u*grad(u_i^old) + u^old*grad(u_i) - u^old*grad(u_i^old)
+ Initfile::get(name + "->non-linear term", nonLinTerm);
+
+ // Parameters for CH-Coupling
+ // 0... u*grad(c), 1... -div(u*c)
+ Initfile::get(name + "->use conservation form", useConservationForm);
+
+ // 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;
+
+ MSG("surface tension: %e\n", surfaceTension);
+
+ 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_TwoPhase::~CahnHilliardNavierStokes_TwoPhase()
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::~CahnHilliardNavierStokes_TwoPhase()");
+
+ if (reinit != NULL) {
+ delete reinit;
+ reinit = NULL;
+ }
+ if (velocity != NULL) {
+ delete velocity;
+ velocity = NULL;
+ }
+ delete fileWriter;
+ fileWriter = NULL;
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::initData()
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::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_TwoPhase::transferInitialSolution(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::transferInitialSolution()");
+
+ calcVelocity();
+
+ for (int i = 0; i < 2+dow; i++)
+ prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
+
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+
+ // initial parameters for detecting mesh changes
+ oldMeshChangeIdx= getMesh()->getChangeIndex();
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::fillOperators()
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::fillOperators()");
+ MSG("CahnHilliardNavierStokes_TwoPhase::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);
+ prob->addMatrixOperator(*uGradC, 0, 0);
+ uGradC->setUhOld(prob->getSolution()->getDOFVector(0));
+ prob->addVectorOperator(*uGradC, 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 *opChWImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChWImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ (doubleWell == 0 ? static_cast<AbstractFunction<double, double>*>(new DoubleWell0Diff(-1.0))
+ : static_cast<AbstractFunction<double, double>*>(new DoubleWell1Diff(-1.0)))));
+ prob->addMatrixOperator(*opChWImpl, 1, 0);
+
+ /// < -2*phi*c_old^3 , psi >
+ // -2*c_old^3 + 3/2*c_old^2
+ Operator *opChWExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opChWExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new Pow3Functor(-2.0)));
+ if (doubleWell == 0) {
+ opChWExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new Pow2Functor(3.0/2.0)));
+ }
+ prob->addVectorOperator(*opChWExpl, 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) {
+ if (abs(density1 - density2) > FLT_TOL)
+ addMaterialDerivativeNonConst(i);
+ else
+ addMaterialDerivativeConst(i);
+
+ /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
+ if (abs(viscosity1 - viscosity2) > FLT_TOL)
+ addLaplaceTermNonConst(i);
+ else
+ addLaplaceTermConst(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
+ // -----------------------------
+ /// < mu_old * grad(c) , psi >
+ Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
+ opMuGradC->addTerm(new VecAndPartialDerivative_FOT(
+ prob->getSolution()->getDOFVector(1), i, -1.0), GRD_PHI);
+ prob->addMatrixOperator(*opMuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
+
+ /// < mu * grad(c_old) , psi >
+ Operator *opMuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
+ opMuGradC2->addTerm(new PartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(0), i, -1.0));
+ prob->addMatrixOperator(*opMuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
+
+ /// < mu_old * grad(c_old) , psi >
+ 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_TwoPhase::addMaterialDerivativeConst(int i)
+{
+ WorldVector<DOFVector<double>* > vel;
+ for (size_t j = 0; j < dow; j++)
+ vel[j] = prob->getSolution()->getDOFVector(2+j);
+
+ /// < (1/tau)*u'_i , psi >
+ Operator *opTime = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opTime->addTerm(new Simple_ZOT(density1));
+ 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));
+ opTimeOld->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(2+i), density1));
+ prob->addVectorOperator(*opTimeOld, 2+i, getInvTau(), getInvTau());
+
+ /// < u^old*grad(u_i^old) , psi >
+ if (nonLinTerm == 0 || nonLinTerm == 3) {
+ Operator *opUGradU0 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradU0->addTerm(new WorldVec_FOT(vel, (nonLinTerm == 0 ? -1.0 : 1.0)*density1),
+ GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(2+i));
+ prob->addVectorOperator(*opUGradU0, 2+i);
+ }
+ if (nonLinTerm == 1 || nonLinTerm == 3) {
+ /// < 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, density1));
+ prob->addMatrixOperator(*opUGradU1, 2+i, 2+j);
+ }
+ }
+ if (nonLinTerm == 2 || nonLinTerm == 3) {
+ /// < 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, density1), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, 2+i, 2+i);
+ }
+ }
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::addMaterialDerivativeNonConst(int i)
+{
+ WorldVector<DOFVector<double>* > vel;
+ for (size_t j = 0; j < dow; j++)
+ vel[j] = prob->getSolution()->getDOFVector(2+j);
+
+ /// < (1/tau)*u'_i , psi >
+ Operator *opTime = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opTime->addTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ new LinearInterpolation(density1, density2)));
+ prob->addMatrixOperator(*opTime, 2+i, 2+i, getInvTau(), getInvTau());
+ /// < (1/tau)*u_i^old , psi >
+ opTime->setUhOld(prob->getSolution()->getDOFVector(2+i));
+ prob->addVectorOperator(*opTime, 2+i, getInvTau(), getInvTau());
+
+ /// < (rho(c^old) + c^old*rho'(c^old)) * u^old*grad(u_i^old) , psi >
+ for(size_t j = 0; j < dow; ++j) {
+ Operator *opCUGradU = new Operator(getFeSpace(2+i),getFeSpace(2+i));
+ opCUGradU->addTerm(new Vec2AndPartialDerivative_FOT(
+ prob->getSolution()->getDOFVector(0),
+ prob->getSolution()->getDOFVector(2+j),
+ j,
+ new TaylorExpansionFactor(density1, density2)), GRD_PHI);
+ opCUGradU->setUhOld(prob->getSolution()->getDOFVector(2+i));
+ prob->addVectorOperator(*opCUGradU, 2+i);
+ }
+
+ /// < rho(c^old)*u*grad(u_i^old) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU1 = new Operator(getFeSpace(2+i),getFeSpace(2+j));
+ opUGradU1->addTerm(new VecAndPartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ prob->getSolution()->getDOFVector(2+i),
+ j,
+ new LinearInterpolation(density1, density2)));
+ prob->addMatrixOperator(*opUGradU1, 2+i, 2+j);
+ }
+
+ /// < rho(c^old)*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 Vec2AndPartialDerivative_FOT(
+ prob->getSolution()->getDOFVector(0),
+ prob->getSolution()->getDOFVector(2+j),
+ j,
+ new LinearInterpolationFactor(density1, density2)), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, 2+i, 2+i);
+ }
+
+ /// < c * rho'(c^old) * u^old*grad(u_i^old) , psi >
+ for(size_t j = 0; j < dow; ++j) {
+ Operator *opCUGradU2 = new Operator(getFeSpace(2+i),getFeSpace(0));
+ opCUGradU2->addTerm(new Vec2AndPartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(0),
+ prob->getSolution()->getDOFVector(2+j),
+ prob->getSolution()->getDOFVector(2+i),
+ j,
+ new LinearInterpolationDiffFactor(density1, density2)));
+ prob->addMatrixOperator(*opCUGradU2, 2+i, 0);
+ }
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermConst(int i)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::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, viscosity1));
+ 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(viscosity1));
+ prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::addLaplaceTermNonConst(int i)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::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 VecAtQP_IJ_SOT(
+ prob->getSolution()->getDOFVector(0),
+ new LinearInterpolation(viscosity1, viscosity2),
+ 1-i, 1-j));
+ 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 VecAtQP_SOT(
+ prob->getSolution()->getDOFVector(0),
+ new LinearInterpolation(viscosity1, viscosity2)));
+ prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase::closeTimestep()");
+
+ calcVelocity();
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h
new file mode 100644
index 0000000000000000000000000000000000000000..4fc241c5a717b7391c819504c0cdb3342cde67de
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase.h
@@ -0,0 +1,96 @@
+/** \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
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
index 38b0df01dc02b6741e6dd35fdc7549c8e6efe58e..d395c5abcc06b27f33aa63765f4ca737a0268efb 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
@@ -146,35 +146,43 @@ void CahnHilliardNavierStokes_TwoPhase_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));
@@ -260,34 +268,40 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
// Navier-Stokes part
// -----------------------------------------------------------------------------------
+ bool variableDensity = abs(density1 - density2) > DBL_TOL;
+ int shift = 2;
+ if (variableDensity)
+ shift = 3+dow;
+
for (size_t i = 0; i < dow; ++i) {
/// < d_t(u_i) , psi >
- prob->addTimeOperator(3+dow+i, 2+i);
-
- /// < z , psi >
- Operator *opZ = new Operator(getFeSpace(3+dow+i), getFeSpace(3+dow+i));
- opZ->addTerm(new Simple_ZOT);
- prob->addMatrixOperator(*opZ, 3+dow+i, 3+dow+i, &minus1, &minus1);
- opZ->setUhOld(prob->getStageSolution(3+dow+i));
- prob->addVectorOperator(*opZ, 3+dow+i);
-
+ prob->addTimeOperator(shift+i, 2+i, (variableDensity ? 1.0 : density1));
+
+ if (variableDensity) {
+ /// < z , psi >
+ Operator *opZ = new Operator(getFeSpace(shift+i), getFeSpace(shift+i));
+ opZ->addTerm(new Simple_ZOT);
+ prob->addMatrixOperator(*opZ, shift+i, shift+i, &minus1, &minus1);
+ opZ->setUhOld(prob->getStageSolution(shift+i));
+ prob->addVectorOperator(*opZ, shift+i);
+ }
/// < u_s*grad(u_s_i) , psi >
- Operator *opUGradU = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opUGradU->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PHI);
+ Operator *opUGradU = new Operator(getFeSpace(shift+i), getFeSpace(2+i));
+ opUGradU->addTerm(new WorldVec_FOT(stage_velocity, (variableDensity ? -1.0 : -density1)), GRD_PHI);
opUGradU->setUhOld(stage_velocity[i]);
- prob->addVectorOperator(*opUGradU, 3+dow+i);
+ prob->addVectorOperator(*opUGradU, shift+i);
/// < u_old*grad(u_i) , psi >
- Operator *opUGradV = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opUGradV->addTerm(new WorldVec_FOT(un_velocity, -1.0), GRD_PHI);
- prob->addMatrixOperator(*opUGradV, 3+dow+i, 2+i, &minus1, &minus1);
+ Operator *opUGradV = new Operator(getFeSpace(shift+i), getFeSpace(2+i));
+ opUGradV->addTerm(new WorldVec_FOT(un_velocity, (variableDensity ? -1.0 : -density1)), GRD_PHI);
+ prob->addMatrixOperator(*opUGradV, shift+i, 2+i, &minus1, &minus1);
/// < u*grad(u_old_i) , psi >
for (size_t j = 0; j < dow; ++j) {
- Operator *opVGradU = new Operator(getFeSpace(2+i), getFeSpace(2+j));
- opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, -1.0));
- prob->addMatrixOperator(*opVGradU, 3+dow+i, 2+j, &minus1, &minus1);
+ Operator *opVGradU = new Operator(getFeSpace(shift+i), getFeSpace(2+j));
+ opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, (variableDensity ? -1.0 : -density1)));
+ prob->addMatrixOperator(*opVGradU, shift+i, 2+j, &minus1, &minus1);
}
}
@@ -305,49 +319,52 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
// forces by Cahn-Hilliard terms
// -----------------------------
- /// < mu_old * grad(c) , psi >
- Operator *JmuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
- JmuGradC->addTerm(new VecAndPartialDerivative_FOT(
- prob->getUnVec(1), i, -1.0), GRD_PHI);
- prob->addMatrixOperator(*JmuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
-
- /// < mu * grad(c_old) , psi >
- Operator *JmuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
- JmuGradC2->addTerm(new PartialDerivative_ZOT(
- prob->getUnVec(0), i, -1.0));
- prob->addMatrixOperator(*JmuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
-
- /// < 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(0), i));
- prob->addVectorOperator(*opMuGradC, 2+i, &surfaceTension, &surfaceTension);
-
- // ---------------------------------------
-
- /// < rho(c_s)*z_s , psi >
- Operator *opDensityExpl = new Operator(prob->getFeSpace(2+i),prob->getFeSpace(3+dow+1));
- opDensityExpl->addZeroOrderTerm(new VecAtQP_ZOT(
- prob->getStageSolution(0),
- new LinearInterpolation(density1, density2, -1.0)));
- opDensityExpl->setUhOld(prob->getStageSolution(3+dow+i));
- prob->addVectorOperator(*opDensityExpl, 2+i);
-
- /// < rho(c_old)*z , psi >
- Operator *opDensityImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
- opDensityImpl->addZeroOrderTerm(new VecAtQP_ZOT(
- prob->getUnVec(0),
- new LinearInterpolation(density1, density2, -1.0)));
- prob->addMatrixOperator(*opDensityImpl, 2+i, 3+dow+i, &minus1, &minus1);
+ if (surfaceTension > DBL_TOL) {
+ /// < mu_old * grad(c) , psi >
+ Operator *JmuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
+ JmuGradC->addTerm(new VecAndPartialDerivative_FOT(
+ prob->getUnVec(1), i, -1.0), GRD_PHI);
+ prob->addMatrixOperator(*JmuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
+
+ /// < mu * grad(c_old) , psi >
+ Operator *JmuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
+ JmuGradC2->addTerm(new PartialDerivative_ZOT(
+ prob->getUnVec(0), i, -1.0));
+ prob->addMatrixOperator(*JmuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
+
+ /// < 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(0), i));
+ prob->addVectorOperator(*opMuGradC, 2+i, &surfaceTension, &surfaceTension);
+ }
- /// < rho'(c_old)*z_old * c , psi >
- Operator *opDensityImpl2 = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
- opDensityImpl2->addZeroOrderTerm(new Vec2AtQP_ZOT(
- prob->getUnVec(0),
- prob->getUnVec(3+dow+i),
- new LinearInterpolationDiffFactor(density1, density2, -1.0)));
- prob->addMatrixOperator(*opDensityImpl2, 2+i, 0, &minus1, &minus1);
+ // ---------------------------------------
+ if (variableDensity) {
+ /// < rho(c_s)*z_s , psi >
+ Operator *opDensityExpl = new Operator(prob->getFeSpace(2+i),prob->getFeSpace(shift+1));
+ opDensityExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getStageSolution(0),
+ new LinearInterpolation(density1, density2, -1.0)));
+ opDensityExpl->setUhOld(prob->getStageSolution(shift+i));
+ prob->addVectorOperator(*opDensityExpl, 2+i);
+
+ /// < rho(c_old)*z , psi >
+ Operator *opDensityImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(shift+i));
+ opDensityImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getUnVec(0),
+ new LinearInterpolation(density1, density2, -1.0)));
+ prob->addMatrixOperator(*opDensityImpl, 2+i, shift+i, &minus1, &minus1);
+
+ /// < rho'(c_old)*z_old * c , psi >
+ Operator *opDensityImpl2 = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opDensityImpl2->addZeroOrderTerm(new Vec2AtQP_ZOT(
+ prob->getUnVec(0),
+ prob->getUnVec(shift+i),
+ new LinearInterpolationDiffFactor(density1, density2, -1.0)));
+ prob->addMatrixOperator(*opDensityImpl2, 2+i, 0, &minus1, &minus1);
+ }
}
// div(u) = 0
@@ -411,53 +428,53 @@ void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
/// < nu*[grad(u)+grad(u)^t] , grad(psi) >
if (laplaceType == 1) {
for (unsigned j = 0; j < dow; ++j) {
- /// < nu(c_s)*d_j(u_s_i) , d_i(psi) >
+ /// < nu(c_s)*d_j(u_s_i) , d_i(psi) >
Operator *opLaplaceUijExpl = new Operator(getFeSpace(2+i), getFeSpace(2+j));
opLaplaceUijExpl->addTerm(new VecAtQP_IJ_SOT(
- prob->getStageSolution(0),
- new LinearInterpolation(viscosity1, viscosity2),
- 1-i, 1-j));
- opLaplaceUijExpl->setUhOld(prob->getStageSolution(2+j));
- prob->addVectorOperator(*opLaplaceUijExpl, 2+i, &minus1, &minus1);
+ prob->getStageSolution(0),
+ new LinearInterpolation(viscosity1, viscosity2),
+ 1-i, 1-j));
+ opLaplaceUijExpl->setUhOld(prob->getStageSolution(2+j));
+ prob->addVectorOperator(*opLaplaceUijExpl, 2+i, &minus1, &minus1);
- /// < nu(c_old)*d_j(u_i) , d_i(psi) >
+ /// < nu(c_old)*d_j(u_i) , d_i(psi) >
Operator *opLaplaceUijImpl = new Operator(getFeSpace(2+i), getFeSpace(2+j));
opLaplaceUijImpl->addTerm(new VecAtQP_IJ_SOT(
- prob->getUnVec(0),
- new LinearInterpolation(viscosity1, viscosity2),
- 1-i, 1-j));
+ prob->getUnVec(0),
+ new LinearInterpolation(viscosity1, viscosity2),
+ 1-i, 1-j));
prob->addMatrixOperator(*opLaplaceUijImpl, 2+i, 2+j);
- /// < c * nu'(c_old)*d_j(u_old_i) , d_i(psi) >
+ /// < c * nu'(c_old)*d_j(u_old_i) , d_i(psi) >
Operator *opLaplaceUijImpl2 = new Operator(getFeSpace(2+i), getFeSpace(2+j));
opLaplaceUijImpl2->addTerm(new VecAndPartialDerivativeIJ_FOT(
- prob->getUnVec(0),
- prob->getUnVec(2+j),
- 1-i, 1-j,
- new LinearInterpolationDiffFactor(viscosity1, viscosity2)), GRD_PSI);
+ prob->getUnVec(0),
+ prob->getUnVec(2+j),
+ 1-i, 1-j,
+ new LinearInterpolationDiffFactor(viscosity1, viscosity2)), GRD_PSI);
prob->addMatrixOperator(*opLaplaceUijImpl2, 2+i, 0);
}
}
/// < nu*grad(u'_i) , grad(psi) >
Operator *opLaplaceUExpl = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opLaplaceUExpl->addTerm(new VecAtQP_SOT(
- prob->getStageSolution(0),
- new LinearInterpolation(viscosity1, viscosity2)));
+ opLaplaceUExpl->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new LinearInterpolation(viscosity1, viscosity2)));
opLaplaceUExpl->setUhOld(prob->getStageSolution(2+i));
prob->addVectorOperator(*opLaplaceUExpl, 2+i, &minus1, &minus1);
Operator *opLaplaceUImpl = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opLaplaceUImpl->addTerm(new VecAtQP_SOT(
- prob->getUnVec(0),
- new LinearInterpolation(viscosity1, viscosity2)));
+ opLaplaceUImpl->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new LinearInterpolation(viscosity1, viscosity2)));
prob->addMatrixOperator(*opLaplaceUImpl, 2+i, 2+i);
Operator *opLaplaceUImpl2 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opLaplaceUImpl2->addTerm(new VecGrad_FOT(
- prob->getUnVec(0),
- prob->getUnVec(2+i),
- new LinearInterpolationDiffVector(viscosity1, viscosity2)), GRD_PSI);
+ opLaplaceUImpl2->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(2+i),
+ new LinearInterpolationDiffVector(viscosity1, viscosity2)), GRD_PSI);
prob->addMatrixOperator(*opLaplaceUImpl2, 2+i, 0);
/// < p , d_i(psi) >
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
index a4ccf45031be5cf4d6a7e4e566c93ff684478705..94030c42c8b9d5e850ba153e5ac4a8369d13399a 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
+++ b/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;
-};
-
-
-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;
-};
-
-struct LinearInterpolation : AbstractFunction<double,double>
-{
- LinearInterpolation(double v0_, double v1_, double factor_ = 1.0) :
- AbstractFunction<double,double>(4),
- v0(v0_), v1(v1_), factor(factor_) { }
-
- double operator()(const double &c) const
- {
- return factor * (v0*(c+1.0) - v1*(c-1.0))/2.0;
- }
-
-private:
- double v0;
- double v1;
- double factor;
-};
-
-struct LinearInterpolationDiffFactor : BinaryAbstractFunction<double,double,double>
-{
- LinearInterpolationDiffFactor(double v0_, double v1_, double factor_ = 1.0) :
- BinaryAbstractFunction<double,double,double>(4),
- v0(v0_), v1(v1_), factor(factor_) { }
-
- double operator()(const double &c, const double& z) const
- {
- return factor * z * (v0 - v1)/2.0;
- }
-
-private:
- double v0;
- double v1;
- double factor;
-};
-
-struct LinearInterpolationDiffVector : BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
-{
- LinearInterpolationDiffVector(double v0_, double v1_, double factor_ = 1.0) :
- BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
- v0(v0_), v1(v1_), factor(factor_) { }
-
- WorldVector<double> operator()(const double &c, const WorldVector<double>& grdU) const
- {
- return (factor * (v0 - v1)/2.0) * grdU;
- }
-
-private:
- double v0;
- double v1;
- double factor;
-};
-
#endif // CAHN_HILLIARD_H
diff --git a/extensions/base_problems/CahnHilliard_RB.cc b/extensions/base_problems/CahnHilliard_RB.cc
index 94cc4798c328fd827562ba340327ee696f99b30c..7b3e839d1443bccc2295c99bcee600c2c1830c01 100644
--- a/extensions/base_problems/CahnHilliard_RB.cc
+++ b/extensions/base_problems/CahnHilliard_RB.cc
@@ -141,35 +141,43 @@ void CahnHilliard_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));
diff --git a/extensions/base_problems/NavierStokes_TaylorHood_RB.h b/extensions/base_problems/NavierStokes_TaylorHood_RB.h
index 56834b4fa24458dcaa8dafd4a226d4f50cabea60..1e471752f8d509af5d2fd8c33074a74fb1188662 100644
--- a/extensions/base_problems/NavierStokes_TaylorHood_RB.h
+++ b/extensions/base_problems/NavierStokes_TaylorHood_RB.h
@@ -7,7 +7,7 @@
#include "Helpers.h"
#include "POperators.h"
#include "BoundaryFunctions.h"
-#include "BaseProblem.h"
+#include "BaseProblem_RB.h"
#include "time/ExtendedRosenbrockStationary.h"
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
@@ -24,11 +24,11 @@ using namespace AMDiS;
* \brief
* Navier-Stokes problem, using Taylor Hood elements
*/
-class NavierStokes_TaylorHood_RB : public BaseProblem<ExtendedRosenbrockStationary>
+class NavierStokes_TaylorHood_RB : public BaseProblem_RB
{
public: // typedefs
- typedef BaseProblem<ExtendedRosenbrockStationary> super;
+ typedef BaseProblem_RB super;
public: // methods
@@ -58,23 +58,6 @@ public: // methods
return velocity;
}
-// Rosenbrock methods
-//_________________________________________________________________________
-
-public:
-
- 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: // methods
virtual void fillOperators();
diff --git a/extensions/base_problems/chns.h b/extensions/base_problems/chns.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc6cf319833f723335bbc4ca7fd88297b39b8a56
--- /dev/null
+++ b/extensions/base_problems/chns.h
@@ -0,0 +1,276 @@
+/** \file chns.h */
+
+#ifndef CHNS_FUNCTORS_H
+#define CHNS_FUNCTORS_H
+
+#include "AMDiS.h"
+
+/** \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;
+};
+
+
+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;
+};
+
+struct LinearInterpolation : AbstractFunction<double,double>
+{
+ LinearInterpolation(double v0_, double v1_, double factor_ = 1.0) :
+ AbstractFunction<double,double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (v0*(c+1.0) - v1*(c-1.0))/2.0;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct LinearInterpolationFactor : BinaryAbstractFunction<double, double, double>
+{
+ LinearInterpolationFactor(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<double, double, double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c, const double& u) const
+ {
+ return factor * u * (v0*(c+1.0) - v1*(c-1.0))/2.0;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct LinearInterpolationDiffFactor : BinaryAbstractFunction<double,double,double>
+{
+ LinearInterpolationDiffFactor(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<double,double,double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c, const double& z) const
+ {
+ return factor * z * (v0 - v1)/2.0;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct LinearInterpolationDiffVector : BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ LinearInterpolationDiffVector(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdU) const
+ {
+ return (factor * (v0 - v1)/2.0) * grdU;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct TaylorExpansionFactor : BinaryAbstractFunction<double,double,double>
+{
+ TaylorExpansionFactor(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<double,double,double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c, const double& u) const
+ {
+ return LinearInterpolationFactor(v0, v1, factor)(c, u) + c * LinearInterpolationDiffFactor(v0, v1, factor)(c, u);
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+#endif // CHNS_FUNCTORS_H
diff --git a/extensions/demo/CMakeLists.txt b/extensions/demo/CMakeLists.txt
index 074280e7a7696d784a9efc0acf00ec8b9eda67b1..d6f1a931ee54376bbf758453905beda9a25b9e42 100644
--- a/extensions/demo/CMakeLists.txt
+++ b/extensions/demo/CMakeLists.txt
@@ -10,47 +10,36 @@ if(AMDIS_FOUND)
SET(BASIS_LIBS ${AMDIS_LIBRARIES})
endif(AMDIS_FOUND)
-
-# set(elliptBase src/elliptBaseProblem.cc)
-# add_executable("elliptBase" ${elliptBase})
-# target_link_libraries("elliptBase" ${BASIS_LIBS})
-#
-# set(elliptImplicit src/elliptImplicit.cc)
-# add_executable("elliptImplicit" ${elliptImplicit})
-# target_link_libraries("elliptImplicit" ${BASIS_LIBS})
-#
-# set(cahn_hilliard src/cahn_hilliard.cc)
-# add_executable("cahn_hilliard" ${cahn_hilliard})
-# target_link_libraries("cahn_hilliard" ${BASIS_LIBS})
-#
-# set(pfc src/pfc.cc)
-# add_executable("pfc" ${pfc})
-# target_link_libraries("pfc" ${BASIS_LIBS})
-#
-# set(drivenCavity src/chns/drivenCavity.cc)
-# add_executable("drivenCavity" ${drivenCavity})
-# target_link_libraries("drivenCavity" ${BASIS_LIBS})
-#
-# include_directories(/home/spraetor/projects/src/common)
-# set(navierStokes src/navierStokes.cc)
-# add_executable("navierStokes" ${navierStokes})
-# target_link_libraries("navierStokes" ${BASIS_LIBS})
-#
-# set(navierStokesDd src/navierStokes_diffuseDomain.cc)
-# add_executable("navierStokesDd" ${navierStokesDd})
-# target_link_libraries("navierStokesDd" ${BASIS_LIBS})
-#
-# set(linearElasticityDd src/linearElasticity.cc)
-# add_executable("linearElasticityDd" ${linearElasticityDd})
-# target_link_libraries("linearElasticityDd" ${BASIS_LIBS})
-
set(fsi src/fsi_explicit/fluidStructureInteraction.cc)
add_executable("fsi" ${fsi})
target_link_libraries("fsi" ${BASIS_LIBS})
+ set_target_properties("fsi" PROPERTIES COMPILE_FLAGS "-DUSE_KD_TREE")
+
+ set(drivenCavity src/drivenCavity.cc)
+ add_executable("drivenCavity" ${drivenCavity})
+ target_link_libraries("drivenCavity" ${BASIS_LIBS})
+
+ include_directories(/opt/amdis/amdis-src/extensions/base_problems/)
+ set(drivenCavity_rb src/drivenCavity_rb.cc /opt/amdis/amdis-src/extensions/base_problems/CahnHilliardNavierStokes_RB.cc)
+ add_executable("drivenCavity_rb" ${drivenCavity_rb})
+ target_link_libraries("drivenCavity_rb" ${BASIS_LIBS})
-# set(ddFsi src/diffuseDomainFsi.cc)
-# add_executable("ddFsi" ${ddFsi})
-# target_link_libraries("ddFsi" ${BASIS_LIBS})
+# include_directories(/opt/amdis/amdis-src/extensions/base_problems/)
+ set(drivenCavity_twophase src/drivenCavity_twophase_rb.cc /opt/amdis/amdis-src/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc)
+ add_executable("drivenCavity_twophase" ${drivenCavity_twophase})
+ target_link_libraries("drivenCavity_twophase" ${BASIS_LIBS})
+
+ set(ddFsi src/diffuseDomainFsi.cc)
+ add_executable("ddFsi" ${ddFsi})
+ target_link_libraries("ddFsi" ${BASIS_LIBS})
+
+ set(movingMesh src/movingMesh.cc)
+ add_executable("movingMesh" ${movingMesh})
+ target_link_libraries("movingMesh" ${BASIS_LIBS})
+
+ set(rosenbrockTest src/rosenbrock_test.cc)
+ add_executable("rosenbrockTest" ${rosenbrockTest})
+ target_link_libraries("rosenbrockTest" ${BASIS_LIBS})
#create the output dir
diff --git a/extensions/demo/init/drivenCavity2.dat.2d b/extensions/demo/init/drivenCavity2.dat.2d
index 951831594610d3e1c3b0b5fa7ee05c5239fdee49..d12fd8e428529dc860819c5df36247405c2fe44a 100644
--- a/extensions/demo/init/drivenCavity2.dat.2d
+++ b/extensions/demo/init/drivenCavity2.dat.2d
@@ -66,12 +66,12 @@ chns->non-linear term: 3 % 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2.
% ====================== USER_PARAMETER - COUPLING ==============
-chns->sigma: 0.0 % surface tension
+chns->sigma: 1 % surface tension
chns->use conservation form: 0
% ====================== TIMESTEPS ========================
-adapt->rosenbrock method: rodasp
+adapt->rosenbrock method: rosI2Pw
adapt->fix first timesteps: 0
adapt->rosenbrock->timestep study: 0
adapt->rosenbrock->timestep study steps: 0
@@ -87,7 +87,7 @@ adapt->timestep: 1.e-2
adapt->max timestep: 1e+10
adapt->min timestep: 1e-6
adapt->start time: 0.0
-adapt->end time: 10
+adapt->end time: 9999999
% =========== OUTPUT ==============================================
diff --git a/extensions/doc/doxyfile b/extensions/doc/doxyfile
index f4e5b24aeed3455ceeebed365f0204a6a2cbf75e..458fde2f83f63d64296727915205bb69da8b1fa4 100644
--- a/extensions/doc/doxyfile
+++ b/extensions/doc/doxyfile
@@ -1103,7 +1103,7 @@ EXT_LINKS_IN_WINDOW = NO
# to manually remove any form_*.png images from the HTML output directory
# to force them to be regenerated.
-FORMULA_FONTSIZE = 10
+FORMULA_FONTSIZE = 11
# Use the FORMULA_TRANPARENT tag to determine whether or not the images
# generated for formulas are transparent PNGs. Transparent PNGs are
diff --git a/extensions/time/ExtendedRosenbrockAdaptInstationary.h b/extensions/time/ExtendedRosenbrockAdaptInstationary.h
index d41c6e5d241f6e778ca32b7e8e4636ebb2e6be65..bb20ab87c6a8d274c98c21581ec114b5a4a053ee 100644
--- a/extensions/time/ExtendedRosenbrockAdaptInstationary.h
+++ b/extensions/time/ExtendedRosenbrockAdaptInstationary.h
@@ -110,6 +110,8 @@ namespace AMDiS {
/// The value tau * gamma, where gamma is a value of the used
/// Rosenbrock method.
double tauGamma, minusTauGamma, invTauGamma, minusInvTauGamma;
+
+ double oldTime;
/// If true, the first timestep is calculated with different timesteps.
/// This is usually used to make a study how the time error estimator
diff --git a/extensions/time/ExtendedRosenbrockAdaptInstationary.hh b/extensions/time/ExtendedRosenbrockAdaptInstationary.hh
index 7a4192fdd84eba96ae1095552193922c600ec603..f338e886f393de6ae1bfa0cd39bb8e0e16c540e9 100644
--- a/extensions/time/ExtendedRosenbrockAdaptInstationary.hh
+++ b/extensions/time/ExtendedRosenbrockAdaptInstationary.hh
@@ -35,6 +35,7 @@ namespace AMDiS {
minusTauGamma(-1.0),
invTauGamma(1.0),
minusInvTauGamma(-1.0),
+ oldTime(adaptInfo->getStartTime()),
dbgTimestepStudy(false)
{
FUNCNAME("ExtendedRosenbrockAdaptInstationary::ExtendedRosenbrockAdaptInstationary()");
@@ -87,7 +88,8 @@ namespace AMDiS {
adaptInfo->setRosenbrockMode(true);
rosenbrockStat->setTauGamma(&tauGamma, &minusTauGamma,
- &invTauGamma, &minusInvTauGamma);
+ &invTauGamma, &minusInvTauGamma,
+ &oldTime);
rosenbrockStat->setTau(&tau);
Parameters::get(name + "->rosenbrock->timestep study", dbgTimestepStudy);
@@ -125,6 +127,7 @@ namespace AMDiS {
minusTauGamma = -1.0;
invTauGamma = 1.0;
minusInvTauGamma = -1.0;
+ oldTime = 0.0;
}
@@ -265,6 +268,8 @@ namespace AMDiS {
(dbgTimestepStudy && (studyTimestep + 1 < dbgTimesteps.size())));
rosenbrockStat->acceptTimestep();
+ oldTime = adaptInfo->getTime();
+
adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep());
adaptInfo->incTimestepNumber();
diff --git a/extensions/time/ExtendedRosenbrockStationary.cc b/extensions/time/ExtendedRosenbrockStationary.cc
index e8c5db1028c6aeebfe93d6aa25b97f88ffd15765..1fb320792197f1313cbe9289a97e3b81159990e8 100644
--- a/extensions/time/ExtendedRosenbrockStationary.cc
+++ b/extensions/time/ExtendedRosenbrockStationary.cc
@@ -16,6 +16,7 @@
#include "SystemVector.h"
#include "OEMSolver.h"
#include "Debug.h"
+#include "POperators_ZOT.h"
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
@@ -65,7 +66,11 @@ using namespace AMDiS;
*newUn = *unVec;
*lowSol = *unVec;
- for (int i = 0; i < rm->getStages(); i++) {
+ for (int i = 0; i < rm->getStages(); i++) {
+
+ stageTime = (*oldTime) + rm->getAlphaI(i) * (*tauPtr);
+ tauGammaI = rm->getGammaI(i) * (*tauPtr);
+
*stageSolution = *unVec;
for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]);
@@ -76,6 +81,13 @@ using namespace AMDiS;
for (unsigned int j = 0; j < boundaries.size(); j++) {
boundaries[j].vec->interpol(boundaries[j].fct);
*(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].row));
+ if (boundaries[j].fctDt != NULL) {
+ DOFVector<double>* tmpDt = new DOFVector<double>(getFeSpace(boundaries[j].row), "tmp");
+ tmpDt->interpol(boundaries[j].fctDt);
+ *tmpDt *= tauGammaI;
+ *(boundaries[j].vec) -= *tmpDt;
+ delete tmpDt;
+ }
}
timeRhsVec->set(0.0);
@@ -100,6 +112,8 @@ using namespace AMDiS;
*lowSol += *tmp;
}
+ stageTime = (*oldTime) + (*tauPtr);
+
for (int i = 0; i < nComponents; i++) {
(*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i+componentShift);
@@ -135,40 +149,41 @@ using namespace AMDiS;
}
- void ExtendedRosenbrockStationary::addTimeOperator(int row, int col)
+ void ExtendedRosenbrockStationary::addTimeOperator(int row, int col, double factor)
{
FUNCNAME("ExtendedRosenbrockStationary::addTimeOperator()");
TEST_EXIT(invTauGamma)("This should not happen!\n");
Operator *op = new Operator(componentSpaces[row], componentSpaces[col]);
- op->addZeroOrderTerm(new Simple_ZOT);
+ op->addZeroOrderTerm(new Simple_ZOT(factor));
super::addMatrixOperator(op, row, col, invTauGamma, invTauGamma);
Operator *opRhs = new Operator(componentSpaces[row]);
- opRhs->addZeroOrderTerm(new VecAtQP_ZOT(timeRhsVec->getDOFVector(col)));
- super::addVectorOperator(opRhs, row, &minusOne, &minusOne);
+ opRhs->addZeroOrderTerm(new Phase_ZOT(timeRhsVec->getDOFVector(col), factor));
+ super::addVectorOperator(opRhs, row);
}
void ExtendedRosenbrockStationary::addDirichletBC(BoundaryType type, int row, int col,
- AbstractFunction<double, WorldVector<double> > *fct)
+ AbstractFunction<double, WorldVector<double> > *fct,
+ AbstractFunction<double, WorldVector<double> > *fctDt)
{
FUNCNAME("ExtendedRosenbrockStationary::addDirichletBC()");
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {fct, vec, row, col};
+ MyRosenbrockBoundary bound(fct, fctDt, vec, row, col);
boundaries.push_back(bound);
super::addDirichletBC(type, row, col, vec);
}
-
+
void ExtendedRosenbrockStationary::addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addSingularDirichletBC(pos, row, col, *vec);
@@ -179,7 +194,7 @@ using namespace AMDiS;
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addSingularDirichletBC(idx, row, col, *vec);
@@ -190,7 +205,7 @@ using namespace AMDiS;
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addImplicitDirichletBC(signedDist, row, col, *vec);
@@ -201,7 +216,7 @@ using namespace AMDiS;
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addImplicitDirichletBC(signedDist, row, col, *vec);
@@ -213,7 +228,7 @@ using namespace AMDiS;
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addManualDirichletBC(nr, meshIndicator, row, col, *vec);
@@ -225,7 +240,7 @@ using namespace AMDiS;
AbstractFunction<double, WorldVector<double> > &fct)
{
DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
- MyRosenbrockBoundary bound = {&fct, vec, row, col};
+ MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
super::addManualDirichletBC(meshIndicator, row, col, *vec);
diff --git a/extensions/time/ExtendedRosenbrockStationary.h b/extensions/time/ExtendedRosenbrockStationary.h
index 181eed911fd2b23bb3b6ad9a3210f49149ec3095..50c1994d689965bf78eec0c38b89f4a44271ddc4 100644
--- a/extensions/time/ExtendedRosenbrockStationary.h
+++ b/extensions/time/ExtendedRosenbrockStationary.h
@@ -40,7 +40,15 @@ using namespace AMDiS;
struct MyRosenbrockBoundary
{
- AbstractFunction<double, WorldVector<double> > *fct;
+ MyRosenbrockBoundary(AbstractFunction<double, WorldVector<double> > * fct_,
+ AbstractFunction<double, WorldVector<double> > * fctDt_,
+ DOFVector<double> *vec_,
+ int row_,
+ int col_
+ ) : fct(fct_), fctDt(fctDt_), vec(vec_), row(row_), col(col_) {}
+
+ AbstractFunction<double, WorldVector<double> > *fct; // f(t,x)
+ AbstractFunction<double, WorldVector<double> > *fctDt; // dt[f](t,x)
DOFVector<double> *vec;
@@ -65,7 +73,8 @@ using namespace AMDiS;
tauGamma(NULL),
minusTauGamma(NULL),
invTauGamma(NULL),
- minusInvTauGamma(NULL)
+ minusInvTauGamma(NULL),
+ stageTime(0.0)
{}
DOFVector<double>* getUnVec(int i)
@@ -111,20 +120,21 @@ using namespace AMDiS;
void addJacobianOperator(Operator &op, int row, int col,
double *factor = NULL, double *estFactor = NULL);
- void addTimeOperator(int i, int j);
+ void addTimeOperator(int i, int j, double factor = 1.0);
void setTau(double *ptr)
{
tauPtr = ptr;
}
- void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3)
+ void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3, double *ptr4)
{
tauGamma = ptr0;
minusTauGamma = ptr1;
invTauGamma = ptr2;
minusInvTauGamma = ptr3;
- }
+ oldTime = ptr4;
+ }
double* getTauGamma()
{
@@ -145,11 +155,33 @@ using namespace AMDiS;
{
return minusInvTauGamma;
}
+
+ double getStageTime()
+ {
+ return stageTime;
+ }
+
+ double* getOldTimePtr()
+ {
+ return oldTime;
+ }
+
+ double* getStageTimePtr()
+ {
+ return &stageTime;
+ }
+
+
+ double* getTauGammaI()
+ {
+ return &tauGammaI;
+ }
/// Adds a Dirichlet boundary condition, where the rhs is given by an
/// abstract function.
void addDirichletBC(BoundaryType type, int row, int col,
- AbstractFunction<double, WorldVector<double> > *b);
+ AbstractFunction<double, WorldVector<double> > *b,
+ AbstractFunction<double, WorldVector<double> > *bDt = NULL);
/// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
/// vector.
@@ -163,33 +195,33 @@ using namespace AMDiS;
/// Adds a Neumann boundary condition, where the rhs is given by an
/// abstract function.
- void addNeumannBC(BoundaryType type, int row, int col,
- AbstractFunction<double, WorldVector<double> > *n)
- {
- FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
-
- ERROR_EXIT("Not yet supported!\n");
- }
+// void addNeumannBC(BoundaryType type, int row, int col,
+// AbstractFunction<double, WorldVector<double> > *n)
+// {
+// FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
+//
+// ERROR_EXIT("Not yet supported!\n");
+// }
/// Adds a Neumann boundary condition, where the rhs is given by an DOF
/// vector.
- void addNeumannBC(BoundaryType type, int row, int col,
- DOFVector<double> *n)
- {
- FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
-
- ERROR_EXIT("Not yet supported!\n");
- }
+// void addNeumannBC(BoundaryType type, int row, int col,
+// DOFVector<double> *n)
+// {
+// FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
+//
+// ERROR_EXIT("Not yet supported!\n");
+// }
/// Adds Robin boundary condition.
- void addRobinBC(BoundaryType type, int row, int col,
- AbstractFunction<double, WorldVector<double> > *n,
- AbstractFunction<double, WorldVector<double> > *r)
- {
- FUNCNAME("ExtendedRosenbrockStationary::addRobinBC()");
-
- ERROR_EXIT("Not yet supported!\n");
- }
+// void addRobinBC(BoundaryType type, int row, int col,
+// AbstractFunction<double, WorldVector<double> > *n,
+// AbstractFunction<double, WorldVector<double> > *r)
+// {
+// FUNCNAME("ExtendedRosenbrockStationary::addRobinBC()");
+//
+// ERROR_EXIT("Not yet supported!\n");
+// }
/// special boundary conditions defined in ExtendedProblemStat
@@ -229,7 +261,7 @@ using namespace AMDiS;
protected:
RosenbrockMethod *rm;
- SystemVector *stageSolution, *unVec, *timeRhsVec, *newUn, *tmp, *lowSol;
+ SystemVector *stageSolution, *unVec, *timeRhsVec, *newUn, *tmp, *lowSol, *timeDrivative;
std::vector<SystemVector*> stageSolutions;
@@ -242,6 +274,10 @@ using namespace AMDiS;
double *tauPtr;
double *tauGamma, *minusTauGamma, *invTauGamma, *minusInvTauGamma;
+
+ double stageTime; // t_n + alpha_i*tau
+ double* oldTime;
+ double tauGammaI; // tau*gamma_i
std::vector<MyRosenbrockBoundary> boundaries;
};
diff --git a/extensions/tutorial.hpp b/extensions/tutorial.hpp
index 4ee32d5a911d93c277cada13fcce37b2ba226e3e..cb161cb88a09bab0729b514a6ddde6fd295005ff 100644
--- a/extensions/tutorial.hpp
+++ b/extensions/tutorial.hpp
@@ -28,7 +28,7 @@ But it should give you enough information to get started.
- \subpage mesh_refinement
- \subpage boundary_conditions
- \subpage tools
-
+- \subpage base_problems
*/
@@ -940,6 +940,116 @@ int main(int argc, char* argv[])
\image html elliptManual.png "Manual Dirichlet condition at the left boundary of the domain" width=300px
+*/
+
+//-----------------------------------------------------------
+
+
+/*! \page base_problems Base Problems
+
+Some common problems are implemented in the directory base_problems:
+- \subpage cahn_hilliard
+- \subpage phase_field_crystal
+- \subpage navier_stokes
+- \subpage chns
+- \subpage chns2
+- \subpage linear_elasticity
+- \subpage fsi
+
+*/
+
+//-----------------------------------------------------------
+
+
+/*! \page fsi Fully coupled fluid structure interaction using diffuse domain approach
+
+The underliing equations are
+
+\f{eqnarray*}{
+ \partial_t u_j + \mathbf{u}\cdot\nabla u_j - \nabla\cdot(\phi\cdot\sigma_f) - \nabla\cdot((1-\phi)\cdot\sigma_s) &=& 0,\quad j=0\ldots d-1 \\
+ \nabla\cdot(\phi\cdot\mathbf{u}) - \nabla\phi \cdot \mathbf{u} + (1-\phi)(p - \alpha\nabla\cdot\eta) &=& 0 \\
+ (1-\phi)\cdot(\partial_t\eta_j + \mathbf{u}\cdot\nabla\eta_j) &=& (1-\phi)\cdot u_j,\quad j=0\ldots d-1
+\f}
+
+with \f$ \alpha = \lambda + \frac{2}{3}\mu,\; \sigma_s = \mu(\nabla\eta+\nabla\eta^\top)+\lambda I\nabla\cdot\eta
+\f$ and \f$ \sigma_f = \nu(\nabla \mathbf{u} + \nabla\mathbf{u}^\top)-Ip \f$.
+
+*/
+
+//-----------------------------------------------------------
+
+
+/*! \page chns Cahn-Hilliard-Navier-Stokes equation
+
+Two phase flow simulations are often implemented using by coupling the Cahn-Hilliard equations and the
+Navier-Stokes equations. This coupling can be done explicitely, i.e. solve the equations one after the other,
+or implicitely as one system of equations. This is implemented in the CahnHillirdNavierStokes base problem:
+
+\f{eqnarray*}{
+ \partial_t u_j + \mathbf{u}\cdot\nabla u_j - \nabla\cdot\sigma &=& \mu\partial_j c,\quad j=0\ldots d-1 \\
+ \nabla\cdot \mathbf{u} &=& 0 \\
+ \partial_t c - \Delta\mu + \mathbf{u}\cdot\nabla c &=& 0 \\
+ -\epsilon^2\Delta c + W'(c) &=& \mu
+\f}
+
+with \f$ \sigma = Ip - \nu\nabla u \f$ and \f$ W(c) = \frac{1}{4}c^2(1-c)^2 \f$.
+
+For the time-discretization the nonlinear terms are linearized around the old timestep:
+
+\f{eqnarray*}{
+ \frac{1}{\tau} (u^{k+1}_j - u^{k}_j) + \mathbf{u}^{k+1}\cdot\nabla u^k_j + \mathbf{u}^{k}\cdot\nabla u^{k+1}_j - \mathbf{u}^{k}\cdot\nabla u^{k}_j - \nabla\cdot(\sigma^{k+1}) &=& \mu^{k+1}\partial_j c^k + \mu^{k}\partial_j c^{k+1} - \mu^{k}\partial_j c^k,\quad j=0\ldots d-1 \\
+ \nabla\cdot \mathbf{u}^{k+1} &=& 0 \\
+ \frac{1}{\tau} (c^{k+1} - c^k) - \Delta\mu^{k+1} &=& - \mathbf{u}^{k+1}\cdot\nabla c^k - \mathbf{u}^{k}\cdot\nabla c^{k+1} + \mathbf{u}^{k}\cdot\nabla c^k \\
+ \mu^{k+1} &=& -\epsilon^2\Delta c^{k+1} + W'(c^k) + W''(c^k)(c^{k+1}-c^k)
+\f}
+
+*/
+
+//-----------------------------------------------------------
+
+
+/*! \page chns2 Cahn-Hilliard-Navier-Stokes equation (variable density/viscosity)
+
+Two phase flow simulations are often implemented using by coupling the Cahn-Hilliard equations and the
+Navier-Stokes equations. This coupling can be done explicitely, i.e. solve the equations one after the other,
+or implicitely as one system of equations. This is implemented in the CahnHillirdNavierStokes base problem:
+
+\f{eqnarray*}{
+ \partial_t c + \mathbf{u}\cdot\nabla c &=& \nabla\cdot \gamma W(c) \nabla \mu \\
+ 0 &=& \mu + \epsilon^2\Delta c - W'(c) \\
+ \rho(c)\partial_t u_j + \rho(c)\mathbf{u}\cdot\nabla u_j &=& \nabla\cdot\nu(c)\nabla u_j - \nabla p + \alpha\mu\partial_j c + F_j,\quad j=0\ldots d-1 \\
+ 0 &=& \nabla\cdot \mathbf{u}
+\f}
+
+with \f$ W(c) = \frac{1}{4}c^2(1-c)^2 \f$ and \f$ \alpha \f$ the surface tension.
+
+\section chns2_euler Semi-implicit Euler discretization
+
+For the time-discretization a semi-implicit euler discretization is chosen and the nonlinear terms are linearized around the old timestep:
+
+\f{eqnarray*}{
+ \frac{1}{\tau} (c^{k+1} - c^k) - \Delta\mu^{k+1} &=& - \mathbf{u}^{k+1}\cdot\nabla c^k - \mathbf{u}^{k}\cdot\nabla c^{k+1} + \mathbf{u}^{k}\cdot\nabla c^k \\
+ \mu^{k+1} &=& -\epsilon^2\Delta c^{k+1} + W'(c^k) + W''(c^k)(c^{k+1}-c^k) \\
+ \rho(c^k)\frac{1}{\tau} (u^{k+1}_j - u^{k}_j) + [\rho(c)\mathbf{u}\cdot\nabla u_j] &=& \nabla\cdot(\nu(c^k)\nabla u_i^{k+1}) - \nabla p^{k+1} + \alpha(\mu^{k+1}\partial_j c^k + \mu^{k}\partial_j c^{k+1} - \mu^{k}\partial_j c^k) + F_i,\quad j=0\ldots d-1 \\
+ \nabla\cdot \mathbf{u}^{k+1} &=& 0
+\f}
+
+with \f[ [\rho(c)\mathbf{u}\cdot\nabla u_j] := \rho(c^k)\cdot\mathbf{u}^{k+1}\cdot\nabla u^k_j + \rho(c^k)\cdot\mathbf{u}^{k}\cdot\nabla u^{k+1}_j + c^{k+1}\cdot\rho'(c^k)\cdot\mathbf{u}^{k}\cdot\nabla u^{k}_j - \big(\rho(c^k) + c^k\cdot\rho'(c^k)\big)\cdot\mathbf{u}^{k}\cdot\nabla u^{k}_j \f].
+
+\section chns2_rb Rosenbrock Discretization
+
+For the time-discretization a Rosenbrock scheme is chosen.
+Therefor we have to introduce a new variable \f$ z_j = \partial_t u_j + \mathbf{u}\cdot\nabla u_j \f$.
+Stage solutions are depicted with supscripts \f$ \cdot^s \f$ and old solution with \f$ \cdot^{old} \f$:
+
+\f{eqnarray*}{
+ \partial^s_t c - \Delta\mu + \mathbf{u}\cdot\nabla c^{old} + \mathbf{u}^{old}\cdot\nabla c &=& \Delta\mu^s - \mathbf{u}^{s}\cdot\nabla c^s \\
+ -\mu - \epsilon^2\Delta c + W''(c^{old})c&=& \mu^s + \epsilon^2\Delta c^{s} - W'(c^s) \\
+ \partial^s_t u_j - z_j + u^{old}\cdot\nabla u_j + u\cdot\nabla u^{old}_j &=& z^s - u^s\cdot\nabla u^s_j,\quad j=0\ldots d-1 \\
+ \rho(c^{old})z + \rho'(c^{old})z^{old}\cdot c - \nabla\cdot\nu(c^{old})\nabla u_j - \nabla\cdot[c\nu'(c^{old})\nabla u^{old}_j] \\
+ + \nabla p - \alpha\mu\partial_j c^{old} - \alpha\mu^{old}\partial_j c &=& \nabla\cdot\nu(c^s)\nabla u^s_j - \nabla p^s + \alpha\mu^{s}\partial_j c^s + F_j,\quad j=0\ldots d-1 \\
+ -\nabla\cdot \mathbf{u} &=& \nabla\cdot \mathbf{u}^s
+\f}
*/