Commit d41f7355 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

CouplingBasePorblems adopted to new namespace structure

parent 508c68fa
......@@ -23,10 +23,10 @@ using namespace AMDiS;
template<class ProblemType=ProblemStat, class BaseProblemType=BaseProblem<ProblemStat> >
class CouplingBaseProblem : public CouplingIterationInterface,
public CouplingTimeInterface,
public CouplingProblemStatImpl<ProblemType>
public detail::CouplingProblemStat<ProblemType>
{
public:
typedef CouplingProblemStatImpl<ProblemType> CProblemStat;
typedef detail::CouplingProblemStat<ProblemType> CProblemStat;
CouplingBaseProblem(std::string name_,
BaseProblemType *prob0_,
......
......@@ -77,10 +77,10 @@ template<class ProblemType=ProblemStat,
BOOST_PP_ENUM_PARAMS_WITH_A_DEFAULT(MAX_NUM_COUPLED_PROBLEMS, class BaseProblemType, void) >
class CouplingBaseProblem : public CouplingIterationInterface,
public CouplingTimeInterface,
public CouplingProblemStatImpl<ProblemType>
public detail::CouplingProblemStat<ProblemType>
{
public:
typedef CouplingProblemStatImpl<ProblemType> CProblemStat;
typedef detail::CouplingProblemStat<ProblemType> CProblemStat;
typedef typename boost::tuples::detail::tie_mapper
< BOOST_PP_ENUM_PARAMS(MAX_NUM_COUPLED_PROBLEMS, BaseProblemType) >::type BaseProblemsTupleType;
......
......@@ -78,10 +78,10 @@ namespace detail {
template<typename ProblemType=ProblemStat, typename... BaseProblemTypes >
class CouplingBaseProblem : public CouplingIterationInterface,
public CouplingTimeInterface,
public CouplingProblemStatImpl<ProblemType>
public detail::CouplingProblemStat<ProblemType>
{
public:
typedef CouplingProblemStatImpl<ProblemType> CProblemStat;
typedef detail::CouplingProblemStat<ProblemType> CProblemStat;
typedef std::tuple<BaseProblemTypes&...> BaseProblemsTupleType;
template<typename... BaseProblemTypes_>
......
......@@ -130,7 +130,7 @@ void NavierStokesCahnHilliard::initData()
if (velocity == NULL)
velocity = new DOFVector<WorldVector<double> >(getFeSpace(0), "velocity");
//fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
super::initData();
}
......@@ -139,8 +139,24 @@ void NavierStokesCahnHilliard::initData()
void NavierStokesCahnHilliard::closeTimestep(AdaptInfo *adaptInfo)
{
super::closeTimestep(adaptInfo);
calcVelocity();
fileWriter->writeFiles(adaptInfo, true);
}
void NavierStokesCahnHilliard::transferInitialSolution(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokesCahnHilliard::transferInitialSolution()");
calcVelocity();
fileWriter->writeFiles(adaptInfo, false);
writeFiles(adaptInfo, false);
// initial parameters for detecting mesh changes
oldMeshChangeIdx= getMesh()->getChangeIndex();
}
void NavierStokesCahnHilliard::initTimestep(AdaptInfo *adaptInfo)
{
super::initTimestep(adaptInfo);
......@@ -186,6 +202,8 @@ void NavierStokesCahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo)
void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
{
Flag initFlag = initDataFromFile(adaptInfo);
int comp = 2+dow;
if (!initFlag.isSet(DATA_ADOPTED)) {
int initialInterface = 0;
......@@ -198,28 +216,28 @@ void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
double a= 0.0, dir= -1.0;
Initfile::get(name + "->line->pos", a);
Initfile::get(name + "->line->direction", dir);
prob->getSolution()->getDOFVector(1+3)->interpol(new Plane(a, dir));
prob->getSolution()->getDOFVector(comp)->interpol(new Plane(a, dir));
}
else if (initialInterface == 1) {
/// schraege Linie
double theta = m_pi/4.0;
prob->getSolution()->getDOFVector(1+3)->interpol(new PlaneRotation(0.0, theta, 1.0));
transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
prob->getSolution()->getDOFVector(comp)->interpol(new PlaneRotation(0.0, theta, 1.0));
transformDOFInterpolation(prob->getSolution()->getDOFVector(comp),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
}
else if (initialInterface == 2) {
/// Ellipse
double a= 1.0, b= 1.0;
Initfile::get(name + "->ellipse->a", a);
Initfile::get(name + "->ellipse->b", b);
prob->getSolution()->getDOFVector(1+3)->interpol(new Ellipse(a,b));
prob->getSolution()->getDOFVector(comp)->interpol(new Ellipse(a,b));
}
else if (initialInterface == 3) {
/// zwei horizontale Linien
double a= 0.75, b= 0.375;
Initfile::get(name + "->lines->pos1", a);
Initfile::get(name + "->lines->pos2", b);
prob->getSolution()->getDOFVector(1+3)->interpol(new Plane(a, -1.0));
transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),new Plane(b, 1.0), new AMDiS::Max<double>);
prob->getSolution()->getDOFVector(comp)->interpol(new Plane(a, -1.0));
transformDOFInterpolation(prob->getSolution()->getDOFVector(comp),new Plane(b, 1.0), new AMDiS::Max<double>);
}
else if (initialInterface == 4) {
/// Kreis
......@@ -229,7 +247,7 @@ void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
Initfile::get(name + "->circle->center_y", y);
WorldVector<double> w;
w[0]=x; w[1]=y+adaptInfo->getTime();
prob->getSolution()->getDOFVector(1+3)->interpol(new Circle(radius, w));
prob->getSolution()->getDOFVector(comp)->interpol(new Circle(radius, w));
} else if (initialInterface == 5) {
/// Rechteck
double width = 0.5;
......@@ -238,15 +256,15 @@ void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
Initfile::get(name + "->rectangle->width", width);
Initfile::get(name + "->rectangle->height", height);
Initfile::get(name + "->rectangle->center", center);
prob->getSolution()->getDOFVector(1+3)->interpol(new Rectangle(width, height, center));
prob->getSolution()->getDOFVector(comp)->interpol(new Rectangle(width, height, center));
}
/// create phase-field from signed-dist-function
if (doubleWell == 0) {
transformDOF(prob->getSolution()->getDOFVector(1+3),
transformDOF(prob->getSolution()->getDOFVector(comp),
new SignedDistToPhaseField(initialEps));
} else {
transformDOF(prob->getSolution()->getDOFVector(1+3),
transformDOF(prob->getSolution()->getDOFVector(comp),
new SignedDistToCh(initialEps));
}
}
......@@ -256,13 +274,14 @@ void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
void NavierStokesCahnHilliard::fillOperators()
{
const FiniteElemSpace* feSpace = prob->getFeSpace(0);
int degree = prob->getFeSpace(dow+1)->getBasisFcts()->getDegree();
int comp = dow+2, comp2 = dow+1;
int degree = prob->getFeSpace(comp)->getBasisFcts()->getDegree();
// c
Operator *opChMnew = new Operator(feSpace,feSpace);
opChMnew->addZeroOrderTerm(new Simple_ZOT);
Operator *opChMold = new Operator(feSpace,feSpace);
opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1+3)));
opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(comp)));
// -nabla*(grad(c))
Operator *opChL = new Operator(feSpace,feSpace);
opChL->addSecondOrderTerm(new Simple_SOT);
......@@ -274,22 +293,22 @@ void NavierStokesCahnHilliard::fillOperators()
// -2*c_old^3 + 3/2*c_old^2
Operator *opChMPowExpl = new Operator(feSpace,feSpace);
opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getSolution()->getDOFVector(1+3),
prob->getSolution()->getDOFVector(comp),
new AMDiS::Pow<3>(-2.0, 3*degree)));
if (doubleWell == 0) {
opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getSolution()->getDOFVector(1+3),
prob->getSolution()->getDOFVector(comp),
new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
}
// -3*c_old^2 * c
Operator *opChMPowImpl = new Operator(feSpace,feSpace);
opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getSolution()->getDOFVector(1+3),
prob->getSolution()->getDOFVector(comp),
new AMDiS::Pow<2>(-3.0, 2*degree)));
if (doubleWell == 0) {
opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getSolution()->getDOFVector(1+3),
prob->getSolution()->getDOFVector(comp),
new AMDiS::Factor<double>(3.0, degree)));
opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(-0.5));
} else {
......@@ -298,21 +317,21 @@ void NavierStokesCahnHilliard::fillOperators()
// mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
// ----------------------------------------------------------------------
prob->addMatrixOperator(*opChMnew,0+3,0+3, &ab); /// < phi*mu , psi >
prob->addMatrixOperator(*opChMnew,comp2,comp2, &ab); /// < phi*mu , psi >
prob->addMatrixOperator(*opChMPowImpl,0+3,1+3, &b); /// < -3*phi*c*c_old^2 , psi >
prob->addMatrixOperator(*opChL,0+3,1+3, &minusEpsSqr); /// < -eps^2*phi*grad(c) , grad(psi) >
prob->addMatrixOperator(*opChMPowImpl,comp2,comp, &b); /// < -3*phi*c*c_old^2 , psi >
prob->addMatrixOperator(*opChL,comp2,comp, &minusEpsSqr); /// < -eps^2*phi*grad(c) , grad(psi) >
// . . . vectorOperators . . . . . . . . . . . . . . .
prob->addVectorOperator(*opChMPowExpl,0+3, &b); /// < -2*phi*c_old^3 , psi >
prob->addVectorOperator(*opChMPowExpl,comp2, &b); /// < -2*phi*c_old^3 , psi >
// setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
// dt(c) = laplace(mu) - u*grad(c)
// -----------------------------------
prob->addMatrixOperator(*opChMnew,1+3,1+3, &b); /// < phi*c/tau , psi >
prob->addMatrixOperator(*opChLM,1+3,0+3, getTau()); /// < phi*grad(mu) , grad(psi) >
prob->addMatrixOperator(*opChMnew,comp,comp, &b); /// < phi*c/tau , psi >
prob->addMatrixOperator(*opChLM,comp,comp2, getTau()); /// < phi*grad(mu) , grad(psi) >
// . . . vectorOperators . . . . . . . . . . . . . . .
prob->addVectorOperator(*opChMold,1+3, &b ); /// < phi*c^old/tau , psi >
prob->addVectorOperator(*opChMold,comp, &b ); /// < phi*c^old/tau , psi >
/**/
......@@ -322,9 +341,8 @@ void NavierStokesCahnHilliard::fillOperators()
// int degree_p = getFeSpace(dow)->getBasisFcts()->getDegree();
WorldVector<DOFVector<double>* > vel;
for (unsigned i = 0; i < dow; i++){
vel[i]= prob->getSolution()->getDOFVector(i+2);
}
for (unsigned i = 0; i < dow; i++)
vel[i]= prob->getSolution()->getDOFVector(i);
// fill operators for prob
......@@ -380,9 +398,9 @@ void NavierStokesCahnHilliard::fillOperators()
prob->addMatrixOperator(*opLaplaceUi, i, i);
/// < p , d_i(psi) >
Operator *opGradP = new Operator(getFeSpace(i),getFeSpace(2));
Operator *opGradP = new Operator(getFeSpace(i),getFeSpace(dow));
opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
prob->addMatrixOperator(*opGradP, i, 2);
prob->addMatrixOperator(*opGradP, i, dow);
/// external force, i.e. gravitational force
if (force[i]*force[i] > DBL_TOL) {
......@@ -392,20 +410,20 @@ void NavierStokesCahnHilliard::fillOperators()
}
/// < d_i(u'_i) , psi >
Operator *opDivU = new Operator(getFeSpace(2),getFeSpace(i));
Operator *opDivU = new Operator(getFeSpace(dow),getFeSpace(i));
opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
prob->addMatrixOperator(*opDivU, 2, i);
prob->addMatrixOperator(*opDivU, dow, i);
///coupling Operators
Operator *opNuGradC = new Operator(getFeSpace(i), getFeSpace(dow+1));
opNuGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i));
prob->addMatrixOperator(opNuGradC, i, dow+1, &surfaceTension, &surfaceTension);
Operator *opNuGradC = new Operator(getFeSpace(i), getFeSpace(comp2));
opNuGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(comp), i));
prob->addMatrixOperator(opNuGradC, i, comp2, &surfaceTension, &surfaceTension);
Operator *opVGradC = new Operator(getFeSpace(dow+2), getFeSpace(i));
opVGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i, b));
prob->addMatrixOperator(opVGradC, dow+2, i, getTau());
Operator *opVGradC = new Operator(getFeSpace(comp), getFeSpace(i));
opVGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(comp), i, b));
prob->addMatrixOperator(opVGradC, comp, i, getTau());
/**/
}
......
......@@ -59,6 +59,7 @@ public: // public methods
virtual void solveInitialProblem(AdaptInfo *adaptInfo);
virtual void solveInitialProblem2(AdaptInfo *adaptInfo);
virtual void transferInitialSolution(AdaptInfo* adaptInfo);
double getEpsilon() { return eps; }
int getDoubleWellType() { return doubleWell; }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment