From 88e2a1793a11b7be34aae5f22074e8a94033e394 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 29 Oct 2013 09:55:48 +0000 Subject: [PATCH] singular dbc corrected --- AMDiS/CMakeLists.txt | 65 ++++++- demo/CMakeLists.txt | 6 + demo/init/ellipt.dat.3d | 6 +- demo/src/ellipt.cc | 10 +- extensions/ExtendedProblemStat.h | 16 +- extensions/ProblemStatMassConserve.h | 8 +- extensions/SingularDirichletBC2.cc | 19 +- extensions/SingularDirichletBC2.h | 41 +++- extensions/base_problems/CahnHilliard.cc | 122 ++++++------ extensions/base_problems/CahnHilliard.h | 67 ++++--- .../CahnHilliardNavierStokes_TwoPhase_RB.cc | 1 - .../CahnHilliardNavierStokes_TwoPhase_RB.h | 12 +- .../base_problems/NavierStokes_TaylorHood.h | 177 ++++++++++-------- extensions/demo/other/CMakeLists.txt | 13 +- .../other/src/cahnHilliard_navierStokes.cc | 20 +- .../src/fsi_explicit/ElasticityNavierStokes.h | 8 +- .../fsi_explicit/fluidStructureInteraction.cc | 2 +- extensions/demo/other/src/movingMesh.h | 8 +- 18 files changed, 367 insertions(+), 234 deletions(-) diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt index f72461f8..3e5d1ec8 100644 --- a/AMDiS/CMakeLists.txt +++ b/AMDiS/CMakeLists.txt @@ -497,7 +497,7 @@ if(ENABLE_EXTENSIONS) # ${EXTENSIONS_DIR}/base_problems/NavierStokes_Chorin.cc ${EXTENSIONS_DIR}/base_problems/NavierStokesCahnHilliard.cc ${EXTENSIONS_DIR}/base_problems/NavierStokesPhase_TaylorHood.cc - ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc +# ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood_RB.cc ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase.cc ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase_RB.cc @@ -708,6 +708,69 @@ list(APPEND deb_add_dirs "share/amdis/doc") list(REMOVE_DUPLICATES deb_add_dirs) +# TESTING + +ENABLE_TESTING() +INCLUDE(CTest) + +set(TEST_FILE amdis_test.cc) +file(WRITE ${TEST_FILE} "#include +int main (int argc, char* argv[]) { + AMDiS::init(argc, argv); + AMDiS::finalize(); + return 0 ; +}") + + +if(WIN32) + set(Boost_USE_STATIC_LIBS ON) + find_package(Boost 1.42 REQUIRED system iostreams filesystem program_options date_time zlib bzip2) +else() + set(Boost_USE_STATIC_LIBS OFF) + find_package(Boost 1.42 REQUIRED system iostreams filesystem program_options date_time) +endif() +if(Boost_FOUND) + list(APPEND AMDIS_LIBRARIES ${Boost_LIBRARIES}) + list(APPEND AMDIS_LIBRARY_DIRS ${Boost_LIBRARY_DIRS}) + list(APPEND AMDIS_INCLUDE_DIRS ${Boost_INCLUDE_DIR}) +else() + message(ERROR "Boost libraries not found") +endif(Boost_FOUND) +find_library(UMFPACK_LIB umfpack + HINTS ENV LIBRARY_PATH + DOC "The UMFPACK library") +if(UMFPACK_LIB) + list(APPEND AMDIS_LIBRARIES ${UMFPACK_LIB}) +else() + message(FATAL_ERROR "Could not find the UMFPACK library") +endif() + +add_executable("amdis_test" ${TEST_FILE}) +target_link_libraries("amdis_test" amdis reinit compositeFEM muparser ${AMDIS_LIBRARIES}) + +add_test(NAME test.amdis + WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} + COMMAND ${CMAKE_COMMAND} --build . --target amdis_test --config $ +) + +# add_test(NAME test.build.amdis +# WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} +# COMMAND ${CMAKE_COMMAND} --build . --target amdis --config $ +# ) +# add_test(NAME test.build.reinit +# WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} +# COMMAND ${CMAKE_COMMAND} --build . --target reinit --config $ +# ) +# add_test(NAME test.build.compositeFEM +# WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} +# COMMAND ${CMAKE_COMMAND} --build . --target compositeFEM --config $ +# ) +# add_test(NAME test.build.muparser +# WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} +# COMMAND ${CMAKE_COMMAND} --build . --target muparser --config $ +# ) + + # PACKAGES # ======================================================== set(CPACK_PACKAGE_NAME "AMDIS") diff --git a/demo/CMakeLists.txt b/demo/CMakeLists.txt index cd6db035..2c382945 100644 --- a/demo/CMakeLists.txt +++ b/demo/CMakeLists.txt @@ -71,6 +71,12 @@ add_executable("vecheat" ${vecheat}) target_link_libraries("vecheat" ${BASIS_LIBS}) install(TARGETS vecheat RUNTIME DESTINATION bin) + +set(test src/test.cc) +add_executable("test" ${test}) +target_link_libraries("test" ${BASIS_LIBS}) +install(TARGETS test RUNTIME DESTINATION bin) + # add_executable("matrix_test" src/matrix_test.cc) # target_link_libraries("matrix_test" ${BASIS_LIBS}) # install(TARGETS vecheat RUNTIME DESTINATION bin) diff --git a/demo/init/ellipt.dat.3d b/demo/init/ellipt.dat.3d index 289cabb3..5b4869b0 100644 --- a/demo/init/ellipt.dat.3d +++ b/demo/init/ellipt.dat.3d @@ -15,17 +15,17 @@ ellipt->solver->info: 2 ellipt->solver->left precon: diag ellipt->solver->right precon: no -ellipt->estimator[0]: residual +%ellipt->estimator[0]: residual ellipt->estimator[0]->error norm: H1_NORM % 1: H1_NORM, 2: L2_NORM ellipt->estimator[0]->C1: 0.01 % constant of jump residual -ellipt->marker[0]->strategy: 2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS +ellipt->marker[0]->strategy: 0 % 2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS ellipt->marker[0]->MSGamma: 0.5 ellipt->adapt[0]->tolerance: 1e-4 ellipt->adapt[0]->refine bisections: 2 -ellipt->adapt->max iteration: 10 +ellipt->adapt->max iteration: 0 % 10 ellipt->output->filename: output/ellipt.3d ellipt->output->ParaView format: 1 diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc index 213a91d0..8bba1235 100644 --- a/demo/src/ellipt.cc +++ b/demo/src/ellipt.cc @@ -62,17 +62,17 @@ int main(int argc, char* argv[]) // ===== create matrix operator ===== Operator matrixOperator(ellipt.getFeSpace()); - matrixOperator.addTerm(new Simple_SOT); + matrixOperator.addTerm(new Simple_ZOT); ellipt.addMatrixOperator(matrixOperator, 0, 0); // ===== create rhs operator ===== - Operator rhsOperator(ellipt.getFeSpace()); - rhsOperator.addTerm(new CoordsAtQP_ZOT(new F)); - ellipt.addVectorOperator(rhsOperator, 0); +// Operator rhsOperator(ellipt.getFeSpace()); +// rhsOperator.addTerm(new CoordsAtQP_ZOT(new F)); +// ellipt.addVectorOperator(rhsOperator, 0); // ===== add boundary conditions ===== - ellipt.addDirichletBC(1, 0, 0, new G); + ellipt.addDirichletBC(1, 0, 0, new AMDiS::Const >(1.0)); // ===== start adaption loop ===== diff --git a/extensions/ExtendedProblemStat.h b/extensions/ExtendedProblemStat.h index 6423a524..46fea460 100644 --- a/extensions/ExtendedProblemStat.h +++ b/extensions/ExtendedProblemStat.h @@ -68,11 +68,15 @@ public: ~ExtendedProblemStat() { - for (int i = 0; i < nComponents; ++i) + for (int i = 0; i < nComponents; ++i) { if (exactSolutions[i]) { delete exactSolutions[i]; exactSolutions[i] = NULL; } + } + + for (size_t i = 0; i < DirichletBcDataList.size(); i++) + delete DirichletBcDataList[i]; } @@ -80,7 +84,7 @@ public: void initialize(Flag initFlag, ProblemStatSeq *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) - { + { ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag); for (int i = 0; i < nComponents; ++i) exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution"); @@ -144,7 +148,7 @@ public: // apply dirichlet boundary conditions for (size_t k = 0; k < singularDirichletBC.size(); k++) applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector); - + // update solverMatrix if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) { solverMatrix.setMatrix(*getSystemMatrix()); @@ -230,7 +234,7 @@ public: { WARNING("Depreciated! Use addDirichletBC instead!!!\n"); DofIndex* dofIndex = new DofIndex(idx); - addDirichletBC(*dofIndex, row, col, values); + addDirichletBC(dofIndex, row, col, &values); } @@ -333,7 +337,9 @@ public: protected: // traverse matrix rows and set unity row where dirichlet condition shall be applied. - void applyDirichletBC(size_t row_, size_t col_, DegreeOfFreedom idx, double value, bool asmMatrix = true, bool asmVector = true) + void applyDirichletBC(size_t row_, size_t col_, + DegreeOfFreedom idx, double value, + bool asmMatrix = true, bool asmVector = true) { using namespace mtl; typedef DOFMatrix::base_matrix_type Matrix; diff --git a/extensions/ProblemStatMassConserve.h b/extensions/ProblemStatMassConserve.h index f02215c7..e1512f25 100644 --- a/extensions/ProblemStatMassConserve.h +++ b/extensions/ProblemStatMassConserve.h @@ -125,7 +125,7 @@ struct ProblemStatMassConserve : public ExtendedProblemStat prob2->setMeshes(meshes2); prob2->initialize(INIT_ALL - INIT_MESH); - for (size_t i = 1; i < meshes2.size(); i++) { + for (size_t i = 0; i < meshes2.size(); i++) { int globalRefinements = 0; // If AMDiS is compiled for parallel computations, the global refinements are @@ -139,12 +139,10 @@ struct ProblemStatMassConserve : public ExtendedProblemStat bool initMesh = true; // Initialize the meshes if there is no serialization file. - if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) + if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) { meshes2[i]->initialize(); - - // do global refinements - if (initMesh && meshes2[i]) prob2->getRefinementManager()->globalRefine(meshes2[i], globalRefinements); + } } diff --git a/extensions/SingularDirichletBC2.cc b/extensions/SingularDirichletBC2.cc index 0177a826..9c6a51c9 100644 --- a/extensions/SingularDirichletBC2.cc +++ b/extensions/SingularDirichletBC2.cc @@ -43,7 +43,7 @@ void getIndicesAux >( // AbstractFunction (signedDistFct) template<> -void getIndicesAux( +void getIndicesAux > >( AbstractFunction >& pos, BoundaryType boundary_nr, const FiniteElemSpace *feSpace, @@ -67,7 +67,7 @@ void getIndicesAux >( // Meshindicator template<> -void getIndicesAux( +void getIndicesAux > >( AbstractFunction >& pos, BoundaryType boundary_nr, const FiniteElemSpace *feSpace, @@ -85,7 +85,8 @@ void getValuesAux( const double &val, const FiniteElemSpace *feSpace, std::map &indicesValues, - boost::false_type is_not_abstract_fct) + boost::false_type is_not_abstract_fct, + int depth) { std::map::iterator it; for (it = indicesValues.begin(); it != indicesValues.end(); it++) @@ -97,13 +98,14 @@ void getValuesAux > >( const AbstractFunction >& val, const FiniteElemSpace *feSpace, std::map &indicesValues, - boost::true_type is_abstract_fct) + boost::true_type is_abstract_fct, + int depth) { - WorldVector x; + DOFVector > x(feSpace, "coords"); + feSpace->getMesh()->getDofIndexCoords(x); std::map::iterator it; for (it = indicesValues.begin(); it != indicesValues.end(); it++) - if (feSpace->getMesh()->getDofIndexCoords(it->first, feSpace, x)) - it->second = val(x); + it->second = val(x[it->first]); } // DOFVector @@ -112,7 +114,8 @@ void getValuesAux >( const DOFVector& val, const FiniteElemSpace *feSpace, std::map &indicesValues, - boost::false_type is_not_abstract_fct) + boost::false_type is_not_abstract_fct, + int depth) { std::map::iterator it; for (it = indicesValues.begin(); it != indicesValues.end(); it++) diff --git a/extensions/SingularDirichletBC2.h b/extensions/SingularDirichletBC2.h index 7321b2d2..26f88c6a 100644 --- a/extensions/SingularDirichletBC2.h +++ b/extensions/SingularDirichletBC2.h @@ -86,6 +86,23 @@ namespace details { const FiniteElemSpace* feSpace, std::map &indices, boost::true_type is_abstract_fct); + + // template specialization for AbstractFunctions + template<> + void getIndicesAux > >( + AbstractFunction > &pos, + BoundaryType b, + const FiniteElemSpace* + feSpace, std::map &indices, + boost::true_type is_abstract_fct); + template<> + void getIndicesAux > >( + AbstractFunction > &pos, + BoundaryType b, + const FiniteElemSpace* + feSpace, std::map &indices, + boost::true_type is_abstract_fct); + template void getIndicesAux(PosType &pos, BoundaryType b, @@ -114,20 +131,34 @@ namespace details { void getValuesAux(const ValueType &val, const FiniteElemSpace* feSpace, std::map &indices, - boost::true_type is_abstract_fct); + boost::true_type is_abstract_fct, + int depth = 0) + { + TEST_EXIT(depth == 0)("Recursion bricht nicht ab!\n"); + details::getValuesAux(dynamic_cast >&>(val), feSpace, indices, is_abstract_fct, depth+1); + } + // template specialization for AbstractFunctions + template<> + void getValuesAux > >( + const AbstractFunction >& val, + const FiniteElemSpace *feSpace, + std::map &indicesValues, + boost::true_type is_abstract_fct, + int depth); + template void getValuesAux(const ValueType &val, const FiniteElemSpace* feSpace, std::map &indices, - boost::false_type is_not_abstract_fct); + boost::false_type is_not_abstract_fct, + int depth = 0); template void getValues(const ValueType &val, const FiniteElemSpace* feSpace, std::map &indices) { - details::getValuesAux(val, feSpace, indices, - typename boost::is_base_of >, ValueType>::type()); + details::getValuesAux(val, feSpace, indices, typename boost::is_base_of >, ValueType>::type()); } /** @@ -300,7 +331,7 @@ struct DirichletBcData : DirichletBcDataBase { DirichletBcData(int i_, int j_, BoundaryType nr_, PosType &meshIndicator_, const ValueType &val_) : DirichletBcDataBase(i_, j_), boundary_nr(nr_), pos(meshIndicator_), val(val_) {} - virtual void addToList(const FiniteElemSpace *feSpace, std::vector &list) + void addToList(const FiniteElemSpace *feSpace, std::vector &list) override { std::map indicesValues; diff --git a/extensions/base_problems/CahnHilliard.cc b/extensions/base_problems/CahnHilliard.cc index 96ad1bf6..5f402708 100644 --- a/extensions/base_problems/CahnHilliard.cc +++ b/extensions/base_problems/CahnHilliard.cc @@ -24,7 +24,10 @@ using namespace AMDiS; -CahnHilliard::CahnHilliard(const std::string &name_) : +namespace details { + +template +CahnHilliard

