diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt index 2342aaa395f498d99942f797665765d6da64466d..adec7b9dced99cf30237f5bd0d9d6987d45b4a7c 100644 --- a/AMDiS/CMakeLists.txt +++ b/AMDiS/CMakeLists.txt @@ -16,12 +16,35 @@ endif () # ------------------------------------------------------------------------------ # some options to control the build process -option(ENABLE_PARALLEL_DOMAIN "Use parallel domain decomposition" false) +option(ENABLE_MPI "Compile with MPI compiler and library" false) option(ENABLE_UMFPACK "Enable support for the UMFPACK solver" true) option(ENABLE_EXTENSIONS "Use extensions for AMDiS" false) option(ENABLE_HYPRE "Use HYPRE AMG solver" false) option(ENABLE_PNG "Use png reader/writer" false) -option(ENABLE_SEQ_PETSC "Use sequential PETSc solvers" false) + +include(CMakeDependentOption) +cmake_dependent_option(ENABLE_PARALLEL_DOMAIN + "Use parallel domain decomposition" true + "ENABLE_MPI" false) + +# additional options for parallel_domain +# { +include(CMakeDependentOption) +cmake_dependent_option(ENABLE_ZOLTAN + "Add support for the Parallel Partitioning suite Zoltan" false + "ENABLE_PARALLEL_DOMAIN" false) + +include(CMakeDependentOption) +cmake_dependent_option(ENABLE_PARALLEL_SOLVERS + "Add some problem dependent solver, e.g. Feti, Navier-Stokes and Cahn-Hilliard" true + "ENABLE_PARALLEL_DOMAIN" false) + +include(CMakeDependentOption) +cmake_dependent_option(ENABLE_BDDCML + "Add support for the BDDCML library" false + "ENABLE_PARALLEL_DOMAIN" false) +# } + # enable/disable some more features option(ENABLE_COMPRESSION "Use output compression for vtu and arh files" true) @@ -33,6 +56,12 @@ mark_as_advanced(ENABLE_REINIT) option(ENABLE_COMPOSITE_FEM "Compile compositeFEM library" true) mark_as_advanced(ENABLE_COMPOSITE_FEM) +option(ENABLE_OPENMP "Compile with OpenMP flags. (Experimental)" false) +mark_as_advanced(ENABLE_OPENMP) + +option(ENABLE_SEQ_PETSC "Use sequential PETSc solvers" false) +mark_as_advanced(ENABLE_SEQ_PETSC) + option(BUILD_SHARED_LIBS "Build all libraries as shared or static, default: shared" true) mark_as_advanced(BUILD_SHARED_LIBS) diff --git a/AMDiS/cmake/amdis_parallel.cmake b/AMDiS/cmake/amdis_parallel.cmake index c5ea733b2d7d8ead3c2eb21941399c3f9f9d7df6..0d325c3bcbfa9317fc5b32b210d1e48f69d683be 100644 --- a/AMDiS/cmake/amdis_parallel.cmake +++ b/AMDiS/cmake/amdis_parallel.cmake @@ -1,10 +1,4 @@ if (ENABLE_PARALLEL_DOMAIN) - option(ENABLE_ZOLTAN "Add support for the Parallel Partitioning suite Zoltan" false) - option(ENABLE_PARALLEL_SOLVERS "Add some problem dependent solver, e.g. Feti, Navier-Stokes and Cahn-Hilliard" true) - option(ENABLE_BDDCML "Add support for the BDDCML library" false) - mark_as_advanced(ENABLE_PARALLEL_SOLVERS) - mark_as_advanced(ENABLE_BDDCML) - set(PARALLEL_DOMAIN_AMDIS_SRC ${SOURCE_DIR}/parallel/DofComm.cc ${SOURCE_DIR}/parallel/CheckerPartitioner.cc diff --git a/AMDiS/cmake/enable_mtl4.cmake b/AMDiS/cmake/enable_mtl4.cmake index ae2b3e3d97e11df2ca814c7080c333475ef5e9a2..3b24ae48165888d5b6d1585a161f98069b2c7e4e 100644 --- a/AMDiS/cmake/enable_mtl4.cmake +++ b/AMDiS/cmake/enable_mtl4.cmake @@ -12,22 +12,14 @@ macro(enable_mtl4 _FLAGS_ _INC_DIRS_ _LIBS_) endif (LINK_EXECUTABLE) list(APPEND ${_FLAGS_} "-DMTL_ASSERT_FOR_THROW=1") - - if (ENABLE_CXX11 OR AMDIS_NEED_CXX11) - set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") - foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) - list(APPEND ${_FLAGS_} "-DMTL_WITH_${feature}") - endforeach () - endif () - + set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") + foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) + list(APPEND ${_FLAGS_} "-DMTL_WITH_${feature}") + endforeach () if (ENABLE_OPENMP) find_package(OpenMP REQUIRED) - if (OPENMP_FOUND) - list(APPEND ${_FLAGS_} "-DMTL_WITH_OPENMP" ${OpenMP_CXX_FLAGS}) - else () - message(FATAL_ERROR "OpenMP not found") - endif (OPENMP_FOUND) + list(APPEND ${_FLAGS_} "-DMTL_WITH_OPENMP" "-DHAVE_OPENMP" ${OpenMP_CXX_FLAGS}) endif (ENABLE_OPENMP) endmacro(enable_mtl4) diff --git a/AMDiS/cmake3/AMDIS.cmake.in b/AMDiS/cmake3/AMDIS.cmake.in index 0d3dfa26c24392590c726d57d011e5a76eda8c60..b98c63c2111eb848cdea932362d79726c83f9ad5 100644 --- a/AMDiS/cmake3/AMDIS.cmake.in +++ b/AMDiS/cmake3/AMDIS.cmake.in @@ -1,4 +1,3 @@ -set(AMDIS_NEED_CXX11 @ENABLE_CXX11@) set(AMDIS_NEED_UMFPACK @ENABLE_UMFPACK@) set(AMDIS_NEED_COMPRESSION @ENABLE_COMPRESSION@) set(AMDIS_NEED_EXTENSIONS @ENABLE_EXTENSIONS@) @@ -12,22 +11,23 @@ set(AMDIS_NEED_PNG @ENABLE_PNG@) set(AMDIS_NEED_SEQ_PETSC @ENABLE_SEQ_PETSC@) set(BUILD_SHARED_LIBS @BUILD_SHARED_LIBS@) -set(CMAKE_BUILD_TYPE @CMAKE_BUILD_TYPE@) add_library(amdis_base INTERFACE) add_library(AMDiS ALIAS amdis_base) target_compile_definitions(amdis_base INTERFACE) -if (AMDIS_NEED_CXX11) - target_enable_cxx11(AMDIS_NEED_CXX11 amdis_base INTERFACE) - if (NOT AMDIS_NEED_CXX11) - message(FATAL_ERROR "AMDiS was compiled with c++11 support, but the current compiler does not support this feature!") - endif (NOT AMDIS_NEED_CXX11) - target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=1) -else () - target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=0) -endif (AMDIS_NEED_CXX11) + +target_enable_cxx14(SUPPORTS_CXX14 amdis_base INTERFACE) +if (NOT SUPPORTS_CXX14) + target_enable_cxx11(SUPPORTS_CXX11 amdis_base INTERFACE) +endif () + +if (NOT SUPPORTS_CXX11 AND NOT SUPPORTS_CXX14) + message(FATAL_ERROR "AMDiS was compiled with c++11 support, but the current compiler does not support this feature!") +endif () +target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=1) + # try to detect the AMDiS include directory diff --git a/AMDiS/cmake3/CMakeLists.txt b/AMDiS/cmake3/CMakeLists.txt index fb308342960b552f210091772440656a8429d487..b9c5002f4f6737de045e1cd0d91802383813a8a8 100644 --- a/AMDiS/cmake3/CMakeLists.txt +++ b/AMDiS/cmake3/CMakeLists.txt @@ -23,12 +23,19 @@ add_library(AMDiS::base ALIAS amdis_base) target_include_directories(amdis_base INTERFACE ${SOURCE_DIR}) target_compile_definitions(amdis_base INTERFACE) -target_enable_cxx11(ENABLE_CXX11 amdis_base INTERFACE) -if (ENABLE_CXX11) - target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=1) + +target_enable_cxx14(SUPPORTS_CXX14 amdis_base INTERFACE) +if (SUPPORTS_CXX14) + target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX14=1) else () - target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=0) -endif (ENABLE_CXX11) + target_enable_cxx11(SUPPORTS_CXX11 amdis_base INTERFACE) +endif () + +if (NOT SUPPORTS_CXX11 AND NOT SUPPORTS_CXX14) + message(FATAL_ERROR "AMDiS needs c++11 support, but the current compiler does not support this feature!") +endif () +target_compile_definitions(amdis_base INTERFACE AMDIS_HAS_CXX11=1) + target_enable_mtl4(amdis_base INTERFACE) target_enable_boost(amdis_base INTERFACE) @@ -65,6 +72,7 @@ add_library(amdis ${SOURCE_DIR}/DOFAdmin.cc ${SOURCE_DIR}/DOFIndexed.cc ${SOURCE_DIR}/DOFMatrix.cc + ${SOURCE_DIR}/DOFSerializer.cc ${SOURCE_DIR}/DOFVector.cc ${SOURCE_DIR}/DirichletBC.cc ${SOURCE_DIR}/DualTraverse.cc @@ -171,6 +179,11 @@ add_library(amdis ${SOURCE_DIR}/ProblemStatDbg.cc ) +if (ENABLE_MPI) + target_enable_mpi(amdis_base INTERFACE) + target_compile_definitions(amdis_base INTERFACE HAVE_MPI=1) +endif (ENABLE_MPI) + include(amdis_parallel) # adds sources to amdis include(amdis_extensions) # -> target AMDiS::extensions include(muparser) # -> target muparser diff --git a/AMDiS/cmake3/amdis_parallel.cmake b/AMDiS/cmake3/amdis_parallel.cmake index ab09f39cb3177bdb8b78cb2a708d5d90bdb18693..3a9ebc8b31063956110a6ecf680213b9558689bb 100644 --- a/AMDiS/cmake3/amdis_parallel.cmake +++ b/AMDiS/cmake3/amdis_parallel.cmake @@ -1,8 +1,4 @@ if (ENABLE_PARALLEL_DOMAIN) - option(ENABLE_ZOLTAN "Add support for the Parallel Partitioning suite Zoltan" false) - option(ENABLE_PARALLEL_SOLVERS "Add some problem dependent solver, e.g. Feti, Navier-Stokes and Cahn-Hilliard" true) - mark_as_advanced(ENABLE_PARALLEL_SOLVERS) - add_library(amdis_parallel INTERFACE) target_sources(amdis PRIVATE ${SOURCE_DIR}/parallel/DofComm.cc diff --git a/AMDiS/cmake3/muparser.cmake b/AMDiS/cmake3/muparser.cmake index 016a16daa29203f7d13c8c04b4f736273dc8fca8..82910f8ce8c4d7d7049b88f4f3ea5acc0efb0773 100644 --- a/AMDiS/cmake3/muparser.cmake +++ b/AMDiS/cmake3/muparser.cmake @@ -17,10 +17,8 @@ add_library(muparser STATIC target_include_directories(muparser PUBLIC ${MUPARSER_INCLUDE_DIR}) set_property(TARGET muparser PROPERTY POSITION_INDEPENDENT_CODE ON) -if (ENABLE_CXX11) - target_enable_cxx11(MUPARSER_ENABLE_CXX11 muparser PRIVATE) - target_compile_definitions(muparser PRIVATE AMDIS_HAS_CXX11=1) -endif (ENABLE_CXX11) +target_enable_cxx11(MUPARSER_ENABLE_CXX11 muparser PRIVATE) +target_compile_definitions(muparser PRIVATE AMDIS_HAS_CXX11=1) # specify how to install this target: # ----------------------------------- diff --git a/AMDiS/cmake3/target_enable_cxx11.cmake b/AMDiS/cmake3/target_enable_cxx11.cmake index 687fb51ecfda7a590d36349271730271d04644fd..13c5819bcd39556bb09c8e922933374ee6758f1c 100644 --- a/AMDiS/cmake3/target_enable_cxx11.cmake +++ b/AMDiS/cmake3/target_enable_cxx11.cmake @@ -1,6 +1,37 @@ include(CheckCXXCompilerFlag) include(CheckCXXSourceCompiles) +macro(target_enable_cxx14 RESULT_VAR _TARGET_ _SCOPE_) + check_cxx_compiler_flag("-std=c++14" COMPILER_SUPPORTS_CXX14_FLAG) + + set(CXX14_CODE " + auto f() { return 1; } + int main(){ + auto f1 = [](auto const& y) { return y; }; + auto y0 = f(); + auto y1 = f1(1); + }") + + if (COMPILER_SUPPORTS_CXX14_FLAG) + set(CMAKE_REQUIRED_FLAGS "-std=c++14") + check_cxx_source_compiles("${CXX14_CODE}" CXX14_COMPILES_WITH_CXX14_FLAG) + set(CMAKE_REQUIRED_FLAGS "") + endif () + + if (COMPILER_SUPPORTS_CXX14_FLAG AND CXX14_COMPILES_WITH_CXX14_FLAG) + target_compile_options(${_TARGET_} ${_SCOPE_} "-std=c++14") + set(${RESULT_VAR} true CACHE BOOL "Enable C++14 compiler features" FORCE) + else () + check_cxx_source_compiles("${CXX14_CODE}" CXX14_COMPILES) + + if (CXX14_COMPILES) + set(${RESULT_VAR} true CACHE BOOL "Enable C++14 compiler features" FORCE) + else () + set(${RESULT_VAR} false CACHE BOOL "Enable C++14 compiler features" FORCE) + endif () + endif () +endmacro(target_enable_cxx14) + macro(target_enable_cxx11 RESULT_VAR _TARGET_ _SCOPE_) check_cxx_compiler_flag("-std=c++11" COMPILER_SUPPORTS_CXX11_FLAG) diff --git a/AMDiS/cmake3/target_enable_mtl4.cmake b/AMDiS/cmake3/target_enable_mtl4.cmake index b718b73b8f95f0684a2e057ccfaae973d12ae6f7..517eae51e412598a03097ef543c1d3b4ea201786 100644 --- a/AMDiS/cmake3/target_enable_mtl4.cmake +++ b/AMDiS/cmake3/target_enable_mtl4.cmake @@ -1,32 +1,25 @@ macro(target_enable_mtl4 _TARGET_ _SCOPE_) - if (${ARGC} GREATER 2) - set(LINK_EXECUTABLE ON) - else () - set(LINK_EXECUTABLE OFF) - endif () + if (${ARGC} GREATER 2) + set(LINK_EXECUTABLE ON) + else () + set(LINK_EXECUTABLE OFF) + endif () - if (LINK_EXECUTABLE) - target_include_directories(${_TARGET_} ${_SCOPE_} ${AMDIS_INCLUDE_DIR}/mtl4) - else () - target_include_directories(${_TARGET_} ${_SCOPE_} ${BASE_DIR}/lib/mtl4) - endif (LINK_EXECUTABLE) - target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_ASSERT_FOR_THROW=1) + if (LINK_EXECUTABLE) + target_include_directories(${_TARGET_} ${_SCOPE_} ${AMDIS_INCLUDE_DIR}/mtl4) + else () + target_include_directories(${_TARGET_} ${_SCOPE_} ${BASE_DIR}/lib/mtl4) + endif (LINK_EXECUTABLE) + target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_ASSERT_FOR_THROW=1) - if (ENABLE_CXX11) - set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") - foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) - target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_WITH_${feature}) - endforeach () - endif (ENABLE_CXX11) - - if (ENABLE_OPENMP) - find_package(OpenMP REQUIRED) - if (OPENMP_FOUND) - target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_WITH_OPENMP) - target_compile_options(${_TARGET_} ${_SCOPE_} ${OpenMP_CXX_FLAGS}) - else () - message(FATAL_ERROR "OpenMP not found") - endif (OPENMP_FOUND) - endif (ENABLE_OPENMP) + set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") + foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) + target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_WITH_${feature}) + endforeach () + if (ENABLE_OPENMP) + find_package(OpenMP REQUIRED) + target_compile_definitions(${_TARGET_} ${_SCOPE_} MTL_WITH_OPENMP HAVE_OPENMP) + target_compile_options(${_TARGET_} ${_SCOPE_} ${OpenMP_CXX_FLAGS}) + endif (ENABLE_OPENMP) endmacro(target_enable_mtl4) diff --git a/AMDiS/src/AMDiS.cc b/AMDiS/src/AMDiS.cc index 897250554df76da57dc5a30df5353445f7d07301..b81d6b1cf020c81d0746f5f80b6743a0c3d1bb3b 100644 --- a/AMDiS/src/AMDiS.cc +++ b/AMDiS/src/AMDiS.cc @@ -115,7 +115,7 @@ namespace AMDiS { Parameters::init(initFileName); } -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_OPENMP) || defined(HAVE_MPI) Parameters::get("parallel->log main rank", Msg::outputMainRank); #endif diff --git a/AMDiS/src/AdaptInfo.h b/AMDiS/src/AdaptInfo.h index b13206d9db86554cbe11f3d655898efb6e02fc06..b7434de7b75c2f7a4fdf88bc86ccefb64f9c82e4 100644 --- a/AMDiS/src/AdaptInfo.h +++ b/AMDiS/src/AdaptInfo.h @@ -156,11 +156,8 @@ namespace AMDiS { rosenbrockMode(false) { init(); - char number[5]; - for (int i = 0; i < size; i++) { - sprintf(number, "[%d]", i); - scalContents[i] = new ScalContent(name + std::string(number)); - } + for (int i = 0; i < size; i++) + scalContents[i] = new ScalContent(name + "[" + std::to_string(i) + "]"); } /// Destructor. diff --git a/AMDiS/src/CouplingIterationInterface.cc b/AMDiS/src/CouplingIterationInterface.cc index 1cda2797181d9b77579fdc4fb8e3bd26092a94dc..7a3b2d92eae0158a600642f9bfaa2e75fae1d1ce 100644 --- a/AMDiS/src/CouplingIterationInterface.cc +++ b/AMDiS/src/CouplingIterationInterface.cc @@ -88,7 +88,7 @@ namespace AMDiS { /// Returns number of managed problems - int CouplingIterationInterface::getNumProblems() + int CouplingIterationInterface::getNumProblems() const { size_t num= 0; for (size_t i = 0; i < problems.size(); ++i) @@ -122,7 +122,7 @@ namespace AMDiS { /// Returns the name of the problem. - std::string CouplingIterationInterface::getName(size_t number) + std::string CouplingIterationInterface::getName(size_t number) const { if (static_cast(problems.size()) <= number) throw(std::runtime_error("Problem number out of range.")); diff --git a/AMDiS/src/CouplingIterationInterface.h b/AMDiS/src/CouplingIterationInterface.h index e141f4ef774fb40a60f414ef376d016b17e517dc..c926cd81f205118ac2ada4ca5ba953c08887a6c2 100644 --- a/AMDiS/src/CouplingIterationInterface.h +++ b/AMDiS/src/CouplingIterationInterface.h @@ -61,7 +61,7 @@ namespace AMDiS { /// Called after each adaption loop iteration. virtual void endIteration(AdaptInfo *adaptInfo); - virtual int getNumProblems(); + virtual int getNumProblems() const; /// Returns number of managed problems virtual size_t getNumIterationInterfaces() { return problems.size(); } @@ -75,8 +75,8 @@ namespace AMDiS { virtual ProblemIterationInterface *getIterationInterface(size_t number = 0) { return problems[number]; } /// Returns the name of the problem with the given number. - virtual std::string getName(size_t number); - virtual std::string getName() { return getName(0); } + virtual std::string getName(size_t number) const; + virtual std::string getName() const { return getName(0); } virtual void setSolveProblem(size_t number, bool flag = true) { solveProblem[number] = flag; } virtual void setSolveProblem(std::string name, bool flag = true); diff --git a/AMDiS/src/CouplingProblemStat.h b/AMDiS/src/CouplingProblemStat.h index a7d5cabc4da3a78a69fdadeee92e2a0a754a5875..c51dfd6c96efc361129e1fb10dd4a26bd46e0a1b 100644 --- a/AMDiS/src/CouplingProblemStat.h +++ b/AMDiS/src/CouplingProblemStat.h @@ -75,8 +75,7 @@ namespace AMDiS { virtual void initialize(Flag initFlag, ProblemStatSeq *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) - { FUNCNAME("CouplingProblemStat::initialize()"); - + { super::initialize(initFlag - INIT_MESH); const Flag DEFAULT_INIT = (INIT_FE_SPACE | INIT_MESH | CREATE_MESH | INIT_SYSTEM | INIT_SOLVER | INIT_ESTIMATOR | INIT_MARKER | INIT_FILEWRITER); @@ -130,7 +129,7 @@ namespace AMDiS { problems[i]->componentMeshes.resize(nComponents + nAddComponents); - for (size_t j = 0; j < nComponents + nAddComponents; j++) { + for (int j = 0; j < nComponents + nAddComponents; j++) { // name of the mesh std::string meshName(""); Parameters::get(problems[i]->getName() + "->mesh", meshName); @@ -186,7 +185,7 @@ namespace AMDiS { problems[p]->componentSpaces.resize(nComponents + nAddComponents, NULL); problems[p]->traverseInfo.resize(nComponents); - for (size_t i = 0; i < nComponents + nAddComponents; i++) { + for (int i = 0; i < nComponents + nAddComponents; i++) { std::string componentString = "[" + to_string(i) + "]"; @@ -335,7 +334,7 @@ namespace AMDiS { /// Returns number of managed problems - virtual int getNumProblems() + virtual int getNumProblems() const { return problems.size(); } @@ -348,17 +347,35 @@ namespace AMDiS { { return problems[number]; } + + virtual ProblemStatType const* getProblem(int number = 0) const + { + return problems[number]; + } /// Returns \ref meshes[i] - inline Mesh* getMesh(int number = 0) + Mesh* getMesh(int number = 0) + { + return meshes[number]; + } + + /// Returns \ref meshes[i] + Mesh const* getMesh(int number = 0) const { return meshes[number]; } - using super::getNumComponents; using super::getRefinementManager; using super::getCoarseningManager; using super::getName; + + virtual int getNumComponents() const override + { + int num = 0; + for (int i = 0; i < getNumProblems(); ++i) + num += getProblem(i)->getNumComponents(); + return num; + } protected: /// unqiue mesh-dimension for all problems diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc index f87dabf3398e21284085f3d1b1330d013135df64..ec7a8ad0feb572707e09d7c8c52edbe785705e63 100644 --- a/AMDiS/src/DOFAdmin.cc +++ b/AMDiS/src/DOFAdmin.cc @@ -215,10 +215,15 @@ namespace AMDiS { TEST_EXIT(dofIndexed)("no dofIndexed\n"); - if (dofIndexed->getSize() < size) - dofIndexed->resize(size); - - dofIndexedList.push_back(dofIndexed); + #ifdef _OPENMP + #pragma omp critical (addDOFIndexed) + #endif + { + if (dofIndexed->getSize() < size) + dofIndexed->resize(size); + + dofIndexedList.push_back(dofIndexed); + } } @@ -226,18 +231,17 @@ namespace AMDiS { { FUNCNAME("DOFAdmin::removeDOFIndexed()"); - bool removed = false; - std::list::iterator it; - std::list::iterator end = dofIndexedList.end(); - for (it = dofIndexedList.begin(); it != end; ++it) { - if (*it == dofIndexed) { - dofIndexedList.erase(it); - removed = true; - break; + #ifdef _OPENMP + #pragma omp critical (removeDOFIndexed) + #endif + { + auto it = std::find(dofIndexedList.begin(), dofIndexedList.end(), dofIndexed); + if (it != dofIndexedList.end()) + dofIndexedList.erase(it); + else { + ERROR_EXIT("DOFIndexed not in list!\n"); } } - - TEST_EXIT(removed)("DOFIndexed not in list!\n"); } @@ -255,16 +259,17 @@ namespace AMDiS { { FUNCNAME("DOFAdmin::removeDOFContainer()"); - std::list::iterator it; - std::list::iterator end = dofContainerList.end(); - for (it = dofContainerList.begin(); it != end; ++it) { - if (*it == cont) { - dofContainerList.erase(it); - return; + #ifdef _OPENMP + #pragma omp critical (removeDOFContainer) + #endif + { + auto it = std::find(dofContainerList.begin(), dofContainerList.end(), cont); + if (it != dofContainerList.end()) + dofContainerList.erase(it); + else { + ERROR_EXIT("Container not in list!\n"); } } - - ERROR("Container not in list!\n"); } diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h index 4583f00339aa0fd45c16a2c4132c3896363d6464..1a148267202738362b15465dc231e70fa2501be3 100644 --- a/AMDiS/src/DOFAdmin.h +++ b/AMDiS/src/DOFAdmin.h @@ -182,6 +182,12 @@ namespace AMDiS { return mesh; } + /// Returns \ref mesh + inline Mesh* getMesh() + { + return mesh; + } + /// Returns \ref dofFree, the array denoting DOFs to be either free or used. inline const std::vector& getDofFree() const { diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc index 02bd1ad546bb41f173e8244f9d920c71fe25bc63..9f85c975baa34e572da9782bf1ad390385497758 100644 --- a/AMDiS/src/DOFMatrix.cc +++ b/AMDiS/src/DOFMatrix.cc @@ -43,7 +43,7 @@ namespace AMDiS { DOFMatrix::DOFMatrix() : rowFeSpace(NULL), colFeSpace(NULL), - elementMatrix(3, 3), + elementMatrix(0, 0), nRow(0), nCol(0), nnzPerRow(0), diff --git a/AMDiS/src/DOFSerializer.cc b/AMDiS/src/DOFSerializer.cc new file mode 100644 index 0000000000000000000000000000000000000000..0ca8de2c678eea2094a9b352078747b000cc3b4f --- /dev/null +++ b/AMDiS/src/DOFSerializer.cc @@ -0,0 +1,135 @@ +#include "DOFSerializer.h" + +#include + +#include "DOFVector.h" +#include "Global.h" +#include "Traverse.h" + +namespace AMDiS +{ + +// satre values of DOFVector on all macro elements to internal container +void DOFSerializer::gather(DOFVector const* vec, std::vector& values) +{ + std::fill(visited_.begin(), visited_.end(), false); + values.clear(); + + for (auto* macroEl : mesh_->getMacroElements()) + gather(macroEl->getIndex(), vec, values, false); +} + + +// satre values of DOFVector on macro element `macroIndex` to internal container +void DOFSerializer::gather(int macroIndex, DOFVector const* vec, std::vector& values, bool reset) +{ + if (reset) { + std::fill(visited_.begin(), visited_.end(), false); + values.clear(); + } + + TEST_EXIT(mesh_ == vec->getFeSpace()->getMesh())("Incompatible meshes!\n"); + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirstOneMacro(mesh_, macroIndex, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + Element *el = elInfo->getElement(); + + if (el->isLeaf()) { + gather(VERTEX, elInfo, vec, values); + if (mesh_->getDim() > 1) + gather(EDGE, elInfo, vec, values); + if (mesh_->getDim() == 3) + gather(FACE, elInfo, vec, values); + gather(CENTER, elInfo, vec, values); + } + + elInfo = stack.traverseNext(elInfo); + } +} + + +void DOFSerializer::gather(GeoIndex geo, ElInfo* elInfo, DOFVector const* vec, std::vector& values) +{ + int nd; + if ((nd = admin_->getNumberOfDofs(geo))) { + int entities = mesh_->getGeo(geo); + int nd0 = admin_->getNumberOfPreDofs(geo); + int n0 = mesh_->getNode(geo); + for (int n = 0; n < entities; n++) { + for(int d = 0; d < nd; d++) { + DegreeOfFreedom globalDof = elInfo->getElement()->getDof(n0 + n, nd0 + d); + TEST_EXIT(globalDof < visited_.size())("visited container not large enough!\n"); + if (!visited_[globalDof]) { + visited_[globalDof] = true; + values.push_back((*vec)[globalDof]); + } + } + } + } +} + + +// assign stored values to DOFVector, for all macro elements +void DOFSerializer::scatter(std::vector const& values, DOFVector* vec) const +{ + std::fill(visited_.begin(), visited_.end(), false); + counter_ = 0u; + + for (auto* macroEl : mesh_->getMacroElements()) + scatter(macroEl->getIndex(), values, vec, false); +} + + +// assign stored values to DOFVector, on macroElement with index `macroIndex` +void DOFSerializer::scatter(int macroIndex, std::vector const& values, DOFVector* vec, bool reset) const +{ + if (reset) { + std::fill(visited_.begin(), visited_.end(), false); + counter_ = 0u; + } + + TEST_EXIT(mesh_ == vec->getFeSpace()->getMesh())("Incompatible meshes!\n"); + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirstOneMacro(mesh_, macroIndex, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + Element *el = elInfo->getElement(); + + if (el->isLeaf()) { + scatter(VERTEX, elInfo, values, vec); + if (mesh_->getDim() > 1) + scatter(EDGE, elInfo, values, vec); + if (mesh_->getDim() == 3) + scatter(FACE, elInfo, values, vec); + scatter(CENTER, elInfo, values, vec); + } + + elInfo = stack.traverseNext(elInfo); + } +} + + +void DOFSerializer::scatter(GeoIndex geo, ElInfo* elInfo, std::vector const& values, DOFVector* vec) const +{ + int nd; + if ((nd = admin_->getNumberOfDofs(geo))) { + int entities = mesh_->getGeo(geo); + int nd0 = admin_->getNumberOfPreDofs(geo); + int n0 = mesh_->getNode(geo); + for (int n = 0; n < entities; n++) { + for(int d = 0; d < nd; d++) { + DegreeOfFreedom globalDof = elInfo->getElement()->getDof(n0 + n, nd0 + d); + TEST_EXIT(globalDof < visited_.size())("visited container not large enough!\n"); + if (!visited_[globalDof]) { + visited_[globalDof] = true; + TEST_EXIT(counter_ < values.size())("Not enough values in value-container!\n"); + (*vec)[globalDof] = values[counter_++]; + } + } + } + } +} + + +} // end namespace AMDiS diff --git a/AMDiS/src/DOFSerializer.h b/AMDiS/src/DOFSerializer.h new file mode 100644 index 0000000000000000000000000000000000000000..3327c4562e3ab1eb349f9dfc9d5e4d2aeea7ac35 --- /dev/null +++ b/AMDiS/src/DOFSerializer.h @@ -0,0 +1,91 @@ +/****************************************************************************** + * + * AMDiS - Adaptive multidimensional simulations + * + * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. + * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis + * + * Authors: + * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al. + * + * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE + * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. + * + * + * This file is part of AMDiS + * + * See also license.opensource.txt in the distribution. + * + ******************************************************************************/ + + + +/** \file DOFSerializer.h */ + +#ifndef AMDIS_DOFSERIALIZER_H +#define AMDIS_DOFSERIALIZER_H + +#include +#include + +#include "DOFAdmin.h" +#include "ElInfo.h" +#include "Global.h" +#include "Mesh.h" +#include "Traverse.h" + +namespace AMDiS +{ + // Collect all DOFValues into a std::vector, by traversing the mesh macro-elenment wise + // Can be used to transfer data on macro elements or to store values in a file, e.g. in ARHWriter + class DOFSerializer + { + public: + + /// Constructor, stores a pointer to the DOFAdmin. If the number of values to store can be estimated, + /// pass a second argument `numValues` + DOFSerializer(DOFAdmin* admin, std::size_t numValues = 0) + : admin_(admin) + , mesh_(admin->getMesh()) + , visited_(admin->getUsedSize(), false) + {} + + + // satre values of DOFVector on all macro elements to internal container + void gather(DOFVector const* vec, std::vector& values); + + // satre values of DOFVector on macro element `macroIndex` to internal container + void gather(int macroIndex, DOFVector const* vec, std::vector& values, bool reset = true); + + + // assign stored values to DOFVector, for all macro elements + void scatter(std::vector const& values, DOFVector* vec) const; + + // assign stored values to DOFVector, on macroElement with index `macroIndex` + void scatter(int macroIndex, std::vector const& values, DOFVector* vec, bool reset = true) const; + + + protected: + + // collect values of one geometry-type + void gather(GeoIndex geo, ElInfo* elInfo, DOFVector const* vec, std::vector& values); + + // assign values of one geometry type + void scatter(GeoIndex geo, ElInfo* elInfo, std::vector const& values, DOFVector* vec) const; + + + protected: + + DOFAdmin* admin_; + Mesh* mesh_; + + // stored which DOFs were already visited during traversal + mutable std::vector visited_; + + // a counter used during assignment of values. To guarantee the same order as during gather + mutable std::size_t counter_ = 0u; + }; + +} // end namespace AMDiS + +#endif // AMDIS_DOFSERIALIZER_H diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc index ac2a3f7a2e5aedbaa009e06a05eb6dc82d6433ca..9ee3948153e05fff98134fccf7ff87b85d0313eb 100644 --- a/AMDiS/src/DOFVector.cc +++ b/AMDiS/src/DOFVector.cc @@ -273,7 +273,7 @@ namespace AMDiS { template<> - void DOFVector::interpol(DOFVector *source, double factor) + void DOFVector::interpol(const DOFVector *source, double factor) { const FiniteElemSpace *sourceFeSpace = source->getFeSpace(); @@ -383,7 +383,7 @@ namespace AMDiS { template<> - void DOFVector >::interpol(DOFVector > *v, + void DOFVector >::interpol(const DOFVector > *v, double factor) { FUNCNAME("DOFVector >::interpol()"); diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h index 75e4d0fc838623c8b6e33d6d915e926a76b84c37..348aa9be94d241598e71428f73c82b0ab3f0ac2f 100644 --- a/AMDiS/src/DOFVector.h +++ b/AMDiS/src/DOFVector.h @@ -574,7 +574,7 @@ namespace AMDiS { /// stores these in the DOFVector void interpol(AbstractFunction > *fct); - void interpol(DOFVector *v, double factor = 1.0); + void interpol(const DOFVector *v, double factor = 1.0); /// Eval DOFVector at given point p. If oldElInfo != NULL the search for diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 4c229e78bc3fed3acf22cf7bb23774cc3f0c734e..61b1c7678340bcc300616a8ba1d1af8edb274a7c 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -124,8 +124,10 @@ namespace AMDiS { { this->name = n; this->feSpace = f; + if (this->feSpace && this->feSpace->getAdmin()) (this->feSpace->getAdmin())->addDOFIndexed(this); + this->boundaryManager = new BoundaryManager(f); #ifdef HAVE_PARALLEL_DOMAIN_AMDIS if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL) @@ -140,6 +142,10 @@ namespace AMDiS { if ( Parallel::MeshDistributor::globalMeshDistributor != NULL) Parallel::MeshDistributor::globalMeshDistributor->removeInterchangeVector(this); #endif + + #ifdef _OPENMP + #pragma omp critical + #endif if (this->feSpace && this->feSpace->getAdmin()) (this->feSpace->getAdmin())->removeDOFIndexed(this); diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc index 8a90e6445a9678acd58ce54474d8800d46adc512..32932afe04263ccce0ca5f26f1f9e56ba77a10ea 100644 --- a/AMDiS/src/DualTraverse.cc +++ b/AMDiS/src/DualTraverse.cc @@ -111,6 +111,93 @@ namespace AMDiS { } + bool DualTraverse::traverseFirstOneMacro(Mesh *mesh1, + Mesh *mesh2, + int macroIndex, + int level1, + int level2, + Flag flag1, + Flag flag2, + ElInfo **elInfo1, + ElInfo **elInfo2, + ElInfo **elInfoSmall, + ElInfo **elInfoLarge) + { + FUNCNAME("DualTraverse::traverseFirst()"); + // replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain) + if (flag1.isSet(Mesh::CALL_EL_LEVEL)) { + flag1 &= ~Mesh::CALL_EL_LEVEL; + flag1 |= Mesh::CALL_MG_LEVEL; + level1_ = level1; + } else { + level1_ = -1; + } + + if (flag2.isSet(Mesh::CALL_EL_LEVEL)) { + flag2 &= ~Mesh::CALL_EL_LEVEL; + flag2 |= Mesh::CALL_MG_LEVEL; + level2_ = level2; + } else { + level2_ = -1; + } + + // replace CALL_LEAF_EL_LEVEL by CALL_MG_LEVEL (covers whole domain) + if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) { + flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL; + flag1 |= Mesh::CALL_MG_LEVEL; + level1_ = level1; + callLeafElLevel1_ = true; + } else { + level1_ = -1; + callLeafElLevel1_ = false; + } + + if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) { + flag2 &= ~Mesh::CALL_LEAF_EL_LEVEL; + flag2 |= Mesh::CALL_MG_LEVEL; + level2_ = level2; + callLeafElLevel2_ = true; + } else { + level2_ = -1; + callLeafElLevel2_ = false; + } + + // call standard traverse + *elInfo1 = stack1.traverseFirstOneMacro(mesh1, macroIndex, level1, flag1); + while (*elInfo1 != NULL && skipEl1(*elInfo1)) { + *elInfo1 = stack1.traverseNext(*elInfo1); + } + + *elInfo2 = stack2.traverseFirstOneMacro(mesh2, macroIndex, level2, flag2); + while (*elInfo2 != NULL && skipEl2(*elInfo2)) { + *elInfo2 = stack2.traverseNext(*elInfo2); + } + + // finished ? + if (*elInfo1 == NULL || *elInfo2 == NULL) { + TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n"); + return false; + } + + rest = 1.0; + + bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge); + + // check for non domain covering level traverse + if (!accepted || + (level1_ != -1 && (*elInfo1)->getLevel() != level1_) || + (callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) || + (level2_ != -1 && (*elInfo2)->getLevel() != level2_) || + (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) { + return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge); + } + + fillSubElInfo(*elInfo1, *elInfo2, *elInfoSmall, *elInfoLarge); + + return true; + } + + bool DualTraverse::traverseNext(ElInfo **elInfo1, ElInfo **elInfo2, ElInfo **elInfoSmall, diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h index f5990540ca9ad82882d04211810aff8f3d0d4c36..3a43e2a864b836a8550b2e2fad075372fb3c8bdf 100644 --- a/AMDiS/src/DualTraverse.h +++ b/AMDiS/src/DualTraverse.h @@ -78,6 +78,18 @@ namespace AMDiS { &(dualElInfo.smallElInfo), &(dualElInfo.largeElInfo)); } + + bool traverseFirstOneMacro(Mesh *mesh1, + Mesh *mesh2, + int macroIndex, + int level1, + int level2, + Flag flag1, + Flag flag2, + ElInfo **elInfo1, + ElInfo **elInfo2, + ElInfo **elInfoSmall, + ElInfo **elInfoLarge); /// Get next ElInfo combination bool traverseNext(ElInfo **elInfoNext1, diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc index e4836a0520c40c825737957ef49cd1f978eb4df4..71a21d1f98114cfd014eb5515d982e70f1f370f2 100644 --- a/AMDiS/src/ElInfo1d.cc +++ b/AMDiS/src/ElInfo1d.cc @@ -161,8 +161,6 @@ namespace AMDiS { double length = (coord[1][0] - a); int dim = mesh->getDim(); - static DimVec vec(dim, NO_INIT); - TEST_EXIT_DBG(lambda)("lambda must not be NULL\n"); TEST_EXIT_DBG(dim == 1)("dim!=1\n"); TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n"); diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc index 7d7e2119205023152506d350940a623b42dbef97..3bc7d64db20b327dcda9339faaa13ad825deaa37 100644 --- a/AMDiS/src/ElInfo2d.cc +++ b/AMDiS/src/ElInfo2d.cc @@ -679,7 +679,6 @@ namespace AMDiS { DimVec > edge(mesh->getDim(), NO_INIT); WorldVector x; - static DimVec vec(mesh->getDim(), NO_INIT); int dim = mesh->getDim(); @@ -875,4 +874,4 @@ namespace AMDiS { { return getSubElemCoordsMat(degree); } -} \ No newline at end of file +} diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc index 87fd26af81825bb39a8b7f223b503bac608e9846..f029d67a6edb63d43c25b17d6d78b4d9c7a85545 100644 --- a/AMDiS/src/ElInfo3d.cc +++ b/AMDiS/src/ElInfo3d.cc @@ -336,8 +336,6 @@ namespace AMDiS { WorldVector x; double x0, det, det0, det1, det2; - static DimVec vec(mesh->getDim(), NO_INIT); - TEST_EXIT_DBG(lambda)("lambda must not be NULL\n"); int dim = mesh->getDim(); diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc index e332728ee07eda1bfb846cbb0641284dbd2dce7f..996863876a55a3ad6eddcfcf885ce179bc94a91a 100644 --- a/AMDiS/src/FirstOrderAssembler.cc +++ b/AMDiS/src/FirstOrderAssembler.cc @@ -204,11 +204,16 @@ namespace AMDiS { void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) { if (firstCall) { - const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); - psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); - basFcts = colFeSpace->getBasisFcts(); - phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); + psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); + basFcts = colFeSpace->getBasisFcts(); + phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); + firstCall = false; + } } int nPoints = quadrature->getNumPoints(); @@ -240,11 +245,16 @@ namespace AMDiS { void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec) { if (firstCall) { - const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); - psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); - basFcts = colFeSpace->getBasisFcts(); - phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); + psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); + basFcts = colFeSpace->getBasisFcts(); + phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); + firstCall = false; + } } int nPoints = quadrature->getNumPoints(); @@ -280,11 +290,16 @@ namespace AMDiS { const double *values; if (firstCall) { - q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(), - colFeSpace->getBasisFcts(), - quadrature); - q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(), + colFeSpace->getBasisFcts(), + quadrature); + q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); + firstCall = false; + } } const int **nEntries = q10->getNumberEntries(); @@ -362,11 +377,16 @@ namespace AMDiS { void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) { if (firstCall) { - const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); - psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); - basFcts = colFeSpace->getBasisFcts(); - phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); + psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); + basFcts = colFeSpace->getBasisFcts(); + phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI); + firstCall = false; + } } int nPoints = quadrature->getNumPoints(); @@ -410,11 +430,16 @@ namespace AMDiS { const double *values; if (firstCall) { - q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(), - colFeSpace->getBasisFcts(), - quadrature); - q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(), + colFeSpace->getBasisFcts(), + quadrature); + q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); + firstCall = false; + } } const int **nEntries = q01->getNumberEntries(); @@ -446,11 +471,16 @@ namespace AMDiS { const double *values; if (firstCall) { - q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(), - colFeSpace->getBasisFcts(), - quadrature); - q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(), + colFeSpace->getBasisFcts(), + quadrature); + q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature); + firstCall = false; + } } const int *nEntries = q1->getNumberEntries(); diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc index bace7aee050dbb5d8f314178504325654d01c6ae..865b82feddbb309f38d05b74f34038533da0bdf5 100644 --- a/AMDiS/src/Global.cc +++ b/AMDiS/src/Global.cc @@ -32,11 +32,14 @@ #ifdef HAVE_PARALLEL_PETSC #include "petsc.h" #endif +#ifdef HAVE_OPENMP +#include +#endif namespace AMDiS { const char *funcName = NULL; -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_OPENMP) || defined(HAVE_MPI) bool Msg::outputMainRank = true; #endif @@ -129,10 +132,14 @@ namespace AMDiS { void Msg::print_funcname(const char *funcName) { -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI) if (outputMainRank && MPI::COMM_WORLD.Get_rank() != 0) return; #endif +#ifdef HAVE_OPENMP + if (outputMainRank && omp_get_thread_num() != 0) + return; +#endif if (!out) out = &std::cout; @@ -220,10 +227,14 @@ namespace AMDiS { const char *file, int line) { -#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(NDEBUG) +#if (defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI)) && defined(NDEBUG) if (outputMainRank && MPI::COMM_WORLD.Get_rank() != 0) return; #endif +#ifdef HAVE_OPENMP + if (outputMainRank && omp_get_thread_num() != 0) + return; +#endif static int old_line = -1; @@ -254,9 +265,13 @@ namespace AMDiS { void Msg::print_warn(const char *format, ...) { -#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(NDEBUG) +#if (defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI)) && defined(NDEBUG) if (outputMainRank && MPI::COMM_WORLD.Get_rank() != 0) return; +#endif +#ifdef HAVE_OPENMP + if (outputMainRank && omp_get_thread_num() != 0) + return; #endif va_list arg; char buff[255]; @@ -274,10 +289,14 @@ namespace AMDiS { void Msg::print(const char *format, ...) { #ifndef SUPPRESS_OUTPUT -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI) if (outputMainRank && MPI::COMM_WORLD.Get_rank() != 0) return; #endif +#ifdef HAVE_OPENMP + if (outputMainRank && omp_get_thread_num() != 0) + return; +#endif va_list arg; char buff[255]; @@ -298,7 +317,6 @@ namespace AMDiS { int d = -1; // get dimension - TEST_EXIT(Parameters::initialized())("Parameters not initialized!\n"); Parameters::get("dimension of world",d,0); TEST_EXIT(d > 0)("Cannot initialize dimension!\n"); TEST_EXIT((d == 1) || (d == 2) || (d == 3))("Invalid world dimension %d!\n",d); diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h index 9c20ce4b368c394b9c9e9ca038d8674e3b159326..1b3e0f6c64a36ad8b47b2b03a4d1b7238e5b6957 100644 --- a/AMDiS/src/Global.h +++ b/AMDiS/src/Global.h @@ -53,7 +53,7 @@ #include #endif -#if HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI) #include #endif @@ -131,8 +131,10 @@ namespace AMDiS { #if SUPPRESS_OUTPUT #define PRINT_LINE(stream, line) #else -#if HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI) #define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line +#elif defined(HAVE_OPENMP) +#define PRINT_LINE(stream, line) stream << "[" << omp_get_thread_num() << "] " << line #else #define PRINT_LINE(stream, line) stream << line #endif @@ -331,7 +333,7 @@ namespace AMDiS { } public: -#if HAVE_PARALLEL_DOMAIN_AMDIS +#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_OPENMP) || defined(HAVE_MPI) /// In parallel computations, when this variable is true, only the 0 rank will /// print messages to the output stream. Error messages and warnings are always /// printed from all ranks. diff --git a/AMDiS/src/Initfile.cc b/AMDiS/src/Initfile.cc index bd7976ba17066869320aeca8250bcdda71fa804d..222f7f6047261270b5cc3da6b546bb466ecc431a 100644 --- a/AMDiS/src/Initfile.cc +++ b/AMDiS/src/Initfile.cc @@ -67,17 +67,13 @@ namespace AMDiS { } } - Initfile* Initfile::singlett = NULL; - std::set Initfile::fn_include_list; - - /// initialize singleton object an global parameters + /// initialize singleton object and global parameters void Initfile::init(std::string in) { - initIntern(); - fn_include_list.clear(); - singlett->read(in); - singlett->getInternalParameters(); + getIncludeList().clear(); + singlett().read(in); + singlett().getInternalParameters(); } @@ -85,13 +81,13 @@ namespace AMDiS { void Initfile::read(std::string fn, bool force) { // read file if its not parsed already - if (fn_include_list.find(fn) == fn_include_list.end() || force) { + if (getIncludeList().find(fn) == getIncludeList().end() || force) { std::ifstream inputFile; inputFile.open(fn.c_str(), std::ios::in); if (!inputFile.is_open()) throw runtime_error("init-file '" + fn + "' cannot be opened for reading"); - fn_include_list.insert(fn); + getIncludeList().insert(fn); read(inputFile); } } @@ -233,8 +229,6 @@ namespace AMDiS { void Initfile::readArgv(std::string parameters, int debugInfo) { - initIntern(); - char seperator = ';'; typedef boost::escaped_list_separator TokenizerFunc; typedef boost::tokenizer Tokenizer; @@ -286,8 +280,7 @@ namespace AMDiS { /// print all parameters to std::cout void Initfile::printParameters() { - initIntern(); - for (Initfile::iterator it = singlett->begin(); it != singlett->end(); it++) + for (Initfile::iterator it = singlett().begin(); it != singlett().end(); it++) std::cout << (*it).first << " => " << (*it).second << std::endl; } diff --git a/AMDiS/src/Initfile.h b/AMDiS/src/Initfile.h index add610cf4912231d9597388678a3fa3b5629af84..260d38d3c1ad647710f646701f11ced35edcd97e 100644 --- a/AMDiS/src/Initfile.h +++ b/AMDiS/src/Initfile.h @@ -37,6 +37,7 @@ #include #include +#include // requires c++14 #include @@ -130,7 +131,7 @@ namespace AMDiS { /// convert string to string - inline void convert(const std::string valStr, std::string& value) + inline void convert(const std::string& valStr, std::string& value) { value = trim(valStr); } @@ -144,7 +145,7 @@ namespace AMDiS { boost::mpl::not_< boost::is_enum > >, void >::type - convert(const std::string valStr, T& value) + convert(const std::string& valStr, T& value) { using boost::lexical_cast; using boost::numeric_cast; @@ -165,7 +166,7 @@ namespace AMDiS { < boost::is_enum, void >::type - convert(const std::string valStr, T& value) + convert(const std::string& valStr, T& value) { int swap = 0; // try { @@ -178,7 +179,7 @@ namespace AMDiS { /// convert special enums - inline void convert(const std::string valStr, Norm& value) + inline void convert(const std::string& valStr, Norm& value) { std::string swapStr = boost::to_upper_copy(valStr); @@ -199,7 +200,7 @@ namespace AMDiS { /// convert value of arbitrary type to string using stringstream and /// operator<< for type template - inline void convert(const T value, std::string& valStr) + inline void convert(const T& value, std::string& valStr) { std::stringstream ss; ss.precision(6); @@ -222,20 +223,20 @@ namespace AMDiS { // forward declarations template< typename T > - inline void convert(const std::string valStr, WorldVector& c); + inline void convert(const std::string& valStr, WorldVector& c); template - inline void convert(const std::string valStr, std::list& value); + inline void convert(const std::string& valStr, std::list& value); template - inline void convert(const std::string valStr, std::vector& value); + inline void convert(const std::string& valStr, std::vector& value); /** parse an container from tag tag. The Container must have the properties: * - type value_type * - member function push_back */ template< typename Container > - inline void getContainer(const std::string val_, Container& c) + inline void getContainer(const std::string& val_, Container& c) { // accepted brackets and delimiters for vector input std::string begBrackets= "{[("; @@ -285,7 +286,7 @@ namespace AMDiS { /// convert string to WorldVector template< typename T > - inline void convert(const std::string valStr, WorldVector& c) + inline void convert(const std::string& valStr, WorldVector& c) { std::vector temp_vec; getContainer(valStr, temp_vec); @@ -298,7 +299,7 @@ namespace AMDiS { /// convert string to std::list using begBrackets, endBrackets and delims template - inline void convert(const std::string valStr, std::list& value) + inline void convert(const std::string& valStr, std::list& value) { getContainer(valStr, value); } @@ -306,7 +307,7 @@ namespace AMDiS { /// convert string to std::vector using begBrackets, endBrackets and delims template - inline void convert(const std::string valStr, std::vector& value) + inline void convert(const std::string& valStr, std::vector& value) { getContainer(valStr, value); } @@ -422,18 +423,19 @@ namespace AMDiS { template static void get(const std::string tag, T& value, int debugInfo = -1) { - initIntern(); +// boost::shared_lock lock(singlett().mutex_); + if (debugInfo == -1) - debugInfo = singlett->getMsgInfo(); + debugInfo = singlett().getMsgInfo(); else { int swap(debugInfo); - debugInfo = singlett->getMsgInfo(); - singlett->msgInfo=swap; + debugInfo = singlett().getMsgInfo(); + singlett().msgInfo=swap; } std::string valStr; try { - int error_code = singlett->checkedGet(tag, valStr); + int error_code = singlett().checkedGet(tag, valStr); if (error_code == 0) { valStr = trim(valStr); detail::convert(valStr, value); @@ -454,7 +456,7 @@ namespace AMDiS { std::cout << "Parameter '" << tag << "'" << " initialized with: " << value << std::endl; } - singlett->msgInfo = debugInfo; + singlett().msgInfo = debugInfo; } @@ -471,8 +473,7 @@ namespace AMDiS { std::map &pm, int debugInfo = -1) { - initIntern(); - for (Initfile::iterator it = singlett->begin(); it != singlett->end(); it++){ + for (Initfile::iterator it = singlett().begin(); it != singlett().end(); it++){ std::string longTag= (*it).first ; std::string value=(*it).second; if(longTag.length()>tag.length() && @@ -497,10 +498,12 @@ namespace AMDiS { /// return InitEntry object for tag tag static InitEntry get(const std::string tag) { +// boost::shared_lock lock(singlett().mutex_); + InitEntry result; std::string valStr; - int error_code = singlett->checkedGet(tag, valStr); + int error_code = singlett().checkedGet(tag, valStr); if (error_code == 0) { valStr = trim(valStr); result = InitEntry(valStr); @@ -522,17 +525,18 @@ else if(error_code == TAG_NOT_FOUND_BREAK) /// update map tag->value_old to tag->value in singleton template - static void set(const std::string tag, T& value, int debugInfo= -1) + static void set(const std::string tag, T const& value, int debugInfo= -1) { - initIntern(); +// boost::unique_lock lock(singlett().mutex_); + if (debugInfo == -1) - debugInfo = singlett->getMsgInfo(); + debugInfo = singlett().getMsgInfo(); std::string swap = ""; detail::convert(value, swap); - (*singlett)[trim(tag)] = swap; + singlett()[trim(tag)] = swap; // update msg parameters msgInfo, msgWait, paramInfo - singlett->getInternalParameters(); + singlett().getInternalParameters(); if (debugInfo == 2) std::cout << "Parameter '" << tag << "'" << " set to: " << value << std::endl; @@ -541,7 +545,7 @@ else if(error_code == TAG_NOT_FOUND_BREAK) /// add map tag->value to data in singleton template< typename T > - static void add(const std::string tag, T& value, int debugInfo = -1) + static void add(const std::string tag, T const& value, int debugInfo = -1) { set(tag, value, debugInfo); } @@ -554,28 +558,14 @@ else if(error_code == TAG_NOT_FOUND_BREAK) /// Returns specified info level static int getMsgInfo() { - return (singlett != NULL) ? singlett->msgInfo : 0; + return singlett().msgInfo; } /// Returns specified wait value static int getMsgWait() { - return (singlett != NULL) ? singlett->msgWait : 0; - } - - - /// Checks whether parameters are initialized. if not, call init() - static bool initialized() - { - return (singlett != NULL); - } - - - /// return pointer to singleton - static Initfile *getSingleton() - { - return singlett; + return singlett().msgWait; } @@ -586,21 +576,24 @@ else if(error_code == TAG_NOT_FOUND_BREAK) /// clear data in singleton static void clearData() { - initIntern(); - singlett->clear(); + singlett().clear(); } - + /// save singlett-data to file with filename fn static void save(std::string fn) { - initIntern(); - singlett->write(fn); + singlett().write(fn); } + - static std::set< std::string > getIncludeList(){ - return fn_include_list; + // list of included filenames + static std::set< std::string >& getIncludeList() + { + static std::set< std::string > include_list; + return include_list; } + protected: Initfile() : msgInfo(0), @@ -609,22 +602,14 @@ protected: breakOnMissingTag(0) {} - - static void initIntern() + /// return reference to singleton + static Initfile& singlett() { - if (singlett == NULL) - singlett = new Initfile; + static Initfile singlett_{}; + return singlett_; } - /// list of processed files - static std::set< std::string > fn_include_list; - - - /// pointer to the singleton that contains the data - static Initfile* singlett; - - /// return the value of the given tag or throws an exception if the tag /// does not exist int checkedGet(const std::string& tag, std::string &valStr) const @@ -671,6 +656,8 @@ protected: void getInternalParameters(); int msgInfo, msgWait, paramInfo, breakOnMissingTag; + +// mutable boost::shared_mutex mutex_; }; typedef Initfile Parameters; diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index 2c9ae1ebb93176d80e3c7cf19907f3b4e1c754a8..122a4019945edd3d0a5446a802af3789ef656d84 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -911,6 +911,7 @@ namespace AMDiS { friend struct io::MacroWriter; friend class MacroElement; friend class Element; + friend class MeshRestriction; }; } diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc index b23b4911d93d3d990b4f35cee85b5c03bdfa785c..059d4356b0707e7dca87ecf80d7ee3a668a4cf55 100644 --- a/AMDiS/src/MeshStructure.cc +++ b/AMDiS/src/MeshStructure.cc @@ -280,25 +280,25 @@ namespace AMDiS { bool cont = true; while (cont) { bool cont1; -#ifndef NDEBUG +// #ifndef NDEBUG bool cont2; -#endif +// #endif if (structure1->isLeafElement() == structure2->isLeafElement()) { cont1 = structure1->nextElement(result); -#ifndef NDEBUG +// #ifndef NDEBUG cont2 = structure2->nextElement(); -#endif +// #endif } else { if (structure1->isLeafElement()) { cont1 = structure1->nextElement(); -#ifndef NDEBUG +// #ifndef NDEBUG cont2 = structure2->skipBranch(result); -#endif +// #endif } else { cont1 = structure1->skipBranch(result); -#ifndef NDEBUG +// #ifndef NDEBUG cont2 = structure2->nextElement(); -#endif +// #endif } } TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n"); @@ -504,7 +504,7 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } } - + void MeshStructure::setMeshStructureValues(int macroElIndex, DOFVector* vec, @@ -519,8 +519,8 @@ namespace AMDiS { Mesh *mesh = feSpace->getMesh(); bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1); int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); - unsigned int counter = 0; - + + std::size_t counter = 0; if (withElIndex) { TEST_EXIT(static_cast(values[0]) == macroElIndex) ("Value structure code was created for macro element index %d, but should be set to macro element index %d\n", @@ -568,5 +568,4 @@ namespace AMDiS { TEST_EXIT_DBG(values.size() == counter) ("Should not happen! values size %d, counter %d\n", values.size(), counter); } - } diff --git a/AMDiS/src/ProblemInstat.h b/AMDiS/src/ProblemInstat.h index e7ea6f5ffe5f95b6fc72ea0974fb30f67bdf7e01..2c22001b7fb58ba8f216d4d7bc6ce825cb04c262 100644 --- a/AMDiS/src/ProblemInstat.h +++ b/AMDiS/src/ProblemInstat.h @@ -101,7 +101,7 @@ namespace AMDiS { {} /// Returns \ref name. - inline std::string getName() + inline std::string getName() const { return name; } diff --git a/AMDiS/src/ProblemIterationInterface.h b/AMDiS/src/ProblemIterationInterface.h index b3240ff654ba67c0c35b58ebd114e6fbfb49dd62..b12bec1187702af7ac874cf752f4f28b16aba380 100644 --- a/AMDiS/src/ProblemIterationInterface.h +++ b/AMDiS/src/ProblemIterationInterface.h @@ -68,13 +68,19 @@ namespace AMDiS { virtual void endIteration(AdaptInfo *adaptInfo) {} /// Returns number of managed problems - virtual int getNumProblems() = 0; + virtual int getNumProblems() const = 0; /** \brief * Returns the problem with the given number. If only one problem * is managed by this master problem, the number hasn't to be given. */ virtual ProblemStatBase *getProblem(int number = 0) = 0; + + // dirty hack + virtual ProblemStatBase const* getProblem(int number = 0) const + { + return const_cast(const_cast(this)->getProblem(number)); + } /// Returns the problem with the given name. virtual ProblemStatBase *getProblem(std::string name) @@ -83,7 +89,7 @@ namespace AMDiS { } /// Returns the name of the problem. - virtual std::string getName() = 0; + virtual std::string getName() const = 0; /// Function that serializes the problem plus information about the iteration. virtual void serialize(std::ostream &out) = 0; diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h index 449a1de3b79d24004aaaa989adb246b8f4790923..dffaa5024732bc1d96ab6756ea22b2e9b80ca7fa 100644 --- a/AMDiS/src/ProblemStat.h +++ b/AMDiS/src/ProblemStat.h @@ -155,13 +155,13 @@ namespace AMDiS { bool asmMatrix = true, bool asmVector = true); /// Returns nr of components \ref nComponents - virtual int getNumComponents() + virtual int getNumComponents() const { return nComponents; } /// Returns nr of additional components \ref nAddComponents - virtual int getNumAddComponents() + virtual int getNumAddComponents() const { return nAddComponents; } @@ -280,12 +280,22 @@ namespace AMDiS { { return solution; } + + inline SystemVector const* getSolution() const + { + return solution; + } inline DOFVector* getSolution(int i) { return solution->getDOFVector(i); } + inline DOFVector const* getSolution(int i) const + { + return solution->getDOFVector(i); + } + /// Returns \ref rhs. inline SystemVector* getRhs() { @@ -317,6 +327,14 @@ namespace AMDiS { ("invalid component number\n"); return componentMeshes[comp]; } + + Mesh const* getMesh(int comp = 0) const + { + FUNCNAME("ProblemStatSeq::getMesh()"); + TEST_EXIT(comp < static_cast(componentMeshes.size()) && comp >= 0) + ("invalid component number\n"); + return componentMeshes[comp]; + } /// Returns \ref meshes inline std::vector getMeshes() @@ -325,7 +343,7 @@ namespace AMDiS { } /// Returns \ref feSpace_. - inline const FiniteElemSpace* getFeSpace(int comp = 0) + inline const FiniteElemSpace* getFeSpace(int comp = 0) const { FUNCNAME("ProblemStatSeq::getFeSpace()"); TEST_EXIT(comp < static_cast(componentSpaces.size()) && comp >= 0) @@ -388,13 +406,13 @@ namespace AMDiS { } /// Returns the name of the problem - inline std::string getName() override + inline std::string getName() const override { return name; } /// Returns the name of the problem - inline std::string getComponentName(int comp = 0) + inline std::string getComponentName(int comp = 0) const { TEST_EXIT(comp < static_cast(componentNames.size()) && comp >= 0) ("invalid component number\n"); @@ -402,19 +420,19 @@ namespace AMDiS { } /// Returns \ref useGetBound. - inline bool getBoundUsed() + inline bool getBoundUsed() const { return useGetBound; } /// Returns \ref info. - int getInfo() + int getInfo() const { return info; } /// Returns \ref deserialized; - bool isDeserialized() + bool isDeserialized() const { return deserialized; } @@ -556,13 +574,13 @@ namespace AMDiS { } /// Returns \ref solutionTime. - double getSolutionTime() + double getSolutionTime() const { return solutionTime; } /// Returns \ref buildTime. - double getBuildTime() + double getBuildTime() const { return buildTime; } @@ -714,7 +732,7 @@ namespace AMDiS { /// Returns number of managed problems // implements StandardProblemIteration::getNumProblems() - virtual int getNumProblems() override + virtual int getNumProblems() const override { return 1; } @@ -726,6 +744,12 @@ namespace AMDiS { { return this; } + + /// Returns the problem with the given number. + virtual ProblemStatBase const* getProblem(int number = 0) const override + { + return this; + } /// Determines the execution order of the single adaption steps. If adapt is /// true, mesh adaption will be performed. This allows to avoid mesh adaption, diff --git a/AMDiS/src/ProblemStatBase.h b/AMDiS/src/ProblemStatBase.h index e1d7c2be9e76043a89ec79a476937d75f30f6d9c..25ec3e856f857f2b68f01b8ab964103d783b4df6 100644 --- a/AMDiS/src/ProblemStatBase.h +++ b/AMDiS/src/ProblemStatBase.h @@ -140,7 +140,7 @@ namespace AMDiS { virtual void estimate(AdaptInfo *adaptInfo) = 0; /// Returns the name of the problem. - virtual std::string getName() = 0; + virtual std::string getName() const = 0; /// Function that serializes the problem plus information about the iteration. virtual void serialize(std::ostream &out) = 0; diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc index dce37e565226af6f8a99177444c88cc3fd9df71b..bbc08352c5621274f8f6d69e1401bbf5a36c6e11 100644 --- a/AMDiS/src/Quadrature.cc +++ b/AMDiS/src/Quadrature.cc @@ -34,21 +34,13 @@ namespace AMDiS { int FastQuadrature::max_points = 0; Quadrature::Quadrature(const Quadrature& q) - { - name = q.name; - degree = q.degree; - dim = q.dim; - n_points = q.n_points; - - // copy barycentric coordinates - lambda = new VectorOfFixVecs >(*(q.lambda)); - - // copy weights - w = new double[n_points]; - - for (int i = 0; i < n_points; i++) - w[i] = q.w[i]; - } + : name(q.name) + , degree(q.degree) + , dim(q.dim) + , n_points(q.n_points) + , lambda(new VectorOfFixVecs >(*(q.lambda))) + , w(q.w) + {} /****************************************************************************/ /* initialize gradient values of a function f in local coordinates at the */ @@ -115,76 +107,7 @@ namespace AMDiS { } - Quadrature **Quadrature::quad_nd[4]; - Quadrature *Quadrature::quad_0d[1]; - Quadrature *Quadrature::quad_1d[20]; - Quadrature *Quadrature::quad_2d[18]; - Quadrature *Quadrature::quad_3d[8]; - - VectorOfFixVecs > *Quadrature::x_0d; - double *Quadrature::w_0d; - - VectorOfFixVecs > *Quadrature::x0_1d = NULL; - VectorOfFixVecs > *Quadrature::x1_1d; - VectorOfFixVecs > *Quadrature::x2_1d; - VectorOfFixVecs > *Quadrature::x3_1d; - VectorOfFixVecs > *Quadrature::x4_1d; - VectorOfFixVecs > *Quadrature::x5_1d; - VectorOfFixVecs > *Quadrature::x6_1d; - VectorOfFixVecs > *Quadrature::x7_1d; - VectorOfFixVecs > *Quadrature::x8_1d; - VectorOfFixVecs > *Quadrature::x9_1d; - double *Quadrature::w0_1d; - double *Quadrature::w1_1d; - double *Quadrature::w2_1d; - double *Quadrature::w3_1d; - double *Quadrature::w4_1d; - double *Quadrature::w5_1d; - double *Quadrature::w6_1d; - double *Quadrature::w7_1d; - double *Quadrature::w8_1d; - double *Quadrature::w9_1d; - - VectorOfFixVecs > *Quadrature::x1_2d; - VectorOfFixVecs > *Quadrature::x2_2d; - VectorOfFixVecs > *Quadrature::x3_2d; - VectorOfFixVecs > *Quadrature::x4_2d; - VectorOfFixVecs > *Quadrature::x5_2d; - VectorOfFixVecs > *Quadrature::x7_2d; - VectorOfFixVecs > *Quadrature::x8_2d; - VectorOfFixVecs > *Quadrature::x9_2d; - VectorOfFixVecs > *Quadrature::x10_2d; - VectorOfFixVecs > *Quadrature::x11_2d; - VectorOfFixVecs > *Quadrature::x12_2d; - VectorOfFixVecs > *Quadrature::x17_2d; - double *Quadrature::w1_2d; - double *Quadrature::w2_2d; - double *Quadrature::w3_2d; - double *Quadrature::w4_2d; - double *Quadrature::w5_2d; - double *Quadrature::w7_2d; - double *Quadrature::w8_2d; - double *Quadrature::w9_2d; - double *Quadrature::w10_2d; - double *Quadrature::w11_2d; - double *Quadrature::w12_2d; - double *Quadrature::w17_2d; - - VectorOfFixVecs > *Quadrature::x1_3d; - VectorOfFixVecs > *Quadrature::x2_3d; - VectorOfFixVecs > *Quadrature::x3_3d; - VectorOfFixVecs > *Quadrature::x4_3d; - VectorOfFixVecs > *Quadrature::x5_3d; - VectorOfFixVecs > *Quadrature::x7_3d; - double *Quadrature::w1_3d; - double *Quadrature::w2_3d; - double *Quadrature::w3_3d; - double *Quadrature::w4_3d; - double *Quadrature::w5_3d; - double *Quadrature::w7_3d; - - - void Quadrature::initStaticQuadratures() + QuadratureFactory::QuadratureFactory() { TEST_EXIT(x0_1d == NULL)("static quadratures already initialized\n"); @@ -201,7 +124,7 @@ namespace AMDiS { w_0d = createAndInitArray(1, StdVol * 1.0); - Quadrature::quad_0d[0] = new Quadrature("0d", 0, 0, 1, Quadrature::x_0d, Quadrature::w_0d); + quad_0d[0].reset(new Quadrature("0d", 0, 0, 1, x_0d, w_0d)); /****************************************************************************/ /* 1d quadrature formulas using 2 barycentric coordinates */ @@ -386,26 +309,26 @@ namespace AMDiS { StdVol * 0.074725674575291, StdVol * 0.033335672154344); - Quadrature::quad_1d[0]= new Quadrature("1d-Gauss: P_1", 1, 1, 1, Quadrature::x0_1d, Quadrature::w0_1d); /* P_0 */ - Quadrature::quad_1d[1]= new Quadrature("1d-Gauss: P_1", 1, 1, 1, Quadrature::x0_1d, Quadrature::w0_1d); /* P_1 */ - Quadrature::quad_1d[2]= new Quadrature("1d-Gauss: P_3", 3, 1, 2, Quadrature::x1_1d, Quadrature::w1_1d); /* P_2 */ - Quadrature::quad_1d[3]= new Quadrature("1d-Gauss: P_3", 3, 1, 2, Quadrature::x1_1d, Quadrature::w1_1d); /* P_3 */ - Quadrature::quad_1d[4]= new Quadrature("1d-Gauss: P_5", 5, 1, 3, Quadrature::x2_1d, Quadrature::w2_1d); /* P_4 */ - Quadrature::quad_1d[5]= new Quadrature("1d-Gauss: P_5", 5, 1, 3, Quadrature::x2_1d, Quadrature::w2_1d); /* P_5 */ - Quadrature::quad_1d[6]= new Quadrature("1d-Gauss: P_7", 7, 1, 4, Quadrature::x3_1d, Quadrature::w3_1d); /* P_6 */ - Quadrature::quad_1d[7]= new Quadrature("1d-Gauss: P_7", 7, 1, 4, Quadrature::x3_1d, Quadrature::w3_1d); /* P_7 */ - Quadrature::quad_1d[8]= new Quadrature("1d-Gauss: P_9", 9, 1, 5, Quadrature::x4_1d, Quadrature::w4_1d); /* P_8 */ - Quadrature::quad_1d[9]= new Quadrature("1d-Gauss: P_9", 9, 1, 5, Quadrature::x4_1d, Quadrature::w4_1d); /* P_9 */ - Quadrature::quad_1d[10]= new Quadrature("1d-Gauss: P_11", 11, 1, 6, Quadrature::x5_1d, Quadrature::w5_1d); /* P_10 */ - Quadrature::quad_1d[11]= new Quadrature("1d-Gauss: P_11", 11, 1, 6, Quadrature::x5_1d, Quadrature::w5_1d); /* P_11 */ - Quadrature::quad_1d[12]= new Quadrature("1d-Gauss: P_13", 13, 1, 7, Quadrature::x6_1d, Quadrature::w6_1d); /* P_12 */ - Quadrature::quad_1d[13]= new Quadrature("1d-Gauss: P_13", 13, 1, 7, Quadrature::x6_1d, Quadrature::w6_1d); /* P_13 */ - Quadrature::quad_1d[14]= new Quadrature("1d-Gauss: P_15", 15, 1, 8, Quadrature::x7_1d, Quadrature::w7_1d); /* P_14 */ - Quadrature::quad_1d[15]= new Quadrature("1d-Gauss: P_15", 15, 1, 8, Quadrature::x7_1d, Quadrature::w7_1d); /* P_15 */ - Quadrature::quad_1d[16]= new Quadrature("1d-Gauss: P_17", 17, 1, 9, Quadrature::x8_1d, Quadrature::w8_1d); /* P_16 */ - Quadrature::quad_1d[17]= new Quadrature("1d-Gauss: P_17", 17, 1, 9, Quadrature::x8_1d, Quadrature::w8_1d); /* P_17 */ - Quadrature::quad_1d[18]= new Quadrature("1d-Gauss: P_19", 19, 1, 10, Quadrature::x9_1d, Quadrature::w9_1d); /* P_18 */ - Quadrature::quad_1d[19]= new Quadrature("1d-Gauss: P_19", 19, 1, 10, Quadrature::x9_1d, Quadrature::w9_1d);/* P_19 */ + quad_1d[0].reset(new Quadrature("1d-Gauss: P_1", 1, 1, 1, x0_1d, w0_1d)); /* P_0 */ + quad_1d[1].reset(new Quadrature("1d-Gauss: P_1", 1, 1, 1, x0_1d, w0_1d)); /* P_1 */ + quad_1d[2].reset(new Quadrature("1d-Gauss: P_3", 3, 1, 2, x1_1d, w1_1d)); /* P_2 */ + quad_1d[3].reset(new Quadrature("1d-Gauss: P_3", 3, 1, 2, x1_1d, w1_1d)); /* P_3 */ + quad_1d[4].reset(new Quadrature("1d-Gauss: P_5", 5, 1, 3, x2_1d, w2_1d)); /* P_4 */ + quad_1d[5].reset(new Quadrature("1d-Gauss: P_5", 5, 1, 3, x2_1d, w2_1d)); /* P_5 */ + quad_1d[6].reset(new Quadrature("1d-Gauss: P_7", 7, 1, 4, x3_1d, w3_1d)); /* P_6 */ + quad_1d[7].reset(new Quadrature("1d-Gauss: P_7", 7, 1, 4, x3_1d, w3_1d)); /* P_7 */ + quad_1d[8].reset(new Quadrature("1d-Gauss: P_9", 9, 1, 5, x4_1d, w4_1d)); /* P_8 */ + quad_1d[9].reset(new Quadrature("1d-Gauss: P_9", 9, 1, 5, x4_1d, w4_1d)); /* P_9 */ + quad_1d[10].reset(new Quadrature("1d-Gauss: P_11", 11, 1, 6, x5_1d, w5_1d)); /* P_10 */ + quad_1d[11].reset(new Quadrature("1d-Gauss: P_11", 11, 1, 6, x5_1d, w5_1d)); /* P_11 */ + quad_1d[12].reset(new Quadrature("1d-Gauss: P_13", 13, 1, 7, x6_1d, w6_1d)); /* P_12 */ + quad_1d[13].reset(new Quadrature("1d-Gauss: P_13", 13, 1, 7, x6_1d, w6_1d)); /* P_13 */ + quad_1d[14].reset(new Quadrature("1d-Gauss: P_15", 15, 1, 8, x7_1d, w7_1d)); /* P_14 */ + quad_1d[15].reset(new Quadrature("1d-Gauss: P_15", 15, 1, 8, x7_1d, w7_1d)); /* P_15 */ + quad_1d[16].reset(new Quadrature("1d-Gauss: P_17", 17, 1, 9, x8_1d, w8_1d)); /* P_16 */ + quad_1d[17].reset(new Quadrature("1d-Gauss: P_17", 17, 1, 9, x8_1d, w8_1d)); /* P_17 */ + quad_1d[18].reset(new Quadrature("1d-Gauss: P_19", 19, 1, 10, x9_1d, w9_1d)); /* P_18 */ + quad_1d[19].reset(new Quadrature("1d-Gauss: P_19", 19, 1, 10, x9_1d, w9_1d));/* P_19 */ #undef StdVol /****************************************************************************/ @@ -1131,24 +1054,24 @@ namespace AMDiS { #undef w14 #undef w15 - Quadrature::quad_2d[0] = new Quadrature("2d-P_1", 1, 2, N1, Quadrature::x1_2d, Quadrature::w1_2d); /* P 0 */ - Quadrature::quad_2d[1] = new Quadrature("2d-P_1", 1, 2, N1, Quadrature::x1_2d, Quadrature::w1_2d); /* P 1 */ - Quadrature::quad_2d[2] = new Quadrature("2d Stroud: P_2", 2, 2, N2, Quadrature::x2_2d, Quadrature::w2_2d); /* P 2 */ - Quadrature::quad_2d[3] = new Quadrature("2d Stroud: P_3", 3, 2, N3, Quadrature::x3_2d, Quadrature::w3_2d); /* P 3 */ - Quadrature::quad_2d[4] = new Quadrature("2d Dunavant: P_4", 4, 2, N4, Quadrature::x4_2d, Quadrature::w4_2d); /* P 4 */ - Quadrature::quad_2d[5] = new Quadrature("2d Dunavant: P_5", 5, 2, N5, Quadrature::x5_2d, Quadrature::w5_2d); /* P 5 */ - Quadrature::quad_2d[6] = new Quadrature("2d Gattermann: P_7", 7, 2, N7, Quadrature::x7_2d, Quadrature::w7_2d); /* P 6 */ - Quadrature::quad_2d[7] = new Quadrature("2d Gattermann: P_7", 7, 2, N7, Quadrature::x7_2d, Quadrature::w7_2d); /* P 7 */ - Quadrature::quad_2d[8] = new Quadrature("2d Dunavant: P_8", 8, 2, N8, Quadrature::x8_2d, Quadrature::w8_2d); /* P 8 */ - Quadrature::quad_2d[9] = new Quadrature("2d Dunavant: P_9", 9, 2, N9, Quadrature::x9_2d, Quadrature::w9_2d); /* P 9 */ - Quadrature::quad_2d[10] = new Quadrature("2d Dunavant: P_10", 10, 2, N10, Quadrature::x10_2d, Quadrature::w10_2d);/* P 10 */ - Quadrature::quad_2d[11] = new Quadrature("2d Dunavant: P_11", 11, 2, N11, Quadrature::x11_2d, Quadrature::w11_2d);/* P 11 */ - Quadrature::quad_2d[12] = new Quadrature("2d Dunavant: P_12", 12, 2, N12, Quadrature::x12_2d, Quadrature::w12_2d);/* P 12 */ - Quadrature::quad_2d[13] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 13 */ - Quadrature::quad_2d[14] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 14 */ - Quadrature::quad_2d[15] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 15 */ - Quadrature::quad_2d[16] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d);/* P 16 */ - Quadrature::quad_2d[17] = new Quadrature("2d Dunavant: P_17", 17, 2, N17, Quadrature::x17_2d, Quadrature::w17_2d); /* P 17 */ + quad_2d[0].reset(new Quadrature("2d-P_1", 1, 2, N1, x1_2d, w1_2d)); /* P 0 */ + quad_2d[1].reset(new Quadrature("2d-P_1", 1, 2, N1, x1_2d, w1_2d)); /* P 1 */ + quad_2d[2].reset(new Quadrature("2d Stroud: P_2", 2, 2, N2, x2_2d, w2_2d)); /* P 2 */ + quad_2d[3].reset(new Quadrature("2d Stroud: P_3", 3, 2, N3, x3_2d, w3_2d)); /* P 3 */ + quad_2d[4].reset(new Quadrature("2d Dunavant: P_4", 4, 2, N4, x4_2d, w4_2d)); /* P 4 */ + quad_2d[5].reset(new Quadrature("2d Dunavant: P_5", 5, 2, N5, x5_2d, w5_2d)); /* P 5 */ + quad_2d[6].reset(new Quadrature("2d Gattermann: P_7", 7, 2, N7, x7_2d, w7_2d)); /* P 6 */ + quad_2d[7].reset(new Quadrature("2d Gattermann: P_7", 7, 2, N7, x7_2d, w7_2d)); /* P 7 */ + quad_2d[8].reset(new Quadrature("2d Dunavant: P_8", 8, 2, N8, x8_2d, w8_2d)); /* P 8 */ + quad_2d[9].reset(new Quadrature("2d Dunavant: P_9", 9, 2, N9, x9_2d, w9_2d)); /* P 9 */ + quad_2d[10].reset(new Quadrature("2d Dunavant: P_10", 10, 2, N10, x10_2d, w10_2d));/* P 10 */ + quad_2d[11].reset(new Quadrature("2d Dunavant: P_11", 11, 2, N11, x11_2d, w11_2d));/* P 11 */ + quad_2d[12].reset(new Quadrature("2d Dunavant: P_12", 12, 2, N12, x12_2d, w12_2d));/* P 12 */ + quad_2d[13].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 13 */ + quad_2d[14].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 14 */ + quad_2d[15].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 15 */ + quad_2d[16].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d));/* P 16 */ + quad_2d[17].reset(new Quadrature("2d Dunavant: P_17", 17, 2, N17, x17_2d, w17_2d)); /* P 17 */ #undef StdVol @@ -1401,14 +1324,14 @@ namespace AMDiS { /* use that of degree (only on function evaluation also) */ /****************************************************************************/ - Quadrature::quad_3d[0] = new Quadrature("3d Stroud: P_1", 1, 3, 1, Quadrature::x1_3d, Quadrature::w1_3d); /* P_0 */ - Quadrature::quad_3d[1] = new Quadrature("3d Stroud: P_1", 1, 3, 1, Quadrature::x1_3d, Quadrature::w1_3d); /* P_1 */ - Quadrature::quad_3d[2] = new Quadrature("3d Stroud: P_2", 2, 3, 4, Quadrature::x2_3d, Quadrature::w2_3d); /* P_2 */ - Quadrature::quad_3d[3] = new Quadrature("3d Stroud: P_3", 3, 3, 8, Quadrature::x3_3d, Quadrature::w3_3d); /* P_3 */ - Quadrature::quad_3d[4] = new Quadrature("3d ???: P_5", 5, 3, 15, Quadrature::x5_3d, Quadrature::w5_3d); /* P_4 */ - Quadrature::quad_3d[5] = new Quadrature("3d ???: P_5", 5, 3, 15, Quadrature::x5_3d, Quadrature::w5_3d); /* P_5 */ - Quadrature::quad_3d[6] = new Quadrature("3d ???: P_7", 7, 3, 64, Quadrature::x7_3d, Quadrature::w7_3d); /* P_6 */ - Quadrature::quad_3d[7] = new Quadrature("3d ???: P_7", 7, 3, 64, Quadrature::x7_3d, Quadrature::w7_3d); /* P_7 */ + quad_3d[0].reset(new Quadrature("3d Stroud: P_1", 1, 3, 1, x1_3d, w1_3d)); /* P_0 */ + quad_3d[1].reset(new Quadrature("3d Stroud: P_1", 1, 3, 1, x1_3d, w1_3d)); /* P_1 */ + quad_3d[2].reset(new Quadrature("3d Stroud: P_2", 2, 3, 4, x2_3d, w2_3d)); /* P_2 */ + quad_3d[3].reset(new Quadrature("3d Stroud: P_3", 3, 3, 8, x3_3d, w3_3d)); /* P_3 */ + quad_3d[4].reset(new Quadrature("3d ???: P_5", 5, 3, 15, x5_3d, w5_3d)); /* P_4 */ + quad_3d[5].reset(new Quadrature("3d ???: P_5", 5, 3, 15, x5_3d, w5_3d)); /* P_5 */ + quad_3d[6].reset(new Quadrature("3d ???: P_7", 7, 3, 64, x7_3d, w7_3d)); /* P_6 */ + quad_3d[7].reset(new Quadrature("3d ???: P_7", 7, 3, 64, x7_3d, w7_3d)); /* P_7 */ #undef StdVol @@ -1417,16 +1340,17 @@ namespace AMDiS { /* integration in different dimensions */ /****************************************************************************/ - Quadrature::quad_nd[0] = Quadrature::quad_0d; - Quadrature::quad_nd[1] = Quadrature::quad_1d; - Quadrature::quad_nd[2] = Quadrature::quad_2d; - Quadrature::quad_nd[3] = Quadrature::quad_3d; +// quad_nd[0] = quad_0d; +// quad_nd[1] = quad_1d; +// quad_nd[2] = quad_2d; +// quad_nd[3] = quad_3d; } - + Quadrature* Quadrature::provideQuadrature(int dim_, int degree_) { FUNCNAME("Quadrature::provideQuadrature()"); + static QuadratureFactory factory{}; switch (dim_) { case 0: @@ -1445,17 +1369,7 @@ namespace AMDiS { ERROR_EXIT("invalid dim\n"); } - if (x0_1d == NULL) - initStaticQuadratures(); - - return (quad_nd[dim_][degree_]); - } - - - Quadrature::~Quadrature() - { - delete lambda; - delete [] w; + return factory.get(dim_, degree_); } @@ -1482,8 +1396,10 @@ namespace AMDiS { Flag init_flag) { FastQuadrature *quad_fast = NULL; - -// #pragma omp critical + + #ifdef _OPENMP + #pragma omp critical (provideFastQuad) + #endif { list::iterator fast = fastQuadList.begin(); for (; fast != fastQuadList.end(); fast++) diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h index d6a7419e96d3d5445bd96d6f78e52ce82bb56884..d2ec5c779d0042091c8f4bf90e3f65cbb777b611 100644 --- a/AMDiS/src/Quadrature.h +++ b/AMDiS/src/Quadrature.h @@ -25,8 +25,10 @@ #ifndef AMDIS_QUADRATURE_H #define AMDIS_QUADRATURE_H +#include #include #include +#include #include #include "AMDiS_fwd.h" #include "BasisFunction.h" @@ -34,7 +36,9 @@ #include "FixVec.h" namespace AMDiS { - + + class QuadratureFactory; + /** * \ingroup Assembler * @@ -50,6 +54,8 @@ namespace AMDiS { */ class Quadrature { + friend class QuadratureFactory; + protected: /// Avoids call of default constructor Quadrature(); @@ -71,7 +77,7 @@ namespace AMDiS { dim(dim_), n_points(n_points_), lambda(lambda_), - w(w_) + w(w_, w_+n_points_) {} public: @@ -79,7 +85,7 @@ namespace AMDiS { Quadrature(const Quadrature&); /// Destructor - ~Quadrature(); +// ~Quadrature(); /// Returns a Quadrature for dimension dim exact for degree degree. static Quadrature *provideQuadrature(int dim, int degree); @@ -97,7 +103,7 @@ namespace AMDiS { double integrateStdSimplex(AbstractFunction > *f); /// Returns \ref name - inline std::string getName() + inline std::string getName() const { return name; } @@ -115,9 +121,15 @@ namespace AMDiS { } /// Returns \ref w. - inline double* getWeight() const + inline double const* getWeight() const + { + return w.data(); + } + + /// Returns \ref w. + inline double* getWeight() { - return w; + return w.data(); } /// Returns \ref dim @@ -199,84 +211,114 @@ namespace AMDiS { VectorOfFixVecs > *lambda; /// Vector of quadrature weights - double *w; - + std::vector w; + }; + + + class QuadratureFactory + { + friend class Quadrature; + + private: + + QuadratureFactory(); + + public: + + Quadrature* get(const int dim, const int degree) const + { + switch (dim) { + case 0: + return quad_0d[degree].get(); break; + case 1: + return quad_1d[degree].get(); break; + case 2: + return quad_2d[degree].get(); break; + case 3: + return quad_3d[degree].get(); break; + default: + ERROR_EXIT("Unknown dimension!"); + return NULL; + } + } + protected: - /// Initialisation of all static Quadrature objects which will be returned - /// by \ref provideQuadrature() - static void initStaticQuadratures(); - /** \name static quadratures, used weights, and barycentric coords + /** \name (static) quadratures, used weights, and barycentric coords * \{ */ - static Quadrature **quad_nd[4]; - static Quadrature *quad_0d[1]; - static Quadrature *quad_1d[20]; - static Quadrature *quad_2d[18]; - static Quadrature *quad_3d[8]; - - static VectorOfFixVecs > *x_0d; - static double *w_0d; - - static VectorOfFixVecs > *x0_1d; - static VectorOfFixVecs > *x1_1d; - static VectorOfFixVecs > *x2_1d; - static VectorOfFixVecs > *x3_1d; - static VectorOfFixVecs > *x4_1d; - static VectorOfFixVecs > *x5_1d; - static VectorOfFixVecs > *x6_1d; - static VectorOfFixVecs > *x7_1d; - static VectorOfFixVecs > *x8_1d; - static VectorOfFixVecs > *x9_1d; - static double *w0_1d; - static double *w1_1d; - static double *w2_1d; - static double *w3_1d; - static double *w4_1d; - static double *w5_1d; - static double *w6_1d; - static double *w7_1d; - static double *w8_1d; - static double *w9_1d; - - static VectorOfFixVecs > *x1_2d; - static VectorOfFixVecs > *x2_2d; - static VectorOfFixVecs > *x3_2d; - static VectorOfFixVecs > *x4_2d; - static VectorOfFixVecs > *x5_2d; - static VectorOfFixVecs > *x7_2d; - static VectorOfFixVecs > *x8_2d; - static VectorOfFixVecs > *x9_2d; - static VectorOfFixVecs > *x10_2d; - static VectorOfFixVecs > *x11_2d; - static VectorOfFixVecs > *x12_2d; - static VectorOfFixVecs > *x17_2d; - static double *w1_2d; - static double *w2_2d; - static double *w3_2d; - static double *w4_2d; - static double *w5_2d; - static double *w7_2d; - static double *w8_2d; - static double *w9_2d; - static double *w10_2d; - static double *w11_2d; - static double *w12_2d; - static double *w17_2d; - - static VectorOfFixVecs > *x1_3d; - static VectorOfFixVecs > *x2_3d; - static VectorOfFixVecs > *x3_3d; - static VectorOfFixVecs > *x4_3d; - static VectorOfFixVecs > *x5_3d; - static VectorOfFixVecs > *x7_3d; - static double *w1_3d; - static double *w2_3d; - static double *w3_3d; - static double *w4_3d; - static double *w5_3d; - static double *w7_3d; - +// Quadrature **quad_nd[4]; + std::array,1> quad_0d; + std::array,20> quad_1d; + std::array,18> quad_2d; + std::array,8> quad_3d; +// Quadrature *quad_0d[1]; +// Quadrature *quad_1d[20]; +// Quadrature *quad_2d[18]; +// Quadrature *quad_3d[8]; + + VectorOfFixVecs > *x_0d; + double *w_0d; + + VectorOfFixVecs > *x0_1d; + VectorOfFixVecs > *x1_1d; + VectorOfFixVecs > *x2_1d; + VectorOfFixVecs > *x3_1d; + VectorOfFixVecs > *x4_1d; + VectorOfFixVecs > *x5_1d; + VectorOfFixVecs > *x6_1d; + VectorOfFixVecs > *x7_1d; + VectorOfFixVecs > *x8_1d; + VectorOfFixVecs > *x9_1d; + double *w0_1d; + double *w1_1d; + double *w2_1d; + double *w3_1d; + double *w4_1d; + double *w5_1d; + double *w6_1d; + double *w7_1d; + double *w8_1d; + double *w9_1d; + + VectorOfFixVecs > *x1_2d; + VectorOfFixVecs > *x2_2d; + VectorOfFixVecs > *x3_2d; + VectorOfFixVecs > *x4_2d; + VectorOfFixVecs > *x5_2d; + VectorOfFixVecs > *x7_2d; + VectorOfFixVecs > *x8_2d; + VectorOfFixVecs > *x9_2d; + VectorOfFixVecs > *x10_2d; + VectorOfFixVecs > *x11_2d; + VectorOfFixVecs > *x12_2d; + VectorOfFixVecs > *x17_2d; + double *w1_2d; + double *w2_2d; + double *w3_2d; + double *w4_2d; + double *w5_2d; + double *w7_2d; + double *w8_2d; + double *w9_2d; + double *w10_2d; + double *w11_2d; + double *w12_2d; + double *w17_2d; + + VectorOfFixVecs > *x1_3d; + VectorOfFixVecs > *x2_3d; + VectorOfFixVecs > *x3_3d; + VectorOfFixVecs > *x4_3d; + VectorOfFixVecs > *x5_3d; + VectorOfFixVecs > *x7_3d; + double *w1_3d; + double *w2_3d; + double *w3_3d; + double *w4_3d; + double *w5_3d; + double *w7_3d; + /** \} */ }; diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc index db720fcfad1a77baea9e631545575357e741c7f2..0c400e8969e55a709b60e92e8ccb7df4078b0ac8 100644 --- a/AMDiS/src/SecondOrderAssembler.cc +++ b/AMDiS/src/SecondOrderAssembler.cc @@ -187,17 +187,22 @@ namespace AMDiS { const int nPoints = quadrature->getNumPoints(); if (firstCall) { - dimVec.change_dim(dim + 1); - LALt.resize(nPoints); - for (int j = 0; j < nPoints; j++) - LALt[j].change_dim(dim + 1, dim + 1); - - psiFast = - updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI); - phiFast = - updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI); - - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + dimVec.change_dim(dim + 1); + LALt.resize(nPoints); + for (int j = 0; j < nPoints; j++) + LALt[j].change_dim(dim + 1, dim + 1); + + psiFast = + updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI); + phiFast = + updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI); + + firstCall = false; + } } for (int i = 0; i < nPoints; i++) diff --git a/AMDiS/src/StandardProblemIteration.cc b/AMDiS/src/StandardProblemIteration.cc index ba792ddaf2e0d4378c73c3539f08cb6a3f317556..ff173623e240cd5191f0eacd78a8c901eaa5a608 100644 --- a/AMDiS/src/StandardProblemIteration.cc +++ b/AMDiS/src/StandardProblemIteration.cc @@ -112,7 +112,7 @@ namespace AMDiS { } - std::string StandardProblemIteration::getName() + std::string StandardProblemIteration::getName() const { return problem->getName(); } diff --git a/AMDiS/src/StandardProblemIteration.h b/AMDiS/src/StandardProblemIteration.h index 95967c70a8cb426784ae533973c17c0ad07aad23..9648c2643b1b31dbe9b7b292b0d9d2930ac37ce2 100644 --- a/AMDiS/src/StandardProblemIteration.h +++ b/AMDiS/src/StandardProblemIteration.h @@ -51,7 +51,7 @@ namespace AMDiS { void endIteration(AdaptInfo *adaptInfo) override; /// Implementation of \ref ProblemIterationInterface::getNumProblems() - int getNumProblems() override + int getNumProblems() const override { return 1; } @@ -61,9 +61,15 @@ namespace AMDiS { { return problem; } + + /// Implementation of \ref ProblemIterationInterface::getProblem(int) + ProblemStatBase const* getProblem(int number = 0) const override + { + return problem; + } /// Returns the name of the problem. - std::string getName() override; + std::string getName() const override; /// Function that serializes the problem plus information about the iteration. void serialize(std::ostream &out) override; diff --git a/AMDiS/src/SystemVector.h b/AMDiS/src/SystemVector.h index 1ac301ea9054d3deed0ffdcbf7236c1a752229fb..afcc5323fc638412a2b2b72442560042fc3424bc 100644 --- a/AMDiS/src/SystemVector.h +++ b/AMDiS/src/SystemVector.h @@ -55,7 +55,8 @@ namespace AMDiS { SystemVector(const SystemVector& rhs) : name(rhs.getName()), componentSpaces(rhs.getFeSpaces()), - vectors(rhs.getNumVectors()) + vectors(rhs.getNumVectors()), + createVec(true) { for (unsigned int i = 0; i < vectors.size(); i++) vectors[i] = new DOFVector(*rhs.getDOFVector(i)); diff --git a/AMDiS/src/Timer.cc b/AMDiS/src/Timer.cc index 10598c806703101a9608d8b618323140d630264a..31718f6854c84e4f1742e83a13cedc6732214bdf 100644 --- a/AMDiS/src/Timer.cc +++ b/AMDiS/src/Timer.cc @@ -42,7 +42,7 @@ namespace AMDiS { #endif } - double Timer::elapsed() + double Timer::elapsed() const { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS return MPI::Wtime() - first_mpi; diff --git a/AMDiS/src/Timer.h b/AMDiS/src/Timer.h index 2620c73f67934a5d4203429fb7831123f0da9999..7b71a85b1fb423a960ffacc993c1cadc4d494e83 100644 --- a/AMDiS/src/Timer.h +++ b/AMDiS/src/Timer.h @@ -42,7 +42,7 @@ namespace AMDiS { void reset(); /// returns the elapsed time (from construction or last reset) to now in seconds - double elapsed(); + double elapsed() const; }; } diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc index 12fa11ede22ccd368eb7af1221e0f2ffb6bdc63f..db100f8d0a3fb836e7b00e910026a740963648b1 100644 --- a/AMDiS/src/ZeroOrderAssembler.cc +++ b/AMDiS/src/ZeroOrderAssembler.cc @@ -195,11 +195,16 @@ namespace AMDiS { int nPoints = quadrature->getNumPoints(); if (firstCall) { - const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); - psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); - basFcts = colFeSpace->getBasisFcts(); - phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); + psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); + basFcts = colFeSpace->getBasisFcts(); + phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); + firstCall = false; + } } if (num_rows(c) < static_cast(nPoints)) @@ -248,11 +253,16 @@ namespace AMDiS { int nPoints = quadrature->getNumPoints(); if (firstCall) { - const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); - psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); - basFcts = colFeSpace->getBasisFcts(); - phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + const BasisFunction *basFcts = rowFeSpace->getBasisFcts(); + psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); + basFcts = colFeSpace->getBasisFcts(); + phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); + firstCall = false; + } } ElementVector c(nPoints, 0.0); @@ -282,11 +292,16 @@ namespace AMDiS { ElementMatrix& mat) { if (firstCall) { - q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), - colFeSpace->getBasisFcts(), - quadrature); - q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), + colFeSpace->getBasisFcts(), + quadrature); + q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature); + firstCall = false; + } } ElementVector c(1, 0.0); @@ -318,11 +333,16 @@ namespace AMDiS { ElementVector& vec) { if (firstCall) { - q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), - colFeSpace->getBasisFcts(), - quadrature); - q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature); - firstCall = false; + #ifdef _OPENMP + #pragma omp critical + #endif + { + q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), + colFeSpace->getBasisFcts(), + quadrature); + q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature); + firstCall = false; + } } vector::iterator termIt; diff --git a/AMDiS/src/expressions/functorN_expr.hpp b/AMDiS/src/expressions/functorN_expr.hpp index 8692faefb745c8f3f9a586b677c16e4c5ffaf21a..e60ba4d2e1431e22df1fd4a80d83239f7f351605 100644 --- a/AMDiS/src/expressions/functorN_expr.hpp +++ b/AMDiS/src/expressions/functorN_expr.hpp @@ -201,8 +201,7 @@ namespace AMDiS { std::tuple term_tuple; - template - LazyOperatorTerms(Terms_... terms_) + LazyOperatorTerms(Terms const&... terms_) : term_tuple(terms_...) { } @@ -249,8 +248,7 @@ namespace AMDiS F f; ///< the functor - template - FunctionN(F const& f_, Terms_... terms_) + FunctionN(F const& f_, Terms const&... terms_) : super(terms_...), f(f_) {} // call f.getDegree() function @@ -404,7 +402,7 @@ namespace AMDiS template inline typename result_of::FunctionN::type - function_(F const& f, Terms... ts) + function_(F const& f, Terms const&... ts) { return expressions::FunctionN::to::type...> (f, traits::to_expr::to::get(ts)...); @@ -412,7 +410,7 @@ namespace AMDiS template inline typename result_of::FunctionN::type - func(F const& f, Terms... ts) + func(F const& f, Terms const&... ts) { return expressions::FunctionN::to::type...> (f, traits::to_expr::to::get(ts)...); @@ -420,7 +418,7 @@ namespace AMDiS template inline typename result_of::FunctionN::type - eval(F const& f, Term0 t0, Terms... ts) + eval(F const& f, Term0 t0, Terms const&... ts) { return expressions::FunctionN::to::type, typename traits::to_expr::to::type...> diff --git a/AMDiS/src/expressions/valueOf.hpp b/AMDiS/src/expressions/valueOf.hpp index cbf17a9a6b0979f5b689232d04157f688a05d542..b7d35738d9c4403fa071c6211006d561c2f7ff91 100644 --- a/AMDiS/src/expressions/valueOf.hpp +++ b/AMDiS/src/expressions/valueOf.hpp @@ -212,7 +212,7 @@ namespace AMDiS /// Expressions that extracts the vector-value of a Vector at QPs - template class Vector, typename T, typename Name> + template class Vector, typename T, typename Name> struct ValueOf*>, Name, typename boost::enable_if >::type>::type > : public LazyOperatorTermBase @@ -224,7 +224,12 @@ namespace AMDiS mutable mtl::dense_vector vec; mutable Vector > coeff; - ValueOf(Vector*>& vector) : vecDV(vector) + ValueOf(Vector*> const& vector) : vecDV(vector) + { + resize(coeff, num_rows(vecDV)); + } + + ValueOf(Vector*>&& vector) : vecDV(std::move(vector)) { resize(coeff, num_rows(vecDV)); } @@ -232,7 +237,7 @@ namespace AMDiS template inline void insertFeSpaces(List& feSpaces) const { - for (size_t i = 0; i < num_rows(vecDV); i++) + for (size_t i = 0; i < size_t(num_rows(vecDV)); i++) feSpaces.insert(at(vecDV, i)->getFeSpace()); } @@ -247,7 +252,7 @@ namespace AMDiS const BasisFunction *basisFct = NULL) { Vector > helper; resize(helper, num_rows(vecDV)); - for (size_t i = 0; i < num_rows(vecDV); i++) { + for (size_t i = 0; i < size_t(num_rows(vecDV)); i++) { if (ot && subAssembler) ot->getVectorAtQPs(at(vecDV, i), elInfo, subAssembler, quad, at(helper, i)); else if (quad) @@ -269,7 +274,7 @@ namespace AMDiS vec.change_dim(num_rows(at(helper, 0))); for (size_t iq = 0; iq < num_rows(at(helper, 0)); iq++) { value_type tmp; resize(tmp, num_rows(vecDV)); - for (size_t i = 0; i < num_rows(vecDV); i++) + for (size_t i = 0; i < size_t(num_rows(vecDV)); i++) at(tmp, i) = at(helper, i)[iq]; vec[iq] = tmp; } @@ -390,13 +395,13 @@ namespace AMDiS template class Matrix, typename T> typename boost::enable_if >::type, expressions::ValueOf*>, Name > >::type - valueOf(Matrix*> &mat) + valueOf(Matrix*> const& mat) { return expressions::ValueOf*>, Name >(mat); } template class Vector, typename T> typename boost::enable_if >::type, expressions::ValueOf*>, Name > >::type - valueOf(Vector*> &vector) + valueOf(Vector*> const& vector) { return expressions::ValueOf*>, Name >(vector); } @@ -412,13 +417,13 @@ namespace AMDiS template class Matrix, typename T> typename boost::enable_if >::type, expressions::ValueOf*>, _unknown > >::type - valueOf(Matrix*> &mat) + valueOf(Matrix*> const& mat) { return expressions::ValueOf*>, _unknown >(mat); } - template class Vector, typename T> + template class Vector, typename T> typename boost::enable_if >::type, expressions::ValueOf*>, _unknown > >::type - valueOf(Vector*> &vector) + valueOf(Vector*> const& vector) { return expressions::ValueOf*>, _unknown >(vector); } diff --git a/AMDiS/src/io/ElementFileWriter.cc b/AMDiS/src/io/ElementFileWriter.cc index 13ba796bc3fa88d9a93ac0858817c8b0b1f44c93..b0e99daff10dd24ed7e7acb93c2a5978f6033752 100644 --- a/AMDiS/src/io/ElementFileWriter.cc +++ b/AMDiS/src/io/ElementFileWriter.cc @@ -192,7 +192,7 @@ namespace AMDiS { namespace io { //boost::iostreams seems not to truncate the file ofstream swapfile(filename.c_str(), ios::out | ios::trunc); TEST_EXIT(swapfile.is_open()) - ("Cannot open file %s for writing!\n", name.c_str()); + ("Cannot open file %s for writing!\n", filename.c_str()); swapfile.close(); } file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc)); @@ -311,7 +311,7 @@ namespace AMDiS { namespace io { //boost::iostreams seems not to truncate the file ofstream swapfile(filename.c_str(), ios::out | ios::trunc); TEST_EXIT(swapfile.is_open()) - ("Cannot open file %s for writing!\n", name.c_str()); + ("Cannot open file %s for writing!\n", filename.c_str()); swapfile.close(); } file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc)); diff --git a/AMDiS/src/io/detail/Arh3Writer.h b/AMDiS/src/io/detail/Arh3Writer.h index 61986009adf2fb904f628d59e9d2f08e34298f0c..aa75baa947ce3984724730bb579fc930e80fb606 100644 --- a/AMDiS/src/io/detail/Arh3Writer.h +++ b/AMDiS/src/io/detail/Arh3Writer.h @@ -20,7 +20,7 @@ namespace AMDiS { namespace io { ZLIB = 1, // zlib compression BZIP2 = 2 // bzip2 compression } Value; - }; + } namespace Macroformat { @@ -29,7 +29,7 @@ namespace AMDiS { namespace io { PT_MACROFILE = 1, // pointer to macro file SELF = 2 // pointer to this file, at the end of this file } Value; - }; + } typedef enum{SI08, SI16, SI32, SI64, UI08, UI16, UI32, UI64, SF32, SF64} Valformat; diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index bc938a06b5c3691d578d433c4c39cdf725674a4d..42c70e72f5e0cd84432ff4d9ecbfec5fa5b3449f 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -151,7 +151,7 @@ namespace AMDiS { namespace Parallel { /// the program if it finds a non macro element in the mesh. void testForMacroMesh(); - inline std::string getName() + inline std::string getName() const { return name; } diff --git a/AMDiS/src/reinit/HL_SignedDistTraverse.cc b/AMDiS/src/reinit/HL_SignedDistTraverse.cc index 29bd4b8825871296017353924054b64e3f297d28..2bb3b49ee70265701c9a06c42c6bc7f0e692c42e 100644 --- a/AMDiS/src/reinit/HL_SignedDistTraverse.cc +++ b/AMDiS/src/reinit/HL_SignedDistTraverse.cc @@ -21,6 +21,7 @@ #include "HL_SignedDistTraverse.h" #include "VelocityExtFromVelocityField.h" +#include "VertexVector.h" namespace reinit { @@ -145,6 +146,11 @@ void HL_SignedDistTraverse::HL_updateIteration() update_DOF = sD_DOF; else update_DOF = sDOld_DOF; + + std::vector associations; + for (int b : boundaryTypes) { + associations.push_back( &(feSpace->getMesh()->getPeriodicAssociations(BoundaryType{b}, feSpace->getAdmin())) ); + } // ===== Iteration loop: proceed until tolerance is reached. ===== TraverseStack stack; @@ -159,8 +165,38 @@ void HL_SignedDistTraverse::HL_updateIteration() elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS); + + int dim = feSpace->getMesh()->getDim(); + GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim); + std::vector localIndices(feSpace->getBasisFcts()->getNumber()); + while (elInfo) { HL_elementUpdate(elInfo); + + // synch parallel DOFs + if (!boundaryTypes.empty()) + feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices); + + for (std::size_t i = 0; i < boundaryTypes.size(); ++i) { + BoundaryType boundaryType = boundaryTypes[i]; + VertexVector& association = *associations[i]; + + auto absmin = [](double v0, double v1) { return std::abs(v0) < std::abs(v1) ? v0 : v1; }; + + for (int side = 0; side <= dim; side++) { + if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) { + for (int vertex = 0; vertex < dim; vertex++) { + DegreeOfFreedom right = localIndices[vertex]; + DegreeOfFreedom left = association[right]; + + + double value = absmin((*sD_DOF)[left], (*sD_DOF)[right]); + (*sD_DOF)[left] = (*sD_DOF)[right] = value; + } + } + } + } + elInfo = stack.traverseNext(elInfo); } diff --git a/AMDiS/src/reinit/HL_SignedDistTraverse.h b/AMDiS/src/reinit/HL_SignedDistTraverse.h index c0ebce30dd389e19e7d4c88410340fb1f67fac75..fbd3979dd66a379720e3984f5d53d6f366c71bd7 100644 --- a/AMDiS/src/reinit/HL_SignedDistTraverse.h +++ b/AMDiS/src/reinit/HL_SignedDistTraverse.h @@ -86,6 +86,7 @@ public: Parameters::get(name + "->tolerance", tol); Parameters::get(name + "->maximal number of iteration steps", maxIt); Parameters::get(name + "->Gauss-Seidel iteration", GaussSeidelFlag); + Parameters::get(name + "->periodic boundaries", boundaryTypes); TEST_EXIT(tol > 0)("illegal tolerance !\n"); @@ -190,6 +191,8 @@ public: int setUpdate_Cntr; // ---> end: for test purposes + // boundary types for periodic mesh boundaries + std::vector boundaryTypes; }; } // end namespace reinit diff --git a/configure.sh b/configure.sh index b3e18f4400e3022fbd7cd06920b73fd97c5e59c2..af4b439bcd5f4903cf39fedd689efa1d00383879 100755 --- a/configure.sh +++ b/configure.sh @@ -15,7 +15,7 @@ ENABLE_COMPRESSION="OFF" ENABLE_PARALLEL="OFF" BOOST_PREFIX="${BOOST_ROOT}" DOWNLOAD_BOOST="0" -BOOST_VERSION="1.62.0" +BOOST_VERSION="1.64.0" PETSC_PREFIX="${PETSC_DIR}" DOWNLOAD_PETSC="0" PETSC_VERSION="3.5.4" @@ -196,11 +196,11 @@ if [ "${DOWNLOAD_BOOST}" -eq "1" ]; then BOOST_PREFIX=${SRC_DIR}/packages/boost/${BOOST_VERSION} cd /tmp/src/boost - ./bootstrap.sh --prefix=${BOOST_PREFIX} --with-libraries=system,iostreams,filesystem,program_options,date_time,test + ./bootstrap.sh --prefix=${BOOST_PREFIX} --with-libraries=system,iostreams,filesystem,program_options,date_time,test,thread,mpi,serialization if [ "${ENABLE_COMPRESSION}" = "ON" ]; then - ./b2 -s cxxflags="-std=c++11" --build-type=minimal variant=release install + ./b2 -s cxxflags="-std=c++14" --build-type=minimal variant=release install else - ./b2 -s NO_BZIP2=1 -s NO_ZLIB=1 cxxflags="-std=c++11" --build-type=minimal variant=release install + ./b2 -s NO_BZIP2=1 -s NO_ZLIB=1 cxxflags="-std=c++14" --build-type=minimal variant=release install fi rm -rf /tmp/src/boost fi diff --git a/extensions/RefinementExpression.h b/extensions/RefinementExpression.h index ccba17ad40353bbeda77c99cddff8fd40fb44b68..502caa3e508fd64794446f6fbad0eb1f217be291 100644 --- a/extensions/RefinementExpression.h +++ b/extensions/RefinementExpression.h @@ -37,8 +37,6 @@ public: RefinementExpression(ProblemStat& prob, int component = 0, bool doCalcMeshSize_ = false) : mesh(prob.getMesh(component)), - adaptInfo(NULL), - refineOperation(NULL), numRefinements0(15), doCalcMeshSize(doCalcMeshSize_), onlyRefine(false), @@ -48,13 +46,13 @@ public: coarseningManager = prob.getCoarseningManager(); refinementManager = prob.getRefinementManager(); numRefinements = numRefinements0; - refineOperation = new StandardRefineOperation; + + standardRefineOperation.reset(new StandardRefineOperation); + refineOperation = standardRefineOperation.get(); } RefinementExpression(Mesh *mesh_, bool doCalcMeshSize_ = false) : mesh(mesh_), - adaptInfo(NULL), - refineOperation(NULL), numRefinements0(15), doCalcMeshSize(doCalcMeshSize_), onlyRefine(false), @@ -81,8 +79,8 @@ public: } numRefinements = numRefinements0; - - refineOperation = new StandardRefineOperation; + standardRefineOperation.reset(new StandardRefineOperation); + refineOperation = standardRefineOperation.get(); } virtual ~RefinementExpression() @@ -93,9 +91,6 @@ public: delete coarseningManager; delete refinementManager; } - - if (refineOperation) - delete refineOperation; } void finalize() @@ -315,9 +310,6 @@ public: void setRefineOperation(AdaptInfo* adaptInfo_, StandardRefineOperation* refineOperation_) { - if (refineOperation) - delete refineOperation; - adaptInfo = adaptInfo_; refineOperation = refineOperation_; } @@ -338,8 +330,10 @@ protected: RefinementManager* refinementManager; CoarseningManager* coarseningManager; - AdaptInfo* adaptInfo; - StandardRefineOperation* refineOperation; + AdaptInfo* adaptInfo = nullptr; + + std::shared_ptr standardRefineOperation; + StandardRefineOperation* refineOperation = nullptr; int numRefinements; int numRefinements0; diff --git a/extensions/VectorOperations.h b/extensions/VectorOperations.h index 71dcc888345ade92a4238ad86b9bf5d62e31043b..5616be83c8a37c5567fe09727e68d10ada5a7fb4 100644 --- a/extensions/VectorOperations.h +++ b/extensions/VectorOperations.h @@ -112,7 +112,7 @@ namespace vector_operations { { typename ValueType::type value; nullify(value); - for (size_t i = 0; i < num_rows(b); ++i) + for (size_t i = 0; i < size_t(num_rows(b)); ++i) value += sqr(b[i]); return std::sqrt(value); } diff --git a/extensions/base_problems/BaseProblem.h b/extensions/base_problems/BaseProblem.h index 6bd710f83cb157b1e4ae60ca0b3dbb52ba186b0e..ab575d95b0ce76733201a486b330b428467d504f 100644 --- a/extensions/base_problems/BaseProblem.h +++ b/extensions/base_problems/BaseProblem.h @@ -126,17 +126,17 @@ public: } /// name of the baseBroblem - std::string getName() + std::string getName() const { return name; } - int getNumProblems() + int getNumProblems() const { return 1; } - int getNumComponents() + int getNumComponents() const { return prob->getNumComponents(); } diff --git a/extensions/base_problems/CouplingBaseProblem.h b/extensions/base_problems/CouplingBaseProblem.h index ab42513b561a9644b5e66b1f56a0eade427ce562..4ee75c209bb03f58b949e868509f3ce1b6a98362 100644 --- a/extensions/base_problems/CouplingBaseProblem.h +++ b/extensions/base_problems/CouplingBaseProblem.h @@ -108,7 +108,7 @@ public: return baseProblems[i]->getFeSpace(j); } - std::string getName() { return name; } + std::string getName() const { return name; } protected: std::vector baseProblems; diff --git a/extensions/base_problems/CouplingBaseProblem2.h b/extensions/base_problems/CouplingBaseProblem2.h index e48512ea6e4f08f53f9db8b96db73e278bc63d7f..e6d2404a1e592060a26b81460a57f51192a7bc3b 100644 --- a/extensions/base_problems/CouplingBaseProblem2.h +++ b/extensions/base_problems/CouplingBaseProblem2.h @@ -204,7 +204,7 @@ public: return _GET_(baseProblems).getFeSpace(j); } - std::string getName() { return name; } + std::string getName() const { return name; } template @@ -254,7 +254,7 @@ public: } // using CProblemStat::getNumProblems - virtual int getNumProblems() override + virtual int getNumProblems() const override { return CProblemStat::getNumProblems(); } diff --git a/extensions/base_problems/CouplingBaseProblem2_cxx11.h b/extensions/base_problems/CouplingBaseProblem2_cxx11.h index 8bb01a021b3093ccd1b71aff55d656563599f41a..d0e482294d026067723f105c867343978439b83f 100644 --- a/extensions/base_problems/CouplingBaseProblem2_cxx11.h +++ b/extensions/base_problems/CouplingBaseProblem2_cxx11.h @@ -103,27 +103,28 @@ namespace extensions { * * \brief Structur to couple BaseProblems of variouse types */ -template -class CouplingBaseProblem : public CouplingIterationInterface, - public CouplingTimeInterface, - public AMDiS::detail::CouplingProblemStat +template +class CouplingBaseProblem + : public CouplingIterationInterface + , public CouplingTimeInterface + , public AMDiS::detail::CouplingProblemStat { -public: typedef AMDiS::detail::CouplingProblemStat CProblemStat; typedef std::tuple BaseProblemsTupleType; - template - CouplingBaseProblem(std::string name_, BaseProblemTypes_&&... baseProblems_) - : ProblemIterationInterface(), // virtual base class constructor - ProblemTimeInterface(), // virtual base class constructor - CProblemStat(name_), - baseProblems(baseProblems_...), - name(name_) +public: + CouplingBaseProblem(std::string name_, BaseProblemTypes&... baseProblems_) + : ProblemIterationInterface() // virtual base class constructor + , ProblemTimeInterface() // virtual base class constructor + , CProblemStat(name_) + , baseProblems(baseProblems_...) + , name(name_) { dow = Global::getGeo(WORLD); + Parameters::get(name_ + "->dim", dim); } - ~CouplingBaseProblem() { } + virtual ~CouplingBaseProblem() { } /** * Add the problems to the iterationInterface, timeInterface and couplingProblemStat. @@ -144,13 +145,17 @@ public: // initialize all ProblemStat CProblemStat::initialize(initFlag, adoptProblem, adoptFlag); + } - // set mesh dimension - dim = CProblemStat::getMesh(0)->getDim(); + virtual void initData() + { + tools::FOR_EACH< detail::InitData >::loop(baseProblems); } - virtual void initData() {} - virtual void finalizeData() {} + virtual void finalizeData() + { + tools::FOR_EACH< detail::FinalizeData >::loop(baseProblems); + } /** * At first the initData method is called for all baseProblems, then @@ -159,16 +164,9 @@ public: **/ virtual void initTimeInterface() { - tools::FOR_EACH< detail::InitData >::loop(baseProblems); initData(); - - tools::FOR_EACH< detail::FillOperators >::loop(baseProblems); - fillCouplingOperators(); - - tools::FOR_EACH< detail::FillBoundaryConditions >::loop(baseProblems); - fillCouplingBoundaryConditions(); - - tools::FOR_EACH< detail::FinalizeData >::loop(baseProblems); + fillOperators(); + fillBoundaryConditions(); finalizeData(); } @@ -176,6 +174,17 @@ public: virtual void fillCouplingOperators() {} virtual void fillCouplingBoundaryConditions() {} + virtual void fillOperators() + { + tools::FOR_EACH< detail::FillOperators >::loop(baseProblems); + fillCouplingOperators(); + } + + virtual void fillBoundaryConditions() + { + tools::FOR_EACH< detail::FillBoundaryConditions >::loop(baseProblems); + fillCouplingBoundaryConditions(); + } /// get the j-th solution-vector of the i-th problem template @@ -198,7 +207,7 @@ public: return _GET_(baseProblems).getFeSpace(j); } - std::string getName() { return name; } + std::string getName() const { return name; } ProblemType *getProblem(std::string name_) @@ -246,7 +255,7 @@ public: } // using CProblemStat::getNumProblems - virtual int getNumProblems() override + virtual int getNumProblems() const override { return CProblemStat::getNumProblems(); } diff --git a/extensions/base_problems/deprecated/NavierStokes_Chorin.h b/extensions/base_problems/deprecated/NavierStokes_Chorin.h index f5bacff65d5341234d7382eafe0ab62ed24ee9a9..27ecb19d080694e849c6563f4efae9d8a5fe14ad 100644 --- a/extensions/base_problems/deprecated/NavierStokes_Chorin.h +++ b/extensions/base_problems/deprecated/NavierStokes_Chorin.h @@ -88,9 +88,9 @@ public: return predictorProb->getMesh(comp); } - std::string getName() { return name; }; + std::string getName() const { return name; }; - int getNumProblems() { return 3; }; + int getNumProblems() const { return 3; }; ProblemStatType *getProblem(int number= 0) { FUNCNAME("NavierStokes_Chorin::getProblem()"); diff --git a/tools/install_boost.sh b/tools/install_boost.sh index 36b56bfa012b52ec5bccdea16d9ff6e7f4a05b53..6a483cca78a4e316820990621c2988a8af817e9c 100755 --- a/tools/install_boost.sh +++ b/tools/install_boost.sh @@ -5,22 +5,23 @@ set -x ROOT=${PWD} -BOOST_VERSION="1.62.0" +BOOST_VERSION="1.64.0" BOOST_FILENAME="boost_${BOOST_VERSION//[.]/_}.tar.gz" -INSTALL_PREFIX=${ROOT}/install -mkdir -p ${INSTALL_PREFIX} +BOOST_PREFIX=${ROOT}/../packages/boost/${BOOST_VERSION} +mkdir -p ${BOOST_PREFIX} BUILD_DIR=/tmp/$USER/boost_build mkdir -p ${BUILD_DIR} +# download boost +curl -SL "http://netcologne.dl.sourceforge.net/project/boost/boost/${BOOST_VERSION}/${BOOST_FILENAME}" \ + | tar --strip-components=1 -xzC ${BUILD_DIR} +cd ${BUILD_DIR} + # install boost -if [ ! -d ${INSTALL_PREFIX}/boost ]; then - curl -SL "http://netcologne.dl.sourceforge.net/project/boost/boost/${BOOST_VERSION}/${BOOST_FILENAME}" \ - | tar --strip-components=1 -xzC ${BUILD_DIR} - cd ${BUILD_DIR} - ./bootstrap.sh --prefix=${INSTALL_PREFIX}/boost \ - --with-libraries=system,iostreams,filesystem,program_options,date_time,test - ./b2 -s NO_BZIP2=1 cxxflags="-std=c++11" --build-type=minimal variant=release -j 4 install - rm -rf ${BUILD_DIR} -fi +./bootstrap.sh --prefix=${BOOST_PREFIX} \ + --with-libraries=system,iostreams,filesystem,program_options,date_time,test,thread,mpi,serialization +echo "using mpi ;" >> project-config.jam +./b2 cxxflags="-std=c++14" --build-type=minimal variant=release install +rm -rf ${BUILD_DIR}