Skip to content

GitLab

Commit bebadadf authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

New code style for the first demos.

parent 645118a8
Hide whitespace changes
Inline Side-by-side
demo/init/vecellipt.dat.2d
View file @ bebadadf
......@@ -11,11 +11,11 @@ vecellipt->components: 2
vecellipt->polynomial degree[0]: 1
vecellipt->polynomial degree[1]: 1
vecellipt->solver: bicgstab
vecellipt->solver: umfpack
vecellipt->solver->max iteration: 1000
vecellipt->solver->tolerance: 1.e-8
vecellipt->solver->info: 2
vecellipt->solver->left precon: diag
vecellipt->solver->left precon: no
vecellipt->solver->right precon: no
vecellipt->estimator[0]: residual
......
demo/init/vecheat.dat.2d
View file @ bebadadf
......@@ -10,11 +10,11 @@ vecheat->space->mesh: vecheatMesh
vecheat->space->components: 2
vecheat->space->solver: cg
vecheat->space->solver: umfpack
vecheat->space->solver->max iteration: 1000
vecheat->space->solver->tolerance: 1.e-8
vecheat->space->solver->info: 2
vecheat->space->solver->left precon: diag
vecheat->space->solver->left precon: no
vecheat->space->solver->right precon: no
vecheat->space->estimator[0]: residual
......
demo/src/ellipt.cc
View file @ bebadadf
......@@ -47,37 +47,43 @@ int main(int argc, char* argv[])
// ===== check for init file =====
TEST_EXIT(argc == 2)("usage: ellipt initfile\n");
// ===== init parameters =====
Parameters::init(true, argv[1]);
// ===== create and init the scalar problem =====
ProblemScal ellipt("ellipt");
ellipt.initialize(INIT_ALL);
// === create adapt info ===
AdaptInfo *adaptInfo = new AdaptInfo("ellipt->adapt");
AdaptInfo adaptInfo("ellipt->adapt");
// === create adapt ===
AdaptStationary *adapt = new AdaptStationary("ellipt->adapt",
&ellipt,
adaptInfo);
AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);
// ===== add boundary conditions =====
ellipt.addDirichletBC(1, new G);
// ===== create matrix operator =====
Operator matrixOperator(Operator::MATRIX_OPERATOR, ellipt.getFESpace());
matrixOperator.addSecondOrderTerm(new Laplace_SOT);
ellipt.addMatrixOperator(&matrixOperator);
ellipt.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
int degree = ellipt.getFESpace()->getBasisFcts()->getDegree();
Operator rhsOperator(Operator::VECTOR_OPERATOR, ellipt.getFESpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(&rhsOperator);
ellipt.addVectorOperator(rhsOperator);
// ===== start adaption loop =====
adapt->adapt();
adapt.adapt();
ellipt.writeFiles(adaptInfo, true);
}
......
demo/src/heat.cc
View file @ bebadadf
......@@ -46,7 +46,7 @@ public:
class Heat : public ProblemInstatScal
{
public:
Heat(ProblemScal *heatSpace)
Heat(ProblemScal &heatSpace)
: ProblemInstatScal("heat", heatSpace)
{
// init theta scheme
......@@ -169,81 +169,74 @@ int main(int argc, char** argv)
Parameters::init(false, argv[1]);
Parameters::readArgv(argc, argv);
// ===== create and init stationary problem =====
ProblemScal *heatSpace = new ProblemScal("heat->space");
heatSpace->initialize(INIT_ALL);
ProblemScal heatSpace("heat->space");
heatSpace.initialize(INIT_ALL);
// ===== create instationary problem =====
Heat *heat = new Heat(heatSpace);
heat->initialize(INIT_ALL);
Heat heat(heatSpace);
heat.initialize(INIT_ALL);
// create adapt info
AdaptInfo *adaptInfo = new AdaptInfo("heat->adapt");
AdaptInfo adaptInfo("heat->adapt");
// create initial adapt info
AdaptInfo *adaptInfoInitial = new AdaptInfo("heat->initial->adapt");
AdaptInfo adaptInfoInitial("heat->initial->adapt");
// create instationary adapt
AdaptInstationary *adaptInstat =
new AdaptInstationary("heat->adapt",
heatSpace,
adaptInfo,
heat,
adaptInfoInitial);
AdaptInstationary adaptInstat("heat->adapt",
heatSpace,
adaptInfo,
heat,
adaptInfoInitial);
// ===== create boundary functions =====
G *boundaryFct = new G;
boundaryFct->setTimePtr(heat->getBoundaryTimePtr());
heat->setExactSolution(boundaryFct);
boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
heat.setExactSolution(boundaryFct);
heatSpace.addDirichletBC(DIRICHLET, boundaryFct);
heatSpace->addDirichletBC(DIRICHLET, boundaryFct);
// ===== create rhs functions =====
int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree();
int degree = heatSpace.getFESpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
rhsFct->setTimePtr(heat->getRHSTimePtr());
rhsFct->setTimePtr(heat.getRHSTimePtr());
// ===== create operators =====
double one = 1.0;
double zero = 0.0;
// create laplace
Operator *A = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
heatSpace->getFESpace());
A->addSecondOrderTerm(new Laplace_SOT);
A->setUhOld(heat->getOldSolution());
if ((*(heat->getThetaPtr())) != 0.0)
heatSpace->addMatrixOperator(A, heat->getThetaPtr(), &one);
if ((*(heat->getTheta1Ptr())) != 0.0)
heatSpace->addVectorOperator(A, heat->getTheta1Ptr(), &zero);
Operator A(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
heatSpace.getFESpace());
A.addSecondOrderTerm(new Laplace_SOT);
A.setUhOld(heat.getOldSolution());
if (*(heat.getThetaPtr()) != 0.0)
heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one);
if (*(heat.getTheta1Ptr()) != 0.0)
heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);
// create zero order operator
Operator *C = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
heatSpace->getFESpace());
C->addZeroOrderTerm(new Simple_ZOT);
C->setUhOld(heat->getOldSolution());
heatSpace->addMatrixOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
Operator C(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
heatSpace.getFESpace());
C.addZeroOrderTerm(new Simple_ZOT);
C.setUhOld(heat.getOldSolution());
heatSpace.addMatrixOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
heatSpace.addVectorOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
// create RHS operator
Operator *F = new Operator(Operator::VECTOR_OPERATOR,
heatSpace->getFESpace());
F->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
Operator F(Operator::VECTOR_OPERATOR, heatSpace.getFESpace());
F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
heatSpace.addVectorOperator(F);
heatSpace->addVectorOperator(F);
// ===== start adaption loop =====
int errorCode = adaptInstat->adapt();
int errorCode = adaptInstat.adapt();
return errorCode;
}
demo/src/vecellipt.cc
View file @ bebadadf
......@@ -45,64 +45,56 @@ int main(int argc, char* argv[])
// ===== check for init file =====
TEST_EXIT(argc == 2)("usage: vecellipt initfile\n");
// ===== init parameters =====
Parameters::init(false, argv[1]);
// ===== create and init the scalar problem =====
ProblemVec vecellipt("vecellipt");
vecellipt.initialize(INIT_ALL);
// === create adapt info ===
AdaptInfo *adaptInfo = new AdaptInfo("vecellipt->adapt",
vecellipt.getNumComponents());
AdaptInfo adaptInfo("vecellipt->adapt", vecellipt.getNumComponents());
// === create adapt ===
AdaptStationary *adapt = new AdaptStationary("vecellipt->adapt",
&vecellipt,
adaptInfo);
AdaptStationary adapt("vecellipt->adapt", vecellipt, adaptInfo);
// ===== add boundary conditions =====
vecellipt.addDirichletBC(1, 0, 0, new G);
// ===== create operators =====
Operator matrixOperator00(Operator::MATRIX_OPERATOR,
vecellipt.getFESpace(0),
vecellipt.getFESpace(0));
matrixOperator00.addSecondOrderTerm(new Laplace_SOT);
vecellipt.addMatrixOperator(&matrixOperator00, 0, 0);
Operator rhsOperator0(Operator::VECTOR_OPERATOR,
vecellipt.getFESpace(0));
int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();
rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
vecellipt.addVectorOperator(&rhsOperator0, 0);
// ===== create matrix operators =====
Operator matrixOperator00(Operator::MATRIX_OPERATOR,
vecellipt.getFESpace(0), vecellipt.getFESpace(0));
matrixOperator00.addSecondOrderTerm(new Laplace_SOT);
vecellipt.addMatrixOperator(matrixOperator00, 0, 0);
Operator matrixOperator10(Operator::MATRIX_OPERATOR,
vecellipt.getFESpace(1),
vecellipt.getFESpace(0));
Operator matrixOperator10(Operator::MATRIX_OPERATOR,
vecellipt.getFESpace(1), vecellipt.getFESpace(0));
matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
vecellipt.addMatrixOperator(matrixOperator10, 1, 0);
Operator matrixOperator11(Operator::MATRIX_OPERATOR,
vecellipt.getFESpace(1),
vecellipt.getFESpace(1));
matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
vecellipt.getFESpace(1), vecellipt.getFESpace(1));
matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0));
vecellipt.addMatrixOperator(matrixOperator11, 1, 1);
vecellipt.addMatrixOperator(&matrixOperator10, 1, 0);
matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0));
// ===== create rhs operator =====
int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();
Operator rhsOperator0(Operator::VECTOR_OPERATOR, vecellipt.getFESpace(0));
rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
vecellipt.addVectorOperator(rhsOperator0, 0);
vecellipt.addMatrixOperator(&matrixOperator11, 1, 1);
// ===== start adaption loop =====
adapt->adapt();
adapt.adapt();
vecellipt.writeFiles(adaptInfo, true);
return 0;
}
demo/src/vecheat.cc
View file @ bebadadf
......@@ -48,7 +48,7 @@ class Vecheat : public ProblemInstatVec
public:
/// Constructor
Vecheat(ProblemVec *vecheatSpace)
Vecheat(ProblemVec &vecheatSpace)
: ProblemInstatVec("vecheat", vecheatSpace)
{
// init theta scheme
......@@ -217,145 +217,106 @@ int main(int argc, char** argv)
Parameters::readArgv(argc, argv);
// ===== create and init stationary problem =====
ProblemVec *vecheatSpace = new ProblemVec("vecheat->space");
vecheatSpace->initialize(INIT_ALL);
ProblemVec vecheatSpace("vecheat->space");
vecheatSpace.initialize(INIT_ALL);
// ===== create instationary problem =====
Vecheat *vecheat = new Vecheat(vecheatSpace);
vecheat->initialize(INIT_ALL);
Vecheat vecheat(vecheatSpace);
vecheat.initialize(INIT_ALL);
// create adapt info
AdaptInfo *adaptInfo =
new AdaptInfo("vecheat->adapt",
vecheatSpace->getNumComponents());
AdaptInfo adaptInfo("vecheat->adapt", vecheatSpace.getNumComponents());
// create initial adapt info
AdaptInfo *adaptInfoInitial =
new AdaptInfo("vecheat->initial->adapt",
vecheatSpace->getNumComponents());
AdaptInfo adaptInfoInitial("vecheat->initial->adapt", vecheatSpace.getNumComponents());
// create instationary adapt
AdaptInstationary *adaptInstat =
new AdaptInstationary("vecheat->adapt",
vecheatSpace,
adaptInfo,
vecheat,
adaptInfoInitial);
double fac = (*(vecheat->getThetaPtr())) != 0 ? 1.0 : 1.0e-3;
int nRefine = 0, dim = vecheatSpace->getMesh(0)->getDim();;
GET_PARAMETER(0, vecheatSpace->getMesh(0)->getName()
AdaptInstationary adaptInstat("vecheat->adapt",
vecheatSpace,
adaptInfo,
vecheat,
adaptInfoInitial);
double fac = *(vecheat.getThetaPtr()) != 0.0 ? 1.0 : 1.0e-3;
int nRefine = 0;
int dim = vecheatSpace.getMesh(0)->getDim();
GET_PARAMETER(0, vecheatSpace.getMesh(0)->getName()
+ "->global refinements", "%d", &nRefine);
if((*(vecheat->getThetaPtr())) == 0.5) {
(*(adaptInfo->getTimestepPtr())) *=
fac * pow(2.0, ((double) -nRefine)/dim);
} else {
(*(adaptInfo->getTimestepPtr())) *=
fac * pow(2.0, -nRefine);
}
if (*(vecheat.getThetaPtr()) == 0.5)
*(adaptInfo.getTimestepPtr()) *= fac * pow(2.0, static_cast<double>(-nRefine) / dim);
else
*(adaptInfo.getTimestepPtr()) *= fac * pow(2.0, -nRefine);
// ===== create boundary functions =====
G *boundaryFct = new G;
boundaryFct->setTimePtr(vecheat->getBoundaryTimePtr());
vecheat->setBoundaryFct(boundaryFct);
boundaryFct->setTimePtr(vecheat.getBoundaryTimePtr());
vecheat.setBoundaryFct(boundaryFct);
vecheatSpace->addDirichletBC(DIRICHLET, 0, 0, boundaryFct);
vecheatSpace->addDirichletBC(DIRICHLET, 1, 1, boundaryFct);
vecheatSpace.addDirichletBC(DIRICHLET, 0, 0, boundaryFct);
vecheatSpace.addDirichletBC(DIRICHLET, 1, 1, boundaryFct);
// ===== create rhs functions =====
F *rhsFct0 = new F(vecheatSpace->getFESpace(0)->getBasisFcts()->getDegree());
rhsFct0->setTimePtr(vecheat->getRHSTimePtr());
F *rhsFct1 = new F(vecheatSpace->getFESpace(1)->getBasisFcts()->getDegree());
rhsFct1->setTimePtr(vecheat->getRHSTimePtr());
F *rhsFct0 = new F(vecheatSpace.getFESpace(0)->getBasisFcts()->getDegree());
rhsFct0->setTimePtr(vecheat.getRHSTimePtr());
F *rhsFct1 = new F(vecheatSpace.getFESpace(1)->getBasisFcts()->getDegree());
rhsFct1->setTimePtr(vecheat.getRHSTimePtr());
// ===== create operators =====
double one = 1.0;
double zero = 0.0;
// create laplace
Operator *A00 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(0),
vecheatSpace->getFESpace(0));
A00->addSecondOrderTerm(new Laplace_SOT);
A00->setUhOld(vecheat->getOldSolution()->getDOFVector(0));
Operator *A11 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(1),
vecheatSpace->getFESpace(1));
A11->addSecondOrderTerm(new Laplace_SOT);
A11->setUhOld(vecheat->getOldSolution()->getDOFVector(1));
if ((*(vecheat->getThetaPtr())) != 0.0) {
vecheatSpace->addMatrixOperator(A00, 0, 0,
vecheat->getThetaPtr(),
&one);
vecheatSpace->addMatrixOperator(A11, 1, 1,
vecheat->getThetaPtr(),
&one);
Operator A00 (Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
vecheatSpace.getFESpace(0), vecheatSpace.getFESpace(0));
A00.addSecondOrderTerm(new Laplace_SOT);
A00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));
Operator A11(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
vecheatSpace.getFESpace(1), vecheatSpace.getFESpace(1));
A11.addSecondOrderTerm(new Laplace_SOT);
A11.setUhOld(vecheat.getOldSolution()->getDOFVector(1));
if (*(vecheat.getThetaPtr()) != 0.0) {
vecheatSpace.addMatrixOperator(A00, 0, 0, vecheat.getThetaPtr(), &one);
vecheatSpace.addMatrixOperator(A11, 1, 1, vecheat.getThetaPtr(), &one);
}
if ((*(vecheat->getTheta1Ptr())) != 0.0) {
vecheatSpace->addVectorOperator(A00, 0, vecheat->getTheta1Ptr(), &zero);
vecheatSpace->addVectorOperator(A11, 1, vecheat->getTheta1Ptr(), &zero);
if (*(vecheat.getTheta1Ptr()) != 0.0) {
vecheatSpace.addVectorOperator(A00, 0, vecheat.getTheta1Ptr(), &zero);
vecheatSpace.addVectorOperator(A11, 1, vecheat.getTheta1Ptr(), &zero);
}
// create zero order operator
Operator *C00 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(0),
vecheatSpace->getFESpace(0));
C00->addZeroOrderTerm(new Simple_ZOT);
C00->setUhOld(vecheat->getOldSolution()->getDOFVector(0));
vecheatSpace->addMatrixOperator(C00, 0, 0,
vecheat->getTau1Ptr(),
vecheat->getTau1Ptr());
vecheatSpace->addVectorOperator(C00, 0,
vecheat->getTau1Ptr(),
vecheat->getTau1Ptr());
Operator *C11 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(1),
vecheatSpace->getFESpace(1));
C11->addZeroOrderTerm(new Simple_ZOT);
C11->setUhOld(vecheat->getOldSolution()->getDOFVector(1));
vecheatSpace->addMatrixOperator(C11, 1, 1,
vecheat->getTau1Ptr(),
vecheat->getTau1Ptr());
vecheatSpace->addVectorOperator(C11, 1,
vecheat->getTau1Ptr(),
vecheat->getTau1Ptr());
Operator C00(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
vecheatSpace.getFESpace(0), vecheatSpace.getFESpace(0));
C00.addZeroOrderTerm(new Simple_ZOT);
C00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));
vecheatSpace.addMatrixOperator(C00, 0, 0,
vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
vecheatSpace.addVectorOperator(C00, 0, vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
Operator C11(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
vecheatSpace.getFESpace(1), vecheatSpace.getFESpace(1));
C11.addZeroOrderTerm(new Simple_ZOT);
C11.setUhOld(vecheat.getOldSolution()->getDOFVector(1));
vecheatSpace.addMatrixOperator(C11, 1, 1,
vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
vecheatSpace.addVectorOperator(C11, 1,
vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
// create RHS operator
Operator *F0 = new Operator(Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(0));
F0->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct0));
vecheatSpace->addVectorOperator(F0, 0);
Operator *F1 = new Operator(Operator::VECTOR_OPERATOR,
vecheatSpace->getFESpace(1));
F1->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct1));
Operator F0(Operator::VECTOR_OPERATOR, vecheatSpace.getFESpace(0));
F0.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct0));
vecheatSpace.addVectorOperator(F0, 0);
vecheatSpace->addVectorOperator(F1, 1);
Operator F1 = Operator(Operator::VECTOR_OPERATOR, vecheatSpace.getFESpace(1));
F1.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct1));
vecheatSpace.addVectorOperator(F1, 1);
// ===== start adaption loop =====
int errorCode = adaptInstat->adapt();
int errorCode = adaptInstat.adapt();
return errorCode;
}
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