::CahnHilliard(const std::string &name_) : super(name_), useMobility(false), useReinit(false), @@ -45,7 +48,7 @@ CahnHilliard::CahnHilliard(const std::string &name_) : // type of double well: 0= [0,1], 1= [-1,1] Parameters::get(name_ + "->double-well type", doubleWell); - Parameters::get(name + "->use reinit", useReinit); + Parameters::get(name_ + "->use reinit", useReinit); // transformation of the parameters minusEps = -eps; @@ -56,106 +59,106 @@ CahnHilliard::CahnHilliard(const std::string &name_) : } -void CahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo) -{ FUNCNAME("CahnHilliard::solveInitialProblem()"); - - Flag initFlag = initDataFromFile(adaptInfo); +template +void CahnHilliard

::solveInitialProblem(AdaptInfo *adaptInfo) +{ + Flag initFlag = self::initDataFromFile(adaptInfo); if (!initFlag.isSet(DATA_ADOPTED)) { int initialInterface = 0; - Initfile::get(name + "->initial interface", initialInterface); + Initfile::get(self::name + "->initial interface", initialInterface); double initialEps = eps; - Initfile::get(name + "->initial epsilon", initialEps); + Initfile::get(self::name + "->initial epsilon", initialEps); if (initialInterface == 0) { /// horizontale Linie double a= 0.0, dir= -1.0; - Initfile::get(name + "->line->pos", a); - Initfile::get(name + "->line->direction", dir); - prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir)); + Initfile::get(self::name + "->line->pos", a); + Initfile::get(self::name + "->line->direction", dir); + self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir)); } else if (initialInterface == 1) { /// schraege Linie double theta = m_pi/4.0; - prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0)); - transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min); + self::prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0)); + transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min); } else if (initialInterface == 2) { /// Ellipse double a= 1.0, b= 1.0; - Initfile::get(name + "->ellipse->a", a); - Initfile::get(name + "->ellipse->b", b); - prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b)); + Initfile::get(self::name + "->ellipse->a", a); + Initfile::get(self::name + "->ellipse->b", b); + self::prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b)); } else if (initialInterface == 3) { /// zwei horizontale Linien double a= 0.75, b= 0.375; - Initfile::get(name + "->lines->pos1", a); - Initfile::get(name + "->lines->pos2", b); - prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0)); - transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max); + Initfile::get(self::name + "->lines->pos1", a); + Initfile::get(self::name + "->lines->pos2", b); + self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0)); + transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max); } else if (initialInterface == 4) { /// Kreis double radius= 1.0; - Initfile::get(name + "->kreis->radius", radius); - prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius)); + Initfile::get(self::name + "->kreis->radius", radius); + self::prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius)); } else if (initialInterface == 5) { /// Rechteck double width = 0.5; double height = 0.3; WorldVector center; center.set(0.5); - Initfile::get(name + "->rectangle->width", width); - Initfile::get(name + "->rectangle->height", height); - Initfile::get(name + "->rectangle->center", center); - prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center)); + Initfile::get(self::name + "->rectangle->width", width); + Initfile::get(self::name + "->rectangle->height", height); + Initfile::get(self::name + "->rectangle->center", center); + self::prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center)); } // TODO: Redistancing einfügen! if (useReinit) { FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace( - const_cast(getMesh()->getVertexAdmin()), - Lagrange::getLagrange(getMesh()->getDim(), 1), - getMesh(), + const_cast(self::getMesh()->getVertexAdmin()), + Lagrange::getLagrange(self::getMesh()->getDim(), 1), + self::getMesh(), "P1"); DOFVector tmp(feSpace, "tmp"); - tmp.interpol(prob->getSolution()->getDOFVector(0)); + tmp.interpol(self::prob->getSolution()->getDOFVector(0)); - HL_SignedDistTraverse reinit("reinit", getMesh()->getDim()); + HL_SignedDistTraverse reinit("reinit", self::getMesh()->getDim()); reinit.calcSignedDistFct(adaptInfo, &tmp); #ifndef HAVE_PARALLEL_DOMAIN_AMDIS Recovery recovery(L2_NORM, 1); - recovery.recoveryUh(&tmp, *prob->getSolution()->getDOFVector(0)); + recovery.recoveryUh(&tmp, *self::prob->getSolution()->getDOFVector(0)); #else - prob->getSolution()->getDOFVector(0)->interpol(&tmp); + self::prob->getSolution()->getDOFVector(0)->interpol(&tmp); #endif } /// create phase-field from signed-dist-function if (doubleWell == 0) { - forEachDOF(prob->getSolution()->getDOFVector(0), + forEachDOF(self::prob->getSolution()->getDOFVector(0), new SignedDistToPhaseField(initialEps)); } else { - forEachDOF(prob->getSolution()->getDOFVector(0), + forEachDOF(self::prob->getSolution()->getDOFVector(0), new SignedDistToCh(initialEps)); } } } -void CahnHilliard::fillOperators() -{ FUNCNAME("CahnHilliard::fillOperators()"); - - const FiniteElemSpace* feSpace = prob->getFeSpace(); +template +void CahnHilliard

::fillOperators() +{ + const FiniteElemSpace* feSpace = self::prob->self::getFeSpace(); int degree = feSpace->getBasisFcts()->getDegree(); // c Operator *opChMnew = new Operator(feSpace, feSpace); opChMnew->addTerm(new Simple_ZOT); Operator *opChMold = new Operator(feSpace, feSpace); - opChMold->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0))); + opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0))); // -nabla*(grad(c)) Operator *opChL = new Operator(feSpace, feSpace); opChL->addTerm(new Simple_SOT); @@ -165,11 +168,11 @@ void CahnHilliard::fillOperators() if (useMobility) { if (doubleWell == 0) opChLM->addTerm(new VecAtQP_SOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), new MobilityCH0(gamma, degree))); else opChLM->addTerm(new VecAtQP_SOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), new MobilityCH1(gamma, degree))); } else opChLM->addTerm(new Simple_SOT(gamma)); @@ -177,22 +180,22 @@ void CahnHilliard::fillOperators() // -2*c_old^3 + 3/2*c_old^2 Operator *opChMPowExpl = new Operator(feSpace, feSpace); opChMPowExpl->addTerm(new VecAtQP_ZOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), new AMDiS::Pow<3>(-2.0, 3*degree))); if (doubleWell == 0) { opChMPowExpl->addTerm(new VecAtQP_ZOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), new AMDiS::Pow<2>(3.0/2.0, 2*degree))); } // -3*c_old^2 * c Operator *opChMPowImpl = new Operator(feSpace, feSpace); opChMPowImpl->addTerm(new VecAtQP_ZOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), new AMDiS::Pow<2>(-3.0, 2*degree))); if (doubleWell == 0) { opChMPowImpl->addTerm(new VecAtQP_ZOT( - prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(0), NULL, 3.0)); opChMPowImpl->addTerm(new Simple_ZOT(-0.5)); } else { @@ -201,29 +204,30 @@ void CahnHilliard::fillOperators() // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC] // ---------------------------------------------------------------------- - prob->addMatrixOperator(*opChMPowImpl,0,0); /// < -3*phi*c*c_old^2 , psi > - prob->addMatrixOperator(*opChL,0,0, &minusEpsSqr); /// < -eps^2*phi*grad(c) , grad(psi) > - prob->addMatrixOperator(*opChMnew,0,1); /// < phi*mu , psi > + self::prob->addMatrixOperator(*opChMPowImpl,0,0); /// < -3*phi*c*c_old^2 , psi > + self::prob->addMatrixOperator(*opChL,0,0, &minusEpsSqr); /// < -eps^2*phi*grad(c) , grad(psi) > + self::prob->addMatrixOperator(*opChMnew,0,1); /// < phi*mu , psi > // . . . vectorOperators . . . . . . . . . . . . . . . - prob->addVectorOperator(*opChMPowExpl,0); /// < -2*phi*c_old^3 , psi > + self::prob->addVectorOperator(*opChMPowExpl,0); /// < -2*phi*c_old^3 , psi > // dt(c) = laplace(mu) - u*grad(c) // ----------------------------------- - prob->addMatrixOperator(*opChMnew,1,0, getInvTau()); /// < phi*c/tau , psi > - prob->addMatrixOperator(*opChLM,1,1); /// < phi*grad(mu) , grad(psi) > + self::prob->addMatrixOperator(*opChMnew,1,0, self::getInvTau()); /// < phi*c/tau , psi > + self::prob->addMatrixOperator(*opChLM,1,1); /// < phi*grad(mu) , grad(psi) > // . . . vectorOperators . . . . . . . . . . . . . . . - prob->addVectorOperator(*opChMold,1, getInvTau()); /// < phi*c^old/tau , psi > + self::prob->addVectorOperator(*opChMold,1, self::getInvTau()); /// < phi*c^old/tau , psi > } -void CahnHilliard::finalizeData() -{ FUNCNAME("CahnHilliard::finalizeData()"); - - setAssembleMatrixOnlyOnce_butTimestepChange(0,1); - setAssembleMatrixOnlyOnce_butTimestepChange(1,0); +template +void CahnHilliard

::finalizeData() +{ + self::setAssembleMatrixOnlyOnce_butTimestepChange(0,1); + self::setAssembleMatrixOnlyOnce_butTimestepChange(1,0); if (!useMobility) - setAssembleMatrixOnlyOnce_butTimestepChange(1,1); + self::setAssembleMatrixOnlyOnce_butTimestepChange(1,1); } +} // end namespace details diff --git a/extensions/base_problems/CahnHilliard.h b/extensions/base_problems/CahnHilliard.h index bcacdf45..668a9a10 100644 --- a/extensions/base_problems/CahnHilliard.h +++ b/extensions/base_problems/CahnHilliard.h @@ -21,48 +21,55 @@ #include "AMDiS.h" #include "BaseProblem.h" #include "chns.h" - -using namespace AMDiS; -class CahnHilliard : public BaseProblem -{ -public: // definition of types +namespace detail { + + template + class CahnHilliard : public BaseProblem + { + public: // definition of types - typedef BaseProblem super; + typedef BaseProblem super; + typedef CahnHilliard self; -public: // public methods + public: // public methods - CahnHilliard(const std::string &name_); - ~CahnHilliard() {}; + CahnHilliard(const std::string &name_); + ~CahnHilliard() {}; - void solveInitialProblem(AdaptInfo *adaptInfo) override; + void solveInitialProblem(AMDiS::AdaptInfo *adaptInfo) override; - double getEpsilon() { return eps; } - int getDoubleWellType() { return doubleWell; } + double getEpsilon() { return eps; } + int getDoubleWellType() { return doubleWell; } - void fillOperators() override; - void fillBoundaryConditions() override {} - - void finalizeData() override; + void fillOperators() override; + void fillBoundaryConditions() override {} + + void finalizeData() override; -protected: // protected variables + protected: // protected variables - bool useMobility; - bool useReinit; + bool useMobility; + bool useReinit; - unsigned dim; + unsigned dim; - int doubleWell; + int doubleWell; - double gamma; - double eps; - double minusEps; - double epsInv; - double minusEpsInv; - double epsSqr; - double minusEpsSqr; -}; - + double gamma; + double eps; + double minusEps; + double epsInv; + double minusEpsInv; + double epsSqr; + double minusEpsSqr; + }; + +} // end namespace detail + +#include "CahnHilliard.hh" + +typedef ::detail::CahnHilliard CahnHilliard; #endif // CAHN_HILLIARD_H diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc index 697f009c..3be87481 100644 --- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc +++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc @@ -44,7 +44,6 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const c(0.0), sigma(1.0), surfaceTension(1.0), - reinit(NULL), velocity(NULL), fileWriter(NULL) { diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h index eb8ca733..c6b91db2 100644 --- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h +++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h @@ -55,10 +55,10 @@ public: // public methods - virtual void addMaterialDerivativeConst(int i, int j,double factor = 1.0); // rho = const - virtual void addMaterialDerivativeNonConst(int i, int j); // rho = rho(x) - virtual void addLaplaceTermConst(int i); // viscosity = const - virtual void addLaplaceTermNonConst(int i); // viscosity = viscosity(x) +// virtual void addMaterialDerivativeConst(int i, int j,double factor = 1.0); // rho = const +// virtual void addMaterialDerivativeNonConst(int i, int j); // rho = rho(x) +// virtual void addLaplaceTermConst(int i); // viscosity = const +// virtual void addLaplaceTermNonConst(int i); // viscosity = viscosity(x) protected: // protected variables @@ -75,13 +75,13 @@ protected: // protected variables } bool useMobility; -// bool useConservationForm; + bool useConservationForm; unsigned dim; int laplaceType; int nonLinTerm; -// int doubleWell; + int doubleWell; // navierStokes parameters double oldTimestep; diff --git a/extensions/base_problems/NavierStokes_TaylorHood.h b/extensions/base_problems/NavierStokes_TaylorHood.h index 3589c936..328585f7 100644 --- a/extensions/base_problems/NavierStokes_TaylorHood.h +++ b/extensions/base_problems/NavierStokes_TaylorHood.h @@ -19,99 +19,110 @@ #define NAVIER_STOKES_TAYLOR_HOOD_H #include "AMDiS.h" -#include "Helpers.h" -#include "POperators.h" #include "BaseProblem.h" #include "ExtendedProblemStat.h" -using namespace AMDiS; -/** \ingroup NavierStokes_TaylorHood -* \brief -* Navier-Stokes problem, using Taylor Hood elements -* and theta-scheme time-integration -*/ -class NavierStokes_TaylorHood : public BaseProblem +namespace detail { -public: // typedefs + using namespace AMDiS; - typedef NavierStokes_TaylorHood self; - typedef BaseProblem super; + /** \ingroup NavierStokes_TaylorHood + * \brief + * Navier-Stokes problem, using Taylor Hood elements + * and theta-scheme time-integration + */ -public: // methods - - NavierStokes_TaylorHood(const std::string &name_, bool createProblem = true); - - ~NavierStokes_TaylorHood(); - - void initData() override; - - void transferInitialSolution(AdaptInfo *adaptInfo) override; - - void closeTimestep(AdaptInfo *adaptInfo) override; - - void writeFiles(AdaptInfo *adaptInfo, bool force = false) override; - - // === getting/setting methods === - - DOFVector >* getVelocity() - { FUNCNAME_DBG("getVelocity()"); - TEST_EXIT_DBG(velocity != NULL) - ("velocity is NULL!"); - return velocity; - } - - DOFVector *getOldSolution(int i) - { FUNCNAME_DBG("getOldSolution()"); - TEST_EXIT_DBG(oldSolution[i] != NULL) - ("Index out of range, or oldSolution not provided for this component!\n"); - return oldSolution[i]; - } - - double getDensity() { return density; } - double getViscosity() { return viscosity; } - - void fillOperators() override; - - virtual void addLaplaceTerm(int i); - -protected: // variables - - DOFVector >* velocity; - - void calcVelocity() + template + class NavierStokes_TaylorHood : public BaseProblem { - if (dow == 1) - transformDOF(prob->getSolution()->getDOFVector(0), - velocity, new AMDiS::Vec1WorldVec); - else if (dow == 2) - transformDOF(prob->getSolution()->getDOFVector(0), - prob->getSolution()->getDOFVector(1), - velocity, new AMDiS::Vec2WorldVec); - else if (dow == 3) - transformDOF(prob->getSolution()->getDOFVector(0), - prob->getSolution()->getDOFVector(1), - prob->getSolution()->getDOFVector(2), - velocity, new AMDiS::Vec3WorldVec); - } - - int laplaceType; - int nonLinTerm; - - double oldTimestep; - double viscosity; - double density; - - double minus1; - double theta; - double theta1; - double minusTheta1; + public: // typedefs + + typedef NavierStokes_TaylorHood self; + typedef BaseProblem super; + + public: // methods + + NavierStokes_TaylorHood(const std::string &name_, bool createProblem = true); + + ~NavierStokes_TaylorHood(); + + void initData() override; + + void transferInitialSolution(AdaptInfo *adaptInfo) override; + + void closeTimestep(AdaptInfo *adaptInfo) override; + + void writeFiles(AdaptInfo *adaptInfo, bool force = false) override; + + // === getting/setting methods === + + DOFVector >* getVelocity() + { FUNCNAME_DBG("getVelocity()"); + TEST_EXIT_DBG(velocity != NULL) + ("velocity is NULL!"); + return velocity; + } + + DOFVector *getOldSolution(int i) + { FUNCNAME_DBG("getOldSolution()"); + TEST_EXIT_DBG(oldSolution[i] != NULL) + ("Index out of range, or oldSolution not provided for this component!\n"); + return oldSolution[i]; + } + + double getDensity() { return density; } + double getViscosity() { return viscosity; } + + void fillOperators() override; + + virtual void addLaplaceTerm(int i); + + protected: // variables + + DOFVector >* velocity; + + void calcVelocity() + { + if (self::dow == 1) + transformDOF(self::prob->getSolution()->getDOFVector(0), + velocity, new AMDiS::Vec1WorldVec); + else if (self::dow == 2) + transformDOF(self::prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(1), + velocity, new AMDiS::Vec2WorldVec); + else if (self::dow == 3) + transformDOF(self::prob->getSolution()->getDOFVector(0), + self::prob->getSolution()->getDOFVector(1), + self::prob->getSolution()->getDOFVector(2), + velocity, new AMDiS::Vec3WorldVec); + } + + int laplaceType; + int nonLinTerm; + + double oldTimestep; + double viscosity; + double density; + + double minus1; + double theta; + double theta1; + double minusTheta1; + + WorldVector force; + + FileVectorWriter *fileWriter; + std::vector*> oldSolution; + }; + +} // end namespace detail + +#include "NavierStokes_TaylorHood.hh" + +typedef ::detail::NavierStokes_TaylorHood NavierStokes_TaylorHood; - WorldVector force; - FileVectorWriter *fileWriter; - std::vector*> oldSolution; -}; #endif // NAVIER_STOKES_TAYLOR_HOOD_H diff --git a/extensions/demo/other/CMakeLists.txt b/extensions/demo/other/CMakeLists.txt index 65d28642..cdb89b62 100644 --- a/extensions/demo/other/CMakeLists.txt +++ b/extensions/demo/other/CMakeLists.txt @@ -10,22 +10,23 @@ if(AMDIS_FOUND) SET(BASIS_LIBS ${AMDIS_LIBRARIES}) endif(AMDIS_FOUND) +set(extensions_dir ../..) +set(base_problems_dir ${extensions_dir}/base_problems) + set(fsi src/fsi_explicit/fluidStructureInteraction.cc) add_executable("fsi" ${fsi}) target_link_libraries("fsi" ${BASIS_LIBS}) - set_target_properties("fsi" PROPERTIES COMPILE_FLAGS "-DUSE_KD_TREE") + set_target_properties("fsi" PROPERTIES COMPILE_FLAGS "-DUSE_KD_TREE") set(drivenCavity src/drivenCavity.cc) add_executable("drivenCavity" ${drivenCavity}) target_link_libraries("drivenCavity" ${BASIS_LIBS}) - include_directories(/opt/amdis/amdis-src/extensions/base_problems/) - set(drivenCavity_rb src/drivenCavity_rb.cc /opt/amdis/amdis-src/extensions/base_problems/CahnHilliardNavierStokes_RB.cc) + set(drivenCavity_rb src/drivenCavity_rb.cc ${base_problems_dir}/CahnHilliardNavierStokes_RB.cc) add_executable("drivenCavity_rb" ${drivenCavity_rb}) target_link_libraries("drivenCavity_rb" ${BASIS_LIBS}) -# include_directories(/opt/amdis/amdis-src/extensions/base_problems/) - set(drivenCavity_twophase src/drivenCavity_twophase_rb.cc /opt/amdis/amdis-src/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc) + set(drivenCavity_twophase src/drivenCavity_twophase_rb.cc ${base_problems_dir}/CahnHilliardNavierStokes_TwoPhase_RB.cc) add_executable("drivenCavity_twophase" ${drivenCavity_twophase}) target_link_libraries("drivenCavity_twophase" ${BASIS_LIBS}) @@ -56,7 +57,7 @@ endif(AMDIS_FOUND) # add_executable("navierStokesDd3" ${navierStokesDd3}) # target_link_libraries("navierStokesDd3" ${BASIS_LIBS}) - set(ch src/cahnHilliard.cc) + set(ch src/cahnHilliard.cc ${extensions_dir}/preconditioner/CahnHilliard_.cc) add_executable("ch" ${ch}) target_link_libraries("ch" ${BASIS_LIBS}) diff --git a/extensions/demo/other/src/cahnHilliard_navierStokes.cc b/extensions/demo/other/src/cahnHilliard_navierStokes.cc index 91089306..7d5da7c2 100644 --- a/extensions/demo/other/src/cahnHilliard_navierStokes.cc +++ b/extensions/demo/other/src/cahnHilliard_navierStokes.cc @@ -22,24 +22,28 @@ #include "CouplingBaseProblem2.h" #endif #include "CahnHilliard.h" +#include "ProblemStatMassConserve.h" #include "NavierStokes_TaylorHood.h" #include "Refinement.h" #include "MeshFunction_Level.h" using namespace AMDiS; -class CahnHilliardNavierStokes : public CouplingBaseProblem +typedef ::details::CahnHilliard CHType; +typedef NavierStokes_TaylorHood NSType; + +class CahnHilliardNavierStokes : public CouplingBaseProblem< ProblemStat, CHType, NSType > { public: // typedefs - typedef CouplingBaseProblem super; + typedef CouplingBaseProblem< ProblemStat, CHType, NSType > super; public: // methods /// Constructor CahnHilliardNavierStokes(std::string name, - CahnHilliard &chProb_, - NavierStokes_TaylorHood &nsProb_) + CHType &chProb_, + NSType &nsProb_) : super(name, chProb_, nsProb_), chProb(&chProb_), nsProb(&nsProb_), @@ -113,8 +117,8 @@ public: // methods protected: // variables - CahnHilliard *chProb; - NavierStokes_TaylorHood *nsProb; + CHType *chProb; + NSType *nsProb; PhaseFieldRefinement* refFunction; RefinementLevelDOF *refinement; }; @@ -126,8 +130,8 @@ int main(int argc, char** argv) AMDiS::init(argc, argv); Timer t; - CahnHilliard chProb("ch"); - NavierStokes_TaylorHood nsProb("ns"); + CHType chProb("ch"); + NSType nsProb("ns"); CahnHilliardNavierStokes chns("chns", chProb, nsProb); chns.initialize(INIT_ALL); diff --git a/extensions/demo/other/src/fsi_explicit/ElasticityNavierStokes.h b/extensions/demo/other/src/fsi_explicit/ElasticityNavierStokes.h index 143bb632..fbc1ae9a 100644 --- a/extensions/demo/other/src/fsi_explicit/ElasticityNavierStokes.h +++ b/extensions/demo/other/src/fsi_explicit/ElasticityNavierStokes.h @@ -215,12 +215,12 @@ public: CouplingTimeInterface::solveInitialProblem(adaptInfo); DOFVector > coords1(parametricCoords1[0]->getFeSpace(), "coords"); - parametricCoords1[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords1[0]->getFeSpace(), coords1); + parametricCoords1[0]->getFeSpace()->getMesh()->getDofIndexCoords(coords1); DOFVector > coords2(parametricCoords2[0]->getFeSpace(), "coords"); - parametricCoords2[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords2[0]->getFeSpace(), coords2); + parametricCoords2[0]->getFeSpace()->getMesh()->getDofIndexCoords(coords2); for (int i = 0; i < Global::getGeo(WORLD); i++) { - transformDOF(&coords1, parametricCoords1[i], new Component(i)); - transformDOF(&coords2, parametricCoords2[i], new Component(i)); + transformDOF(&coords1, parametricCoords1[i], new AMDiS::Component2<>(i)); + transformDOF(&coords2, parametricCoords2[i], new AMDiS::Component2<>(i)); } // ===== enable parametric traverse ===== diff --git a/extensions/demo/other/src/fsi_explicit/fluidStructureInteraction.cc b/extensions/demo/other/src/fsi_explicit/fluidStructureInteraction.cc index ae7bc682..88dbc7ca 100644 --- a/extensions/demo/other/src/fsi_explicit/fluidStructureInteraction.cc +++ b/extensions/demo/other/src/fsi_explicit/fluidStructureInteraction.cc @@ -147,7 +147,7 @@ public: super::initTimestep(adaptInfo); TEST_EXIT(boundaryStress != NULL)("NULL-Pointer problem!\n"); for (size_t i = 0; i < dow; i++) - transformDOF_coords(*boundaryStress, *stress[i], new Component(i)); + transformDOF_coords(*boundaryStress, *stress[i], new AMDiS::Component2<>(i)); fileWriterStress->writeFiles(adaptInfo, false); } diff --git a/extensions/demo/other/src/movingMesh.h b/extensions/demo/other/src/movingMesh.h index 6d7d1ef7..abacd9d3 100644 --- a/extensions/demo/other/src/movingMesh.h +++ b/extensions/demo/other/src/movingMesh.h @@ -282,12 +282,12 @@ public: CouplingTimeInterface::solveInitialProblem(adaptInfo); DOFVector > coords1(parametricCoords1[0]->getFeSpace(), "coords"); - parametricCoords1[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords1[0]->getFeSpace(), coords1); + parametricCoords1[0]->getFeSpace()->getMesh()->getDofIndexCoords(coords1); DOFVector > coords2(parametricCoords2[0]->getFeSpace(), "coords"); - parametricCoords2[0]->getFeSpace()->getMesh()->getDofIndexCoords(parametricCoords2[0]->getFeSpace(), coords2); + parametricCoords2[0]->getFeSpace()->getMesh()->getDofIndexCoords(coords2); for (int i = 0; i < Global::getGeo(WORLD); i++) { - transformDOF(&coords1, parametricCoords1[i], new Component(i)); - transformDOF(&coords2, parametricCoords2[i], new Component(i)); + transformDOF(&coords1, parametricCoords1[i], new AMDiS::Component2<>(i)); + transformDOF(&coords2, parametricCoords2[i], new AMDiS::Component2<>(i)); } // ===== enable parametric traverse ===== -- GitLab