From 83b6de2bc9d8eb197e7bafc54237795acc625689 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski <thomas.witkowski@gmx.de> Date: Mon, 25 May 2009 11:39:21 +0000 Subject: [PATCH] Adapted demos to new AMDiS (DELETE/delete and dirichlet boundary conditions). --- demo/src/ball.cc | 41 ++- demo/src/bunny.cc | 39 ++- demo/src/couple.cc | 35 ++- demo/src/ellipt.cc | 44 ++-- demo/src/heat.cc | 157 +++++------- demo/src/mpcci1.cc | 50 ++-- demo/src/mpcci2.cc | 2 +- demo/src/mpccitest.cc | 54 ++-- demo/src/multigrid.cc | 49 ++-- demo/src/navierstokes.cc | 355 -------------------------- demo/src/neumann.cc | 46 ++-- demo/src/nonlin.cc | 122 ++++----- demo/src/nonlin2.cc | 100 +++----- demo/src/nonlin3.cc | 102 ++++---- demo/src/parallelellipt.cc | 47 ++-- demo/src/parallelheat.cc | 132 ++++------ demo/src/parametric.cc | 76 +++--- demo/src/periodic.cc | 42 ++-- demo/src/smitest.cc | 5 +- demo/src/sphere.cc | 25 +- demo/src/stokesnonlin.cc | 493 ------------------------------------- demo/src/torus.cc | 61 ++--- demo/src/vecellipt.cc | 42 ++-- demo/src/vecheat.cc | 245 +++++++----------- demo/src/vecmultigrid.cc | 45 ++-- demo/src/vecnonlin.cc | 128 ++++------ 26 files changed, 662 insertions(+), 1875 deletions(-) delete mode 100644 demo/src/navierstokes.cc delete mode 100644 demo/src/stokesnonlin.cc diff --git a/demo/src/ball.cc b/demo/src/ball.cc index 610e923f..fc024502 100644 --- a/demo/src/ball.cc +++ b/demo/src/ball.cc @@ -7,36 +7,27 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = x * x; double ux = exp(-10.0*r2); @@ -61,7 +52,7 @@ int main(int argc, char* argv[]) // ===== create projection ===== WorldVector<double> ballCenter; ballCenter.set(0.0); - NEW BallProject(1, + new BallProject(1, BOUNDARY_PROJECTION, ballCenter, 1.0); @@ -71,25 +62,25 @@ int main(int argc, char* argv[]) ball.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("ball->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("ball->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("ball->adapt", + AdaptStationary *adapt = new AdaptStationary("ball->adapt", &ball, adaptInfo); // ===== add boundary conditions ===== - ball.addDirichletBC(1, NEW G); + ball.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, ball.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); ball.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = ball.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, ball.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); ball.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/bunny.cc b/demo/src/bunny.cc index 3e6506b7..f52ea411 100644 --- a/demo/src/bunny.cc +++ b/demo/src/bunny.cc @@ -7,36 +7,27 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return -2 * x[0]; } }; -/** \brief - * boundary - */ +/// boundary class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return 10000.0; } }; @@ -58,7 +49,7 @@ int main(int argc, char* argv[]) // ===== create projection ===== WorldVector<double> ballCenter; ballCenter.set(0.0); - NEW BallProject(1, + new BallProject(1, VOLUME_PROJECTION, ballCenter, 1.0); @@ -68,20 +59,20 @@ int main(int argc, char* argv[]) bunny.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("bunny->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("bunny->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("bunny->adapt", + AdaptStationary *adapt = new AdaptStationary("bunny->adapt", &bunny, adaptInfo); // ===== add boundary conditions ===== - //bunny.addDirichletBC(1111, NEW G); + //bunny.addDirichletBC(1111, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, bunny.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); bunny.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== @@ -90,7 +81,7 @@ int main(int argc, char* argv[]) int degree = bunny.getFESpace()->getBasisFcts()->getDegree(); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); bunny.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/couple.cc b/demo/src/couple.cc index 91b15d75..9556ad72 100644 --- a/demo/src/couple.cc +++ b/demo/src/couple.cc @@ -6,9 +6,8 @@ using namespace AMDiS; class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - double operator()(const WorldVector<double>& x) const { + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; @@ -16,11 +15,10 @@ public: class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - double operator()(const WorldVector<double>& x) const { + double operator()(const WorldVector<double>& x) const + { int dow = x.getSize(); double r2 = (x * x); double ux = exp(-10.0 * r2); @@ -111,9 +109,10 @@ class Identity : public AbstractFunction<double, double> public: MEMORY_MANAGED(Identity); - Identity(int degree) : AbstractFunction<double, double>(degree) {}; + Identity(int degree) : AbstractFunction<double, double>(degree) {} - double operator()(const double& x) const { + double operator()(const double& x) const + { return x; } }; @@ -145,35 +144,35 @@ int main(int argc, char* argv[]) adoptFlag); // ===== add boundary conditions for problem1 ===== - problem1.addDirichletBC(1, NEW G); + problem1.addDirichletBC(1, new G); // ===== add boundary conditions for problem1 ===== - //problem2.addDirichletBC(1, NEW G); + //problem2.addDirichletBC(1, new G); // ===== create operators for problem1 ===== Operator matrixOperator1(Operator::MATRIX_OPERATOR, problem1.getFESpace()); - matrixOperator1.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator1.addSecondOrderTerm(new Laplace_SOT); problem1.addMatrixOperator(&matrixOperator1); int degree = problem1.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator1(Operator::VECTOR_OPERATOR, problem1.getFESpace()); - rhsOperator1.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator1.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); problem1.addVectorOperator(&rhsOperator1); // ===== create operators for problem2 ===== Operator matrixOperator2(Operator::MATRIX_OPERATOR, problem2.getFESpace()); - matrixOperator2.addZeroOrderTerm(NEW Simple_ZOT); + matrixOperator2.addZeroOrderTerm(new Simple_ZOT); problem2.addMatrixOperator(&matrixOperator2); Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFESpace()); - rhsOperator2.addZeroOrderTerm(NEW VecAtQP_ZOT(problem1.getSolution(), - NEW Identity(degree))); + rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(), + new Identity(degree))); problem2.addVectorOperator(&rhsOperator2); // ===== create adaptation loop and iteration interface ===== - AdaptInfo *adaptInfo = NEW AdaptInfo("couple->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("couple->adapt", 1); MyCoupledIteration coupledIteration(&problem1, &problem2); - AdaptStationary *adapt = NEW AdaptStationary("couple->adapt", + AdaptStationary *adapt = new AdaptStationary("couple->adapt", &coupledIteration, adaptInfo); diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc index 9d303c63..d924ab4c 100644 --- a/demo/src/ellipt.cc +++ b/demo/src/ellipt.cc @@ -3,35 +3,30 @@ using namespace AMDiS; using namespace std; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +// =========================================================================== +// ===== function definitions ================================================ +// =========================================================================== + +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dow = x.getSize(); double r2 = (x * x); double ux = exp(-10.0 * r2); @@ -39,7 +34,10 @@ public: } }; -// // ===== main program // +// =========================================================================== +// ===== main program ======================================================== +// =========================================================================== + int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -55,24 +53,24 @@ int main(int argc, char* argv[]) ellipt.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("ellipt->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("ellipt->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("ellipt->adapt", + AdaptStationary *adapt = new AdaptStationary("ellipt->adapt", &ellipt, adaptInfo); // ===== add boundary conditions ===== - ellipt.addDirichletBC(1, NEW G); + ellipt.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, ellipt.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); ellipt.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== Operator rhsOperator(Operator::VECTOR_OPERATOR, ellipt.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F())); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F())); ellipt.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/heat.cc b/demo/src/heat.cc index 72aacec8..c2477f5f 100644 --- a/demo/src/heat.cc +++ b/demo/src/heat.cc @@ -7,38 +7,29 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = x * x; double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2); @@ -51,17 +42,10 @@ public: // ===== instationary problem ================================================ // =========================================================================== -/** \brief - * Instationary problem - */ +/// Instationary problem class Heat : public ProblemInstatScal { public: - MEMORY_MANAGED(Heat); - - /** \brief - * Constructor - */ Heat(ProblemScal *heatSpace) : ProblemInstatScal("heat", heatSpace) { @@ -77,121 +61,98 @@ public: theta1 = theta - 1; } - // ===== ProblemInstatBase methods =================================== - /** \brief - * set the time in all needed functions! - */ - void setTime(AdaptInfo *adaptInfo) { + /// set the time in all needed functions! + void setTime(AdaptInfo *adaptInfo) + { rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep(); boundaryTime = adaptInfo->getTime(); tau1 = 1.0 / adaptInfo->getTimestep(); } - void closeTimestep(AdaptInfo *adaptInfo) { + void closeTimestep(AdaptInfo *adaptInfo) + { ProblemInstatScal::closeTimestep(adaptInfo); WAIT; } // ===== initial problem methods ===================================== - /** \brief - * Used by \ref problemInitial to solve the system of the initial problem - */ - void solve(AdaptInfo *adaptInfo) { + /// Used by \ref problemInitial to solve the system of the initial problem + void solve(AdaptInfo *adaptInfo) + { problemStat->getSolution()->interpol(exactSolution); } - /** \brief - * Used by \ref problemInitial to do error estimation for the initial - * problem. - */ - void estimate(AdaptInfo *adaptInfo) { + /// Used by \ref problemInitial to do error estimation for the initial problem. + void estimate(AdaptInfo *adaptInfo) + { double errMax, errSum; errSum = Error<double>::L2Err(*exactSolution, *(problemStat->getSolution()), 0, &errMax, false); - adaptInfo->setEstSum(errSum, 0); adaptInfo->setEstMax(errMax, 0); } // ===== setting methods =============================================== - /** \brief - * Sets \ref exactSolution; - */ - void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct) { + /// Sets \ref exactSolution; + void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct) + { exactSolution = fct; } // ===== getting methods =============================================== - /** \brief - * Returns pointer to \ref theta. - */ - double *getThetaPtr() { + /// Returns pointer to \ref theta. + double *getThetaPtr() + { return θ } - /** \brief - * Returns pointer to \ref theta1. - */ - double *getTheta1Ptr() { + /// Returns pointer to \ref theta1. + double *getTheta1Ptr() + { return &theta1; } - /** \brief - * Returns pointer to \ref tau1 - */ - double *getTau1Ptr() { + /// Returns pointer to \ref tau1 + double *getTau1Ptr() + { return &tau1; } - /** \brief - * Returns pointer to \ref rhsTime. - */ - double *getRHSTimePtr() { + /// Returns pointer to \ref rhsTime. + double *getRHSTimePtr() + { return &rhsTime; } - /** \brief - * Returns pointer to \ref theta1. - */ - double *getBoundaryTimePtr() { + /// Returns pointer to \ref theta1. + double *getBoundaryTimePtr() + { return &boundaryTime; } private: - /** \brief - * Used for theta scheme. - */ + /// Used for theta scheme. double theta; - /** \brief - * theta - 1 - */ + /// theta - 1 double theta1; - /** \brief - * 1.0 / timestep - */ + /// 1.0 / timestep double tau1; - /** \brief - * time for right hand side functions. - */ + /// time for right hand side functions. double rhsTime; - /** \brief - * time for boundary functions. - */ + /// time for boundary functions. double boundaryTime; - /** \brief - * Pointer to boundary function. Needed for initial problem. - */ + /// Pointer to boundary function. Needed for initial problem. AbstractFunction<double, WorldVector<double> > *exactSolution; }; @@ -209,7 +170,7 @@ int main(int argc, char** argv) Parameters::readArgv(argc, argv); // ===== create and init stationary problem ===== - ProblemScal *heatSpace = NEW ProblemScal("heat->space"); + ProblemScal *heatSpace = new ProblemScal("heat->space"); heatSpace->initialize(INIT_ALL); // ===== create instationary problem ===== @@ -217,21 +178,21 @@ int main(int argc, char** argv) heat->initialize(INIT_ALL); // create adapt info - AdaptInfo *adaptInfo = NEW AdaptInfo("heat->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("heat->adapt"); // create initial adapt info - AdaptInfo *adaptInfoInitial = NEW AdaptInfo("heat->initial->adapt"); + AdaptInfo *adaptInfoInitial = new AdaptInfo("heat->initial->adapt"); // create instationary adapt AdaptInstationary *adaptInstat = - NEW AdaptInstationary("heat->adapt", + new AdaptInstationary("heat->adapt", heatSpace, adaptInfo, heat, adaptInfoInitial); // ===== create boundary functions ===== - G *boundaryFct = NEW G; + G *boundaryFct = new G; boundaryFct->setTimePtr(heat->getBoundaryTimePtr()); heat->setExactSolution(boundaryFct); @@ -239,7 +200,7 @@ int main(int argc, char** argv) // ===== create rhs functions ===== int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree(); - F *rhsFct = NEW F(degree); + F *rhsFct = new F(degree); rhsFct->setTimePtr(heat->getRHSTimePtr()); // ===== create operators ===== @@ -247,7 +208,7 @@ int main(int argc, char** argv) double zero = 0.0; // create laplace - Operator *A = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *A = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); @@ -255,18 +216,18 @@ int main(int argc, char** argv) A->setUhOld(heat->getOldSolution()); - if((*(heat->getThetaPtr())) != 0.0) + if ((*(heat->getThetaPtr())) != 0.0) heatSpace->addMatrixOperator(A, heat->getThetaPtr(), &one); - if((*(heat->getTheta1Ptr())) != 0.0) + if ((*(heat->getTheta1Ptr())) != 0.0) heatSpace->addVectorOperator(A, heat->getTheta1Ptr(), &zero); // create zero order operator - Operator *C = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *C = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); - C->addZeroOrderTerm(NEW Simple_ZOT); + C->addZeroOrderTerm(new Simple_ZOT); C->setUhOld(heat->getOldSolution()); @@ -274,10 +235,10 @@ int main(int argc, char** argv) heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr()); // create RHS operator - Operator *F = NEW Operator(Operator::VECTOR_OPERATOR, + Operator *F = new Operator(Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); - F->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct)); + F->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct)); heatSpace->addVectorOperator(F); diff --git a/demo/src/mpcci1.cc b/demo/src/mpcci1.cc index 77ea2776..231e515d 100644 --- a/demo/src/mpcci1.cc +++ b/demo/src/mpcci1.cc @@ -6,37 +6,30 @@ using namespace std; using namespace AMDiS; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +// =========================================================================== +// ===== function definitions ================================================ +// =========================================================================== + +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ + /// Implementation of AbstractFunction::operator(). double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = x * x; double ux = exp(-10.0 * r2); @@ -57,8 +50,8 @@ public: static MultiGridSolver *mgSolver = NULL; if (!mgSolver) - mgSolver = NEW MultiGridSolver(feSpace_, - NEW GSSmoother(1.0), + mgSolver = new MultiGridSolver(feSpace_, + new GSSmoother(1.0), systemMatrix_, solution_, rhs_); @@ -67,7 +60,10 @@ public: } }; -// // ===== main program // +// =========================================================================== +// ===== main program ======================================================== +// =========================================================================== + int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -84,27 +80,27 @@ int main(int argc, char* argv[]) mpcci1.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("mpcci1->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("mpcci1->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("mpcci1->adapt", + AdaptStationary *adapt = new AdaptStationary("mpcci1->adapt", &mpcci1, adaptInfo); //mpcci1.setAdapt(adapt); // ===== add boundary conditions ===== - mpcci1.addDirichletBC(1, NEW G); + mpcci1.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, mpcci1.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); mpcci1.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = mpcci1.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, mpcci1.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); mpcci1.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/mpcci2.cc b/demo/src/mpcci2.cc index c3933e6b..f6197e7f 100644 --- a/demo/src/mpcci2.cc +++ b/demo/src/mpcci2.cc @@ -31,7 +31,7 @@ int main(int argc, char* argv[]) 1, &quantityID, vec, 1, &syncPointID); - while(convergence != CCI_STOP) { + while (convergence != CCI_STOP) { mpcciAdapter.checkConvergence(CCI_CONTINUE, &convergence); mpcciAdapter.reachSyncPoint(1); mpcciAdapter.getNodes(1); diff --git a/demo/src/mpccitest.cc b/demo/src/mpccitest.cc index a09e79f0..895b677c 100644 --- a/demo/src/mpccitest.cc +++ b/demo/src/mpccitest.cc @@ -4,37 +4,31 @@ using namespace std; using namespace AMDiS; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +// =========================================================================== +// ===== function definitions ================================================ +// =========================================================================== + +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x * x); double ux = exp(-10.0 * r2); @@ -42,7 +36,10 @@ public: } }; -// // ===== main program // + +// =========================================================================== +// ===== main program ======================================================== +// =========================================================================== int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -61,27 +58,27 @@ int main(int argc, char* argv[]) mpccitest.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("mpccitest->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("mpccitest->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("mpccitest->adapt", + AdaptStationary *adapt = new AdaptStationary("mpccitest->adapt", &mpccitest, adaptInfo); //mpccitest.setAdapt(adapt); // ===== add boundary conditions ===== - mpccitest.addDirichletBC(1, NEW G); + mpccitest.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, mpccitest.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); mpccitest.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = mpccitest.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, mpccitest.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); mpccitest.addVectorOperator(&rhsOperator); mpccitest.getMesh()->setPreserveCoarseDOFs(true); @@ -101,7 +98,7 @@ int main(int argc, char* argv[]) int quantityIDs[2] = {1, 2}; DOFVector<double> *vectors[2] = {&vec1, &vec2}; - if(codeID == 1) { + if (codeID == 1) { MpCCIAdapter mpcciAdapter(123, 1, 1, mpccitest.getFESpace(), 1, -1, 2, quantityIDs, vectors, @@ -132,7 +129,8 @@ int main(int argc, char* argv[]) vec1.print(); vec2.print(); } - if(codeID == 2) { + + if (codeID == 2) { MpCCIAdapter mpcciAdapter(124, 1, 1, mpccitest.getFESpace(), 1, -1, 2, quantityIDs, vectors, diff --git a/demo/src/multigrid.cc b/demo/src/multigrid.cc index 80d001aa..be05e06d 100644 --- a/demo/src/multigrid.cc +++ b/demo/src/multigrid.cc @@ -2,41 +2,33 @@ #include "AMDiS.h" #include "MultiGridWrapper.h" -//using namespace ost; using namespace std; using namespace AMDiS; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +// =========================================================================== +// ===== function definitions ================================================ +// =========================================================================== + +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x*x); double ux = exp(-10.0*r2); @@ -44,7 +36,10 @@ public: } }; -// ===== main program // + +// =========================================================================== +// ===== main program ======================================================== +// =========================================================================== int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -60,27 +55,27 @@ int main(int argc, char* argv[]) multigrid.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("multigrid->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("multigrid->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("multigrid->adapt", + AdaptStationary *adapt = new AdaptStationary("multigrid->adapt", &multigrid, adaptInfo); //multigrid.setAdapt(adapt); // ===== add boundary conditions ===== - multigrid.addDirichletBC(1, NEW G); + multigrid.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, multigrid.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); multigrid.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = multigrid.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, multigrid.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); multigrid.addVectorOperator(&rhsOperator); adapt->adapt(); diff --git a/demo/src/navierstokes.cc b/demo/src/navierstokes.cc deleted file mode 100644 index 73831a44..00000000 --- a/demo/src/navierstokes.cc +++ /dev/null @@ -1,355 +0,0 @@ -#include "AMDiS.h" - -using namespace std; -using namespace AMDiS; - -// =========================================================================== -// ===== function definitions ================================================ -// =========================================================================== - -class Project : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(Project); - - /** \brief - * Constructor - */ - Project(int i) : comp(i) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return x[comp]; - } - -protected: - /** \brief - * coordinate plane to project at - */ - int comp; -}; - - -/** \brief - * RHS function - */ -class ConstantFct : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(ConstantFct); - - ConstantFct(double c) - : AbstractFunction<double, WorldVector<double> >(0) , - constant(c) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return constant; - } - -protected: - double constant; -}; - -class F : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(F); - - F() - : AbstractFunction<double, WorldVector<double> >(0) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return (x[0] == 0.0) || (x[0] == 1.0) ? 0.0 : 1.0; - } - -protected: - double constant; -}; - - - -// =========================================================================== -// ===== instationary problem ================================================ -// =========================================================================== - -/** \brief - * Instationary problem - */ -class Navierstokes : public ProblemInstatVec -{ -public: - MEMORY_MANAGED(Navierstokes); - - /** \brief - * Constructor - */ - Navierstokes(ProblemVec *navierstokesSpace) - : ProblemInstatVec("navierstokes", navierstokesSpace) - { - int i, dow = Global::getGeo(WORLD); - TEST_EXIT(dow == navierstokesSpace->getNumComponents() - 1) - ("number of components != dow + 1\n"); - boundaryFct.resize(dow+1, NULL); - - epsilon = -1.0; - GET_PARAMETER(0, "epsilon", "%f", &epsilon); - TEST_EXIT(epsilon != -1.0)("invalid epsilon\n"); - }; - - // ===== initial problem methods ===================================== - - /** \brief - * Used by \ref problemInitial to solve the system of the initial problem - */ - void solve(AdaptInfo *adaptInfo) - { - int i, size = problemStat->getNumComponents()-1; - for(i = 0; i < size; i++) { - //problemStat->getMesh(i)->dofCompress(); - problemStat->getSolution()->getDOFVector(i)->interpol(boundaryFct[i]); - } - }; - - /** \brief - * Used by \ref problemInitial to do error estimation for the initial - * problem. - */ - void estimate(AdaptInfo *adaptInfo) - { - int i, size = problemStat->getNumComponents()-1; - double errMax, errSum; - - for(i = 0; i < size; i++) { - errSum = - Error<double>::L2Err(*boundaryFct[i], - *(problemStat->getSolution()->getDOFVector(i)), - 0, &errMax, false); - - adaptInfo->setEstSum(errSum, i); - } - }; - - /** \brief - * Sets \ref boundaryFct; - */ - void setBoundaryFct(AbstractFunction<double, WorldVector<double> > *fct, int i) - { - boundaryFct[i] = fct; - }; - - /** \brief - * Set the time in all needed functions! Called by adaption loop. - */ - void setTime(AdaptInfo *adaptInfo) - { - time = adaptInfo->getTime(); - tau1 = 1.0 / adaptInfo->getTimestep(); - }; - - /** \brief - * Returns address of \ref time. - */ - double *getTimePtr() { - return &time; - }; - - /** \brief - * Returns address of \ref tau1. - */ - double *getTau1Ptr() { - return &tau1; - }; - - /** \brief - * Implementation of ProblemInstatBase::closeTimestep() - * pressure correction step - */ - void closeTimestep(AdaptInfo *adaptInfo) { - problemStat->writeFiles(adaptInfo, true); - - WAIT; - - int i, numComponents = problemStat->getNumComponents(); - int dow = Global::getGeo(WORLD); - - double timestep = adaptInfo->getTimestep(); - - static DOFVector<double> pInterpol(problemStat->getFESpace(0), "p interpolated"); - static DOFVector<WorldVector<double> > gradPInterpol(problemStat->getFESpace(0), "gradient p interpolated"); - DOFVector<double> *p = problemStat->getSolution()->getDOFVector(numComponents-1); - - pInterpol.interpol(p, timestep * epsilon); - pInterpol.getRecoveryGradient(&gradPInterpol); - - for(i = 0; i < dow; i++) { - DOFVector<double>::Iterator - vIt(problemStat->getSolution()->getDOFVector(i), USED_DOFS); - DOFVector<WorldVector<double> >::Iterator - gradIt(&gradPInterpol, USED_DOFS); - - for(vIt.reset(), gradIt.reset(); !vIt.end(); ++vIt, ++gradIt) { - *vIt -= (*gradIt)[i]; - } - } - - problemStat->writeFiles(adaptInfo, true); - }; - -private: - /** \brief - * boundary function - */ - std::vector<AbstractFunction<double, WorldVector<double> >*> boundaryFct; - - /** \brief - * time - */ - double time; - - /** \brief - * 1.0/timestep - */ - double tau1; - - /** \brief - */ - double epsilon; -}; - -// =========================================================================== -// ===== main program ======================================================== -// =========================================================================== - -int main(int argc, char** argv) -{ - FUNCNAME("main"); - - // ===== check for init file ===== - TEST_EXIT(argc > 1)("usage: navierstokes initfile\n"); - - // ===== init parameters ===== - Parameters::init(false, argv[1]); - Parameters::readArgv(argc, argv); - - // read nu (1/Re) - double nu = -1.0; - GET_PARAMETER(0, "nu", "%f", &nu); - TEST_EXIT(nu != -1.0)("invalid nu\n"); - - // ===== create and init stationary problem ===== - ProblemVec *navierstokesSpace = NEW ProblemVec("navierstokes->space"); - navierstokesSpace->initialize(INIT_ALL); - - int numComponents = navierstokesSpace->getNumComponents(); - - // ===== create instationary problem ===== - Navierstokes *navierstokes = NEW Navierstokes(navierstokesSpace); - navierstokes->initialize(INIT_ALL); - - // ===== create adapt Info ===== - AdaptInfo *adaptInfo = NEW AdaptInfo("navierstokes->adapt", - numComponents); - - // create initial adapt - AdaptInfo *adaptInfoInitial = NEW AdaptInfo("navierstokes->initial->adapt", - numComponents); - - // create instationary adapt - AdaptInstationary *adaptInstat = - NEW AdaptInstationary("navierstokes->adapt", - navierstokesSpace, - adaptInfo, - navierstokes, - adaptInfoInitial); - - // ===== create boundary functions ===== - - double *tau1 = navierstokes->getTau1Ptr(); - int i; - - // === boundary functions === - ConstantFct *boundaryFct = NEW ConstantFct(0.0); - - for(i = 0; i < numComponents-1; i++) { - navierstokes->setBoundaryFct(boundaryFct, i); - navierstokesSpace->addDirichletBC(1, i, boundaryFct); - } - navierstokesSpace->addDirichletBC(2, 0, NEW F()); - navierstokes->setBoundaryFct(boundaryFct, numComponents-1); - - - // ===== create operators ===== - - // ===== velocity ===== - for(i = 0; i < numComponents-1; i++) { - // create laplace - Operator *laplaceV = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(i), - navierstokesSpace->getFESpace(i)); - - laplaceV->addSecondOrderTerm(NEW FactorLaplace_SOT(nu)); - navierstokesSpace->addMatrixOperator(laplaceV, i, i); - - // create zero order operator - Operator *dtv = NEW Operator(Operator::MATRIX_OPERATOR | - Operator::VECTOR_OPERATOR, - navierstokesSpace->getFESpace(i), - navierstokesSpace->getFESpace(i)); - - dtv->setUhOld(navierstokes->getOldSolution()->getDOFVector(i)); - dtv->addZeroOrderTerm(NEW Simple_ZOT); - - navierstokesSpace->addMatrixOperator(dtv, i, i, tau1, tau1); - navierstokesSpace->addVectorOperator(dtv, i, tau1, tau1); - - - Operator *gradP = NEW Operator(Operator::VECTOR_OPERATOR, - navierstokesSpace->getFESpace(i)); - gradP->addZeroOrderTerm(NEW FctGradient_ZOT(navierstokes->getOldSolution()-> - getDOFVector(numComponents-1), - NEW Project(i))); - double minusOne = 1.0; - navierstokesSpace->addVectorOperator(gradP, i, &minusOne, &minusOne); - } - - // ===== pressure ===== - Operator *laplaceP = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(numComponents-1), - navierstokesSpace->getFESpace(numComponents-1)); - - laplaceP->addSecondOrderTerm(NEW Laplace_SOT()); - - navierstokesSpace->addMatrixOperator(laplaceP, - numComponents-1, - numComponents-1); - - WorldVector<double> b; - - for(i = 0; i < numComponents-1; i++) { - Operator *divV = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(numComponents-1), - navierstokesSpace->getFESpace(i)); - - b.set(0.0); - b[i] = 1.0; - - divV->addFirstOrderTerm(NEW Vector_FOT(b), GRD_PHI); - - navierstokesSpace->addMatrixOperator(divV, numComponents-1, i, tau1, tau1); - } - - // ===== start adaption loop ===== - int errorCode = adaptInstat->adapt(); - - return errorCode; -} diff --git a/demo/src/neumann.cc b/demo/src/neumann.cc index 0259a1aa..8de2f500 100644 --- a/demo/src/neumann.cc +++ b/demo/src/neumann.cc @@ -6,44 +6,32 @@ using namespace std; class N : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(N); - - double operator()(const WorldVector<double>& x) const { + double operator()(const WorldVector<double>& x) const + { return 1.0; } }; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dow = x.getSize(); double r2 = (x*x); double ux = exp(-10.0*r2); @@ -67,26 +55,26 @@ int main(int argc, char* argv[]) neumann.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("neumann->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("neumann->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("neumann->adapt", + AdaptStationary *adapt = new AdaptStationary("neumann->adapt", &neumann, adaptInfo); // ===== add boundary conditions ===== - neumann.addNeumannBC(1, NEW N); - neumann.addDirichletBC(2, NEW G); + neumann.addNeumannBC(1, new N); + neumann.addDirichletBC(2, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, neumann.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); neumann.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = neumann.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, neumann.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); neumann.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/nonlin.cc b/demo/src/nonlin.cc index 00907933..4e22f24f 100644 --- a/demo/src/nonlin.cc +++ b/demo/src/nonlin.cc @@ -7,45 +7,36 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return exp(-10.0*(x*x)); + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { + return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - /** \brief - * Constructor - */ + /// Constructor F(int degree, double sigma_, double k_, int dim_) - : AbstractFunction<double, WorldVector<double> >(degree), sigma(sigma_), k(k_), dim(dim_) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - double r2 = x*x; - double ux = exp(-10.0*r2); - double ux4 = ux*ux*ux*ux; - return sigma*ux4 - k*(400.0*r2 - 20.0*dim)*ux; + : AbstractFunction<double, WorldVector<double> >(degree), + sigma(sigma_), + k(k_), + dim(dim_) + {} + + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { + double r2 = x * x; + double ux = exp(-10.0 * r2); + double ux4 = ux * ux * ux * ux; + return sigma * ux4 - k * (400.0 * r2 - 20.0 * dim) * ux; } private: @@ -54,25 +45,18 @@ private: int dim; }; -/** \brief - * Needed for zero order term. - */ +/// Needed for zero order term. class CQPFct : public AbstractFunction<double, double> { public: - MEMORY_MANAGED(CQPFct); + CQPFct() : AbstractFunction<double, double>(3) {} - CQPFct() : AbstractFunction<double, double>(3) {}; + /// Constructor. + CQPFct(double sigma_) : sigma(sigma_) {} - /** \brief - * Constructor. - */ - CQPFct(double sigma_) : sigma(sigma_) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const double& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const double& x) const + { return sigma * x * x * x; } @@ -85,41 +69,37 @@ private: // ===== class NonLin ======================================================== // =========================================================================== -/** \brief - * Non linear problem. - */ +/// Non linear problem. class NonLin : public ProblemNonLinScal { public: - /** \brief - * Constructor. - */ + /// Constructor. NonLin() : ProblemNonLinScal("nonlin"), sigma(1.0) - {}; + {} - /** \brief - * Returns \ref sigma. - */ - double getSigma() { return sigma; }; + /// Returns \ref sigma. + double getSigma() + { + return sigma; + } - /** \brief - * Sets \ref sigma. - */ - void setSigma(double sig) { sigma = sig; }; + /// Sets \ref sigma. + void setSigma(double sig) + { + sigma = sig; + } private: - /** \brief - * Stefan-Bolzmann constant. - */ + /// Stefan-Bolzmann constant. double sigma; }; + // =========================================================================== // ===== main program ======================================================== // =========================================================================== - int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -140,10 +120,10 @@ int main(int argc, char* argv[]) nonlin.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("nonlin->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("nonlin->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("nonlin->adapt", + AdaptStationary *adapt = new AdaptStationary("nonlin->adapt", &nonlin, adaptInfo); @@ -151,25 +131,25 @@ int main(int argc, char* argv[]) nonlin.setSigma(1.0); // ===== create operators ===== - Operator *nonLinOperator0 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator0 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonLinOperator0->setUhOld(nonlin.getSolution()); nonLinOperator0->addZeroOrderTerm - (NEW VecAtQP_ZOT(NULL, NEW CQPFct(nonlin.getSigma()))); + (new VecAtQP_ZOT(NULL, new CQPFct(nonlin.getSigma()))); nonlin.addMatrixOperator(nonLinOperator0, &factor1, &factor1); nonlin.addVectorOperator(nonLinOperator0); - Operator *nonLinOperator2 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator2 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonLinOperator2->setUhOld(nonlin.getSolution()); - nonLinOperator2->addSecondOrderTerm(NEW FactorLaplace_SOT(&k)); + nonLinOperator2->addSecondOrderTerm(new FactorLaplace_SOT(&k)); nonlin.addMatrixOperator(nonLinOperator2); nonlin.addVectorOperator(nonLinOperator2); @@ -178,9 +158,9 @@ int main(int argc, char* argv[]) int degree = nonlin.getFESpace()->getBasisFcts()->getDegree(); Operator* rhsFunctionOperator = - NEW Operator(Operator::VECTOR_OPERATOR, nonlin.getFESpace()); + new Operator(Operator::VECTOR_OPERATOR, nonlin.getFESpace()); rhsFunctionOperator->addZeroOrderTerm - (NEW CoordsAtQP_ZOT(NEW F(degree, + (new CoordsAtQP_ZOT(new F(degree, nonlin.getSigma(), k, nonlin.getMesh()->getDim()) @@ -190,7 +170,7 @@ int main(int argc, char* argv[]) nonlin.addVectorOperator(rhsFunctionOperator, &factor2, &factor2); // ===== add boundary conditions ===== - nonlin.addDirichletBC(DIRICHLET, NEW G); + nonlin.addDirichletBC(DIRICHLET, new G); // ===== start adaption loop ===== adapt->adapt(); diff --git a/demo/src/nonlin2.cc b/demo/src/nonlin2.cc index fc84f16c..9a590287 100644 --- a/demo/src/nonlin2.cc +++ b/demo/src/nonlin2.cc @@ -1,4 +1,4 @@ -#include "AMDiS.h" +#include "AMDi.h" using namespace std; using namespace AMDiS; @@ -7,18 +7,13 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; @@ -26,32 +21,24 @@ public: class Zero : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(Zero); - - const double& operator()(const WorldVector<double>& x) const { + const double& operator()(const WorldVector<double>& x) const + { return 0.0; - }; + } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - /** \brief - * Constructor - */ + /// Constructor F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) - {}; + {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dow = x.getSize(); double r2 = x*x; double ux = exp(-10.0*r2); @@ -60,20 +47,15 @@ public: } }; -/** \brief - * Needed for zero order term. - */ +/// Needed for zero order term. class X3 : public AbstractFunction<double, double> { public: - MEMORY_MANAGED(X3); + X3() : AbstractFunction<double, double>(3) {} - X3() : AbstractFunction<double, double>(3) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const double& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const double& x) const + { return x * x * x; } }; @@ -82,9 +64,7 @@ public: // ===== class NonLin ======================================================== // =========================================================================== -/** \brief - * Non linear problem. - */ +/// Non linear problem. class NonLin : public ProblemScal { public: @@ -97,33 +77,33 @@ public: newtonMaxIter = 50; GET_PARAMETER(0, std::string(name) + "->newton->max iteration", "%d", &newtonMaxIter); - }; + } void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) { ProblemScal::initialize(initFlag, adoptProblem, adoptFlag); - correction = NEW DOFVector<double>(this->getFESpace(), "old solution"); + correction = new DOFVector<double>(this->getFESpace(), "old solution"); correction->set(0.0); - dirichletZero = NEW DirichletBC(1, &zero, feSpace_); - dirichletG = NEW DirichletBC(1, &g, feSpace_); + dirichletZero = new DirichletBC(1, &zero, feSpace_); + dirichletG = new DirichletBC(1, &g, feSpace_); systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichletZero); solution_->getBoundaryManager()->addBoundaryCondition(dirichletG); rhs_->getBoundaryManager()->addBoundaryCondition(dirichletZero); correction->getBoundaryManager()->addBoundaryCondition(dirichletZero); - }; + } ~NonLin() { - DELETE correction; - DELETE dirichletZero; - DELETE dirichletG; - }; + delete correction; + delete dirichletZero; + delete dirichletG; + } - void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag) {}; + void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag) {} void solve(AdaptInfo *adaptInfo) { @@ -160,7 +140,7 @@ public: newtonIteration, res, newtonTolerance); } while((res > newtonTolerance) && (newtonIteration < newtonMaxIter)); solution_ = sol; - }; + } private: double newtonTolerance; @@ -190,10 +170,10 @@ int main(int argc, char* argv[]) nonlin.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("nonlin->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("nonlin->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("nonlin->adapt", + AdaptStationary *adapt = new AdaptStationary("nonlin->adapt", &nonlin, adaptInfo); @@ -203,32 +183,32 @@ int main(int argc, char* argv[]) double zero = 0.0; double minusOne = -1.0; - Operator *nonlinOperator0 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonlinOperator0 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonlinOperator0->setUhOld(nonlin.getSolution()); - nonlinOperator0->addZeroOrderTerm(NEW VecAtQP_ZOT(nonlin.getSolution(), - NEW X3)); + nonlinOperator0->addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(), + new X3)); nonlin.addMatrixOperator(nonlinOperator0, &four, &one); nonlin.addVectorOperator(nonlinOperator0, &one, &zero); - Operator *nonlinOperator2 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonlinOperator2 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonlinOperator2->setUhOld(nonlin.getSolution()); - nonlinOperator2->addSecondOrderTerm(NEW Laplace_SOT); + nonlinOperator2->addSecondOrderTerm(new Laplace_SOT); nonlin.addMatrixOperator(nonlinOperator2, &one, &one); nonlin.addVectorOperator(nonlinOperator2, &one, &zero); int degree = nonlin.getFESpace()->getBasisFcts()->getDegree(); - Operator* rhsFunctionOperator = NEW Operator(Operator::VECTOR_OPERATOR, + Operator* rhsFunctionOperator = new Operator(Operator::VECTOR_OPERATOR, nonlin.getFESpace()); - rhsFunctionOperator->addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsFunctionOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one); diff --git a/demo/src/nonlin3.cc b/demo/src/nonlin3.cc index b27beea7..7e0e1da7 100644 --- a/demo/src/nonlin3.cc +++ b/demo/src/nonlin3.cc @@ -7,18 +7,13 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; @@ -26,32 +21,24 @@ public: class Zero : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(Zero); - - double operator()(const WorldVector<double>& x) const { + double operator()(const WorldVector<double>& x) const + { return 0.0; } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - /** \brief - * Constructor - */ + /// Constructor F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) - {}; + {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dow = x.getSize(); double r2 = x *x; double ux = exp(-10.0*r2); @@ -60,20 +47,15 @@ public: } }; -/** \brief - * Needed for zero order term. - */ +/// Needed for zero order term. class X3 : public AbstractFunction<double, double> { public: - MEMORY_MANAGED(X3); - - X3() : AbstractFunction<double, double>(3) {}; + X3() : AbstractFunction<double, double>(3) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const double& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const double& x) const + { return x * x * x; } }; @@ -103,7 +85,7 @@ public: GET_PARAMETER(0, std::string(name) + "->max iteration", "%d", &newtonMaxIter); solution = problemNonlin->getSolution(); correction = newtonStep->getCorrection(); - }; + } void beginIteration(AdaptInfo *adaptInfo) @@ -143,7 +125,7 @@ public: if(toDo.isSet(ESTIMATE)) problemNonlin->estimate(adaptInfo); return flag; - }; + } void endIteration(AdaptInfo *adaptInfo) { @@ -162,7 +144,7 @@ public: if(number == 0) return problemNonlin; ERROR_EXIT("invalid problem number\n"); return NULL; - }; + } private: ProblemScal *problemNonlin; @@ -186,31 +168,31 @@ public: // newtonMaxIter = 50; // GET_PARAMETER(0, std::string(name) + "->newton->max iteration", "%d", // &newtonMaxIter); - }; + } void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) { ProblemScal::initialize(initFlag, adoptProblem, adoptFlag); - correction = NEW DOFVector<double>(this->getFESpace(), "old solution"); + correction = new DOFVector<double>(this->getFESpace(), "old solution"); correction->set(0.0); - dirichletZero = NEW DirichletBC(1, &zero, feSpace_); - dirichletG = NEW DirichletBC(1, &g, feSpace_); + dirichletZero = new DirichletBC(1, &zero, feSpace_); + dirichletG = new DirichletBC(1, &g, feSpace_); systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichletZero); solution_->getBoundaryManager()->addBoundaryCondition(dirichletG); rhs_->getBoundaryManager()->addBoundaryCondition(dirichletZero); correction->getBoundaryManager()->addBoundaryCondition(dirichletZero); - }; + } ~Nonlin() { - DELETE correction; - DELETE dirichletZero; - DELETE dirichletG; - }; + delete correction; + delete dirichletZero; + delete dirichletG; + } void initNewtonStep(AdaptInfo *adaptInfo) { solution_->getBoundaryManager()->initVector(solution_); @@ -227,19 +209,19 @@ public: tmp = solution_; solution_ = correction; - }; + } void exitNewtonStep(AdaptInfo *adaptInfo) { solution_ = tmp; - }; + } void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) { ProblemScal::buildAfterCoarsen(adaptInfo, flag); - }; + } void solveNewtonStep(AdaptInfo *adaptInfo) { ProblemScal::solve(adaptInfo); - }; + } DOFVector<double> *getCorrection() { return correction; }; @@ -272,12 +254,12 @@ int main(int argc, char* argv[]) nonlin.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("nonlin->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("nonlin->adapt", 1); NewtonMethod newtonMethod("nonlin->newton", &nonlin, &nonlin); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("nonlin->adapt", + AdaptStationary *adapt = new AdaptStationary("nonlin->adapt", &newtonMethod, adaptInfo); @@ -287,32 +269,32 @@ int main(int argc, char* argv[]) double zero = 0.0; double minusOne = -1.0; - Operator *nonlinOperator0 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonlinOperator0 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonlinOperator0->setUhOld(nonlin.getSolution()); - nonlinOperator0->addZeroOrderTerm(NEW VecAtQP_ZOT(nonlin.getSolution(), - NEW X3)); + nonlinOperator0->addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(), + new X3)); nonlin.addMatrixOperator(nonlinOperator0, &four, &one); nonlin.addVectorOperator(nonlinOperator0, &one, &zero); - Operator *nonlinOperator2 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonlinOperator2 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, nonlin.getFESpace()); nonlinOperator2->setUhOld(nonlin.getSolution()); - nonlinOperator2->addSecondOrderTerm(NEW Laplace_SOT); + nonlinOperator2->addSecondOrderTerm(new Laplace_SOT); nonlin.addMatrixOperator(nonlinOperator2, &one, &one); nonlin.addVectorOperator(nonlinOperator2, &one, &zero); int degree = nonlin.getFESpace()->getBasisFcts()->getDegree(); - Operator* rhsFunctionOperator = NEW Operator(Operator::VECTOR_OPERATOR, + Operator* rhsFunctionOperator = new Operator(Operator::VECTOR_OPERATOR, nonlin.getFESpace()); - rhsFunctionOperator->addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsFunctionOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one); diff --git a/demo/src/parallelellipt.cc b/demo/src/parallelellipt.cc index dfd85111..bbe369e4 100644 --- a/demo/src/parallelellipt.cc +++ b/demo/src/parallelellipt.cc @@ -5,49 +5,36 @@ using namespace std; using namespace AMDiS; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0*(x*x)); } }; class GrdG : public AbstractFunction<WorldVector<double>, WorldVector<double> > { - MEMORY_MANAGED(GrdG); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - WorldVector<double> operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + WorldVector<double> operator()(const WorldVector<double>& x) const + { return x * -20.0 * exp(-10.0*(x*x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x*x); double ux = exp(-10.0*r2); @@ -77,25 +64,25 @@ int main(int argc, char* argv[]) ParallelProblemScal parallelellipt("ellipt->parallel", &ellipt, NULL, vectors); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("ellipt->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("ellipt->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("ellipt->adapt", + AdaptStationary *adapt = new AdaptStationary("ellipt->adapt", ¶llelellipt, adaptInfo); // ===== add boundary conditions ===== - ellipt.addDirichletBC(1, NEW G); + ellipt.addDirichletBC(1, new G); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, ellipt.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); 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))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); ellipt.addVectorOperator(&rhsOperator); parallelellipt.initParallelization(adaptInfo); diff --git a/demo/src/parallelheat.cc b/demo/src/parallelheat.cc index 85acd3d3..868337f6 100644 --- a/demo/src/parallelheat.cc +++ b/demo/src/parallelheat.cc @@ -9,49 +9,36 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& argX) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& argX) const + { WorldVector<double> x = argX; int dim = x.getSize(); - for (int i = 0; i < dim; i++) { + for (int i = 0; i < dim; i++) x[i] -= *timePtr; - } return sin(M_PI*(*timePtr)) * exp(-10.0*(x*x)); } -}; +} -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& argX) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& argX) const + { WorldVector<double> x = argX; int dim = x.getSize(); - for (int i = 0; i < dim; i++) { + for (int i = 0; i < dim; i++) x[i] -= *timePtr; - } - double r2 = (x*x); double ux = sin(M_PI * (*timePtr)) * exp(-10.0*r2); double ut = M_PI * cos(M_PI*(*timePtr)) * exp(-10.0*r2); @@ -63,65 +50,47 @@ public: // ===== instationary problem ================================================ // =========================================================================== -/** \brief - * Instationary problem - */ +/// Instationary problem class Heat : public ProblemInstatScal { public: - MEMORY_MANAGED(Heat); - - /** \brief - * Constructor - */ + /// Constructor Heat(ProblemScal *heatSpace) : ProblemInstatScal("heat", heatSpace) - {}; - + {} // ===== ProblemInstatBase methods =================================== - /** \brief - * set the time in all needed functions! - */ - void setTime(AdaptInfo *adaptInfo) { + /// set the time in all needed functions! + void setTime(AdaptInfo *adaptInfo) + { time = adaptInfo->getTime(); tau1 = 1.0 / adaptInfo->getTimestep(); - }; + } // ===== initial problem methods ===================================== - /** \brief - * Used by \ref problemInitial to solve the system of the initial problem - */ + /// Used by \ref problemInitial to solve the system of the initial problem void solve(AdaptInfo *adaptInfo) { problemStat->getMesh()->dofCompress(); time = adaptInfo->getStartTime(); problemStat->getSolution()->interpol(boundaryFct); - }; + } - /** \brief - * Used by \ref problemInitial to do error estimation for the initial - * problem. - */ + /// Used by \ref problemInitial to do error estimation for the initial problem. void estimate(AdaptInfo *adaptInfo) { double errMax, errSum; - time = adaptInfo->getStartTime(); - errSum = Error<double>::L2Err(*boundaryFct, *(problemStat->getSolution()), 0, &errMax, false); - adaptInfo->setEstSum(errSum, 0); - }; + } // ===== setting methods =============================================== - /** \brief - * Sets \ref boundaryFct; - */ + /// Sets \ref boundaryFct; void setBoundaryFct(AbstractFunction<double, WorldVector<double> > *fct) { boundaryFct = fct; @@ -129,30 +98,23 @@ public: // ===== getting methods =============================================== - - /** \brief - * Returns pointer to \ref tau1 - */ + /// Returns pointer to \ref tau1 double *getTau1Ptr() { return &tau1; - }; + } - /** \brief - * Returns pointer to \ref time. - */ + /// Returns pointer to \ref time. double *getTimePtr() { return &time; - }; + } - /** \brief - * Returns \ref boundaryFct; - */ + /// Returns \ref boundaryFct; AbstractFunction<double, WorldVector<double> > *getBoundaryFct() { return boundaryFct; - }; + } // void closeTimestep(AdaptInfo *adaptInfo) { @@ -162,19 +124,13 @@ public: // }; private: - /** \brief - * 1.0 / timestep - */ + /// 1.0 / timestep double tau1; - /** \brief - * time - */ + /// time double time; - /** \brief - * Pointer to boundary function. Needed for initial problem. - */ + /// Pointer to boundary function. Needed for initial problem. AbstractFunction<double, WorldVector<double> > *boundaryFct; }; @@ -193,7 +149,7 @@ int main(int argc, char** argv) Parameters::init(false, argv[1]); // ===== create and init stationary problem ===== - ProblemScal *heatSpace = NEW ProblemScal("heat->space"); + ProblemScal *heatSpace = new ProblemScal("heat->space"); heatSpace->initialize(INIT_ALL); // ===== create instationary problem ===== @@ -201,10 +157,10 @@ int main(int argc, char** argv) heat->initialize(INIT_ALL); // create adapt info - AdaptInfo *adaptInfo = NEW AdaptInfo("heat->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("heat->adapt"); // create initial adapt info - AdaptInfo *adaptInfoInitial = NEW AdaptInfo("heat->initial->adapt"); + AdaptInfo *adaptInfoInitial = new AdaptInfo("heat->initial->adapt"); // create parallel problem std::vector<DOFVector<double>*> vectors; @@ -212,14 +168,14 @@ int main(int argc, char** argv) // create instationary adapt AdaptInstationary *adaptInstat = - NEW AdaptInstationary("heat->adapt", + new AdaptInstationary("heat->adapt", ¶llelheat, adaptInfo, ¶llelheat, adaptInfoInitial); // ===== create boundary functions ===== - G *boundaryFct = NEW G; + G *boundaryFct = new G; boundaryFct->setTimePtr(heat->getTimePtr()); heat->setBoundaryFct(boundaryFct); @@ -227,7 +183,7 @@ int main(int argc, char** argv) // ===== create rhs functions ===== int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree(); - F *rhsFct = NEW F(degree); + F *rhsFct = new F(degree); rhsFct->setTimePtr(heat->getTimePtr()); // ===== create operators ===== @@ -235,7 +191,7 @@ int main(int argc, char** argv) double zero = 0.0; // create laplace - Operator *A = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *A = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); @@ -246,11 +202,11 @@ int main(int argc, char** argv) heatSpace->addMatrixOperator(A, &one, &one); // create zero order operator - Operator *C = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *C = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); - C->addZeroOrderTerm(NEW Simple_ZOT); + C->addZeroOrderTerm(new Simple_ZOT); C->setUhOld(heat->getOldSolution()); @@ -258,10 +214,10 @@ int main(int argc, char** argv) heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr()); // create RHS operator - Operator *F = NEW Operator(Operator::VECTOR_OPERATOR, + Operator *F = new Operator(Operator::VECTOR_OPERATOR, heatSpace->getFESpace()); - F->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct)); + F->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct)); heatSpace->addVectorOperator(F); diff --git a/demo/src/parametric.cc b/demo/src/parametric.cc index 1808b6d0..90ead584 100644 --- a/demo/src/parametric.cc +++ b/demo/src/parametric.cc @@ -7,20 +7,15 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return -2.0 * x[0]; } }; @@ -33,27 +28,22 @@ public: class ParametricSphere : public ProblemScal { public: - /** \brief - * constructor - */ + /// constructor ParametricSphere(const char *name) : ProblemScal(name), parametric(NULL) - {}; + {} - /** \brief - * destructor - */ - ~ParametricSphere() { - for (int i = 0; i < Global::getGeo(WORLD); i++) { - DELETE parametricCoords[i]; - } - DELETE parametric; + /// destructor + ~ParametricSphere() + { + for (int i = 0; i < Global::getGeo(WORLD); i++) + delete parametricCoords[i]; + + delete parametric; } - /** \brief - * initialization of the base class and creation of \ref parametric. - */ + /// initialization of the base class and creation of \ref parametric. void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) @@ -63,20 +53,18 @@ public: // ===== create projection ===== WorldVector<double> ballCenter; ballCenter.set(0.0); - NEW BallProject(1, + new BallProject(1, VOLUME_PROJECTION, ballCenter, 1.0); // ===== create coords vector ===== - int i; - for(i = 0; i < Global::getGeo(WORLD); i++) { - parametricCoords[i] = NEW DOFVector<double>(this->getFESpace(), + for (itn i = 0; i < Global::getGeo(WORLD); i++) { + parametricCoords[i] = new DOFVector<double>(this->getFESpace(), "parametric coords"); - } // ===== create parametric object ===== - parametric = NEW ParametricFirstOrder(¶metricCoords); + parametric = new ParametricFirstOrder(¶metricCoords); // ===== enable parametric traverse ===== this->getMesh()->setParametric(parametric); @@ -88,7 +76,6 @@ public: void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) { FUNCNAME("ParametricSphere::buildAfterCoarsen()"); MSG("calculation of parametric coordinates\n"); - int i, j; int preDOFs = this->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX); int dim = this->getMesh()->getDim(); int dow = Global::getGeo(WORLD); @@ -103,14 +90,13 @@ public: ElInfo *elInfo = NULL; elInfo = stack.traverseFirst(this->getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); - while(elInfo) { + while (elInfo) { dof = elInfo->getElement()->getDOF(); - for(i = 0; i < dim + 1; i++) { + for (int i = 0; i < dim + 1; i++) { vertexCoords = elInfo->getCoord(i); Projection::getProjection(1)->project(vertexCoords); - for(j = 0; j < dow; j++) { + for (int j = 0; j < dow; j++) (*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j]; - } } elInfo = stack.traverseNext(elInfo); } @@ -118,17 +104,13 @@ public: // ===== enable parametric traverse ===== this->getMesh()->setParametric(parametric); - }; + } protected: - /** \brief - * DOFVector storing parametric coordinates. - */ + /// DOFVector storing parametric coordinates. WorldVector<DOFVector<double>*> parametricCoords; - /** \brief - * Parametrizes vertex coordinates while mesh traversal. - */ + /// Parametrizes vertex coordinates while mesh traversal. ParametricFirstOrder *parametric; }; @@ -151,17 +133,17 @@ int main(int argc, char* argv[]) parametric.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("parametric->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("parametric->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("parametric->adapt", + AdaptStationary *adapt = new AdaptStationary("parametric->adapt", ¶metric, adaptInfo); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, parametric.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); parametric.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== @@ -170,7 +152,7 @@ int main(int argc, char* argv[]) int degree = parametric.getFESpace()->getBasisFcts()->getDegree(); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); parametric.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/periodic.cc b/demo/src/periodic.cc index 8219a71b..99a1f164 100644 --- a/demo/src/periodic.cc +++ b/demo/src/periodic.cc @@ -3,45 +3,33 @@ using namespace std; using namespace AMDiS; -// ===== function definitions // -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return 0.0; } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x*x); double ux = exp(-10.0*r2); return -(400.0*r2 - 20.0*dim)*ux; - }; + } }; -// // ===== main program // int main(int argc, char* argv[]) { FUNCNAME("main"); @@ -57,28 +45,28 @@ int main(int argc, char* argv[]) periodic.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("periodic->adapt", 1); + AdaptInfo *adaptInfo = new AdaptInfo("periodic->adapt", 1); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("periodic->adapt", + AdaptStationary *adapt = new AdaptStationary("periodic->adapt", &periodic, adaptInfo); // ===== add boundary conditions ===== - periodic.addDirichletBC(1, NEW G); + periodic.addDirichletBC(1, new G); // ===== add periodic conditions ===== periodic.addPeriodicBC(-1); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, periodic.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); periodic.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== int degree = periodic.getFESpace()->getBasisFcts()->getDegree(); Operator rhsOperator(Operator::VECTOR_OPERATOR, periodic.getFESpace()); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); periodic.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/smitest.cc b/demo/src/smitest.cc index ba939365..348e5565 100644 --- a/demo/src/smitest.cc +++ b/demo/src/smitest.cc @@ -5,8 +5,6 @@ using namespace ost; using namespace std; - -// // ===== main program // int main(int argc, char* argv[]) { SMI_Connect_to_server("10.1.13.3", 5432); @@ -21,9 +19,8 @@ int main(int argc, char* argv[]) SMI_Get_all_mesh_ids(applicationID, &numMeshes, &meshIDs); std::cout << numMeshes << std::endl; - for(int i = 0; i < numMeshes; i++) { + for (int i = 0; i < numMeshes; i++) std::cout << meshIDs[i] << " "; - } std::cout << std::endl; SMI_Begin_read_transaction(applicationID, 1); diff --git a/demo/src/sphere.cc b/demo/src/sphere.cc index 28d7eb0e..cc70c3d5 100644 --- a/demo/src/sphere.cc +++ b/demo/src/sphere.cc @@ -7,20 +7,15 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return -2.0 * x[0]; } }; @@ -42,7 +37,7 @@ int main(int argc, char* argv[]) // ===== create projection ===== WorldVector<double> ballCenter; ballCenter.set(0.0); - NEW BallProject(1, + new BallProject(1, BOUNDARY_PROJECTION, ballCenter, 1.0); @@ -52,17 +47,17 @@ int main(int argc, char* argv[]) sphere.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("sphere->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("sphere->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("sphere->adapt", + AdaptStationary *adapt = new AdaptStationary("sphere->adapt", &sphere, adaptInfo); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, sphere.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); sphere.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== @@ -71,7 +66,7 @@ int main(int argc, char* argv[]) int degree = sphere.getFESpace()->getBasisFcts()->getDegree(); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); sphere.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== diff --git a/demo/src/stokesnonlin.cc b/demo/src/stokesnonlin.cc deleted file mode 100644 index 220bf1a5..00000000 --- a/demo/src/stokesnonlin.cc +++ /dev/null @@ -1,493 +0,0 @@ -#include "AMDiS.h" - -using namespace std; -using namespace AMDiS; - -// =========================================================================== -// ===== function definitions ================================================ -// =========================================================================== - -class Project : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(Project); - - /** \brief - * Constructor - */ - Project(int i) : comp(i) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return x[comp]; - } - -protected: - /** \brief - * coordinate plane to project at - */ - int comp; -}; - - -/** \brief - * RHS function - */ -class ConstantFct : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(ConstantFct); - - ConstantFct(double c) - : AbstractFunction<double, WorldVector<double> >(0) , - constant(c) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return constant; - } - -protected: - double constant; -}; - -class F : public AbstractFunction<double, WorldVector<double> > -{ -public: - MEMORY_MANAGED(F); - - F() - : AbstractFunction<double, WorldVector<double> >(0) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return (x[0] == 0.0) || (x[0] == 1.0) ? 0.0 : 1.0; - } - -protected: - double constant; -}; - -// =========================================================================== -// ===== Updater ============================================================= -// =========================================================================== - -/** \brief - * Updater for the Stokes operator - */ - -class StokesUpdater: public NonLinUpdaterVec -{ -public: - /** \brief - * Constructor. - */ - StokesUpdater(Matrix<DOFMatrix*> *m, SystemVector *sol) : - NonLinUpdaterVec(m),sol_(sol) - {}; - - - virtual void update(bool updateDF, SystemVector *F) { - - FUNCNAME("StokesUpdater::update()"); - - int i, j, size = df->getNumRows(); - Matrix<DOFMatrix*>* lhsMatrix=NULL; - static Matrix<DOFMatrix*>* lhsMat=NULL; // Storage for lhs matrix if DF is not updated - - TraverseStack stack; - ElInfo *el_info; - Flag fill_flag; - - TEST_EXIT(F || df)("neither F nor df are set\n"); - - if ((!updateDF) && (F==NULL)) { - WARNING("No RHS vector and no update for system matrix! Updater does nothing!\n"); - return; - } - - const FiniteElemSpace *feSpace = - F ? - F->getDOFVector(0)->getFESpace() : - (*df)[0][0]->getFESpace(); - - if (updateDF) { - TEST_EXIT(df)("df not set but update tried!\n"); - for(i = 0; i < size; i++) { - for(j = 0; j < size; j++) { - if((*df)[i][j]) { - (*df)[i][j]->clear(); - } - } - } - lhsMatrix=df; // lhs matrix is DF, which has to be updated. - } - else { // DF must not be changed, so we provide storage for lhs matrix. - if (lhsMat==NULL) { - lhsMat=NEW Matrix<DOFMatrix*>(size,size); - for(i = 0; i < size; i++) { - for(j = 0; j < size; j++) { - if ((!df)||(*df)[i][j]) { - (*lhsMat)[i][j]=NEW DOFMatrix(*((*df)[i][j])); - (*lhsMat)[i][j]->clear(); - } - else { - (*lhsMat)[i][j]=NULL; - } - } - } - } - lhsMatrix=lhsMat; - } - - BasisFunction *basFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts()); - - if (F) { // clear RHS vector if neccessary. - for(i = 0; i < size; i++) { - F->getDOFVector(i)->set(0.0); - } - } - - fill_flag = - Mesh::CALL_LEAF_EL| - Mesh::FILL_COORDS| - Mesh::FILL_BOUND| - Mesh::FILL_DET| - Mesh::FILL_GRD_LAMBDA; - - el_info = stack.traverseFirst(feSpace->getMesh(), -1, fill_flag); - while (el_info) { - const BoundaryType *bound = basFcts->getBound(el_info, NULL); - // assemble LHS matrix. - for(i = 0; i < size; i++) { - for(j = 0; j < size; j++) { - if((*lhsMatrix)[i][j]) { - (*lhsMatrix)[i][j]->assemble(1.0, el_info, bound); - if((*lhsMatrix)[i][j]->getBoundaryManager()) - (*lhsMatrix)[i][j]->getBoundaryManager()->fillLocalBCs(el_info, - (*lhsMatrix)[i][j]); - } - } - } - - - if (F) { // assemble RHS vector if neccessary - for(i = 0; i < size; i++) { - F->getDOFVector(i)->assemble(1.0, el_info, bound); - if(F->getDOFVector(i)->getBoundaryManager()) - F->getDOFVector(i)->getBoundaryManager()->fillLocalBCs(el_info, - F->getDOFVector(i)); - } - } - - el_info = stack.traverseNext(el_info); - } - - if (F) { - SystemVector lonR(*F); - - // final RHS vector is LHS matrix times old solution minus assembled RHS vector - // i.e. the residual of the linear equation. - - StandardMatVec<Matrix<DOFMatrix*>,SystemVector> matVec(lhsMatrix); - - matVec.matVec(NoTranspose,*sol_,lonR); - *F *= -1.0; - *F +=lonR; - - } - }; - -private: - /** \brief - * Solution of the nonlinear Problem. Needed to calculate RHS contribution of Matrix. - */ - SystemVector* sol_; - -}; - -// =========================================================================== -// ===== instationary problem ================================================ -// =========================================================================== - -/** \brief - * Instationary problem - */ -class Navierstokes : public ProblemInstatVec -{ -public: - MEMORY_MANAGED(Navierstokes); - - /** \brief - * Constructor - */ - Navierstokes(ProblemVec *navierstokesSpace) - : ProblemInstatVec("navierstokes", navierstokesSpace) - { - int i, dow = Global::getGeo(WORLD); - TEST_EXIT(dow == navierstokesSpace->getNumComponents() - 1) - ("number of components != dow + 1\n"); - boundaryFct.resize(dow+1, NULL); - - epsilon = -1.0; - GET_PARAMETER(0, "epsilon", "%f", &epsilon); - TEST_EXIT(epsilon != -1.0)("invalid epsilon\n"); - }; - - // ===== initial problem methods ===================================== - - /** \brief - * Used by \ref problemInitial to solve the system of the initial problem - */ - void solve(AdaptInfo *adaptInfo) - { - int i, size = problemStat->getNumComponents()-1; - for(i = 0; i < size; i++) { - //problemStat->getMesh(i)->dofCompress(); - problemStat->getSolution()->getDOFVector(i)->interpol(boundaryFct[i]); - } - }; - - /** \brief - * Used by \ref problemInitial to do error estimation for the initial - * problem. - */ - void estimate(AdaptInfo *adaptInfo) - { - int i, size = problemStat->getNumComponents()-1; - double errMax, errSum; - - for(i = 0; i < size; i++) { - errSum = - Error<double>::L2Err(*boundaryFct[i], - *(problemStat->getSolution()->getDOFVector(i)), - 0, &errMax, false); - - adaptInfo->setEstSum(errSum, i); - } - }; - - /** \brief - * Sets \ref boundaryFct; - */ - void setBoundaryFct(AbstractFunction<double, WorldVector<double> > *fct, int i) - { - boundaryFct[i] = fct; - }; - - /** \brief - * Set the time in all needed functions! Called by adaption loop. - */ - void setTime(AdaptInfo *adaptInfo) - { - time = adaptInfo->getTime(); - tau1 = 1.0 / adaptInfo->getTimestep(); - }; - - /** \brief - * Returns address of \ref time. - */ - double *getTimePtr() { - return &time; - }; - - /** \brief - * Returns address of \ref tau1. - */ - double *getTau1Ptr() { - return &tau1; - }; - - /** \brief - * Implementation of ProblemInstatBase::closeTimestep() - * pressure correction step - */ - void closeTimestep(AdaptInfo *adaptInfo) { - problemStat->writeFiles(adaptInfo, true); - - WAIT; - - int i, numComponents = problemStat->getNumComponents(); - int dow = Global::getGeo(WORLD); - - double timestep = adaptInfo->getTimestep(); - - static DOFVector<double> pInterpol(problemStat->getFESpace(0), "p interpolated"); - static DOFVector<WorldVector<double> > gradPInterpol(problemStat->getFESpace(0), "gradient p interpolated"); - DOFVector<double> *p = problemStat->getSolution()->getDOFVector(numComponents-1); - - pInterpol.interpol(p, timestep * epsilon); - pInterpol.getRecoveryGradient(&gradPInterpol); - - for(i = 0; i < dow; i++) { - DOFVector<double>::Iterator - vIt(problemStat->getSolution()->getDOFVector(i), USED_DOFS); - DOFVector<WorldVector<double> >::Iterator - gradIt(&gradPInterpol, USED_DOFS); - - for(vIt.reset(), gradIt.reset(); !vIt.end(); ++vIt, ++gradIt) { - *vIt -= (*gradIt)[i]; - } - } - - problemStat->writeFiles(adaptInfo, true); - }; - -private: - /** \brief - * boundary function - */ - std::vector<AbstractFunction<double, WorldVector<double> >*> boundaryFct; - - /** \brief - * time - */ - double time; - - /** \brief - * 1.0/timestep - */ - double tau1; - - /** \brief - */ - double epsilon; -}; - -// =========================================================================== -// ===== main program ======================================================== -// =========================================================================== - -int main(int argc, char** argv) -{ - FUNCNAME("main"); - - // ===== check for init file ===== - TEST_EXIT(argc == 2)("usage: navierstokes initfile\n"); - - // ===== init parameters ===== - Parameters::init(false, argv[1]); - - // read nu (1/Re) - double nu = -1.0; - GET_PARAMETER(0, "nu", "%f", &nu); - TEST_EXIT(nu != -1.0)("invalid nu\n"); - - // ===== create and init stationary problem ===== - ProblemNonLinVec *navierstokesSpace = NEW ProblemNonLinVec("navierstokes->space"); - navierstokesSpace->initialize(INIT_ALL&~INIT_UPDATER&~INIT_NONLIN_SOLVER); - StokesUpdater* stokesUpdater=new StokesUpdater(navierstokesSpace->getSystemMatrix(), - navierstokesSpace->getSolution()); - navierstokesSpace->setUpdater(stokesUpdater); - - //navierstokesSpace->setUpdater(NEW NonLinUpdaterVec(navierstokesSpace->getSystemMatrix())); - navierstokesSpace->initialize(INIT_NONLIN_SOLVER); - - int numComponents = navierstokesSpace->getNumComponents(); - - // ===== create instationary problem ===== - Navierstokes *navierstokes = NEW Navierstokes(navierstokesSpace); - navierstokes->initialize(INIT_ALL); - - AdaptInfo *adaptInfo = NEW AdaptInfo("navierstokes->adapt", - numComponents); - - // create initial adapt - AdaptInfo *adaptInfoInitial = NEW AdaptInfo("navierstokes->initial->adapt", - numComponents); - - // create instationary adapt - AdaptInstationary *adaptInstat = - NEW AdaptInstationary("navierstokes->adapt", - navierstokesSpace, - adaptInfo, - navierstokes, - adaptInfoInitial); - - // ===== create boundary functions ===== - - double *tau1 = navierstokes->getTau1Ptr(); - int i; - - // === boundary functions === - ConstantFct *boundaryFct = NEW ConstantFct(0.0); - - for(i = 0; i < numComponents-1; i++) { - navierstokes->setBoundaryFct(boundaryFct, i); - navierstokesSpace->addDirichletBC(1, i, boundaryFct); - } - navierstokesSpace->addDirichletBC(2, 0, NEW F()); - navierstokesSpace->addDirichletBC(2, 1, boundaryFct); - navierstokes->setBoundaryFct(boundaryFct, numComponents-1); - - - // ===== create operators ===== - - // ===== velocity ===== - for(i = 0; i < numComponents-1; i++) { - // create laplace - Operator *laplaceV = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(i), - navierstokesSpace->getFESpace(i)); - - laplaceV->addSecondOrderTerm(NEW FactorLaplace_SOT(nu)); - navierstokesSpace->addMatrixOperator(laplaceV, i, i); - - // create zero order operator - Operator *dtv = NEW Operator(Operator::MATRIX_OPERATOR | - Operator::VECTOR_OPERATOR, - navierstokesSpace->getFESpace(i), - navierstokesSpace->getFESpace(i)); - - dtv->setUhOld(navierstokes->getOldSolution()->getDOFVector(i)); - dtv->addZeroOrderTerm(NEW Simple_ZOT); - - navierstokesSpace->addMatrixOperator(dtv, i, i, tau1, tau1); - navierstokesSpace->addVectorOperator(dtv, i, tau1, tau1); - - Operator *gradP = NEW Operator(Operator::VECTOR_OPERATOR, - navierstokesSpace->getFESpace(i)); - gradP->addZeroOrderTerm(NEW FctGradient_ZOT(navierstokes->getOldSolution()-> - getDOFVector(numComponents-1), - NEW Project(i))); - double minusOne = 1.0; - navierstokesSpace->addVectorOperator(gradP, i, &minusOne, &minusOne); - } - - // ===== pressure ===== - Operator *laplaceP = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(numComponents-1), - navierstokesSpace->getFESpace(numComponents-1)); - - laplaceP->addSecondOrderTerm(NEW Laplace_SOT()); - - navierstokesSpace->addMatrixOperator(laplaceP, - numComponents-1, - numComponents-1); - - WorldVector<double> b; - - for(i = 0; i < numComponents-1; i++) { - Operator *divV = NEW Operator(Operator::MATRIX_OPERATOR, - navierstokesSpace->getFESpace(numComponents-1), - navierstokesSpace->getFESpace(i)); - - b.set(0.0); - b[i] = 1.0; - - divV->addFirstOrderTerm(NEW Vector_FOT(b), GRD_PHI); - - navierstokesSpace->addMatrixOperator(divV, numComponents-1, i, tau1, tau1); - } - - // ===== start adaption loop ===== - adaptInstat->adapt(); -} diff --git a/demo/src/torus.cc b/demo/src/torus.cc index 964ba77c..580fb7ef 100644 --- a/demo/src/torus.cc +++ b/demo/src/torus.cc @@ -18,31 +18,27 @@ public: // ===== function definitions ================================================ // =========================================================================== -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree), rotation(0.0) - {}; + {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { WorldVector<double> myX = x; YRotation::rotate(myX, -rotation); return -2.0 * myX[0]; } - void rotate(double r) { + void rotate(double r) + { rotation += r; - }; + } private: double rotation; @@ -51,9 +47,7 @@ private: class TorusProject : public Projection { public: - /** \brief - * Constructor. - */ + /// Constructor. TorusProject(int id, ProjectionType type, double radius1, @@ -61,17 +55,14 @@ public: : Projection(id, type), radius1_(radius1), radius2_(radius2) - {}; + {} - /** \brief - * Destructor. - */ - virtual ~TorusProject() {}; + /// Destructor. + virtual ~TorusProject() {} - /** \brief - * Implementation of Projection::project(); - */ - void project(WorldVector<double> &x) { + /// Implementation of Projection::project(); + void project(WorldVector<double> &x) + { WorldVector<double> xPlane = x; xPlane[2] = 0.0; @@ -92,12 +83,8 @@ public: }; protected: - /** \brief - */ double radius1_; - /** \brief - */ double radius2_; }; @@ -112,10 +99,10 @@ int main(int argc, char* argv[]) Parameters::init(false, argv[1]); // ===== create projection ===== - double r2 = (1.5 - 1.0/sqrt(2.0)) / 2.0; - double r1 = 1.0/sqrt(2.0) + r2; + double r2 = (1.5 - 1.0 / sqrt(2.0)) / 2.0; + double r1 = 1.0 / sqrt(2.0) + r2; - NEW TorusProject(1, + new TorusProject(1, VOLUME_PROJECTION, r1, r2); @@ -125,17 +112,17 @@ int main(int argc, char* argv[]) torus.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("torus->adapt"); + AdaptInfo *adaptInfo = new AdaptInfo("torus->adapt"); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("torus->adapt", + AdaptStationary *adapt = new AdaptStationary("torus->adapt", &torus, adaptInfo); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, torus.getFESpace()); - matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator.addSecondOrderTerm(new Laplace_SOT); torus.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== @@ -145,7 +132,7 @@ int main(int argc, char* argv[]) int degree = torus.getFESpace()->getBasisFcts()->getDegree(); F f(degree); - rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(&f)); + rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(&f)); torus.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== @@ -170,7 +157,7 @@ int main(int argc, char* argv[]) WorldVector<DOFVector<double>*> parametricCoords; for(i = 0; i < dow; i++) { - parametricCoords[i] = NEW DOFVector<double>(feSpace, + parametricCoords[i] = new DOFVector<double>(feSpace, "parametric coords"); } @@ -226,7 +213,7 @@ int main(int argc, char* argv[]) adapt->adapt(); for(i = 0; i < dow; i++) - DELETE parametricCoords[i]; + delete parametricCoords[i]; FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts); } diff --git a/demo/src/vecellipt.cc b/demo/src/vecellipt.cc index 02981630..f8e2103b 100644 --- a/demo/src/vecellipt.cc +++ b/demo/src/vecellipt.cc @@ -7,36 +7,26 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x * x); double ux = exp(-10.0 * r2); @@ -63,23 +53,23 @@ int main(int argc, char* argv[]) vecellipt.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("vecellipt->adapt", + AdaptInfo *adaptInfo = new AdaptInfo("vecellipt->adapt", vecellipt.getNumComponents()); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("vecellipt->adapt", + AdaptStationary *adapt = new AdaptStationary("vecellipt->adapt", &vecellipt, adaptInfo); // ===== add boundary conditions ===== - vecellipt.addDirichletBC(1, 0, NEW G); + 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); + matrixOperator00.addSecondOrderTerm(new Laplace_SOT); vecellipt.addMatrixOperator(&matrixOperator00, 0, 0); Operator rhsOperator0(Operator::VECTOR_OPERATOR, @@ -87,7 +77,7 @@ int main(int argc, char* argv[]) int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree(); - rhsOperator0.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); vecellipt.addVectorOperator(&rhsOperator0, 0); @@ -99,11 +89,11 @@ int main(int argc, char* argv[]) vecellipt.getFESpace(1), vecellipt.getFESpace(1)); - matrixOperator10.addZeroOrderTerm(NEW Simple_ZOT); + matrixOperator10.addZeroOrderTerm(new Simple_ZOT); vecellipt.addMatrixOperator(&matrixOperator10, 1, 0); - matrixOperator11.addZeroOrderTerm(NEW Simple_ZOT(-1.0)); + matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0)); vecellipt.addMatrixOperator(&matrixOperator11, 1, 1); diff --git a/demo/src/vecheat.cc b/demo/src/vecheat.cc index 86093642..8fad0bcd 100644 --- a/demo/src/vecheat.cc +++ b/demo/src/vecheat.cc @@ -7,38 +7,29 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(G); - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x)); - }; + } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> >, public TimedObject { public: - MEMORY_MANAGED(F); + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = x * x; double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2); @@ -51,17 +42,12 @@ public: // ===== instationary problem ================================================ // =========================================================================== -/** \brief - * Instationary problem - */ +/// Instationary problem class Vecheat : public ProblemInstatVec { public: - MEMORY_MANAGED(Vecheat); - /** \brief - * Constructor - */ + /// Constructor Vecheat(ProblemVec *vecheatSpace) : ProblemInstatVec("vecheat", vecheatSpace) { @@ -75,116 +61,83 @@ public: } MSG("theta = %f\n", theta); theta1 = theta - 1; - }; + } // ===== ProblemInstatBase methods =================================== - /** \brief - * set the time in all needed functions! - */ + /// set the time in all needed functions! void setTime(AdaptInfo *adaptInfo) { rhsTime = adaptInfo->getTime() - (1-theta) * adaptInfo->getTimestep(); boundaryTime = adaptInfo->getTime(); tau1 = 1.0 / adaptInfo->getTimestep(); - }; + } // ===== initial problem methods ===================================== - /** \brief - * Used by \ref problemInitial to solve the system of the initial problem - */ + /// Used by \ref problemInitial to solve the system of the initial problem void solveInitialProblem(AdaptInfo *adaptInfo) { - int i, size = problemStat->getNumComponents(); - + int size = problemStat->getNumComponents(); boundaryTime = rhsTime = 0.0; - for(i = 0; i < size; i++) { + for (int i = 0; i < size; i++) { problemStat->getMesh(i)->dofCompress(); problemStat->getSolution()->getDOFVector(i)->interpol(boundaryFct); } - }; + } - /** \brief - * Used by \ref problemInitial to do error estimation for the initial - * problem. - */ + /// Used by \ref problemInitial to do error estimation for the initial problem. void estimateInitial(AdaptInfo *adaptInfo) { - int i, size = problemStat->getNumComponents(); + int size = problemStat->getNumComponents(); double errMax, errSum; - boundaryTime = 0.0; - for(i = 0; i < size; i++) { + for (int i = 0; i < size; i++) { errSum = Error<double>::L2Err(*boundaryFct, *(problemStat->getSolution()-> getDOFVector(i)), 0, &errMax, false); - adaptInfo->setEstSum(errSum, i); } - }; + } + + /// Used by \ref problemInitial to build before refinement. + void buildBeforeRefineInitial(Flag) {} + + /// Used by \ref problemInitial to build before coarsening. + void buildBeforeCoarsenInitial(Flag) {} + + /// Used by \ref problemInitial to build after coarsening. + void buildAfterCoarsenInitial(Flag) {} - /** \brief - * Used by \ref problemInitial to build before refinement. - */ - void buildBeforeRefineInitial(Flag) - {}; - - /** \brief - * Used by \ref problemInitial to build before coarsening. - */ - void buildBeforeCoarsenInitial(Flag) - {}; - - /** \brief - * Used by \ref problemInitial to build after coarsening. - */ - void buildAfterCoarsenInitial(Flag) - {}; - - /** \brief - * Used by \ref problemInitial to mark elements. - */ + /// Used by \ref problemInitial to mark elements. Flag markElementsInitial() { return 0; - }; + } - /** \brief - * Used by \ref problemInitial - */ + /// Used by \ref problemInitial Flag refineMeshInitial() { return 0; - }; + } - /** \brief - * Used by \ref problemInitial - */ + /// Used by \ref problemInitial Flag coarsenMeshInitial() { return 0; - }; + } - /** \brief - * Used by \ref problemInitial - */ - void beginIterationInitial(int ) - {}; + /// Used by \ref problemInitial + void beginIterationInitial(int) {} - /** \brief - * Used by \ref problemInitial - */ - void endIterationInitial(int ) - {}; + /// Used by \ref problemInitial + void endIterationInitial(int) {} // ===== setting methods =============================================== - /** \brief - * Sets \ref boundaryFct; - */ + /// Sets \ref boundaryFct; void setBoundaryFct(AbstractFunction<double, WorldVector<double> > *fct) { boundaryFct = fct; @@ -192,83 +145,59 @@ public: // ===== getting methods =============================================== - /** \brief - * Returns pointer to \ref theta. - */ + /// Returns pointer to \ref theta. double *getThetaPtr() { return θ - }; + } - /** \brief - * Returns pointer to \ref theta1. - */ + /// Returns pointer to \ref theta1. double *getTheta1Ptr() { return &theta1; - }; + } - /** \brief - * Returns pointer to \ref tau1 - */ + /// Returns pointer to \ref tau1 double *getTau1Ptr() { return &tau1; - }; + } - /** \brief - * Returns pointer to \ref rhsTime. - */ + /// Returns pointer to \ref rhsTime. double *getRHSTimePtr() { return &rhsTime; - }; + } - /** \brief - * Returns pointer to \ref theta1. - */ + /// Returns pointer to \ref theta1. double *getBoundaryTimePtr() { return &boundaryTime; - }; + } - /** \brief - * Returns \ref boundaryFct; - */ + /// Returns \ref boundaryFct; AbstractFunction<double, WorldVector<double> > *getBoundaryFct() { return boundaryFct; - }; + } protected: - /** \brief - * Used for theta scheme. - */ + /// Used for theta scheme. double theta; - /** \brief - * theta - 1 - */ + /// theta - 1 double theta1; - /** \brief - * 1.0 / timestep - */ + /// 1.0 / timestep double tau1; - /** \brief - * time for right hand side functions. - */ + /// time for right hand side functions. double rhsTime; - /** \brief - * time for boundary functions. - */ + /// time for boundary functions. double boundaryTime; - /** \brief - * Pointer to boundary function. Needed for initial problem. - */ + /// Pointer to boundary function. Needed for initial problem. AbstractFunction<double, WorldVector<double> > *boundaryFct; }; @@ -288,26 +217,26 @@ int main(int argc, char** argv) Parameters::readArgv(argc, argv); // ===== create and init stationary problem ===== - ProblemVec *vecheatSpace = NEW ProblemVec("vecheat->space"); + ProblemVec *vecheatSpace = new ProblemVec("vecheat->space"); vecheatSpace->initialize(INIT_ALL); // ===== create instationary problem ===== - Vecheat *vecheat = NEW Vecheat(vecheatSpace); + Vecheat *vecheat = new Vecheat(vecheatSpace); vecheat->initialize(INIT_ALL); // create adapt info AdaptInfo *adaptInfo = - NEW AdaptInfo("vecheat->adapt", + new AdaptInfo("vecheat->adapt", vecheatSpace->getNumComponents()); // create initial adapt info AdaptInfo *adaptInfoInitial = - NEW AdaptInfo("vecheat->initial->adapt", + new AdaptInfo("vecheat->initial->adapt", vecheatSpace->getNumComponents()); // create instationary adapt AdaptInstationary *adaptInstat = - NEW AdaptInstationary("vecheat->adapt", + new AdaptInstationary("vecheat->adapt", vecheatSpace, adaptInfo, vecheat, @@ -327,17 +256,17 @@ int main(int argc, char** argv) // ===== create boundary functions ===== - G *boundaryFct = NEW G; + G *boundaryFct = new G; boundaryFct->setTimePtr(vecheat->getBoundaryTimePtr()); vecheat->setBoundaryFct(boundaryFct); - vecheatSpace->addDirichletBC(DIRICHLET, 0, boundaryFct); - vecheatSpace->addDirichletBC(DIRICHLET, 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()); + F *rhsFct0 = new F(vecheatSpace->getFESpace(0)->getBasisFcts()->getDegree()); rhsFct0->setTimePtr(vecheat->getRHSTimePtr()); - F *rhsFct1 = NEW F(vecheatSpace->getFESpace(1)->getBasisFcts()->getDegree()); + F *rhsFct1 = new F(vecheatSpace->getFESpace(1)->getBasisFcts()->getDegree()); rhsFct1->setTimePtr(vecheat->getRHSTimePtr()); // ===== create operators ===== @@ -345,7 +274,7 @@ int main(int argc, char** argv) double zero = 0.0; // create laplace - Operator *A00 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *A00 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(0), vecheatSpace->getFESpace(0)); @@ -353,7 +282,7 @@ int main(int argc, char** argv) A00->addSecondOrderTerm(new Laplace_SOT); A00->setUhOld(vecheat->getOldSolution()->getDOFVector(0)); - Operator *A11 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *A11 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(1), vecheatSpace->getFESpace(1)); @@ -361,7 +290,7 @@ int main(int argc, char** argv) A11->addSecondOrderTerm(new Laplace_SOT); A11->setUhOld(vecheat->getOldSolution()->getDOFVector(1)); - if((*(vecheat->getThetaPtr())) != 0.0) { + if ((*(vecheat->getThetaPtr())) != 0.0) { vecheatSpace->addMatrixOperator(A00, 0, 0, vecheat->getThetaPtr(), &one); @@ -370,22 +299,18 @@ int main(int argc, char** argv) &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 *C00 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(0), vecheatSpace->getFESpace(0)); - C00->addZeroOrderTerm(NEW Simple_ZOT); + C00->addZeroOrderTerm(new Simple_ZOT); C00->setUhOld(vecheat->getOldSolution()->getDOFVector(0)); @@ -397,12 +322,12 @@ int main(int argc, char** argv) vecheat->getTau1Ptr(), vecheat->getTau1Ptr()); - Operator *C11 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *C11 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(1), vecheatSpace->getFESpace(1)); - C11->addZeroOrderTerm(NEW Simple_ZOT); + C11->addZeroOrderTerm(new Simple_ZOT); C11->setUhOld(vecheat->getOldSolution()->getDOFVector(1)); @@ -415,17 +340,17 @@ int main(int argc, char** argv) vecheat->getTau1Ptr()); // create RHS operator - Operator *F0 = NEW Operator(Operator::VECTOR_OPERATOR, + Operator *F0 = new Operator(Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(0)); - F0->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct0)); + F0->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct0)); vecheatSpace->addVectorOperator(F0, 0); - Operator *F1 = NEW Operator(Operator::VECTOR_OPERATOR, + Operator *F1 = new Operator(Operator::VECTOR_OPERATOR, vecheatSpace->getFESpace(1)); - F1->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct1)); + F1->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct1)); vecheatSpace->addVectorOperator(F1, 1); diff --git a/demo/src/vecmultigrid.cc b/demo/src/vecmultigrid.cc index 44b6970b..b57b219b 100644 --- a/demo/src/vecmultigrid.cc +++ b/demo/src/vecmultigrid.cc @@ -8,36 +8,26 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { - return exp(-10.0*(x*x)); + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { + return exp(-10.0 * (x * x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}; + F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {} - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { int dim = x.getSize(); double r2 = (x*x); double ux = exp(-10.0*r2); @@ -64,24 +54,23 @@ int main(int argc, char* argv[]) vecmultigrid.initialize(INIT_ALL); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("vecmultigrid->adapt", + AdaptInfo *adaptInfo = new AdaptInfo("vecmultigrid->adapt", vecmultigrid.getNumComponents()); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("vecmultigrid->adapt", + AdaptStationary *adapt = new AdaptStationary("vecmultigrid->adapt", &vecmultigrid, adaptInfo); // ===== add boundary conditions ===== - vecmultigrid.addDirichletBC(1, 0, NEW G); - //vecmultigrid.addDirichletBC(1, 1, NEW G); + vecmultigrid.addDirichletBC(1, 0, 0, new G); // ===== create operators ===== Operator matrixOperator00(Operator::MATRIX_OPERATOR, vecmultigrid.getFESpace(0), vecmultigrid.getFESpace(0)); - matrixOperator00.addSecondOrderTerm(NEW Laplace_SOT); + matrixOperator00.addSecondOrderTerm(new Laplace_SOT); vecmultigrid.addMatrixOperator(&matrixOperator00, 0, 0); Operator rhsOperator0(Operator::VECTOR_OPERATOR, @@ -89,7 +78,7 @@ int main(int argc, char* argv[]) int degree = vecmultigrid.getFESpace(0)->getBasisFcts()->getDegree(); - rhsOperator0.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); + rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree))); vecmultigrid.addVectorOperator(&rhsOperator0, 0); @@ -101,11 +90,11 @@ int main(int argc, char* argv[]) vecmultigrid.getFESpace(1), vecmultigrid.getFESpace(1)); - matrixOperator10.addZeroOrderTerm(NEW Simple_ZOT); + matrixOperator10.addZeroOrderTerm(new Simple_ZOT); vecmultigrid.addMatrixOperator(&matrixOperator10, 1, 0); - matrixOperator11.addZeroOrderTerm(NEW Simple_ZOT(-1.0)); + matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0)); vecmultigrid.addMatrixOperator(&matrixOperator11, 1, 1); diff --git a/demo/src/vecnonlin.cc b/demo/src/vecnonlin.cc index 162c16ba..64d085c8 100644 --- a/demo/src/vecnonlin.cc +++ b/demo/src/vecnonlin.cc @@ -7,41 +7,32 @@ using namespace AMDiS; // ===== function definitions ================================================ // =========================================================================== -/** \brief - * Dirichlet boundary function - */ +/// Dirichlet boundary function class G : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(G); - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { return exp(-10.0*(x*x)); } }; -/** \brief - * RHS function - */ +/// RHS function class F : public AbstractFunction<double, WorldVector<double> > { public: - MEMORY_MANAGED(F); - - /** \brief - * Constructor - */ + /// Constructor F(int degree, double sigma_, double k_, int dim_) - : AbstractFunction<double, WorldVector<double> >(degree), sigma(sigma_), k(k_), dim(dim_) - {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const WorldVector<double>& x) const { + : AbstractFunction<double, WorldVector<double> >(degree), + sigma(sigma_), + k(k_), + dim(dim_) + {} + + /// Implementation of AbstractFunction::operator(). + double operator()(const WorldVector<double>& x) const + { double r2 = x*x; double ux = exp(-10.0*r2); double ux4 = ux*ux*ux*ux; @@ -54,25 +45,18 @@ private: int dim; }; -/** \brief - * Needed for zero order term. - */ +/// Needed for zero order term. class CQPFct : public AbstractFunction<double, double> { public: - MEMORY_MANAGED(CQPFct); + CQPFct() : AbstractFunction<double, double>(3) {} - CQPFct() : AbstractFunction<double, double>(3) {}; + /// Constructor. + CQPFct(double sigma_) : sigma(sigma_) {} - /** \brief - * Constructor. - */ - CQPFct(double sigma_) : sigma(sigma_) {}; - - /** \brief - * Implementation of AbstractFunction::operator(). - */ - double operator()(const double& x) const { + /// Implementation of AbstractFunction::operator(). + double operator()(const double& x) const + { return sigma * x * x * x; } @@ -85,34 +69,30 @@ private: // ===== class NonLin ======================================================== // =========================================================================== -/** \brief - * Non linear problem. - */ +/// Non linear problem. class NonLinVec : public ProblemNonLinVec { public: - /** \brief - * Constructor. - */ + /// Constructor. NonLinVec(const std::string& name_) : ProblemNonLinVec(name_), sigma(1.0) - {}; + {} - /** \brief - * Returns \ref sigma. - */ - double getSigma() { return sigma; }; + /// Returns \ref sigma. + double getSigma() + { + return sigma; + } - /** \brief - * Sets \ref sigma. - */ - void setSigma(double sig) { sigma = sig; }; + /// Sets \ref sigma. + void setSigma(double sig) + { + sigma = sig; + } private: - /** \brief - * Stefan-Bolzmann constant. - */ + /// Stefan-Bolzmann constant. double sigma; }; @@ -141,20 +121,20 @@ int main(int argc, char* argv[]) vecnonlin.setSigma(1.0); // === create adapt info === - AdaptInfo *adaptInfo = NEW AdaptInfo("vecnonlin->adapt", + AdaptInfo *adaptInfo = new AdaptInfo("vecnonlin->adapt", vecnonlin.getNumComponents()); // === create adapt === - AdaptStationary *adapt = NEW AdaptStationary("vecnonlin->adapt", + AdaptStationary *adapt = new AdaptStationary("vecnonlin->adapt", &vecnonlin, adaptInfo); // ===== create operators ===== - Operator *nonLinOperator0_0 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator0_0 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(0), vecnonlin.getFESpace(0)); - Operator *nonLinOperator0_1 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator0_1 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(1), vecnonlin.getFESpace(1)); @@ -163,16 +143,16 @@ int main(int argc, char* argv[]) nonLinOperator0_1->setUhOld(vecnonlin.getSolution()->getDOFVector(1)); nonLinOperator0_0->addZeroOrderTerm - (NEW VecAtQP_ZOT( + (new VecAtQP_ZOT( NULL, // use uhOld - NEW CQPFct(vecnonlin.getSigma()) + new CQPFct(vecnonlin.getSigma()) ) ); nonLinOperator0_1->addZeroOrderTerm - (NEW VecAtQP_ZOT( + (new VecAtQP_ZOT( NULL, // use uhOld - NEW CQPFct(vecnonlin.getSigma()) + new CQPFct(vecnonlin.getSigma()) ) ); @@ -181,11 +161,11 @@ int main(int argc, char* argv[]) vecnonlin.addVectorOperator(nonLinOperator0_0, 0); vecnonlin.addVectorOperator(nonLinOperator0_1, 1); - Operator *nonLinOperator2_0 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator2_0 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(0), vecnonlin.getFESpace(0)); - Operator *nonLinOperator2_1 = NEW Operator(Operator::MATRIX_OPERATOR | + Operator *nonLinOperator2_1 = new Operator(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(1), vecnonlin.getFESpace(1)); @@ -193,8 +173,8 @@ int main(int argc, char* argv[]) nonLinOperator2_0->setUhOld(vecnonlin.getSolution()->getDOFVector(0)); nonLinOperator2_1->setUhOld(vecnonlin.getSolution()->getDOFVector(1)); - nonLinOperator2_0->addSecondOrderTerm(NEW FactorLaplace_SOT(&k)); - nonLinOperator2_1->addSecondOrderTerm(NEW FactorLaplace_SOT(&k)); + nonLinOperator2_0->addSecondOrderTerm(new FactorLaplace_SOT(&k)); + nonLinOperator2_1->addSecondOrderTerm(new FactorLaplace_SOT(&k)); vecnonlin.addMatrixOperator(nonLinOperator2_0, 0, 0); vecnonlin.addMatrixOperator(nonLinOperator2_1, 1, 1); @@ -207,23 +187,23 @@ int main(int argc, char* argv[]) degree = vecnonlin.getFESpace(0)->getBasisFcts()->getDegree(); Operator* rhsFunctionOperator_0 = - NEW Operator(Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(0)); + new Operator(Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(0)); rhsFunctionOperator_0->addZeroOrderTerm - (NEW CoordsAtQP_ZOT(NEW F(degree, vecnonlin.getSigma(), k, dim))); + (new CoordsAtQP_ZOT(new F(degree, vecnonlin.getSigma(), k, dim))); degree = vecnonlin.getFESpace(1)->getBasisFcts()->getDegree(); Operator* rhsFunctionOperator_1 = - NEW Operator(Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(1)); + new Operator(Operator::VECTOR_OPERATOR, vecnonlin.getFESpace(1)); rhsFunctionOperator_1->addZeroOrderTerm - (NEW CoordsAtQP_ZOT(NEW F(degree, vecnonlin.getSigma(), k, dim))); + (new CoordsAtQP_ZOT(new F(degree, vecnonlin.getSigma(), k, dim))); vecnonlin.addVectorOperator(rhsFunctionOperator_0, 0, &factor2, &factor2); vecnonlin.addVectorOperator(rhsFunctionOperator_1, 1, &factor2, &factor2); // ===== add boundary conditions ===== - vecnonlin.addDirichletBC(DIRICHLET, 0, NEW G); - vecnonlin.addDirichletBC(DIRICHLET, 1, NEW G); + vecnonlin.addDirichletBC(DIRICHLET, 0, 0, new G); + vecnonlin.addDirichletBC(DIRICHLET, 1, 1, new G); // ===== start adaption loop ===== adapt->adapt(); -- GitLab