diff --git a/demo/init/vecellipt.dat.2d b/demo/init/vecellipt.dat.2d index 1930ed5100ef6aaa3a6f40cc010e56720be9c138..8ea2462bb3db0f8ef19f9ed74b0b8ceea3df868a 100644 --- a/demo/init/vecellipt.dat.2d +++ b/demo/init/vecellipt.dat.2d @@ -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 diff --git a/demo/init/vecheat.dat.2d b/demo/init/vecheat.dat.2d index fdddf318e3d54166fe9270d2d52c9272b2bb913e..1a9a95a663b15a1f78c2e00034c7c0adb1583d8b 100644 --- a/demo/init/vecheat.dat.2d +++ b/demo/init/vecheat.dat.2d @@ -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 diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc index 6989f42108cecb2523e1c4468e4fd12fa86eba04..1ba61c89c8e595c122fa2d66d881e1f5d0e624c3 100644 --- a/demo/src/ellipt.cc +++ b/demo/src/ellipt.cc @@ -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); } diff --git a/demo/src/heat.cc b/demo/src/heat.cc index c2477f5ffe58715d87b5928594a4270c7ff28ff9..595f0c09f8ed097f59c5ab9bfe14967408dfa021 100644 --- a/demo/src/heat.cc +++ b/demo/src/heat.cc @@ -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; } diff --git a/demo/src/vecellipt.cc b/demo/src/vecellipt.cc index f8e2103b639af89aa0e0a0420f5913dd6bbfaac9..8b2e5e2d29c5a8e19b5008b95bdb8db84342d22c 100644 --- a/demo/src/vecellipt.cc +++ b/demo/src/vecellipt.cc @@ -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; } diff --git a/demo/src/vecheat.cc b/demo/src/vecheat.cc index 8fad0bcd98c57b567a9853ed187648f3dfb9c125..b19d2a05a3470d1be6c8f0a3b8e14db721644e01 100644 --- a/demo/src/vecheat.cc +++ b/demo/src/vecheat.cc @@ -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; }