diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8047520391cb75d3f46acad5c87a720414f3c040..8f484b0393a938d9995e19666903f3fc8f20b254 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -22,12 +22,20 @@ if("${PUGS_SHORT_VERSION}" STREQUAL "")
message(FATAL_ERROR "Unable to compute short version from PUGS_VERSION=${PUGS_VERSION}")
endif()
-
-set(CMAKE_CONFIGURATION_TYPES "Release;Debug;Coverage" CACHE STRING INTERNAL FORCE )
-
# set project version as PUGS_SHORT_VERSION
project (Pugs VERSION ${PUGS_SHORT_VERSION})
+# -----------------------------------------------------
+# dynamic libraries
+
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+link_libraries("-rdynamic")
+set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+
+#------------------------------------------------------
+
+set(CMAKE_CONFIGURATION_TYPES "Release;Debug;Coverage" CACHE STRING INTERNAL FORCE )
+
#------------------------------------------------------
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
@@ -133,8 +141,16 @@ endif()
#------------------------------------------------------
# Search for ParMETIS
+find_package(ParMETIS)
+if(PARMETIS_LIBRARIES)
+ add_library(PkgConfig::ParMETIS STATIC IMPORTED)
+ set_property(TARGET PkgConfig::ParMETIS PROPERTY
+ IMPORTED_LOCATION "${PARMETIS_LIBRARIES}")
+
+ set(PARMETIS_TARGET PkgConfig::ParMETIS)
+endif()
+
if(${MPI_FOUND})
- find_package(ParMETIS)
if (NOT PARMETIS_LIBRARIES)
if(PUGS_ENABLE_MPI MATCHES "^AUTO$")
message(STATUS "MPI support deactivated: ParMETIS cannot be found!")
@@ -160,7 +176,16 @@ set(PUGS_ENABLE_PETSC AUTO CACHE STRING
if (PUGS_ENABLE_PETSC MATCHES "^(AUTO|ON)$")
if (MPI_FOUND)
# PETSc support is deactivated if MPI is not found
- pkg_check_modules(PETSC PETSc)
+ pkg_check_modules(PETSC IMPORTED_TARGET GLOBAL PETSc)
+
+ set_property(TARGET PkgConfig::PETSC PROPERTY
+ IMPORTED_LOCATION "${PETSC_LIBRARIES}"
+ )
+ set_property(TARGET PkgConfig::PETSC PROPERTY
+ INTERFACE_INCLUDE_DIRECTORIES "${PETSC_INCLUDE_DIRS}"
+ )
+
+ set(PETSC_TARGET PkgConfig::PETSC)
else()
message(STATUS "PETSc support is deactivated since pugs will not be build with MPI support")
set(PETSC_FOUND FALSE)
@@ -188,7 +213,16 @@ set(PUGS_ENABLE_SLEPC AUTO CACHE STRING
if (PUGS_ENABLE_SLEPC MATCHES "^(AUTO|ON)$")
if (PETSC_FOUND)
# SLEPc support is deactivated if PETSc is not found
- pkg_check_modules(SLEPC SLEPc)
+ pkg_check_modules(SLEPC IMPORTED_TARGET GLOBAL SLEPc)
+
+ set_property(TARGET PkgConfig::SLEPC PROPERTY
+ IMPORTED_LOCATION "${SLEPC_LIBRARIES}"
+ )
+ set_property(TARGET PkgConfig::SLEPC PROPERTY
+ INTERFACE_INCLUDE_DIRECTORIES "${SLEPC_INCLUDE_DIRS}"
+ )
+
+ set(SLEPC_TARGET PkgConfig::SLEPC)
else()
message(STATUS "SLEPc support is deactivated since pugs will not be build with PETSc support")
set(SLEPC_FOUND FALSE)
@@ -225,22 +259,37 @@ set(PUGS_ENABLE_HDF5 AUTO CACHE STRING
"Choose one of: AUTO ON OFF")
if (PUGS_ENABLE_HDF5 MATCHES "^(AUTO|ON)$")
- # May be risky. (make to show pugs build options)
- set(HDF5_PREFER_PARALLEL TRUE)
- find_package(HDF5)
- if (HDF5_FOUND)
- # HighFive
- set(HIGHFIVE_USE_BOOST OFF) # no Boost
- set(HIGHFIVE_BUILD_DOCS OFF) # no doc
- set(HIGHFIVE_UNIT_TESTS OFF) # no unit tests
- set(HIGHFIVE_UNIT_TESTS OFF) # no unit tests
- set(HIGHFIVE_EXAMPLES OFF) # no examples
- add_subdirectory(${PUGS_SOURCE_DIR}/packages/HighFive/)
- set(HIGHFIVE_TARGET HighFive)
+ if (MPI_FOUND)
+ # May be risky. (make to show pugs build options)
+ set(HDF5_PREFER_PARALLEL TRUE)
+ find_package(HDF5)
+ if (HDF5_FOUND)
+ # HighFive
+ set(HIGHFIVE_USE_BOOST OFF) # no Boost
+ set(HIGHFIVE_BUILD_DOCS OFF) # no doc
+ set(HIGHFIVE_UNIT_TESTS OFF) # no unit tests
+ set(HIGHFIVE_UNIT_TESTS OFF) # no unit tests
+ set(HIGHFIVE_EXAMPLES OFF) # no examples
+ add_subdirectory(${PUGS_SOURCE_DIR}/packages/HighFive/)
+ set(HIGHFIVE_TARGET HighFive::HighFive)
+
+ add_library(PkgConfig::HDF5 STATIC IMPORTED)
+ set_property(TARGET PkgConfig::HDF5 PROPERTY
+ IMPORTED_LOCATION "${HDF5_LIBRARIES}")
+
+ set_property(TARGET PkgConfig::HDF5 PROPERTY
+ INTERFACE_INCLUDE_DIRECTORIES "${HDF5_INCLUDE_DIRS}")
+
+ set(HDF5_TARGET PkgConfig::HDF5)
+ endif()
+ else()
+ message(STATUS "HDF5 support is deactivated since pugs will not be build with MPI support")
+ set(HDF5_FOUND FALSE)
endif()
set(PUGS_HAS_HDF5 ${HDF5_FOUND})
else()
unset(HIGHFIVE_TARGET)
+ unset(HDF5_TARGET PkgConfig::HDF5)
unset(PUGS_HAS_HDF5)
endif()
@@ -251,7 +300,14 @@ find_package(Slurm)
set(PUGS_HAS_SLURM ${SLURM_FOUND})
-if (${SLURM_FOUND})
+if (SLURM_FOUND)
+ add_library(PkgConfig::Slurm STATIC IMPORTED)
+ set_property(TARGET PkgConfig::Slurm PROPERTY
+ IMPORTED_LOCATION "${SLURM_LIBRARY}")
+ set_property(TARGET PkgConfig::Slurm PROPERTY
+ INTERFACE_INCLUDE_DIRECTORIES "${SLURM_INCLUDE_DIR}")
+
+ set(SLURM_TARGET PkgConfig::Slurm)
include_directories(SYSTEM "${SLURM_INCLUDE_DIR}")
else()
set(SLURM_LIBRARY "")
@@ -364,6 +420,10 @@ endif()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${PUGS_CXX_FLAGS}")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${PUGS_CXX_FLAGS}")
+
+message(STATUS "CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}")
+
+
# Add debug mode for Standard C++ library (not for AppleClang since it is broken)
if (NOT "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -D_GLIBCXX_DEBUG -D_LIBCPP_DEBUG=1")
@@ -372,12 +432,14 @@ endif()
#------------------------------------------------------
# Rang (colors? Useless thus necessary!)
+add_subdirectory(${PUGS_SOURCE_DIR}/packages/rang/)
include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/rang/include")
# CLI11
include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/CLI11/include")
# PEGTL
+add_subdirectory(${PUGS_SOURCE_DIR}/packages/PEGTL/)
include_directories(SYSTEM "${PUGS_SOURCE_DIR}/packages/PEGTL/include/tao")
# Pugs src
@@ -400,7 +462,7 @@ if(${PUGS_HAS_MPI})
endif()
set(MPIEXEC_NUMPROC 3)
set(MPIEXEC_PATH_FLAG --path)
- set(MPI_LAUNCHER ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_NUMPROC} ${MPIEXEC_PATH_FLAG} ${PUGS_BINARY_DIR} ${MPIEXEC_OPTION_FLAGS})
+ set(MPI_LAUNCHER ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_NUMPROC} ${MPIEXEC_OPTION_FLAGS})
endif()
add_custom_target(all_unit_tests
@@ -427,11 +489,11 @@ endif()
add_custom_target(run_mpi_unit_tests
COMMAND ${MPI_LAUNCHER} "${PUGS_BINARY_DIR}/mpi_unit_tests"
+ WORKING_DIRECTORY ${PUGS_BINARY_DIR}
DEPENDS run_unit_tests
COMMENT ${RUN_MPI_UNIT_TESTS_COMMENT}
)
-
add_custom_target(run_unit_tests
COMMAND "${PUGS_BINARY_DIR}/unit_tests"
DEPENDS all_unit_tests
@@ -594,10 +656,6 @@ if("${CMAKE_BUILD_TYPE}" STREQUAL "Coverage")
endif()
endif()
-# -----------------------------------------------------
-
-link_libraries("-rdynamic")
-
# ------------------- Source files --------------------
# Pugs binary
add_executable(
@@ -625,17 +683,20 @@ target_link_libraries(
PugsLanguageUtils
PugsCheckpointing
Kokkos::kokkos
- ${PETSC_LIBRARIES}
- ${SLEPC_LIBRARIES}
- ${PARMETIS_LIBRARIES}
+ ${PETSC_TARGET}
+ ${SLEPC_TARGET}
+ ${PARMETIS_TARGET}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
${KOKKOS_CXX_FLAGS}
${OPENMP_LINK_FLAGS}
${PUGS_STD_LINK_FLAGS}
${HIGHFIVE_TARGET}
- ${SLURM_LIBRARY}
+ ${SLURM_TARGET}
stdc++fs
- )
+)
+
+target_include_directories(pugs PUBLIC ${PETSC_INCLUDE_DIRS})
+target_include_directories(pugs PUBLIC ${HDF5_INCLUDE_DIRS})
# Checkpoint management tool
add_executable(
@@ -660,19 +721,18 @@ target_link_libraries(
PugsMesh
PugsOutput
Kokkos::kokkos
- ${PETSC_LIBRARIES}
- ${SLEPC_LIBRARIES}
- ${PARMETIS_LIBRARIES}
+ ${PETSC_TARGET}
+ ${SLEPC_TARGET}
+ ${PARMETIS_TARGET}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
${KOKKOS_CXX_FLAGS}
${OPENMP_LINK_FLAGS}
${PUGS_STD_LINK_FLAGS}
${HIGHFIVE_TARGET}
- ${SLURM_LIBRARY}
+ ${SLURM_TARGET}
stdc++fs
)
-
# -------------------- Documentation --------------------
include(PugsDoc)
@@ -685,7 +745,6 @@ include(PugsDoxygen)
install(TARGETS
pugs
pugs_checkpoint
- PugsMesh
PugsAlgebra
PugsAnalysis
PugsCheckpointing
@@ -694,18 +753,33 @@ install(TARGETS
PugsLanguage
PugsLanguageAST
PugsLanguageModules
- PugsLanguageAlgorithms
+ PugsLanguageUtils
PugsMesh
- PugsAlgebra
- PugsScheme
- PugsUtils
PugsOutput
- PugsLanguageUtils
+ PugsScheme
kokkos
+ Catch2
+
+ EXPORT "${PROJECT_NAME}Targets"
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
- ARCHIVE DESTINATION lib)
+ ARCHIVE DESTINATION lib
+)
+
+include(CMakePackageConfigHelpers)
+
+write_basic_package_version_file(
+ PugsConfigVersion.cmake
+ VERSION ${PACKAGE_VERSION}
+ COMPATIBILITY AnyNewerVersion
+)
+
+install(EXPORT PugsTargets
+ FILE PugsTargets.cmake
+ NAMESPACE Pugs::
+ DESTINATION lib/cmake/pugs
+)
# ------------------- Build options -------------------
message("")
@@ -840,3 +914,24 @@ endif()
message("================================")
message("")
+
+configure_file(
+ ${PUGS_SOURCE_DIR}/cmake/PugsCompileFlags.cmake.in
+ ${PUGS_BINARY_DIR}/cmake/PugsCompileFlags.cmake
+ @ONLY
+)
+
+install(
+ FILES ${PUGS_BINARY_DIR}/cmake/PugsCompileFlags.cmake
+ DESTINATION lib/cmake/pugs
+)
+
+# Ugly patch to install user headers for Catch2 (for plugins)
+install(
+ DIRECTORY
+ "${PUGS_BINARY_DIR}/packages/Catch2/generated-includes/catch2" # Also install the generated header
+ DESTINATION
+ "include"
+ FILES_MATCHING
+ PATTERN "*.hpp"
+ )
diff --git a/README.md b/README.md
index 36e460a851a53d40b4ab442460594091fa22b20f..c859bc694043fc05a83ed47dae3f515eca0bbb42 100644
--- a/README.md
+++ b/README.md
@@ -11,10 +11,10 @@ that are assembled together by a user friendly language.
## Requirements
For the most basic build, `pugs` only requires
-- a C++-17 compiler.
+- a C++-20 compiler.
- g++-10 or higher versions
- - clang-10 or higher versions
-- `CMake` (at least 3.16)
+ - clang-11 or higher versions
+- `CMake` (at least 3.19)
> **Warning**:<br>
> Building `pugs` in its source directory is **forbidden**. Trying to
diff --git a/cmake/FindParMETIS.cmake b/cmake/FindParMETIS.cmake
index d9b91d33d92617f3282d4eac2184549b62f62418..5fb1bb961c53897ab8f9cd33a2bf6bedfa41a783 100644
--- a/cmake/FindParMETIS.cmake
+++ b/cmake/FindParMETIS.cmake
@@ -1,6 +1,5 @@
# Looking for ParMETIS
-
find_path(PARMETIS_INCLUDE_DIR parmetis.h
PATH_SUFFIX include parmetis $ENV{PARMETIS_INCDIR})
@@ -17,7 +16,8 @@ if(EXISTS "${PARMETIS_INCLUDE_DIR}/parmetis.h")
find_path(METIS_INCLUDE_DIR metis.h $ENV{METIS_INCDIR})
if(EXISTS "${METIS_INCLUDE_DIR}/metis.h")
message(STATUS "Found metis.h in ${METIS_INCLUDE_DIR}")
- set(PARMETIS_LIBRARIES ${LIB_PARMETIS} ${LIB_METIS})
+ set(PARMETIS_LIBRARIES ${LIB_PARMETIS})
+ set(METIS_LIBRARIES ${LIB_METIS})
message(STATUS "Found parmetis/metis libraries ${PARMETIS_LIBRARIES}")
else()
message(WARNING "** Could not find metis.h.\n** Is METIS_INCDIR correctly set (Actual: \"$ENV{METIS_INCDIR}\")?")
@@ -27,3 +27,4 @@ else()
endif()
mark_as_advanced(PARMETIS_INCLUDE_DIR PARMETIS_LIBRARIES)
+mark_as_advanced(METIS_INCLUDE_DIR METIS_LIBRARIES)
diff --git a/cmake/PugsCompileFlags.cmake.in b/cmake/PugsCompileFlags.cmake.in
new file mode 100644
index 0000000000000000000000000000000000000000..83ecd3e6036a82c3e8696d2b75de2f35d724051c
--- /dev/null
+++ b/cmake/PugsCompileFlags.cmake.in
@@ -0,0 +1,10 @@
+@PACKAGE_INIT@
+
+set(PUGS_CMAKE_CXX_FLAGS "@CMAKE_CXX_FLAGS@")
+set(PUGS_CMAKE_CXX_STANDARD "@CMAKE_CXX_STANDARD@")
+
+set(PUGS_CMAKE_BUILD_TYPE "@CMAKE_BUILD_TYPE@")
+set(PUGS_CMAKE_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
+set(PUGS_CMAKE_C_COMPILER "@CMAKE_C_COMPILER@")
+
+set(PUGS_HAS_MPI "@PUGS_HAS_MPI@")
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index eebccac12c3c1c7b80b6b1e13ff81d0e802a8def..30d920bb0c111049ae655f67e3951424134dbb01 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,6 +3,21 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR})
#------------------------------------------------------
+install(
+ DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/"
+ DESTINATION "include"
+ FILES_MATCHING
+ PATTERN "*.hpp"
+)
+
+install(
+ DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/"
+ DESTINATION "include"
+ FILES_MATCHING
+ PATTERN "*.hpp"
+ PATTERN "CMakeFiles" EXCLUDE
+)
+
# Pugs utils
add_subdirectory(utils)
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index 74f40290a2f9250cc2296ec62358b44c366ffc85..d215b09b50d9933f84aab02855388bb4bd25dced 100644
--- a/src/algebra/CMakeLists.txt
+++ b/src/algebra/CMakeLists.txt
@@ -9,7 +9,10 @@ add_library(
target_link_libraries(
PugsAlgebra
- ${PETSC_LIBRARIES}
- ${SLEPC_LIBRARIES}
+ ${PETSC_TARGET}
+ ${SLEPC_TARGET}
${HIGHFIVE_TARGET}
)
+
+target_include_directories(PugsAlgebra PUBLIC ${PETSC_INCLUDE_DIRS})
+target_include_directories(PugsAlgebra PUBLIC ${SLEPC_INCLUDE_DIRS})
diff --git a/src/algebra/ShrinkVectorView.hpp b/src/algebra/ShrinkVectorView.hpp
index cdf786c5a30970b747b0f540642d5637059198b8..83aed7cb2295ee4d20a044f7ba1039f34e7b5ffa 100644
--- a/src/algebra/ShrinkVectorView.hpp
+++ b/src/algebra/ShrinkVectorView.hpp
@@ -23,10 +23,10 @@ class ShrinkVectorView
friend std::ostream&
operator<<(std::ostream& os, const ShrinkVectorView& x)
{
- if (x.size() > 0) {
+ if (x.dimension() > 0) {
os << 0 << ':' << NaNHelper(x[0]);
}
- for (size_t i = 1; i < x.size(); ++i) {
+ for (size_t i = 1; i < x.dimension(); ++i) {
os << ' ' << i << ':' << NaNHelper(x[i]);
}
return os;
diff --git a/src/language/ast/ASTModulesImporter.cpp b/src/language/ast/ASTModulesImporter.cpp
index 4cae0a921c139b3ae68a3decb0ce4f83876c4701..e1dc0043e689b506776540bb57e4851dc10883ae 100644
--- a/src/language/ast/ASTModulesImporter.cpp
+++ b/src/language/ast/ASTModulesImporter.cpp
@@ -26,10 +26,10 @@ ASTModulesImporter::_importModule(ASTNode& import_node)
std::cout << " * importing '" << rang::fgB::green << module_name << rang::style::reset << "' module\n";
}
- m_module_repository.populateSymbolTable(module_name_node, m_symbol_table);
- m_module_repository.registerOperators(module_name);
+ ModuleRepository::getInstance().populateSymbolTable(module_name_node, m_symbol_table);
+ ModuleRepository::getInstance().registerOperators(module_name);
if (CheckpointResumeRepository::isCreated()) {
- m_module_repository.registerCheckpointResume(module_name);
+ ModuleRepository::getInstance().registerCheckpointResume(module_name);
}
}
@@ -49,7 +49,7 @@ ASTModulesImporter::ASTModulesImporter(ASTNode& root_node) : m_symbol_table{*roo
{
Assert(root_node.is_root());
OperatorRepository::instance().reset();
- m_module_repository.populateMandatoryData(root_node, m_symbol_table);
+ ModuleRepository::getInstance().populateMandatoryData(root_node, m_symbol_table);
this->_importAllModules(root_node);
diff --git a/src/language/ast/ASTModulesImporter.hpp b/src/language/ast/ASTModulesImporter.hpp
index 566155c4a19ecd00f0b049a29030c0036f212d8a..15ba8a7dc6cecbbb39bfb11c4d202d39200061a4 100644
--- a/src/language/ast/ASTModulesImporter.hpp
+++ b/src/language/ast/ASTModulesImporter.hpp
@@ -14,8 +14,6 @@ class ASTModulesImporter
std::set<std::string> m_imported_modules;
SymbolTable& m_symbol_table;
- ModuleRepository m_module_repository;
-
void _importModule(ASTNode& import_node);
void _importAllModules(ASTNode& node);
@@ -23,7 +21,7 @@ class ASTModulesImporter
const ModuleRepository&
moduleRepository() const
{
- return m_module_repository;
+ return ModuleRepository::getInstance();
}
ASTModulesImporter(ASTNode& root_node);
diff --git a/src/language/modules/BinaryOperatorRegisterForVh.cpp b/src/language/modules/BinaryOperatorRegisterForVh.cpp
index 3eb83537323d91bbf6266cbefebf4d004d74c2be..d1c20cb78718568a3535e30f0b5c9ce255be1db6 100644
--- a/src/language/modules/BinaryOperatorRegisterForVh.cpp
+++ b/src/language/modules/BinaryOperatorRegisterForVh.cpp
@@ -1,6 +1,6 @@
#include <language/modules/BinaryOperatorRegisterForVh.hpp>
-#include <language/modules/SchemeModule.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/DataHandler.hpp>
#include <language/utils/DataVariant.hpp>
diff --git a/src/language/modules/CoreModule.cpp b/src/language/modules/CoreModule.cpp
index d7b331a1dc606a83cf77dfe405b184b9ff143308..47e0183256413b7579165b26364dff68f0614a74 100644
--- a/src/language/modules/CoreModule.cpp
+++ b/src/language/modules/CoreModule.cpp
@@ -42,7 +42,7 @@
#include <utils/checkpointing/ResumingManager.hpp>
#include <utils/checkpointing/WriteOStream.hpp>
-#include <random>
+#include <language/modules/ModuleRepository.hpp>
CoreModule::CoreModule() : BuiltinModule(true)
{
@@ -219,3 +219,5 @@ CoreModule::registerCheckpointResume() const
#endif // PUGS_HAS_HDF5
}
+
+ModuleRepository::Subscribe<CoreModule> core_module;
diff --git a/src/language/modules/DevUtilsModule.cpp b/src/language/modules/DevUtilsModule.cpp
index 361b4357eda0ebbfe4350cc6c64bbe72cfaa0829..4595ea6c9d9a45431b303c00b4de6defdda38b7b 100644
--- a/src/language/modules/DevUtilsModule.cpp
+++ b/src/language/modules/DevUtilsModule.cpp
@@ -1,5 +1,7 @@
#include <language/modules/DevUtilsModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+
#include <dev/ParallelChecker.hpp>
#include <language/modules/MeshModuleTypes.hpp>
#include <language/modules/SchemeModuleTypes.hpp>
@@ -126,3 +128,5 @@ DevUtilsModule::registerOperators() const
void
DevUtilsModule::registerCheckpointResume() const
{}
+
+ModuleRepository::Subscribe<DevUtilsModule> dev_utils_module;
diff --git a/src/language/modules/LinearSolverModule.cpp b/src/language/modules/LinearSolverModule.cpp
index aef6470d2455909f5ad2d307aa4e5cb859b8bcbb..d3e2de0faeb52c213f372e3a751f075ea2d735d7 100644
--- a/src/language/modules/LinearSolverModule.cpp
+++ b/src/language/modules/LinearSolverModule.cpp
@@ -1,5 +1,7 @@
#include <language/modules/LinearSolverModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+
#include <algebra/LinearSolver.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
@@ -98,3 +100,5 @@ LinearSolverModule::registerOperators() const
void
LinearSolverModule::registerCheckpointResume() const
{}
+
+ModuleRepository::Subscribe<LinearSolverModule> linear_solver_module;
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index b31fe47d0eba87ff2e7f0823a530bf170e410d6c..d9a30e86996a9cd6b99af63d64a8ab33264c7fe3 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -1,6 +1,7 @@
#include <language/modules/MathFunctionRegisterForVh.hpp>
#include <language/modules/SchemeModule.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index 19a35a1529b1b65310c3c51d69fe4c5d02dc4d8a..1caea452428125e7a1c16aa88cb8278a5685fb35 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -1,5 +1,6 @@
#include <language/modules/MathModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <cmath>
@@ -136,3 +137,5 @@ MathModule::registerOperators() const
void
MathModule::registerCheckpointResume() const
{}
+
+ModuleRepository::Subscribe<MathModule> math_module;
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index c76260ce58c75c130e4638054fc394458fedbf63..c060da50d6da17915937edda6f5403a9e604355e 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -1,5 +1,7 @@
#include <language/modules/MeshModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+
#include <algebra/TinyVector.hpp>
#include <language/node_processor/ExecutionPolicy.hpp>
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
@@ -160,27 +162,28 @@ MeshModule::MeshModule()
));
- this
- ->_addBuiltinFunction("interpolate",
- std::function(
+ this->_addBuiltinFunction("interpolate",
+ std::function(
- [](std::shared_ptr<const MeshVariant> mesh_v, std::shared_ptr<const ItemType> item_type,
- const FunctionSymbolId& function_id) -> std::shared_ptr<const ItemValueVariant> {
- return ItemValueVariantFunctionInterpoler{mesh_v, *item_type, function_id}.interpolate();
- }
+ [](std::shared_ptr<const MeshVariant> mesh_v, std::shared_ptr<const ItemType> item_type,
+ const FunctionSymbolId& function_id) -> std::shared_ptr<const ItemValueVariant> {
+ return ItemValueVariantFunctionInterpoler{mesh_v, *item_type, function_id}
+ .interpolate();
+ }
- ));
+ ));
- this->_addBuiltinFunction(
- "interpolate_array",
- std::function(
+ this->_addBuiltinFunction("interpolate_array",
+ std::function(
- [](std::shared_ptr<const MeshVariant> mesh_v, std::shared_ptr<const ItemType> item_type,
- const std::vector<FunctionSymbolId>& function_id_list) -> std::shared_ptr<const ItemArrayVariant> {
- return ItemArrayVariantFunctionInterpoler{mesh_v, *item_type, function_id_list}.interpolate();
- }
+ [](std::shared_ptr<const MeshVariant> mesh_v, std::shared_ptr<const ItemType> item_type,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ -> std::shared_ptr<const ItemArrayVariant> {
+ return ItemArrayVariantFunctionInterpoler{mesh_v, *item_type, function_id_list}
+ .interpolate();
+ }
- ));
+ ));
this->_addBuiltinFunction("transform",
std::function(
@@ -442,3 +445,5 @@ MeshModule::registerCheckpointResume() const
#endif // PUGS_HAS_HDF5
}
+
+ModuleRepository::Subscribe<MeshModule> mesh_module;
diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp
index 59689fb6fee512e45cd76c7c9d3cd5c11d543ee8..c9c487b67b5fdbf63d2b307baf1497519666a5d8 100644
--- a/src/language/modules/ModuleRepository.cpp
+++ b/src/language/modules/ModuleRepository.cpp
@@ -1,14 +1,6 @@
#include <language/modules/ModuleRepository.hpp>
#include <language/ast/ASTNode.hpp>
-#include <language/modules/CoreModule.hpp>
-#include <language/modules/DevUtilsModule.hpp>
-#include <language/modules/LinearSolverModule.hpp>
-#include <language/modules/MathModule.hpp>
-#include <language/modules/MeshModule.hpp>
-#include <language/modules/SchemeModule.hpp>
-#include <language/modules/SocketModule.hpp>
-#include <language/modules/WriterModule.hpp>
#include <language/utils/BasicAffectationRegistrerFor.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/CheckpointResumeRepository.hpp>
@@ -19,9 +11,29 @@
#include <utils/PugsAssert.hpp>
-#include <algorithm>
#include <set>
+ModuleRepository* ModuleRepository::m_instance = nullptr;
+
+ModuleRepository&
+ModuleRepository::getInstance()
+{
+ if (m_instance == nullptr) {
+ m_instance = new ModuleRepository;
+ }
+
+ return *m_instance;
+}
+
+void
+ModuleRepository::destroy()
+{
+ if (m_instance != nullptr) {
+ delete m_instance;
+ m_instance = nullptr;
+ }
+}
+
void
ModuleRepository::_subscribe(std::unique_ptr<IModule> m)
{
@@ -52,18 +64,6 @@ ModuleRepository::_subscribe(std::unique_ptr<IModule> m)
Assert(success, "module has already been subscribed");
}
-ModuleRepository::ModuleRepository()
-{
- this->_subscribe(std::make_unique<CoreModule>());
- this->_subscribe(std::make_unique<LinearSolverModule>());
- this->_subscribe(std::make_unique<MathModule>());
- this->_subscribe(std::make_unique<MeshModule>());
- this->_subscribe(std::make_unique<SchemeModule>());
- this->_subscribe(std::make_unique<SocketModule>());
- this->_subscribe(std::make_unique<DevUtilsModule>());
- this->_subscribe(std::make_unique<WriterModule>());
-}
-
template <typename NameEmbedderMapT, typename EmbedderTableT>
void
ModuleRepository::_populateEmbedderTableT(const ASTNode& module_node,
diff --git a/src/language/modules/ModuleRepository.hpp b/src/language/modules/ModuleRepository.hpp
index 407d7b55bd9b7c8454159a315b9c2f7a256eecb6..35bf6f7aa6a7704f2804125c488ac0bf115ddf39 100644
--- a/src/language/modules/ModuleRepository.hpp
+++ b/src/language/modules/ModuleRepository.hpp
@@ -31,7 +31,23 @@ class ModuleRepository
const IModule::NameValueMap& name_value_map,
SymbolTable& symbol_table);
+ static ModuleRepository* m_instance;
+
public:
+ static ModuleRepository& getInstance();
+
+ static void destroy();
+
+ template <typename ModuleT>
+ struct Subscribe
+ {
+ Subscribe()
+ {
+ static_assert(std::is_base_of_v<IModule, ModuleT>, "module must inherit IModule interface");
+ ModuleRepository::getInstance()._subscribe(std::make_unique<ModuleT>());
+ }
+ };
+
void populateSymbolTable(const ASTNode& module_name_node, SymbolTable& symbol_table);
void populateMandatoryData(const ASTNode& root_node, SymbolTable& symbol_table);
void registerOperators(const std::string& module_name);
@@ -46,7 +62,10 @@ class ModuleRepository
ModuleRepository(const ModuleRepository&) = delete;
ModuleRepository(ModuleRepository&&) = delete;
- ModuleRepository();
+ private:
+ ModuleRepository() = default;
+
+ public:
~ModuleRepository() = default;
};
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e1729dace44e70d88cd35aca7034fde3936c523b..11cb53bc1ee9d79d10edab6a1d3d8ca21c07e0a3 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,4 +1,7 @@
#include <language/modules/SchemeModule.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
+
+#include <language/modules/ModuleRepository.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
@@ -1716,3 +1719,5 @@ SchemeModule::registerCheckpointResume() const
#endif // PUGS_HAS_HDF5
}
+
+ModuleRepository::Subscribe<SchemeModule> scheme_module;
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index c61b3084406dbd8ae9889c79b4172195daaa6fe8..cf169175c4b46f7364799e16bdb4b84633d1f61c 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -2,7 +2,6 @@
#define SCHEME_MODULE_HPP
#include <language/modules/BuiltinModule.hpp>
-#include <language/modules/SchemeModuleTypes.hpp>
class SchemeModule : public BuiltinModule
{
diff --git a/src/language/modules/SocketModule.cpp b/src/language/modules/SocketModule.cpp
index 79321d088f960e0bb7450c82820cf6568e79d174..9c2ecdbe1b40730d0181d7f4ee6f3bf396f70a34 100644
--- a/src/language/modules/SocketModule.cpp
+++ b/src/language/modules/SocketModule.cpp
@@ -1,5 +1,7 @@
#include <language/modules/SocketModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/CheckpointResumeRepository.hpp>
@@ -271,3 +273,5 @@ SocketModule::registerCheckpointResume() const
throw NotImplementedError("checkpoint/resume with sockets");
}));
}
+
+ModuleRepository::Subscribe<SocketModule> socket_module;
diff --git a/src/language/modules/UnaryOperatorRegisterForVh.cpp b/src/language/modules/UnaryOperatorRegisterForVh.cpp
index e2e9c2db9b72a9a687806d24dee836cfd37a22e8..3d52ab54687eabcd6162536d88325b3285d45e1f 100644
--- a/src/language/modules/UnaryOperatorRegisterForVh.cpp
+++ b/src/language/modules/UnaryOperatorRegisterForVh.cpp
@@ -1,6 +1,6 @@
#include <language/modules/UnaryOperatorRegisterForVh.hpp>
-#include <language/modules/SchemeModule.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
#include <language/utils/DataHandler.hpp>
#include <language/utils/DataVariant.hpp>
#include <language/utils/EmbeddedDiscreteFunctionOperators.hpp>
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index 6d84506501de277f1970c13a7c24cd587ed79f20..e524a885639ac9a2a2da6330422565447bcb347d 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -1,7 +1,9 @@
#include <language/modules/WriterModule.hpp>
-#include <language/modules/MeshModule.hpp>
-#include <language/modules/SchemeModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+
+#include <language/modules/MeshModuleTypes.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/CheckpointResumeRepository.hpp>
#include <language/utils/TypeDescriptor.hpp>
@@ -242,3 +244,5 @@ WriterModule::registerCheckpointResume() const
#endif // PUGS_HAS_HDF5
}
+
+ModuleRepository::Subscribe<WriterModule> writer_module{};
diff --git a/src/main.cpp b/src/main.cpp
index c1eeef6a633d5074fef20545a2348a77e5a6888a..cbaf6277e5ee0c1c7d689403e8d78cae54c3ed5f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,7 @@
#include <analysis/QuadratureManager.hpp>
#include <dev/ParallelChecker.hpp>
#include <language/PugsParser.hpp>
+#include <language/modules/ModuleRepository.hpp>
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/DualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
@@ -12,6 +13,8 @@
#include <utils/RandomEngine.hpp>
#include <utils/checkpointing/ResumingManager.hpp>
+#include <utils/PluginsLoader.hpp>
+
int
main(int argc, char* argv[])
{
@@ -23,6 +26,8 @@ main(int argc, char* argv[])
std::string filename = initialize(argc, argv);
+ PluginsLoader plugins_loader;
+
SynchronizerManager::create();
RandomEngine::create();
QuadratureManager::create();
@@ -34,6 +39,8 @@ main(int argc, char* argv[])
parser(filename);
ExecutionStatManager::printInfo();
+ ModuleRepository::destroy();
+
StencilManager::destroy();
DualMeshManager::destroy();
DualConnectivityManager::destroy();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index b4d70359a5a4a7f8a4d2e352ab702ecc26761dbc..781b7e583595eab942d2f924de705cf332841fad 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -50,5 +50,6 @@ add_library(
target_link_libraries(
PugsMesh
+ ${PARMETIS_TARGET}
${HIGHFIVE_TARGET}
)
diff --git a/src/mesh/StencilArray.hpp b/src/mesh/StencilArray.hpp
index 9571c346bae6e4e345a16b09188dce15ab73b456..ffe8638bf0d8fdb472780150603c7fab21ca1eb9 100644
--- a/src/mesh/StencilArray.hpp
+++ b/src/mesh/StencilArray.hpp
@@ -12,6 +12,9 @@ class StencilArray
public:
using ItemToItemMatrixT = ItemToItemMatrix<source_item_type, target_item_type>;
+ using SourceItemId = ItemIdT<source_item_type>;
+ using TargetItemId = ItemIdT<target_item_type>;
+
class BoundaryDescriptorStencilArray
{
private:
@@ -71,9 +74,9 @@ class StencilArray
PUGS_INLINE
auto
- operator[](CellId cell_id) const
+ operator[](SourceItemId source_item_id) const
{
- return m_stencil_array[cell_id];
+ return m_stencil_array[source_item_id];
}
StencilArray(const ConnectivityMatrix& connectivity_matrix,
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 31db762a1002e3c633afeaabd4f06d6a98dff634..20fc3bf61dd6a59ef82ff7bdb25dc5e32c7eef24 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -5,16 +5,23 @@
#include <utils/GlobalVariableManager.hpp>
#include <utils/Messenger.hpp>
-#include <set>
-
template <ItemType item_type>
class StencilBuilder::Layer
{
using ItemId = ItemIdT<item_type>;
+ ItemValue<const int, item_type> m_number_of;
+
std::vector<ItemId> m_item_id_vector;
std::vector<int> m_item_number_vector;
public:
+ void
+ clear()
+ {
+ m_item_id_vector.clear();
+ m_item_number_vector.clear();
+ }
+
size_t
size() const
{
@@ -22,9 +29,17 @@ class StencilBuilder::Layer
return m_item_id_vector.size();
}
+ const std::vector<ItemId>&
+ itemIdList() const
+ {
+ return m_item_id_vector;
+ }
+
bool
- hasItemNumber(const int item_number) const
+ hasItemNumber(const ItemId item_id) const
{
+ const int item_number = m_number_of[item_id];
+
ssize_t begin = 0;
ssize_t end = m_item_number_vector.size();
@@ -37,26 +52,22 @@ class StencilBuilder::Layer
if (begin == mid) {
break;
}
- begin = mid - 1;
+ begin = mid;
} else {
if (end == mid) {
break;
}
- end = mid + 1;
+ end = mid;
}
}
return false;
}
- const std::vector<ItemId>&
- itemIdList() const
- {
- return m_item_id_vector;
- }
-
void
- add(const ItemId item_id, const int item_number)
+ add(const ItemId item_id)
{
+ const int item_number = m_number_of[item_id];
+
ssize_t begin = 0;
ssize_t end = m_item_number_vector.size();
@@ -97,7 +108,14 @@ class StencilBuilder::Layer
}
}
- Layer() = default;
+ Layer& operator=(const Layer&) = default;
+ Layer& operator=(Layer&&) = default;
+
+ template <size_t Dimension>
+ Layer(const Connectivity<Dimension>& connectivity) : m_number_of{connectivity.template number<item_type>()}
+ {}
+
+ // Layer() = default;
Layer(const Layer&) = default;
Layer(Layer&&) = default;
@@ -117,27 +135,25 @@ StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connect
symmetry_boundary_descriptor_list.size()};
symmetry_connecting_item_list.fill(false);
- if constexpr (ConnectivityType::Dimension > 1) {
- auto face_to_connecting_item_matrix =
- connectivity.template getItemToItemMatrix<ItemType::face, connecting_item_type>();
+ if constexpr (ItemTypeId<ConnectivityType::Dimension>::itemTId(connecting_item_type) ==
+ ItemTypeId<ConnectivityType::Dimension>::itemTId(ItemType::face)) {
size_t i_symmetry_boundary = 0;
for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) {
const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
bool found = false;
- for (size_t i_ref_face_list = 0; i_ref_face_list < connectivity.template numberOfRefItemList<ItemType::face>();
- ++i_ref_face_list) {
- const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_face_list);
- if (ref_face_list.refId() == boundary_descriptor) {
- found = true;
- for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) {
- const FaceId face_id = ref_face_list.list()[i_face];
- auto connecting_item_list = face_to_connecting_item_matrix[face_id];
- for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
- const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
+ for (size_t i_ref_connecting_item_list = 0;
+ i_ref_connecting_item_list < connectivity.template numberOfRefItemList<connecting_item_type>();
+ ++i_ref_connecting_item_list) {
+ const auto& ref_connecting_item_list =
+ connectivity.template refItemList<connecting_item_type>(i_ref_connecting_item_list);
+ if (ref_connecting_item_list.refId() == boundary_descriptor) {
+ found = true;
+ auto connecting_item_list = ref_connecting_item_list.list();
+ for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
- symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
- }
+ symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
}
break;
}
@@ -149,26 +165,27 @@ StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connect
throw NormalError(error_msg.str());
}
}
-
- return symmetry_connecting_item_list;
} else {
+ auto face_to_connecting_item_matrix =
+ connectivity.template getItemToItemMatrix<ItemType::face, connecting_item_type>();
size_t i_symmetry_boundary = 0;
for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) {
const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
bool found = false;
- for (size_t i_ref_connecting_item_list = 0;
- i_ref_connecting_item_list < connectivity.template numberOfRefItemList<connecting_item_type>();
- ++i_ref_connecting_item_list) {
- const auto& ref_connecting_item_list =
- connectivity.template refItemList<connecting_item_type>(i_ref_connecting_item_list);
- if (ref_connecting_item_list.refId() == boundary_descriptor) {
- found = true;
- auto connecting_item_list = ref_connecting_item_list.list();
- for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
- const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
+ for (size_t i_ref_face_list = 0; i_ref_face_list < connectivity.template numberOfRefItemList<ItemType::face>();
+ ++i_ref_face_list) {
+ const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_face_list);
+ if (ref_face_list.refId() == boundary_descriptor) {
+ found = true;
+ for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) {
+ const FaceId face_id = ref_face_list.list()[i_face];
+ auto connecting_item_list = face_to_connecting_item_matrix[face_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
- symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
+ symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
+ }
}
break;
}
@@ -180,357 +197,572 @@ StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connect
throw NormalError(error_msg.str());
}
}
-
- return symmetry_connecting_item_list;
}
+
+ return symmetry_connecting_item_list;
}
-template <typename ConnectivityType>
-Array<const uint32_t>
-StencilBuilder::_getRowMap(const ConnectivityType& connectivity) const
+template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType>
+StencilArray<item_type, item_type>
+StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connectivity,
+ const size_t& number_of_layers,
+ const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
{
- auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
- auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ if (number_of_layers == 0) {
+ throw NormalError("number of layers must be greater than 0 to build stencils");
+ }
+
+ auto item_to_connecting_item_matrix = connectivity.template getItemToItemMatrix<item_type, connecting_item_type>();
+ auto connecting_item_to_item_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, item_type>();
+
+ auto item_is_owned = connectivity.template isOwned<item_type>();
+ auto item_number = connectivity.template number<item_type>();
+ auto connecting_item_number = connectivity.template number<connecting_item_type>();
+
+ using ItemId = ItemIdT<item_type>;
+ using ConnectingItemId = ItemIdT<connecting_item_type>;
- auto cell_is_owned = connectivity.cellIsOwned();
+ ItemArray<bool, connecting_item_type> symmetry_item_list =
+ this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list);
- Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+ Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
row_map[0] = 0;
- std::vector<CellId> neighbors;
- for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
- neighbors.resize(0);
- // The stencil is not built for ghost cells
- if (cell_is_owned[cell_id]) {
- auto cell_nodes = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- const NodeId node_id = cell_nodes[i_node];
- auto node_cells = node_to_cell_matrix[node_id];
- for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
- const CellId node_cell_id = node_cells[i_node_cell];
- if (node_cell_id != cell_id) {
- neighbors.push_back(node_cells[i_node_cell]);
+ std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_descriptor_list.size());
+ for (auto&& symmetry_row_map : symmetry_row_map_list) {
+ symmetry_row_map = Array<uint32_t>{connectivity.template numberOf<item_type>() + 1};
+ symmetry_row_map[0] = 0;
+ }
+
+ std::vector<ItemId> column_indices_vector;
+ std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size());
+
+ std::vector<Layer<item_type>> item_layer_list;
+ std::vector<Layer<connecting_item_type>> connecting_layer_list;
+
+ std::vector<Layer<item_type>> symmetry_item_layer_list;
+ std::vector<Layer<connecting_item_type>> symmetry_connecting_layer_list;
+
+ for (size_t i = 0; i < number_of_layers; ++i) {
+ item_layer_list.emplace_back(Layer<item_type>{connectivity});
+ connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity});
+
+ if (symmetry_boundary_descriptor_list.size() > 0) {
+ symmetry_item_layer_list.emplace_back(Layer<item_type>{connectivity});
+ symmetry_connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity});
+ }
+ }
+
+ for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+ if (item_is_owned[item_id]) {
+ for (auto&& item_layer : item_layer_list) {
+ item_layer.clear();
+ }
+ for (auto&& connecting_layer : connecting_layer_list) {
+ connecting_layer.clear();
+ }
+
+ // First layer is a special case
+ {
+ auto& item_layer = item_layer_list[0];
+ auto& connecting_layer = connecting_layer_list[0];
+
+ for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_matrix[item_id].size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = item_to_connecting_item_matrix[item_id][i_connecting_item];
+
+ connecting_layer.add(connecting_item_id);
+ }
+
+ for (auto connecting_item_id : connecting_layer.itemIdList()) {
+ for (size_t i_item = 0; i_item < connecting_item_to_item_matrix[connecting_item_id].size(); ++i_item) {
+ const ItemId layer_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item];
+ if (layer_item_id != item_id) {
+ item_layer.add(layer_item_id);
+ }
}
}
}
- std::sort(neighbors.begin(), neighbors.end());
- neighbors.erase(std::unique(neighbors.begin(), neighbors.end()), neighbors.end());
- }
- // The cell itself is not counted
- row_map[cell_id + 1] = row_map[cell_id] + neighbors.size();
- }
- return row_map;
-}
+ for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) {
+ auto& connecting_layer = connecting_layer_list[i_layer];
-template <typename ConnectivityType>
-Array<const uint32_t>
-StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Array<const uint32_t>& row_map) const
-{
- auto cell_number = connectivity.cellNumber();
+ auto has_connecting_item = [&i_layer, &connecting_layer_list](ConnectingItemId connecting_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (connecting_layer_list[i].hasItemNumber(connecting_item_id)) {
+ return true;
+ }
+ }
- Array<uint32_t> max_index(row_map.size() - 1);
- parallel_for(
- max_index.size(), PUGS_LAMBDA(size_t i) { max_index[i] = row_map[i]; });
-
- auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
- auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
-
- auto cell_is_owned = connectivity.cellIsOwned();
-
- Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
- column_indices.fill(std::numeric_limits<uint32_t>::max());
-
- for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
- // The stencil is not built for ghost cells
- if (cell_is_owned[cell_id]) {
- auto cell_nodes = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- const NodeId node_id = cell_nodes[i_node];
- auto node_cells = node_to_cell_matrix[node_id];
- for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
- const CellId node_cell_id = node_cells[i_node_cell];
- if (node_cell_id != cell_id) {
- bool found = false;
- for (size_t i_index = row_map[cell_id]; i_index < max_index[cell_id]; ++i_index) {
- if (column_indices[i_index] == node_cell_id) {
- found = true;
- break;
+ return false;
+ };
+
+ for (auto&& previous_layer_item_id_list : item_layer_list[i_layer - 1].itemIdList()) {
+ const auto item_to_connecting_item_list = item_to_connecting_item_matrix[previous_layer_item_id_list];
+ for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = item_to_connecting_item_list[i_connecting_item];
+
+ if (not has_connecting_item(connecting_item_id)) {
+ connecting_layer.add(connecting_item_id);
+ }
+ }
+ }
+
+ auto& item_layer = item_layer_list[i_layer];
+
+ auto has_layer_item = [&i_layer, &item_layer_list](ItemId layer_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (item_layer_list[i].hasItemNumber(layer_item_id)) {
+ return true;
+ }
+ }
+
+ return false;
+ };
+
+ for (auto connecting_item_id : connecting_layer.itemIdList()) {
+ const auto& connecting_item_to_item_list = connecting_item_to_item_matrix[connecting_item_id];
+ for (size_t i_item = 0; i_item < connecting_item_to_item_list.size(); ++i_item) {
+ const ItemId layer_item_id = connecting_item_to_item_list[i_item];
+ if ((layer_item_id != item_id) and (not has_layer_item(layer_item_id))) {
+ item_layer.add(layer_item_id);
+ }
+ }
+ }
+ }
+
+ for (size_t i_symmetry = 0; i_symmetry < symmetry_boundary_descriptor_list.size(); ++i_symmetry) {
+ for (auto&& symmetry_item_layer : symmetry_item_layer_list) {
+ symmetry_item_layer.clear();
+ }
+ for (auto&& symmetry_connecting_layer : symmetry_connecting_layer_list) {
+ symmetry_connecting_layer.clear();
+ }
+
+ // First layer is a special case
+ for (auto&& connecting_item_id : connecting_layer_list[0].itemIdList()) {
+ if (symmetry_item_list[connecting_item_id][i_symmetry]) {
+ symmetry_connecting_layer_list[0].add(connecting_item_id);
+ }
+ }
+
+ for (auto&& connecting_item_id : symmetry_connecting_layer_list[0].itemIdList()) {
+ auto item_of_connecting_item_list = connecting_item_to_item_matrix[connecting_item_id];
+ for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < item_of_connecting_item_list.size();
+ ++i_item_of_connecting_item) {
+ const ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item];
+ symmetry_item_layer_list[0].add(item_id_of_connecting_item);
+ }
+ }
+
+ for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) {
+ auto has_connecting_item = [&i_layer,
+ &symmetry_connecting_layer_list](ConnectingItemId connecting_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (symmetry_connecting_layer_list[i].hasItemNumber(connecting_item_id)) {
+ return true;
}
}
- if (not found) {
- const auto node_cell_number = cell_number[node_cell_id];
- size_t i_index = row_map[cell_id];
- // search for position for index
- while ((i_index < max_index[cell_id])) {
- if (node_cell_number > cell_number[CellId(column_indices[i_index])]) {
- ++i_index;
- } else {
- break;
- }
+
+ return false;
+ };
+
+ for (auto&& symmetry_connecting_item_id : connecting_layer_list[i_layer].itemIdList()) {
+ if (symmetry_item_list[symmetry_connecting_item_id][i_symmetry]) {
+ if (not has_connecting_item(symmetry_connecting_item_id)) {
+ symmetry_connecting_layer_list[i_layer].add(symmetry_connecting_item_id);
}
+ }
+ }
- for (size_t i_destination = max_index[cell_id]; i_destination > i_index; --i_destination) {
- const size_t i_source = i_destination - 1;
+ for (auto&& previous_layer_item_id_list : symmetry_item_layer_list[i_layer - 1].itemIdList()) {
+ const auto item_to_connecting_item_list = item_to_connecting_item_matrix[previous_layer_item_id_list];
+ for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = item_to_connecting_item_list[i_connecting_item];
- column_indices[i_destination] = column_indices[i_source];
+ if (not has_connecting_item(connecting_item_id)) {
+ symmetry_connecting_layer_list[i_layer].add(connecting_item_id);
+ }
+ }
+ }
+
+ auto has_symmetry_layer_item = [&i_layer, &symmetry_item_layer_list](ItemId layer_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (symmetry_item_layer_list[i].hasItemNumber(layer_item_id)) {
+ return true;
+ }
+ }
+
+ return false;
+ };
+
+ for (auto&& connecting_item_id : symmetry_connecting_layer_list[i_layer].itemIdList()) {
+ auto item_of_connecting_item_list = connecting_item_to_item_matrix[connecting_item_id];
+ for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < item_of_connecting_item_list.size();
+ ++i_item_of_connecting_item) {
+ const ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item];
+ if (not has_symmetry_layer_item(item_id_of_connecting_item)) {
+ symmetry_item_layer_list[i_layer].add(item_id_of_connecting_item);
}
- ++max_index[cell_id];
- column_indices[i_index] = node_cell_id;
}
}
}
+
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ for (auto symmetry_layer_item_id : symmetry_item_layer_list[i_layer].itemIdList()) {
+ symmetry_column_indices_vector[i_symmetry].push_back(symmetry_layer_item_id);
+ }
+ }
+
+ {
+ size_t stencil_size = 0;
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ auto& item_layer = symmetry_item_layer_list[i_layer];
+ stencil_size += item_layer.size();
+ }
+
+ symmetry_row_map_list[i_symmetry][item_id + 1] = symmetry_row_map_list[i_symmetry][item_id] + stencil_size;
+ }
+ }
+
+ {
+ size_t stencil_size = 0;
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ auto& item_layer = item_layer_list[i_layer];
+ for (auto stencil_item_id : item_layer.itemIdList()) {
+ column_indices_vector.push_back(stencil_item_id);
+ }
+ stencil_size += item_layer.size();
+ }
+ row_map[item_id + 1] = row_map[item_id] + stencil_size;
+ }
+
+ } else {
+ row_map[item_id + 1] = row_map[item_id];
+ for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) {
+ symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id];
}
}
}
- return column_indices;
+ Array<uint32_t> column_indices{column_indices_vector.size()};
+ parallel_for(
+ column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
+
+ ConnectivityMatrix primal_stencil{row_map, column_indices};
+
+ typename StencilArray<item_type, item_type>::BoundaryDescriptorStencilArrayList symmetry_boundary_stencil_list;
+ {
+ size_t i = 0;
+ for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
+ symmetry_boundary_stencil_list.emplace_back(
+ typename StencilArray<item_type, item_type>::
+ BoundaryDescriptorStencilArray{p_boundary_descriptor,
+ ConnectivityMatrix{symmetry_row_map_list[i],
+ convert_to_array(symmetry_column_indices_vector[i])}});
+ ++i;
+ }
+ }
+
+ return {{primal_stencil}, {symmetry_boundary_stencil_list}};
}
-template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType>
-StencilArray<item_type, item_type>
-StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connectivity,
- const size_t& number_of_layers,
- const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+template <ItemType source_item_type,
+ ItemType connecting_item_type,
+ ItemType target_item_type,
+ typename ConnectivityType>
+StencilArray<source_item_type, target_item_type>
+StencilBuilder::_build_for_different_source_and_target(
+ const ConnectivityType& connectivity,
+ const size_t& number_of_layers,
+ const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
{
+ constexpr size_t Dimension = ConnectivityType::Dimension;
+
if (number_of_layers == 0) {
throw NormalError("number of layers must be greater than 0 to build stencils");
}
- if (number_of_layers > 2) {
- throw NotImplementedError("number of layers too large");
- }
- auto item_to_connecting_item_matrix = connectivity.template getItemToItemMatrix<item_type, connecting_item_type>();
- auto connecting_item_to_item_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, item_type>();
- auto item_is_owned = connectivity.template isOwned<item_type>();
- auto item_number = connectivity.template number<item_type>();
+ auto connecting_item_to_target_item_matrix =
+ connectivity.template getItemToItemMatrix<connecting_item_type, target_item_type>();
+ auto target_item_to_connecting_item_matrix =
+ connectivity.template getItemToItemMatrix<target_item_type, connecting_item_type>();
+
+ auto source_item_to_connecting_item_matrix =
+ [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> {
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ return {ConnectivityMatrix{}};
+ } else {
+ return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>();
+ }
+ }(connectivity);
+
+ auto source_item_is_owned = connectivity.template isOwned<source_item_type>();
+ auto target_item_number = connectivity.template number<target_item_type>();
auto connecting_item_number = connectivity.template number<connecting_item_type>();
- using ItemId = ItemIdT<item_type>;
+ using SourceItemId = ItemIdT<source_item_type>;
+ using TargetItemId = ItemIdT<target_item_type>;
using ConnectingItemId = ItemIdT<connecting_item_type>;
- if (symmetry_boundary_descriptor_list.size() == 0) {
- if (number_of_layers == 1) {
- Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
- row_map[0] = 0;
+ ItemArray<bool, connecting_item_type> symmetry_item_list =
+ this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list);
+
+ Array<uint32_t> row_map{connectivity.template numberOf<source_item_type>() + 1};
+ row_map[0] = 0;
+ std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_descriptor_list.size());
+ for (auto&& symmetry_row_map : symmetry_row_map_list) {
+ symmetry_row_map = Array<uint32_t>{connectivity.template numberOf<source_item_type>() + 1};
+ symmetry_row_map[0] = 0;
+ }
+
+ std::vector<TargetItemId> column_indices_vector;
+ std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size());
- std::vector<ItemId> column_indices_vector;
+ std::vector<Layer<target_item_type>> target_item_layer_list;
+ std::vector<Layer<connecting_item_type>> connecting_layer_list;
- for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
- // First layer is a special case
+ std::vector<Layer<target_item_type>> symmetry_target_item_layer_list;
+ std::vector<Layer<connecting_item_type>> symmetry_connecting_layer_list;
- Layer<item_type> item_layer;
- Layer<connecting_item_type> connecting_layer;
+ for (size_t i = 0; i < number_of_layers; ++i) {
+ target_item_layer_list.emplace_back(Layer<target_item_type>{connectivity});
+ connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity});
- if (item_is_owned[item_id]) {
- for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size();
- ++i_connecting_item_1) {
- const ConnectingItemId layer_1_connecting_item_id =
- item_to_connecting_item_matrix[item_id][i_connecting_item_1];
- connecting_layer.add(layer_1_connecting_item_id, connecting_item_number[layer_1_connecting_item_id]);
- }
+ if (symmetry_boundary_descriptor_list.size() > 0) {
+ symmetry_target_item_layer_list.emplace_back(Layer<target_item_type>{connectivity});
+ symmetry_connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity});
+ }
+ }
- for (auto connecting_item_id : connecting_layer.itemIdList()) {
- for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[connecting_item_id].size();
- ++i_item_1) {
- const ItemId layer_1_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item_1];
- if (layer_1_item_id != item_id) {
- item_layer.add(layer_1_item_id, item_number[layer_1_item_id]);
- }
- }
+ for (SourceItemId source_item_id = 0; source_item_id < connectivity.template numberOf<source_item_type>();
+ ++source_item_id) {
+ if (source_item_is_owned[source_item_id]) {
+ for (auto&& target_item_layer : target_item_layer_list) {
+ target_item_layer.clear();
+ }
+ for (auto&& connecting_layer : connecting_layer_list) {
+ connecting_layer.clear();
+ }
+
+ // First layer is a special case
+ {
+ auto& target_item_layer = target_item_layer_list[0];
+ auto& connecting_layer = connecting_layer_list[0];
+
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ connecting_layer.add(ConnectingItemId{static_cast<typename ConnectingItemId::base_type>(source_item_id)});
+ } else {
+ for (size_t i_connecting_item = 0;
+ i_connecting_item < source_item_to_connecting_item_matrix[source_item_id].size(); ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id =
+ source_item_to_connecting_item_matrix[source_item_id][i_connecting_item];
+
+ connecting_layer.add(connecting_item_id);
}
}
- for (auto layer_item_id : item_layer.itemIdList()) {
- column_indices_vector.push_back(layer_item_id);
+ for (auto connecting_item_id : connecting_layer.itemIdList()) {
+ for (size_t i_item = 0; i_item < connecting_item_to_target_item_matrix[connecting_item_id].size(); ++i_item) {
+ const TargetItemId layer_item_id = connecting_item_to_target_item_matrix[connecting_item_id][i_item];
+ target_item_layer.add(layer_item_id);
+ }
}
- row_map[item_id + 1] = row_map[item_id] + item_layer.itemIdList().size();
}
- if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
- throw UnexpectedError("incorrect stencil size");
- }
- Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
- column_indices.fill(std::numeric_limits<uint32_t>::max());
-
- for (size_t i = 0; i < column_indices.size(); ++i) {
- column_indices[i] = column_indices_vector[i];
- }
+ for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) {
+ auto& connecting_layer = connecting_layer_list[i_layer];
- return {ConnectivityMatrix{row_map, column_indices}, {}};
- } else {
- Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
- row_map[0] = 0;
-
- std::vector<ItemId> column_indices_vector;
-
- for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
- if (item_is_owned[item_id]) {
- std::set<ItemId, std::function<bool(ItemId, ItemId)>> item_set(
- [=](ItemId item_0, ItemId item_1) { return item_number[item_0] < item_number[item_1]; });
-
- for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size();
- ++i_connecting_item_1) {
- const ConnectingItemId layer_1_connecting_item_id =
- item_to_connecting_item_matrix[item_id][i_connecting_item_1];
-
- for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[layer_1_connecting_item_id].size();
- ++i_item_1) {
- ItemId item_1_id = connecting_item_to_item_matrix[layer_1_connecting_item_id][i_item_1];
-
- for (size_t i_connecting_item_2 = 0;
- i_connecting_item_2 < item_to_connecting_item_matrix[item_1_id].size(); ++i_connecting_item_2) {
- const ConnectingItemId layer_2_connecting_item_id =
- item_to_connecting_item_matrix[item_1_id][i_connecting_item_2];
-
- for (size_t i_item_2 = 0; i_item_2 < connecting_item_to_item_matrix[layer_2_connecting_item_id].size();
- ++i_item_2) {
- ItemId item_2_id = connecting_item_to_item_matrix[layer_2_connecting_item_id][i_item_2];
-
- if (item_2_id != item_id) {
- item_set.insert(item_2_id);
- }
- }
- }
+ auto has_connecting_item = [&i_layer, &connecting_layer_list](ConnectingItemId connecting_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (connecting_layer_list[i].hasItemNumber(connecting_item_id)) {
+ return true;
}
}
- for (auto stencil_item_id : item_set) {
- column_indices_vector.push_back(stencil_item_id);
+ return false;
+ };
+
+ for (auto&& previous_layer_item_id_list : target_item_layer_list[i_layer - 1].itemIdList()) {
+ const auto target_item_to_connecting_item_list =
+ target_item_to_connecting_item_matrix[previous_layer_item_id_list];
+ for (size_t i_connecting_item = 0; i_connecting_item < target_item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = target_item_to_connecting_item_list[i_connecting_item];
+
+ if (not has_connecting_item(connecting_item_id)) {
+ connecting_layer.add(connecting_item_id);
+ }
}
- row_map[item_id + 1] = row_map[item_id] + item_set.size();
}
- }
- if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
- throw UnexpectedError("incorrect stencil size");
- }
+ auto& target_item_layer = target_item_layer_list[i_layer];
- Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
- column_indices.fill(std::numeric_limits<uint32_t>::max());
+ auto has_layer_item = [&i_layer, &target_item_layer_list](TargetItemId layer_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (target_item_layer_list[i].hasItemNumber(layer_item_id)) {
+ return true;
+ }
+ }
- for (size_t i = 0; i < column_indices.size(); ++i) {
- column_indices[i] = column_indices_vector[i];
- }
- ConnectivityMatrix primal_stencil{row_map, column_indices};
+ return false;
+ };
- return {primal_stencil, {}};
- }
- } else {
- ItemArray<bool, connecting_item_type> symmetry_item_list =
- this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list);
-
- Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
- row_map[0] = 0;
- std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_descriptor_list.size());
- for (auto&& symmetry_row_map : symmetry_row_map_list) {
- symmetry_row_map = Array<uint32_t>{connectivity.template numberOf<item_type>() + 1};
- symmetry_row_map[0] = 0;
- }
+ for (auto connecting_item_id : connecting_layer.itemIdList()) {
+ const auto& connecting_item_to_target_item_list = connecting_item_to_target_item_matrix[connecting_item_id];
+ for (size_t i_item = 0; i_item < connecting_item_to_target_item_list.size(); ++i_item) {
+ const TargetItemId layer_item_id = connecting_item_to_target_item_list[i_item];
+ if (not has_layer_item(layer_item_id)) {
+ target_item_layer.add(layer_item_id);
+ }
+ }
+ }
+ }
- std::vector<uint32_t> column_indices_vector;
- std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size());
+ for (size_t i_symmetry = 0; i_symmetry < symmetry_boundary_descriptor_list.size(); ++i_symmetry) {
+ for (auto&& symmetry_target_item_layer : symmetry_target_item_layer_list) {
+ symmetry_target_item_layer.clear();
+ }
+ for (auto&& symmetry_connecting_layer : symmetry_connecting_layer_list) {
+ symmetry_connecting_layer.clear();
+ }
- for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
- std::set<ItemId> item_set;
- std::vector<std::set<ItemId>> by_boundary_symmetry_item(symmetry_boundary_descriptor_list.size());
+ // First layer is a special case
+ for (auto&& connecting_item_id : connecting_layer_list[0].itemIdList()) {
+ if (symmetry_item_list[connecting_item_id][i_symmetry]) {
+ symmetry_connecting_layer_list[0].add(connecting_item_id);
+ }
+ }
- if (item_is_owned[item_id]) {
- auto item_to_connecting_item_list = item_to_connecting_item_matrix[item_id];
- for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size();
- ++i_connecting_item_of_item) {
- const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item];
- auto connecting_item_to_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item];
- for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < connecting_item_to_item_list.size();
+ for (auto&& connecting_item_id : symmetry_connecting_layer_list[0].itemIdList()) {
+ auto item_of_connecting_item_list = connecting_item_to_target_item_matrix[connecting_item_id];
+ for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < item_of_connecting_item_list.size();
++i_item_of_connecting_item) {
- const ItemId item_id_of_connecting_item = connecting_item_to_item_list[i_item_of_connecting_item];
- if (item_id != item_id_of_connecting_item) {
- item_set.insert(item_id_of_connecting_item);
- }
+ const TargetItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item];
+ symmetry_target_item_layer_list[0].add(item_id_of_connecting_item);
}
}
- {
- std::vector<ItemId> item_vector;
- for (auto&& set_item_id : item_set) {
- item_vector.push_back(set_item_id);
+ for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) {
+ auto has_connecting_item = [&i_layer,
+ &symmetry_connecting_layer_list](ConnectingItemId connecting_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (symmetry_connecting_layer_list[i].hasItemNumber(connecting_item_id)) {
+ return true;
+ }
+ }
+
+ return false;
+ };
+
+ for (auto&& symmetry_connecting_item_id : connecting_layer_list[i_layer].itemIdList()) {
+ if (symmetry_item_list[symmetry_connecting_item_id][i_symmetry]) {
+ if (not has_connecting_item(symmetry_connecting_item_id)) {
+ symmetry_connecting_layer_list[i_layer].add(symmetry_connecting_item_id);
+ }
+ }
}
- std::sort(item_vector.begin(), item_vector.end(),
- [&item_number](const ItemId& item0_id, const ItemId& item1_id) {
- return item_number[item0_id] < item_number[item1_id];
- });
- for (auto&& vector_item_id : item_vector) {
- column_indices_vector.push_back(vector_item_id);
+ for (auto&& previous_layer_target_item_id_list : symmetry_target_item_layer_list[i_layer - 1].itemIdList()) {
+ const auto item_to_connecting_item_list =
+ target_item_to_connecting_item_matrix[previous_layer_target_item_id_list];
+ for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = item_to_connecting_item_list[i_connecting_item];
+
+ if (not has_connecting_item(connecting_item_id)) {
+ symmetry_connecting_layer_list[i_layer].add(connecting_item_id);
+ }
+ }
}
- }
- for (size_t i = 0; i < symmetry_boundary_descriptor_list.size(); ++i) {
- std::set<ItemId> symmetry_item_set;
- for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size();
- ++i_connecting_item_of_item) {
- const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item];
- if (symmetry_item_list[connecting_item_id_of_item][i]) {
- auto item_of_connecting_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item];
- for (size_t i_item_of_connecting_item = 0;
- i_item_of_connecting_item < item_of_connecting_item_list.size(); ++i_item_of_connecting_item) {
- const ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item];
- symmetry_item_set.insert(item_id_of_connecting_item);
+ auto has_symmetry_layer_item = [&i_layer,
+ &symmetry_target_item_layer_list](TargetItemId layer_target_item_id) -> bool {
+ for (ssize_t i = i_layer - 1; i >= 0; --i) {
+ if (symmetry_target_item_layer_list[i].hasItemNumber(layer_target_item_id)) {
+ return true;
+ }
+ }
+
+ return false;
+ };
+
+ for (auto&& connecting_item_id : symmetry_connecting_layer_list[i_layer].itemIdList()) {
+ auto item_of_connecting_item_list = connecting_item_to_target_item_matrix[connecting_item_id];
+ for (size_t i_target_item_of_connecting_item = 0;
+ i_target_item_of_connecting_item < item_of_connecting_item_list.size();
+ ++i_target_item_of_connecting_item) {
+ const TargetItemId target_item_id_of_connecting_item =
+ item_of_connecting_item_list[i_target_item_of_connecting_item];
+ if (not has_symmetry_layer_item(target_item_id_of_connecting_item)) {
+ symmetry_target_item_layer_list[i_layer].add(target_item_id_of_connecting_item);
}
}
}
- by_boundary_symmetry_item[i] = symmetry_item_set;
+ }
- std::vector<ItemId> item_vector;
- for (auto&& set_item_id : symmetry_item_set) {
- item_vector.push_back(set_item_id);
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ for (auto symmetry_layer_target_item_id : symmetry_target_item_layer_list[i_layer].itemIdList()) {
+ symmetry_column_indices_vector[i_symmetry].push_back(symmetry_layer_target_item_id);
}
- std::sort(item_vector.begin(), item_vector.end(),
- [&item_number](const ItemId& item0_id, const ItemId& item1_id) {
- return item_number[item0_id] < item_number[item1_id];
- });
+ }
- for (auto&& vector_item_id : item_vector) {
- symmetry_column_indices_vector[i].push_back(vector_item_id);
+ {
+ size_t stencil_size = 0;
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ auto& target_item_layer = symmetry_target_item_layer_list[i_layer];
+ stencil_size += target_item_layer.size();
}
+
+ symmetry_row_map_list[i_symmetry][source_item_id + 1] =
+ symmetry_row_map_list[i_symmetry][source_item_id] + stencil_size;
}
}
- row_map[item_id + 1] = row_map[item_id] + item_set.size();
- for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) {
- symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id] + by_boundary_symmetry_item[i].size();
+ {
+ size_t stencil_size = 0;
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ auto& item_layer = target_item_layer_list[i_layer];
+ for (auto stencil_item_id : item_layer.itemIdList()) {
+ column_indices_vector.push_back(stencil_item_id);
+ }
+ stencil_size += item_layer.size();
+ }
+ row_map[source_item_id + 1] = row_map[source_item_id] + stencil_size;
}
- }
- ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
-
- typename StencilArray<item_type, item_type>::BoundaryDescriptorStencilArrayList symmetry_boundary_stencil_list;
- {
- size_t i = 0;
- for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
- symmetry_boundary_stencil_list.emplace_back(
- typename StencilArray<item_type, item_type>::
- BoundaryDescriptorStencilArray{p_boundary_descriptor,
- ConnectivityMatrix{symmetry_row_map_list[i],
- convert_to_array(symmetry_column_indices_vector[i])}});
- ++i;
+
+ } else {
+ row_map[source_item_id + 1] = row_map[source_item_id];
+ for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) {
+ symmetry_row_map_list[i][source_item_id + 1] = symmetry_row_map_list[i][source_item_id];
}
}
+ }
+
+ Array<uint32_t> column_indices{column_indices_vector.size()};
+ parallel_for(
+ column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; });
+
+ ConnectivityMatrix primal_stencil{row_map, column_indices};
- return {{primal_stencil}, {symmetry_boundary_stencil_list}};
+ typename StencilArray<source_item_type, target_item_type>::BoundaryDescriptorStencilArrayList
+ symmetry_boundary_stencil_list;
+ {
+ size_t i = 0;
+ for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
+ symmetry_boundary_stencil_list.emplace_back(
+ typename StencilArray<source_item_type, target_item_type>::
+ BoundaryDescriptorStencilArray{p_boundary_descriptor,
+ ConnectivityMatrix{symmetry_row_map_list[i],
+ convert_to_array(symmetry_column_indices_vector[i])}});
+ ++i;
+ }
}
-}
-template <ItemType source_item_type,
- ItemType connecting_item_type,
- ItemType target_item_type,
- typename ConnectivityType>
-StencilArray<source_item_type, target_item_type>
-StencilBuilder::_build_for_different_source_and_target(
- const ConnectivityType& connectivity,
- const size_t& number_of_layers,
- const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
-{
- static_assert(source_item_type != target_item_type);
- throw NotImplementedError("different source target");
+ return {{primal_stencil}, {symmetry_boundary_stencil_list}};
}
template <ItemType source_item_type,
@@ -631,7 +863,35 @@ StencilBuilder::buildC2C(const IConnectivity& connectivity,
}
NodeToCellStencilArray
-StencilBuilder::buildN2C(const IConnectivity&, const StencilDescriptor&, const BoundaryDescriptorList&) const
+StencilBuilder::buildN2C(const IConnectivity& connectivity,
+ const StencilDescriptor& stencil_descriptor,
+ const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
{
- throw NotImplementedError("node to cell stencil");
+ if ((parallel::size() > 1) and
+ (stencil_descriptor.numberOfLayers() > GlobalVariableManager::instance().getNumberOfGhostLayers())) {
+ std::ostringstream error_msg;
+ error_msg << "Stencil builder requires" << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
+ << rang::fg::reset << " layers while parallel number of ghost layer is "
+ << GlobalVariableManager::instance().getNumberOfGhostLayers() << ".\n";
+ error_msg << "Increase the number of ghost layers (using the '--number-of-ghost-layers' option).";
+ throw NormalError(error_msg.str());
+ }
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<1>&>(connectivity),
+ stencil_descriptor, symmetry_boundary_descriptor_list);
+ }
+ case 2: {
+ return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<2>&>(connectivity),
+ stencil_descriptor, symmetry_boundary_descriptor_list);
+ }
+ case 3: {
+ return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<3>&>(connectivity),
+ stencil_descriptor, symmetry_boundary_descriptor_list);
+ }
+ default: {
+ throw UnexpectedError("invalid connectivity dimension");
+ }
+ }
}
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 4f7ead00b618db999dd5386258ac058b664c65a7..2ed1d3a6c45fdac651ce62ba05690532c506cf26 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -19,13 +19,6 @@ class StencilBuilder
template <ItemType item_type>
class Layer;
- template <typename ConnectivityType>
- Array<const uint32_t> _getRowMap(const ConnectivityType& connectivity) const;
-
- template <typename ConnectivityType>
- Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity,
- const Array<const uint32_t>& row_map) const;
-
template <ItemType connecting_item, typename ConnectivityType>
auto _buildSymmetryConnectingItemList(const ConnectivityType& connectivity,
const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
@@ -74,11 +67,12 @@ class StencilBuilder
const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
friend class StencilManager;
+
template <ItemType source_item_type, ItemType target_item_type>
StencilArray<source_item_type, target_item_type>
- build(const IConnectivity& connectivity,
- const StencilDescriptor& stencil_descriptor,
- const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+ _build(const IConnectivity& connectivity,
+ const StencilDescriptor& stencil_descriptor,
+ const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
{
if constexpr ((source_item_type == ItemType::cell) and (target_item_type == ItemType::cell)) {
return buildC2C(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list);
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
index 842f8858f9d1a571e28d3f0397f70e97682da0af..04a8e32a1d50d872d95896e8e7a58e76e32e2189 100644
--- a/src/mesh/StencilManager.cpp
+++ b/src/mesh/StencilManager.cpp
@@ -49,8 +49,8 @@ StencilManager::_getStencilArray(
Key{connectivity.id(), stencil_descriptor, symmetry_boundary_descriptor_list})) {
stored_source_to_target_stencil_map[Key{connectivity.id(), stencil_descriptor, symmetry_boundary_descriptor_list}] =
std::make_shared<StencilArray<source_item_type, target_item_type>>(
- StencilBuilder{}.template build<source_item_type, target_item_type>(connectivity, stencil_descriptor,
- symmetry_boundary_descriptor_list));
+ StencilBuilder{}.template _build<source_item_type, target_item_type>(connectivity, stencil_descriptor,
+ symmetry_boundary_descriptor_list));
}
return *stored_source_to_target_stencil_map.at(
diff --git a/src/scheme/DirichletVectorBoundaryConditionDescriptor.hpp b/src/scheme/DirichletVectorBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c1f1dcc3fadeced1dbed1f560a6e72c86eebf0ed
--- /dev/null
+++ b/src/scheme/DirichletVectorBoundaryConditionDescriptor.hpp
@@ -0,0 +1,58 @@
+#ifndef DIRICHLET_VECTOR_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define DIRICHLET_VECTOR_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
+
+#include <memory>
+#include <vector>
+
+class DirichletVectorBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "dirichlet_vector(" << m_name << ',' << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ const std::string m_name;
+
+ const std::vector<FunctionSymbolId> m_rhs_symbol_id_list;
+
+ public:
+ const std::string&
+ name() const
+ {
+ return m_name;
+ }
+
+ const std::vector<FunctionSymbolId>&
+ rhsSymbolIdList() const
+ {
+ return m_rhs_symbol_id_list;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::dirichlet_vector;
+ }
+
+ DirichletVectorBoundaryConditionDescriptor(const std::string_view name,
+ std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::vector<FunctionSymbolId>& rhs_symbol_id_list)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}, m_name{name}, m_rhs_symbol_id_list{rhs_symbol_id_list}
+ {
+ ;
+ }
+
+ DirichletVectorBoundaryConditionDescriptor(const DirichletVectorBoundaryConditionDescriptor&) = delete;
+ DirichletVectorBoundaryConditionDescriptor(DirichletVectorBoundaryConditionDescriptor&&) = delete;
+
+ ~DirichletVectorBoundaryConditionDescriptor() = default;
+};
+
+#endif // DIRICHLET_VECTOR_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index fc9ff381b97758e79642647a36680364093c193b..63f698e16eb611f9cb03bc82e6d737372924d36d 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -13,6 +13,7 @@ class IBoundaryConditionDescriptor
{
axis,
dirichlet,
+ dirichlet_vector,
external,
fourier,
fixed,
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 0e0b2a5afe28330dad7ad0a82b7dc9b277f027c8..0bc65fc85c5e8f9a650a212c0517a2c0d21b8e6b 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -32,6 +32,14 @@ symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimensio
return u - 2 * dot(u, normal) * normal;
}
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
+{
+ const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
+ return S * A * S;
+}
+
template <size_t Dimension>
PUGS_INLINE auto
symmetrize_coordinates(const TinyVector<Dimension>& origin,
@@ -790,6 +798,44 @@ PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
return mutable_discrete_function_dpk_variant_list;
}
+template <MeshConcept MeshType>
+void
+PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+ for (auto&& discrete_function_variant : discrete_function_variant_list) {
+ std::visit(
+ [&](auto&& discrete_function) {
+ using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunctionT> or
+ is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+ if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (DataType::Dimension != MeshType::Dimension) {
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+ << " using a mesh of dimension " << MeshType::Dimension;
+ throw NormalError(error_msg.str());
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ if constexpr ((DataType::NumberOfRows != MeshType::Dimension) or
+ (DataType::NumberOfColumns != MeshType::Dimension)) {
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+ << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+ throw NormalError(error_msg.str());
+ }
+ }
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("invalid discrete function type");
+ // LCOV_EXCL_STOP
+ }
+ },
+ discrete_function_variant->discreteFunction());
+ }
+}
+
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
@@ -800,6 +846,10 @@ PolynomialReconstruction::_build(
using Rd = TinyVector<MeshType::Dimension>;
+ if (m_descriptor.symmetryBoundaryDescriptorList().size() > 0) {
+ this->_checkDataAndSymmetriesCompatibility<MeshType>(discrete_function_variant_list);
+ }
+
const size_t number_of_columns = this->_getNumberOfColumns(discrete_function_variant_list);
const size_t basis_dimension =
@@ -897,7 +947,7 @@ PolynomialReconstruction::_build(
}
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+ mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
const int32_t t = tokens.acquire();
@@ -954,21 +1004,32 @@ PolynomialReconstruction::_build(
B(index, kB) = qi_qj[k];
}
} else {
- const DataType& qi_qj = discrete_function[cell_i_id] - qj;
- for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
- B(index, kB) = qi_qj[k];
- }
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+ << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
(DataType::NumberOfColumns == MeshType::Dimension)) {
- throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0");
- }
- const DataType& qi_qj = discrete_function[cell_i_id] - qj;
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ const DataType& qi = discrete_function[cell_i_id];
+ const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ }
}
+ } else {
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+ << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
}
}
}
@@ -1031,8 +1092,6 @@ PolynomialReconstruction::_build(
}
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
- throw NotImplementedError("NIY TinyMatrix data");
-
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index b2b855271d1ed4140ebb51d1b21314602caa0dc6..bf6013639f732fd52f9c145ba4982701c0e4bcfc 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -17,6 +17,10 @@ class PolynomialReconstruction
size_t _getNumberOfColumns(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+ template <MeshConcept MeshType>
+ void _checkDataAndSymmetriesCompatibility(
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
template <MeshConcept MeshType>
std::vector<MutableDiscreteFunctionDPkVariant> _createMutableDiscreteFunctionDPKVariantList(
const std::shared_ptr<const MeshType>& p_mesh,
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 13b0db80862c2b28e80dc18ff22951b08157db28..4d74c5e195fe0361819315cd3732661ba7589e56 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -16,6 +16,7 @@ add_library(
Messenger.cpp
Partitioner.cpp
PETScWrapper.cpp
+ PluginsLoader.cpp
PugsUtils.cpp
RandomEngine.cpp
ReproducibleSumManager.cpp
@@ -32,11 +33,15 @@ endif()
target_link_libraries(
PugsUtils
- ${PETSC_LIBRARIES}
- ${SLEPC_LIBRARIES}
+ ${PETSC_TARGET}
+ ${SLEPC_TARGET}
+ ${SLURM_TARGET}
${HIGHFIVE_TARGET}
)
+target_include_directories(PugsUtils PUBLIC ${PETSC_INCLUDE_DIRS})
+target_include_directories(PugsUtils PUBLIC ${SLURM_INCLUDE_DIRS})
+
# --------------- get git revision info ---------------
# Generates revision header file
@@ -58,12 +63,11 @@ set_source_files_properties(
add_custom_command(TARGET PugsGitRevison
COMMAND ${CMAKE_COMMAND} -E remove -f
${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision
- DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision.hpp
+ POST_BUILD
COMMENT ""
)
add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision.hpp
- PRE_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_if_different
${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision
${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision.hpp
@@ -73,7 +77,6 @@ add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision.hpp
)
add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/pugs_git_revision
- PRE_BUILD
COMMAND ${CMAKE_COMMAND} -DPUGS_VERSION=${PUGS_VERSION} -DPUGS_SOURCE_DIR=${PUGS_SOURCE_DIR} -P ${PUGS_SOURCE_DIR}/cmake/GetPugsGitRevision.cmake
COMMENT "Check pugs git status"
VERBATIM
diff --git a/src/utils/ConsoleManager.hpp b/src/utils/ConsoleManager.hpp
index a2afe46fbc0ad16e263d1f5a9b3f8a84d5202d94..693219e30285d1b0acae26d4b43d310882069784 100644
--- a/src/utils/ConsoleManager.hpp
+++ b/src/utils/ConsoleManager.hpp
@@ -1,7 +1,7 @@
#ifndef CONSOLE_MANAGER_HPP
#define CONSOLE_MANAGER_HPP
-#include <string>
+#include <ostream>
class ConsoleManager
{
diff --git a/src/utils/PluginsLoader.cpp b/src/utils/PluginsLoader.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..35e24605e8a01b4d350aadb18bcacda313e9004c
--- /dev/null
+++ b/src/utils/PluginsLoader.cpp
@@ -0,0 +1,104 @@
+#include <utils/PluginsLoader.hpp>
+
+#include <utils/ConsoleManager.hpp>
+
+#include <dlfcn.h>
+
+#include <rang.hpp>
+
+#include <filesystem>
+#include <iostream>
+#include <sstream>
+#include <unistd.h>
+#include <vector>
+
+std::vector<std::string>
+split(const std::string& full_string)
+{
+ std::vector<std::string> split_string;
+
+ std::stringstream is(full_string);
+
+ std::string segment;
+ while (std::getline(is, segment, ';')) {
+ if (segment.size() > 0) {
+ split_string.push_back(segment);
+ }
+ }
+ return split_string;
+}
+
+void
+PluginsLoader::_open(const std::string& plugin)
+{
+ auto handle = dlopen(plugin.c_str(), RTLD_NOW);
+ if (handle != nullptr) {
+ m_dl_handler_stack.push(handle);
+ if (ConsoleManager::showPreamble()) {
+ std::cout << " * \"" << rang::fgB::green << plugin << rang::fg::reset << "\"\n";
+ }
+ } else {
+ std::cerr << " " << rang::fgB::red << "cannot load " << rang::fg::reset << '\"' << rang::fgB::yellow << plugin
+ << rang::fg::reset << "\"\n";
+ }
+}
+
+PluginsLoader::PluginsLoader()
+{
+ std::vector<std::string> plugin_vector;
+
+ {
+ char* env = getenv("PUGS_PLUGIN");
+ if (env != nullptr) {
+ std::string plugins = env;
+ plugin_vector = split(plugins);
+ }
+ }
+
+ {
+ char* env = getenv("PUGS_PLUGIN_DIR");
+ if (env != nullptr) {
+ std::string paths = env;
+ std::vector<std::string> path_vector = split(paths);
+ for (auto&& path : path_vector) {
+ if (access(path.c_str(), R_OK) == -1) {
+ std::cerr << ' ' << rang::fgB::red << 'X' << rang::fg::reset << " cannot access plugin dir \""
+ << rang::fgB::yellow << path << rang::fg::reset << "\"\n";
+ } else {
+ for (auto&& entry :
+ std::filesystem::directory_iterator(path,
+ (std::filesystem::directory_options::follow_directory_symlink |
+ std::filesystem::directory_options::skip_permission_denied))) {
+ if (entry.path().extension() == ".so") {
+ plugin_vector.push_back(entry.path().string());
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // keep unique entries
+ std::sort(plugin_vector.begin(), plugin_vector.end());
+ plugin_vector.resize(std::distance(plugin_vector.begin(), std::unique(plugin_vector.begin(), plugin_vector.end())));
+
+ if (plugin_vector.size() > 0) {
+ if (ConsoleManager::showPreamble()) {
+ std::cout << rang::style::bold << "Loading plugins" << rang::style::reset << '\n';
+ }
+ for (auto&& plugin : plugin_vector) {
+ this->_open(plugin);
+ }
+ if (ConsoleManager::showPreamble()) {
+ std::cout << "-------------------------------------------------------\n";
+ }
+ }
+}
+
+PluginsLoader::~PluginsLoader()
+{
+ while (not m_dl_handler_stack.empty()) {
+ dlclose(m_dl_handler_stack.top());
+ m_dl_handler_stack.pop();
+ }
+}
diff --git a/src/utils/PluginsLoader.hpp b/src/utils/PluginsLoader.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..76be25c6934c5f336360e938d9a40d75f075e6d3
--- /dev/null
+++ b/src/utils/PluginsLoader.hpp
@@ -0,0 +1,23 @@
+#ifndef PLUGINS_LOADER_HPP
+#define PLUGINS_LOADER_HPP
+
+#include <stack>
+#include <string>
+
+class PluginsLoader
+{
+ private:
+ std::stack<void*> m_dl_handler_stack;
+
+ void _open(const std::string& filename);
+
+ public:
+ PluginsLoader();
+
+ PluginsLoader(const PluginsLoader&) = delete;
+ PluginsLoader(PluginsLoader&&) = delete;
+
+ ~PluginsLoader();
+};
+
+#endif // PLUGINS_LOADER_HPP
diff --git a/src/utils/checkpointing/CMakeLists.txt b/src/utils/checkpointing/CMakeLists.txt
index b051975e8538df40a2b9b8a8d7f0486874cda938..06d34a349ab63a6e0de206908e50a82d418a857a 100644
--- a/src/utils/checkpointing/CMakeLists.txt
+++ b/src/utils/checkpointing/CMakeLists.txt
@@ -66,4 +66,8 @@ add_dependencies(PugsCheckpointing
target_link_libraries(
PugsCheckpointing
${HIGHFIVE_TARGET}
+# ${HDF5_TARGET}
)
+
+#target_include_directories(PugsCheckpointing PUBLIC ${HDF5_INCLUDE_DIRS})
+#target_include_directories(PugsAlgebra PUBLIC ${SLEPC_INCLUDE_DIRS})
diff --git a/src/utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp b/src/utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp
index 210b521164075f275829865397ea1771a460c207..4a6e5c92cded1a614d43a17ace96f211a3f48e2d 100644
--- a/src/utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp
+++ b/src/utils/checkpointing/IBoundaryConditionDescriptorHFType.hpp
@@ -10,6 +10,7 @@ create_enum_i_boundary_condition_descriptor_type()
{
return {{"axis", IBoundaryConditionDescriptor::Type::axis},
{"dirichlet", IBoundaryConditionDescriptor::Type::dirichlet},
+ {"dirichlet_vector", IBoundaryConditionDescriptor::Type::dirichlet_vector},
{"external", IBoundaryConditionDescriptor::Type::external},
{"fixed", IBoundaryConditionDescriptor::Type::fixed},
{"fourier", IBoundaryConditionDescriptor::Type::fourier},
diff --git a/src/utils/checkpointing/ReadIBoundaryConditionDescriptor.cpp b/src/utils/checkpointing/ReadIBoundaryConditionDescriptor.cpp
index 55b43c9f04bf1b3b585df208bd7a3ea86a3d2850..1f405b4d8a72a8cd15e7799732a16fc4a71ac12b 100644
--- a/src/utils/checkpointing/ReadIBoundaryConditionDescriptor.cpp
+++ b/src/utils/checkpointing/ReadIBoundaryConditionDescriptor.cpp
@@ -5,6 +5,7 @@
#include <language/utils/EmbeddedData.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DirichletVectorBoundaryConditionDescriptor.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
@@ -48,6 +49,21 @@ readIBoundaryConditionDescriptor(const HighFive::Group& iboundaryconditiondecrip
std::make_shared<const DirichletBoundaryConditionDescriptor>(name, i_boundary_descriptor,
*ResumingData::instance().functionSymbolId(rhs_id));
break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet_vector: {
+ const std::string name = iboundaryconditiondecriptor_group.getAttribute("name").read<std::string>();
+
+ const std::vector function_id_list =
+ iboundaryconditiondecriptor_group.getAttribute("function_id_list").read<std::vector<size_t>>();
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+ for (auto function_id : function_id_list) {
+ function_symbol_id_list.push_back(*ResumingData::instance().functionSymbolId(function_id));
+ }
+
+ bc_descriptor = std::make_shared<const DirichletVectorBoundaryConditionDescriptor>(name, i_boundary_descriptor,
+ function_symbol_id_list);
+ break;
}
// LCOV_EXCL_START
case IBoundaryConditionDescriptor::Type::external: {
diff --git a/src/utils/checkpointing/WriteIBoundaryConditionDescriptor.cpp b/src/utils/checkpointing/WriteIBoundaryConditionDescriptor.cpp
index 66a8670583ae63442313da42c847bbeee758a2d5..8a3c5dd062aca586bd800a19f9b0e27f4fc30e93 100644
--- a/src/utils/checkpointing/WriteIBoundaryConditionDescriptor.cpp
+++ b/src/utils/checkpointing/WriteIBoundaryConditionDescriptor.cpp
@@ -5,6 +5,7 @@
#include <language/utils/DataHandler.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DirichletVectorBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
@@ -47,6 +48,18 @@ writeIBoundaryConditionDescriptor(HighFive::Group& variable_group,
variable_group.createAttribute("name", dirichlet_bc_descriptor.name());
variable_group.createAttribute("rhs_function_id", dirichlet_bc_descriptor.rhsSymbolId().id());
break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet_vector: {
+ const DirichletVectorBoundaryConditionDescriptor& dirichlet_vector_bc_descriptor =
+ dynamic_cast<const DirichletVectorBoundaryConditionDescriptor&>(iboundary_condition_descriptor);
+ variable_group.createAttribute("name", dirichlet_vector_bc_descriptor.name());
+ writeIBoundaryDescriptor(boundary_group, dirichlet_vector_bc_descriptor.boundaryDescriptor());
+ std::vector<size_t> function_id_list;
+ for (auto&& function_symbol_id : dirichlet_vector_bc_descriptor.rhsSymbolIdList()) {
+ function_id_list.push_back(function_symbol_id.id());
+ }
+ variable_group.createAttribute("function_id_list", function_id_list);
+ break;
}
// LCOV_EXCL_START
case IBoundaryConditionDescriptor::Type::external: {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 874df648b6b6719d8f5b1e79cfcd6a46683e16ea..690030bb8aaafde19d6a7f5129a5d676f4d599c7 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -4,6 +4,20 @@ include_directories(${PUGS_SOURCE_DIR}/src)
include_directories(${PUGS_BINARY_DIR}/src)
include_directories(${PUGS_SOURCE_DIR}/tests)
+install(
+ DIRECTORY "${PUGS_SOURCE_DIR}/packages/Catch2/src/catch2"
+ DESTINATION "include"
+ FILES_MATCHING
+ PATTERN "*.hpp"
+)
+
+install(
+ DIRECTORY "${PUGS_BINARY_DIR}/generated-includes/catch2"
+ DESTINATION "include"
+ FILES_MATCHING
+ PATTERN "*.hpp"
+)
+
set(checkpointing_sequential_TESTS
# this one should enventually integrate parallel tests
test_checkpointing_Checkpoint_sequential.cpp
@@ -264,7 +278,8 @@ add_executable (mpi_unit_tests
test_PolynomialReconstruction.cpp
test_PolynomialReconstructionDescriptor.cpp
test_RandomEngine.cpp
- test_StencilBuilder.cpp
+ test_StencilBuilder_cell2cell.cpp
+ test_StencilBuilder_node2cell.cpp
test_SubItemArrayPerItemVariant.cpp
test_SubItemValuePerItem.cpp
test_SubItemValuePerItemVariant.cpp
diff --git a/tests/NbGhostLayersTester.hpp b/tests/NbGhostLayersTester.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0468fe18305367834c9cbd0959ba5b3f443fe870
--- /dev/null
+++ b/tests/NbGhostLayersTester.hpp
@@ -0,0 +1,27 @@
+#ifndef NB_GHOST_LAYERS_TESTER_HPP
+#define NB_GHOST_LAYERS_TESTER_HPP
+
+#include <cstddef>
+#include <utils/GlobalVariableManager.hpp>
+
+class NbGhostLayersTester
+{
+ private:
+ const size_t m_original_number_of_ghost_layers;
+
+ public:
+ PUGS_INLINE
+ NbGhostLayersTester(const size_t number_of_ghost_layers)
+ : m_original_number_of_ghost_layers{GlobalVariableManager::instance().getNumberOfGhostLayers()}
+ {
+ GlobalVariableManager::instance().m_number_of_ghost_layers = number_of_ghost_layers;
+ }
+
+ PUGS_INLINE
+ ~NbGhostLayersTester()
+ {
+ GlobalVariableManager::instance().m_number_of_ghost_layers = m_original_number_of_ghost_layers;
+ }
+};
+
+#endif // NB_GHOST_LAYERS_TESTER_HPP
diff --git a/tests/test_ConnectivityDispatcher.cpp b/tests/test_ConnectivityDispatcher.cpp
index abe07cc21e21a1b34ae3f20dfe13ff7f52ac83b9..09ec9982e10e08d84c808f157fe4e11add3d756a 100644
--- a/tests/test_ConnectivityDispatcher.cpp
+++ b/tests/test_ConnectivityDispatcher.cpp
@@ -9,29 +9,12 @@
#include <utils/Messenger.hpp>
#include <MeshDataBaseForTests.hpp>
+#include <NbGhostLayersTester.hpp>
#include <filesystem>
// clazy:excludeall=non-pod-global-static
-class NbGhostLayersTester
-{
- private:
- const size_t m_original_number_of_ghost_layers;
-
- public:
- NbGhostLayersTester(const size_t number_of_ghost_layers)
- : m_original_number_of_ghost_layers{GlobalVariableManager::instance().getNumberOfGhostLayers()}
- {
- GlobalVariableManager::instance().m_number_of_ghost_layers = number_of_ghost_layers;
- }
-
- ~NbGhostLayersTester()
- {
- GlobalVariableManager::instance().m_number_of_ghost_layers = m_original_number_of_ghost_layers;
- }
-};
-
TEST_CASE("ConnectivityDispatcher", "[mesh]")
{
auto check_number_of_ghost_layers = [](const auto& connectivity, const size_t number_of_layers) {
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
deleted file mode 100644
index b99ba218aa1e3eb5d8ec884e974cfad2f4a157db..0000000000000000000000000000000000000000
--- a/tests/test_StencilBuilder.cpp
+++ /dev/null
@@ -1,349 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <MeshDataBaseForTests.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/ConnectivityUtils.hpp>
-#include <mesh/ItemValue.hpp>
-#include <mesh/ItemValueUtils.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshFlatNodeBoundary.hpp>
-#include <mesh/MeshVariant.hpp>
-#include <mesh/NamedBoundaryDescriptor.hpp>
-#include <mesh/StencilManager.hpp>
-#include <utils/Messenger.hpp>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("StencilBuilder", "[mesh]")
-{
- SECTION("inner stencil")
- {
- auto is_valid = [](const auto& connectivity, const auto& stencil_array) {
- auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
- auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
-
- auto cell_is_owned = connectivity.cellIsOwned();
- auto cell_number = connectivity.cellNumber();
-
- for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
- if (cell_is_owned[cell_id]) {
- std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
- [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
- auto cell_nodes = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- const NodeId node_id = cell_nodes[i_node];
- auto node_cells = node_to_cell_matrix[node_id];
- for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
- const CellId node_cell_id = node_cells[i_node_cell];
- if (node_cell_id != cell_id) {
- cell_set.insert(node_cell_id);
- }
- }
- }
-
- auto cell_stencil = stencil_array[cell_id];
-
- auto i_set_cell = cell_set.begin();
- for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
- if (*i_set_cell != cell_stencil[index]) {
- return false;
- }
- }
- }
- }
- return true;
- };
-
- SECTION("1D")
- {
- SECTION("cartesian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
-
- const Connectivity<1>& connectivity = mesh.connectivity();
-
- auto stencil_array =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes});
-
- REQUIRE(is_valid(connectivity, stencil_array));
- }
-
- SECTION("unordered")
- {
- const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-
- const Connectivity<1>& connectivity = mesh.connectivity();
- auto stencil_array =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes});
-
- REQUIRE(is_valid(connectivity, stencil_array));
- }
- }
-
- SECTION("2D")
- {
- SECTION("cartesian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
-
- const Connectivity<2>& connectivity = mesh.connectivity();
- REQUIRE(
- is_valid(connectivity,
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity,
- StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
- }
-
- SECTION("hybrid")
- {
- const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
-
- const Connectivity<2>& connectivity = mesh.connectivity();
- REQUIRE(
- is_valid(connectivity,
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity,
- StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
- }
- }
-
- SECTION("3D")
- {
- SECTION("carteian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
-
- const Connectivity<3>& connectivity = mesh.connectivity();
- REQUIRE(
- is_valid(connectivity,
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity,
- StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
- }
-
- SECTION("hybrid")
- {
- const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
-
- const Connectivity<3>& connectivity = mesh.connectivity();
- REQUIRE(
- is_valid(connectivity,
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity,
- StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
- }
- }
- }
-
- SECTION("Stencil using symmetries")
- {
- auto check_ghost_cells_have_empty_stencils = [](const auto& stencil_array, const auto& connecticity) {
- bool is_empty = true;
-
- auto cell_is_owned = connecticity.cellIsOwned();
-
- for (CellId cell_id = 0; cell_id < connecticity.numberOfCells(); ++cell_id) {
- if (not cell_is_owned[cell_id]) {
- is_empty &= (stencil_array[cell_id].size() == 0);
- for (size_t i_symmetry_stencil = 0;
- i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) {
- is_empty &=
- (stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray()[cell_id].size() ==
- 0);
- }
- }
- }
-
- return is_empty;
- };
-
- auto symmetry_stencils_are_valid = [](const auto& stencil_array, const auto& mesh) {
- bool is_valid = true;
-
- auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
- auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- auto cell_is_owned = mesh.connectivity().cellIsOwned();
- auto cell_number = mesh.connectivity().cellNumber();
-
- for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry_stencil) {
- const IBoundaryDescriptor& boundary_descriptor =
- stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor();
-
- auto boundary_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray();
- auto boundary_node_list = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
-
- CellValue<bool> boundary_cell{mesh.connectivity()};
- boundary_cell.fill(false);
- auto node_list = boundary_node_list.nodeList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId node_id = node_list[i_node];
- auto node_cell_list = node_to_cell_matrix[node_id];
- for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
- const CellId cell_id = node_cell_list[i_cell];
- boundary_cell[cell_id] = true;
- }
- }
-
- std::set<NodeId> symmetry_node;
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId node_id = node_list[i_node];
- symmetry_node.insert(node_id);
- }
-
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- if ((not boundary_cell[cell_id]) or (not cell_is_owned[cell_id])) {
- is_valid &= (boundary_stencil[cell_id].size() == 0);
- } else {
- auto cell_node_list = cell_to_node_matrix[cell_id];
- std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
- [&](CellId cell0_id, CellId cell1_id) { return cell_number[cell0_id] < cell_number[cell1_id]; });
- for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
- const NodeId node_id = cell_node_list[i_node];
- if (symmetry_node.contains(node_id)) {
- auto node_cell_list = node_to_cell_matrix[node_id];
- for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
- const CellId node_cell_id = node_cell_list[i_node_cell];
- cell_set.insert(node_cell_id);
- }
- }
- }
-
- if (cell_set.size() == boundary_stencil[cell_id].size()) {
- size_t i = 0;
- for (auto&& id : cell_set) {
- is_valid &= (id == boundary_stencil[cell_id][i++]);
- }
- } else {
- is_valid = false;
- }
- }
- }
- }
-
- return is_valid;
- };
-
- SECTION("1D")
- {
- StencilManager::BoundaryDescriptorList list;
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
-
- SECTION("cartesian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
-
- const Connectivity<1>& connectivity = mesh.connectivity();
-
- auto stencil_array =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh));
- }
-
- SECTION("hybrid")
- {
- const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-
- const Connectivity<1>& connectivity = mesh.connectivity();
- auto stencil =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 2);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
- }
- }
-
- SECTION("2D")
- {
- StencilManager::BoundaryDescriptorList list;
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
-
- SECTION("cartesian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
-
- const Connectivity<2>& connectivity = mesh.connectivity();
-
- auto stencil_array =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh));
- }
-
- SECTION("hybrid")
- {
- const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
-
- const Connectivity<2>& connectivity = mesh.connectivity();
- auto stencil =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 4);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
- }
- }
-
- SECTION("3D")
- {
- StencilManager::BoundaryDescriptorList list;
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN"));
- list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX"));
-
- SECTION("cartesian")
- {
- const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
-
- const Connectivity<3>& connectivity = mesh.connectivity();
- auto stencil =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
- }
-
- SECTION("hybrid")
- {
- const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
-
- const Connectivity<3>& connectivity = mesh.connectivity();
- auto stencil =
- StencilManager::instance()
- .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
- list);
-
- REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6);
- REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
- REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
- }
- }
- }
-}
diff --git a/tests/test_StencilBuilder_cell2cell.cpp b/tests/test_StencilBuilder_cell2cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..42ceb6a5da583a9869035db09f5239d008c6f780
--- /dev/null
+++ b/tests/test_StencilBuilder_cell2cell.cpp
@@ -0,0 +1,1081 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/CartesianMeshBuilder.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/StencilManager.hpp>
+#include <utils/Messenger.hpp>
+
+#include <NbGhostLayersTester.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("StencilBuilder cell2cell", "[mesh]")
+{
+ auto is_valid = []<ItemType connecting_item_type>(const auto& connectivity, const auto& stencil_array,
+ const size_t number_of_layers) {
+ auto cell_to_connecting_item_matrix =
+ connectivity.template getItemToItemMatrix<ItemType::cell, connecting_item_type>();
+ auto connecting_to_cell_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, ItemType::cell>();
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
+
+ using ConnectingItemId = ItemIdT<connecting_item_type>;
+
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ if (cell_is_owned[cell_id]) {
+ std::vector<CellId> expected_stencil;
+
+ std::set<CellId> marked_cell_set;
+ marked_cell_set.insert(cell_id);
+
+ std::set<ConnectingItemId> marked_connecting_item_set;
+ std::set<ConnectingItemId> layer_connecting_item_set;
+ {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0; i_connecting_item_of_cell < connecting_item_to_cell_list.size();
+ ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_cell_set.contains(connecting_item_id_of_cell)) {
+ cell_set.insert(connecting_item_id_of_cell);
+ marked_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ layer_connecting_item_set.clear();
+ for (auto layer_cell_id : cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_connecting_item_set.contains(connecting_item_id)) {
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ for (auto&& set_cell_id : cell_set) {
+ expected_stencil.push_back(set_cell_id);
+ }
+ }
+
+ auto cell_stencil = stencil_array[cell_id];
+
+ auto i_set_cell = expected_stencil.begin();
+ if (cell_stencil.size() != expected_stencil.size()) {
+ return false;
+ }
+
+ for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
+ if (*i_set_cell != cell_stencil[index]) {
+ return false;
+ }
+ }
+ }
+ }
+ return true;
+ };
+
+ SECTION("inner stencil")
+ {
+ SECTION("1 layer")
+ {
+ SECTION("1D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+
+ SECTION("2D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+
+ SECTION("3D")
+ {
+ SECTION("carteian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+ }
+
+ SECTION("2 layers")
+ {
+ NbGhostLayersTester nb_ghost_layers_tester(2);
+
+ SECTION("1D")
+ {
+ SECTION("cartesian")
+ {
+ auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh();
+
+ const auto& mesh = *mesh_v->get<Mesh<1>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+
+ SECTION("2D")
+ {
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<2>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+
+ SECTION("3D")
+ {
+ SECTION("carteian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5))
+ .mesh();
+ const auto& mesh = *mesh_v->get<Mesh<3>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+ }
+ }
+
+ SECTION("Stencil using symmetries")
+ {
+ auto check_ghost_cells_have_empty_stencils = [](const auto& stencil_array, const auto& connecticity) {
+ bool is_empty = true;
+
+ auto cell_is_owned = connecticity.cellIsOwned();
+
+ for (CellId cell_id = 0; cell_id < connecticity.numberOfCells(); ++cell_id) {
+ if (not cell_is_owned[cell_id]) {
+ is_empty &= (stencil_array[cell_id].size() == 0);
+ for (size_t i_symmetry_stencil = 0;
+ i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) {
+ is_empty &=
+ (stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray()[cell_id].size() ==
+ 0);
+ }
+ }
+ }
+
+ return is_empty;
+ };
+
+ auto are_symmetry_stencils_valid = []<ItemType connecting_item_type>(const auto& stencil_array, const auto& mesh,
+ const size_t number_of_layers) {
+ bool are_valid_symmetries = true;
+
+ auto connecting_to_cell_matrix =
+ mesh.connectivity().template getItemToItemMatrix<connecting_item_type, ItemType::cell>();
+ auto cell_to_connecting_item_matrix =
+ mesh.connectivity().template getItemToItemMatrix<ItemType::cell, connecting_item_type>();
+ auto cell_is_owned = mesh.connectivity().cellIsOwned();
+ auto cell_number = mesh.connectivity().cellNumber();
+
+ auto connecting_number = mesh.connectivity().template number<connecting_item_type>();
+
+ using ConnectingItemId = ItemIdT<connecting_item_type>;
+ using MeshType = std::decay_t<decltype(mesh)>;
+
+ for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size();
+ ++i_symmetry_stencil) {
+ const IBoundaryDescriptor& boundary_descriptor =
+ stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor();
+
+ auto boundary_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray();
+ auto boundary_face_list = getMeshFaceBoundary(mesh, boundary_descriptor);
+
+ std::set<ConnectingItemId> sym_connecting_item_set;
+
+ if constexpr (ItemTypeId<MeshType::Dimension>::itemTId(connecting_item_type) ==
+ ItemTypeId<MeshType::Dimension>::itemTId(ItemType::face)) {
+ for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) {
+ const FaceId face_id = boundary_face_list.faceList()[i_face];
+ sym_connecting_item_set.insert(ConnectingItemId{FaceId::base_type{face_id}});
+ }
+
+ } else {
+ auto face_to_connecting_item_matrix =
+ mesh.connectivity().template getItemToItemMatrix<ItemType::face, connecting_item_type>();
+
+ for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) {
+ const FaceId face_id = boundary_face_list.faceList()[i_face];
+ auto connecting_item_list = face_to_connecting_item_matrix[face_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+ sym_connecting_item_set.insert(connecting_item_list[i_connecting_item]);
+ }
+ }
+ }
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ if (not cell_is_owned[cell_id]) {
+ are_valid_symmetries &= (boundary_stencil[cell_id].size() == 0);
+ } else {
+ std::vector<CellId> expected_stencil;
+
+ std::set<CellId> marked_cell_set;
+ marked_cell_set.insert(cell_id);
+
+ std::set<ConnectingItemId> marked_connecting_item_set;
+ std::set<ConnectingItemId> layer_connecting_item_set;
+ {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+
+ std::set<ConnectingItemId> marked_sym_connecting_item_set;
+ std::set<ConnectingItemId> layer_sym_connecting_item_set;
+
+ for (auto&& connecting_item_id : marked_connecting_item_set) {
+ if (sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+
+ std::set<CellId> marked_sym_cell_set;
+
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0;
+ i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_cell_set.contains(connecting_item_id_of_cell)) {
+ cell_set.insert(connecting_item_id_of_cell);
+ marked_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ std::set<CellId, std::function<bool(CellId, CellId)>> sym_cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_sym_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0;
+ i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_sym_cell_set.contains(connecting_item_id_of_cell)) {
+ sym_cell_set.insert(connecting_item_id_of_cell);
+ marked_sym_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ for (auto&& set_sym_cell_id : sym_cell_set) {
+ expected_stencil.push_back(set_sym_cell_id);
+ }
+
+ layer_connecting_item_set.clear();
+ for (auto&& layer_cell_id : cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_connecting_item_set.contains(connecting_item_id)) {
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ layer_sym_connecting_item_set.clear();
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ if (sym_connecting_item_set.contains(connecting_item_id)) {
+ if (not marked_sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ for (auto layer_sym_cell_id : sym_cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_sym_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+ }
+
+ auto cell_stencil = boundary_stencil[cell_id];
+
+ if (cell_stencil.size() != expected_stencil.size()) {
+ are_valid_symmetries = false;
+ }
+
+ auto i_set_cell = expected_stencil.begin();
+ for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
+ if (*i_set_cell != cell_stencil[index]) {
+ are_valid_symmetries = false;
+ }
+ }
+ }
+ }
+ }
+
+ return are_valid_symmetries;
+ };
+
+ SECTION("1 layer")
+ {
+ SECTION("1D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(not(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1)));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(not(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1)));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+ }
+
+ SECTION("2 layers")
+ {
+ NbGhostLayersTester nb_ghost_layers_tester(2);
+
+ SECTION("1D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(not(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2)));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5))
+ .mesh();
+ const auto& mesh = *mesh_v->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getCellToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::face>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+ }
+ }
+}
diff --git a/tests/test_StencilBuilder_node2cell.cpp b/tests/test_StencilBuilder_node2cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d9c429c6d4ba903185c90a1e6c37d63bf9c9feb
--- /dev/null
+++ b/tests/test_StencilBuilder_node2cell.cpp
@@ -0,0 +1,1184 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/CartesianMeshBuilder.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/StencilManager.hpp>
+#include <utils/Messenger.hpp>
+
+#include <NbGhostLayersTester.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("StencilBuilder node2cell", "[mesh]")
+{
+ auto is_valid = []<ItemType source_item_type, ItemType connecting_item_type>(const auto& connectivity,
+ const auto& stencil_array,
+ const size_t number_of_layers) {
+ constexpr size_t Dimension = std::decay_t<decltype(connectivity)>::Dimension;
+
+ auto cell_to_connecting_item_matrix =
+ connectivity.template getItemToItemMatrix<ItemType::cell, connecting_item_type>();
+
+ auto source_item_to_connecting_item_matrix =
+ [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> {
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ return ConnectivityMatrix{};
+ } else {
+ return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>();
+ }
+ }(connectivity);
+
+ auto connecting_to_cell_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, ItemType::cell>();
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
+ auto source_item_is_owned = connectivity.template isOwned<source_item_type>();
+
+ using SourceItemId = ItemIdT<source_item_type>;
+ using ConnectingItemId = ItemIdT<connecting_item_type>;
+
+ for (SourceItemId source_item_id = 0; source_item_id < connectivity.template numberOf<source_item_type>();
+ ++source_item_id) {
+ if (source_item_is_owned[source_item_id]) {
+ std::vector<CellId> expected_stencil;
+
+ std::set<CellId> marked_cell_set;
+ if constexpr (source_item_type == ItemType::cell) {
+ marked_cell_set.insert(source_item_id);
+ }
+
+ std::set<ConnectingItemId> marked_connecting_item_set;
+ std::set<ConnectingItemId> layer_connecting_item_set;
+ {
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ ConnectingItemId connecting_id =
+ ConnectingItemId{static_cast<typename SourceItemId::base_type>(source_item_id)};
+ layer_connecting_item_set.insert(connecting_id);
+ marked_connecting_item_set.insert(connecting_id);
+ } else {
+ auto source_item_to_connecting_item_list = source_item_to_connecting_item_matrix[source_item_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < source_item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = source_item_to_connecting_item_list[i_connecting_item];
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0; i_connecting_item_of_cell < connecting_item_to_cell_list.size();
+ ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_cell_set.contains(connecting_item_id_of_cell)) {
+ cell_set.insert(connecting_item_id_of_cell);
+ marked_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ layer_connecting_item_set.clear();
+ for (auto layer_cell_id : cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_connecting_item_set.contains(connecting_item_id)) {
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ for (auto&& set_cell_id : cell_set) {
+ expected_stencil.push_back(set_cell_id);
+ }
+ }
+
+ auto cell_stencil = stencil_array[source_item_id];
+
+ auto i_set_cell = expected_stencil.begin();
+ if (cell_stencil.size() != expected_stencil.size()) {
+ return false;
+ }
+
+ for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
+ if (*i_set_cell != cell_stencil[index]) {
+ return false;
+ }
+ }
+ }
+ }
+ return true;
+ };
+
+ SECTION("inner stencil")
+ {
+ SECTION("1 layer")
+ {
+ SECTION("1D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+
+ SECTION("2D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+
+ SECTION("3D")
+ {
+ SECTION("carteian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 1));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(connectivity,
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 1));
+ }
+ }
+ }
+
+ SECTION("2 layers")
+ {
+ NbGhostLayersTester nb_ghost_layers_tester(2);
+
+ SECTION("1D")
+ {
+ SECTION("cartesian")
+ {
+ auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh();
+
+ const auto& mesh = *mesh_v->get<Mesh<1>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+
+ SECTION("2D")
+ {
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<2>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+
+ SECTION("3D")
+ {
+ SECTION("carteian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5))
+ .mesh();
+ const auto& mesh = *mesh_v->get<Mesh<3>>();
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::node>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_nodes}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::edge>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_edges}),
+ 2));
+
+ REQUIRE(
+ is_valid.template
+ operator()<ItemType::node,
+ ItemType::face>(mesh.connectivity(),
+ StencilManager::instance()
+ .getNodeToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{2, StencilDescriptor::
+ ConnectionType::by_faces}),
+ 2));
+ }
+ }
+ }
+ }
+
+ SECTION("Stencil using symmetries")
+ {
+ auto check_ghost_nodes_have_empty_stencils = [](const auto& stencil_array, const auto& connecticity) {
+ bool is_empty = true;
+
+ auto node_is_owned = connecticity.nodeIsOwned();
+
+ for (NodeId node_id = 0; node_id < connecticity.numberOfNodes(); ++node_id) {
+ if (not node_is_owned[node_id]) {
+ is_empty &= (stencil_array[node_id].size() == 0);
+ for (size_t i_symmetry_stencil = 0;
+ i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) {
+ is_empty &=
+ (stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray()[node_id].size() ==
+ 0);
+ }
+ }
+ }
+
+ return is_empty;
+ };
+
+ auto are_symmetry_stencils_valid =
+ []<ItemType source_item_type, ItemType connecting_item_type>(const auto& stencil_array, const auto& mesh,
+ const size_t number_of_layers) {
+ bool are_valid_symmetries = true;
+
+ constexpr const size_t Dimension = std::decay_t<decltype(mesh)>::Dimension;
+
+ auto connecting_to_cell_matrix =
+ mesh.connectivity().template getItemToItemMatrix<connecting_item_type, ItemType::cell>();
+ auto cell_to_connecting_item_matrix =
+ mesh.connectivity().template getItemToItemMatrix<ItemType::cell, connecting_item_type>();
+ // auto source_item_is_owned = mesh.connectivity().cellIsOwned();
+ auto cell_number = mesh.connectivity().cellNumber();
+ auto source_item_is_owned = mesh.connectivity().template isOwned<source_item_type>();
+
+ auto connecting_number = mesh.connectivity().template number<connecting_item_type>();
+
+ auto source_item_to_connecting_item_matrix =
+ [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> {
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ return ConnectivityMatrix{};
+ } else {
+ return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>();
+ }
+ }(mesh.connectivity());
+
+ using SourceItemId = ItemIdT<source_item_type>;
+ using ConnectingItemId = ItemIdT<connecting_item_type>;
+ using MeshType = std::decay_t<decltype(mesh)>;
+
+ for (size_t i_symmetry_stencil = 0;
+ i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) {
+ const IBoundaryDescriptor& boundary_descriptor =
+ stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor();
+
+ auto boundary_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray();
+ auto boundary_face_list = getMeshFaceBoundary(mesh, boundary_descriptor);
+
+ std::set<ConnectingItemId> sym_connecting_item_set;
+
+ if constexpr (ItemTypeId<MeshType::Dimension>::itemTId(connecting_item_type) ==
+ ItemTypeId<MeshType::Dimension>::itemTId(ItemType::face)) {
+ for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) {
+ const FaceId face_id = boundary_face_list.faceList()[i_face];
+ sym_connecting_item_set.insert(ConnectingItemId{FaceId::base_type{face_id}});
+ }
+
+ } else {
+ auto face_to_connecting_item_matrix =
+ mesh.connectivity().template getItemToItemMatrix<ItemType::face, connecting_item_type>();
+
+ for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) {
+ const FaceId face_id = boundary_face_list.faceList()[i_face];
+ auto connecting_item_list = face_to_connecting_item_matrix[face_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+ sym_connecting_item_set.insert(connecting_item_list[i_connecting_item]);
+ }
+ }
+ }
+
+ for (SourceItemId source_item_id = 0; source_item_id < mesh.template numberOf<source_item_type>();
+ ++source_item_id) {
+ if (not source_item_is_owned[source_item_id]) {
+ are_valid_symmetries &= (boundary_stencil[source_item_id].size() == 0);
+ } else {
+ std::vector<CellId> expected_stencil;
+
+ std::set<CellId> marked_cell_set;
+ if constexpr (source_item_type == ItemType::cell) {
+ marked_cell_set.insert(source_item_id);
+ }
+
+ std::set<ConnectingItemId> marked_connecting_item_set;
+ std::set<ConnectingItemId> layer_connecting_item_set;
+ if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) ==
+ ItemTypeId<Dimension>::itemTId(connecting_item_type)) {
+ const ConnectingItemId connecting_item_id =
+ static_cast<typename SourceItemId::base_type>(source_item_id);
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ } else {
+ auto source_item_to_connecting_item_list = source_item_to_connecting_item_matrix[source_item_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < source_item_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = source_item_to_connecting_item_list[i_connecting_item];
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+
+ std::set<ConnectingItemId> marked_sym_connecting_item_set;
+ std::set<ConnectingItemId> layer_sym_connecting_item_set;
+
+ for (auto&& connecting_item_id : marked_connecting_item_set) {
+ if (sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+
+ std::set<CellId> marked_sym_cell_set;
+
+ for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0;
+ i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_cell_set.contains(connecting_item_id_of_cell)) {
+ cell_set.insert(connecting_item_id_of_cell);
+ marked_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ std::set<CellId, std::function<bool(CellId, CellId)>> sym_cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+
+ for (auto&& connecting_item_id : layer_sym_connecting_item_set) {
+ auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id];
+ for (size_t i_connecting_item_of_cell = 0;
+ i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) {
+ const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell];
+ if (not marked_sym_cell_set.contains(connecting_item_id_of_cell)) {
+ sym_cell_set.insert(connecting_item_id_of_cell);
+ marked_sym_cell_set.insert(connecting_item_id_of_cell);
+ }
+ }
+ }
+
+ for (auto&& set_sym_cell_id : sym_cell_set) {
+ expected_stencil.push_back(set_sym_cell_id);
+ }
+
+ layer_connecting_item_set.clear();
+ for (auto&& layer_cell_id : cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_connecting_item_set.contains(connecting_item_id)) {
+ layer_connecting_item_set.insert(connecting_item_id);
+ marked_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ layer_sym_connecting_item_set.clear();
+
+ for (auto&& connecting_item_id : layer_connecting_item_set) {
+ if (sym_connecting_item_set.contains(connecting_item_id)) {
+ if (not marked_sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+
+ for (auto layer_sym_cell_id : sym_cell_set) {
+ auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_sym_cell_id];
+ for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size();
+ ++i_connecting_item) {
+ const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item];
+ if (not marked_sym_connecting_item_set.contains(connecting_item_id)) {
+ marked_sym_connecting_item_set.insert(connecting_item_id);
+ layer_sym_connecting_item_set.insert(connecting_item_id);
+ }
+ }
+ }
+ }
+
+ auto cell_stencil = boundary_stencil[source_item_id];
+
+ if (cell_stencil.size() != expected_stencil.size()) {
+ are_valid_symmetries = false;
+ }
+
+ auto i_set_cell = expected_stencil.begin();
+ for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
+ if (*i_set_cell != cell_stencil[index]) {
+ are_valid_symmetries = false;
+ }
+ }
+ }
+ }
+ }
+
+ return are_valid_symmetries;
+ };
+
+ SECTION("1 layer")
+ {
+ SECTION("1D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX"));
+
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1));
+ }
+ }
+ }
+ }
+
+ SECTION("2 layers")
+ {
+ NbGhostLayersTester nb_ghost_layers_tester(2);
+
+ SECTION("1D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh();
+ const auto& mesh = *mesh_v->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(not(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ StencilManager::BoundaryDescriptorList list;
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN"));
+ list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX"));
+
+ SECTION("cartesian")
+ {
+ auto mesh_v =
+ CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5))
+ .mesh();
+ const auto& mesh = *mesh_v->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 2));
+ }
+
+ {
+ auto stencil_array =
+ StencilManager::instance()
+ .getNodeToCellStencilArray(connectivity,
+ StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list);
+
+ REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6);
+ REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity));
+ REQUIRE(
+ are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2));
+ REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 2));
+ }
+ }
+ }
+ }
+ }
+}
diff --git a/tests/test_checkpointing_IBoundaryConditionDescriptor.cpp b/tests/test_checkpointing_IBoundaryConditionDescriptor.cpp
index 34505db62eff529b480c912ce5ca36714c7bc17b..ab192f387231dba8c5a1802f8d503303852c7b13 100644
--- a/tests/test_checkpointing_IBoundaryConditionDescriptor.cpp
+++ b/tests/test_checkpointing_IBoundaryConditionDescriptor.cpp
@@ -9,6 +9,7 @@
#include <mesh/NumberedBoundaryDescriptor.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DirichletVectorBoundaryConditionDescriptor.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
@@ -179,6 +180,17 @@ let i: R -> R, x -> x+3;
p_dirichlet_bc_descriptor)},
file, useless_group, symbol_table_group);
+ const std::vector<FunctionSymbolId> dirichlet_vector_function_id_list{FunctionSymbolId{2, symbol_table},
+ FunctionSymbolId{3, symbol_table}};
+ auto p_dirichlet_vector_bc_descriptor =
+ std::make_shared<const DirichletVectorBoundaryConditionDescriptor>("dirichlet_vector_name", p_boundary_1,
+ dirichlet_vector_function_id_list);
+ checkpointing::writeIBoundaryConditionDescriptor("dirichlet_vector_bc_descriptor",
+ EmbeddedData{std::make_shared<
+ DataHandler<const IBoundaryConditionDescriptor>>(
+ p_dirichlet_vector_bc_descriptor)},
+ file, useless_group, symbol_table_group);
+
const FunctionSymbolId neumann_function_id{1, symbol_table};
auto p_neumann_bc_descriptor =
std::make_shared<const NeumannBoundaryConditionDescriptor>("neumann_name", p_boundary_1, neumann_function_id);
@@ -245,6 +257,9 @@ let i: R -> R, x -> x+3;
EmbeddedData read_dirichlet_bc_descriptor =
checkpointing::readIBoundaryConditionDescriptor("dirichlet_bc_descriptor", symbol_table_group);
+ EmbeddedData read_dirichlet_vector_bc_descriptor =
+ checkpointing::readIBoundaryConditionDescriptor("dirichlet_vector_bc_descriptor", symbol_table_group);
+
EmbeddedData read_neumann_bc_descriptor =
checkpointing::readIBoundaryConditionDescriptor("neumann_bc_descriptor", symbol_table_group);
@@ -268,6 +283,7 @@ let i: R -> R, x -> x+3;
REQUIRE_NOTHROW(get_value(read_free_bc_descriptor));
REQUIRE_NOTHROW(get_value(read_fixed_bc_descriptor));
REQUIRE_NOTHROW(get_value(read_dirichlet_bc_descriptor));
+ REQUIRE_NOTHROW(get_value(read_dirichlet_vector_bc_descriptor));
REQUIRE_NOTHROW(get_value(read_neumann_bc_descriptor));
REQUIRE_NOTHROW(get_value(read_fourier_bc_descriptor));
REQUIRE_NOTHROW(get_value(read_inflow_bc_descriptor));
@@ -280,6 +296,8 @@ let i: R -> R, x -> x+3;
REQUIRE(get_value(read_free_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::free);
REQUIRE(get_value(read_fixed_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::fixed);
REQUIRE(get_value(read_dirichlet_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::dirichlet);
+ REQUIRE(get_value(read_dirichlet_vector_bc_descriptor).type() ==
+ IBoundaryConditionDescriptor::Type::dirichlet_vector);
REQUIRE(get_value(read_neumann_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::neumann);
REQUIRE(get_value(read_fourier_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::fourier);
REQUIRE(get_value(read_inflow_bc_descriptor).type() == IBoundaryConditionDescriptor::Type::inflow);
@@ -293,6 +311,8 @@ let i: R -> R, x -> x+3;
REQUIRE_NOTHROW(dynamic_cast<const FixedBoundaryConditionDescriptor&>(get_value(read_fixed_bc_descriptor)));
REQUIRE_NOTHROW(
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(get_value(read_dirichlet_bc_descriptor)));
+ REQUIRE_NOTHROW(dynamic_cast<const DirichletVectorBoundaryConditionDescriptor&>(
+ get_value(read_dirichlet_vector_bc_descriptor)));
REQUIRE_NOTHROW(dynamic_cast<const NeumannBoundaryConditionDescriptor&>(get_value(read_neumann_bc_descriptor)));
REQUIRE_NOTHROW(dynamic_cast<const FourierBoundaryConditionDescriptor&>(get_value(read_fourier_bc_descriptor)));
REQUIRE_NOTHROW(dynamic_cast<const InflowBoundaryConditionDescriptor&>(get_value(read_inflow_bc_descriptor)));
@@ -308,6 +328,8 @@ let i: R -> R, x -> x+3;
auto& read_fixed_bc = dynamic_cast<const FixedBoundaryConditionDescriptor&>(get_value(read_fixed_bc_descriptor));
auto& read_dirichlet_bc =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(get_value(read_dirichlet_bc_descriptor));
+ auto& read_dirichlet_vector_bc =
+ dynamic_cast<const DirichletVectorBoundaryConditionDescriptor&>(get_value(read_dirichlet_vector_bc_descriptor));
auto& read_neumann_bc =
dynamic_cast<const NeumannBoundaryConditionDescriptor&>(get_value(read_neumann_bc_descriptor));
auto& read_fourier_bc =
@@ -326,6 +348,15 @@ let i: R -> R, x -> x+3;
REQUIRE(read_dirichlet_bc.boundaryDescriptor().type() == p_dirichlet_bc_descriptor->boundaryDescriptor().type());
REQUIRE(read_dirichlet_bc.name() == p_dirichlet_bc_descriptor->name());
REQUIRE(read_dirichlet_bc.rhsSymbolId().id() == p_dirichlet_bc_descriptor->rhsSymbolId().id());
+ REQUIRE(read_dirichlet_vector_bc.boundaryDescriptor().type() ==
+ p_dirichlet_vector_bc_descriptor->boundaryDescriptor().type());
+ REQUIRE(read_dirichlet_vector_bc.name() == p_dirichlet_vector_bc_descriptor->name());
+ REQUIRE(read_dirichlet_vector_bc.rhsSymbolIdList().size() ==
+ p_dirichlet_vector_bc_descriptor->rhsSymbolIdList().size());
+ for (size_t i = 0; i < read_dirichlet_vector_bc.rhsSymbolIdList().size(); ++i) {
+ REQUIRE(read_dirichlet_vector_bc.rhsSymbolIdList()[i].id() ==
+ p_dirichlet_vector_bc_descriptor->rhsSymbolIdList()[i].id());
+ }
REQUIRE(read_neumann_bc.boundaryDescriptor().type() == p_neumann_bc_descriptor->boundaryDescriptor().type());
REQUIRE(read_neumann_bc.name() == p_neumann_bc_descriptor->name());
REQUIRE(read_neumann_bc.rhsSymbolId().id() == p_neumann_bc_descriptor->rhsSymbolId().id());
@@ -337,6 +368,8 @@ let i: R -> R, x -> x+3;
REQUIRE(read_inflow_bc.functionSymbolId().id() == p_inflow_bc_descriptor->functionSymbolId().id());
REQUIRE(read_inflow_list_bc.boundaryDescriptor().type() ==
p_inflow_list_bc_descriptor->boundaryDescriptor().type());
+ REQUIRE(read_inflow_list_bc.functionSymbolIdList().size() ==
+ p_inflow_list_bc_descriptor->functionSymbolIdList().size());
for (size_t i = 0; i < read_inflow_list_bc.functionSymbolIdList().size(); ++i) {
REQUIRE(read_inflow_list_bc.functionSymbolIdList()[i].id() ==
p_inflow_list_bc_descriptor->functionSymbolIdList()[i].id());
diff --git a/tools/generate-plugin.sh b/tools/generate-plugin.sh
new file mode 100755
index 0000000000000000000000000000000000000000..275d1795f79753fbd259a2079119d977c89b7d82
--- /dev/null
+++ b/tools/generate-plugin.sh
@@ -0,0 +1,95 @@
+#!/bin/bash
+
+BOLD='\e[1m'
+RESET='\e[0m'
+
+GREEN='\e[92m'
+RED='\e[91m'
+YELLOW='\e[93m'
+
+echo -ne ${BOLD}
+echo -e "---------------------"
+echo -e "pugs plugin generator"
+echo -e "---------------------"
+echo -e ${RESET}
+
+CURRENT_DIR="$(pwd -P)"
+SCRIPT_DIR="$( cd -- "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"
+PUGS_DIR="$(dirname ${SCRIPT_DIR})"
+
+if [[ "${CURRENT_DIR}" =~ "${PUGS_DIR}" ]]
+then
+ echo -e ${RED}"Aborting..."${RESET}
+ echo -e "run this script outside of pugs sources"
+ exit 1
+fi
+
+NAME_RE='^[A-Z][a-zA-Z0-9]*$'
+
+echo " Plugin name must fulfill the following constrains:"
+echo " - be a single word that starts by an upper case,"
+echo " - contains only letters or numbers,"
+echo " and preferably separate words with caps."
+echo
+echo " ex.: MyFirstPlugin"
+echo
+
+while [[ ! "${PLUGIN_NAME}" =~ $NAME_RE ]]
+do
+ echo -n "Give plugin name: "
+ read -r PLUGIN_NAME
+
+ if [[ ! "${PLUGIN_NAME}" =~ $NAME_RE ]]
+ then
+ echo -e ${RED}" invalid name!"${RESET}
+ echo
+ unset PLUGIN_NAME
+ fi
+
+done
+
+PLUGIN_UP="${PLUGIN_NAME^^}"
+PLUGIN_LOW="${PLUGIN_NAME,,}"
+echo
+echo -e "creating plugin ${YELLOW}${PLUGIN_NAME}${RESET} in directory ${YELLOW}${PLUGIN_LOW}${RESET}"
+echo
+
+if [[ -e ${PLUGIN_LOW} ]]
+then
+ echo -e ${RED}"Aborting..."${RESET}
+ echo -e "directory \"${PLUGIN_LOW}\" ${YELLOW}already exists${RESET}!"
+ exit 1
+fi
+
+function substitute()
+{
+ sed s/_PLUGIN_NAME_/${PLUGIN_NAME}/g | sed s/_PLUGIN_LOW_/${PLUGIN_LOW}/g | sed s/_PLUGIN_UP_/${PLUGIN_UP}/g
+}
+
+mkdir -p "${PLUGIN_LOW}/cmake"
+mkdir -p "${PLUGIN_LOW}/tests"
+
+cp "${PUGS_DIR}"/tests/MeshDataBaseForTests.hpp "${PLUGIN_LOW}"/tests/
+cp "${PUGS_DIR}"/tests/MeshDataBaseForTests.cpp "${PLUGIN_LOW}"/tests/
+cp "${PUGS_DIR}"/tests/ParallelCheckerTester.hpp "${PLUGIN_LOW}"/tests/
+cp "${PUGS_DIR}"/tests/ParallelCheckerTester.cpp "${PLUGIN_LOW}"/tests/
+cp "${PUGS_DIR}"/tests/test_main.cpp "${PLUGIN_LOW}"/tests/
+cp "${PUGS_DIR}"/tests/mpi_test_main.cpp "${PLUGIN_LOW}"/tests/
+
+cp "${PUGS_DIR}"/cmake/CheckNotInSources.cmake "${PLUGIN_LOW}"/cmake/
+cp "${PUGS_DIR}"/tools/plugin-template/FindPugs.cmake "${PLUGIN_LOW}"/cmake/
+cp "${PUGS_DIR}"/.gitignore "${PLUGIN_LOW}"
+cp "${PUGS_DIR}"/.clang-format "${PLUGIN_LOW}"
+
+cat "${PUGS_DIR}"/tools/plugin-template/CMakeLists.txt-template | substitute > "${PLUGIN_LOW}"/CMakeLists.txt
+cat "${PUGS_DIR}"/tools/plugin-template/Module.hpp-template | substitute > "${PLUGIN_LOW}"/${PLUGIN_NAME}Module.hpp
+cat "${PUGS_DIR}"/tools/plugin-template/Module.cpp-template | substitute > "${PLUGIN_LOW}"/${PLUGIN_NAME}Module.cpp
+cat "${PUGS_DIR}"/tools/plugin-template/README.md-template | substitute > "${PLUGIN_LOW}"/README.md
+
+cat "${PUGS_DIR}"/tools/plugin-template/tests-CMakeLists.txt-template | substitute > "${PLUGIN_LOW}"/tests/CMakeLists.txt
+
+(cd "${PLUGIN_LOW}"; git init -q)
+(cd "${PLUGIN_LOW}"; git add .)
+(cd "${PLUGIN_LOW}"; git commit -m "init" -q)
+
+echo -e ${GREEN}"Creation finished successfully!"${RESET}
diff --git a/tools/plugin-template/CMakeLists.txt-template b/tools/plugin-template/CMakeLists.txt-template
new file mode 100644
index 0000000000000000000000000000000000000000..9c6478b0711cbbd5e38a41f3eaa6c4736a300e67
--- /dev/null
+++ b/tools/plugin-template/CMakeLists.txt-template
@@ -0,0 +1,152 @@
+cmake_minimum_required (VERSION 3.19)
+
+project("_PLUGIN_LOW_")
+
+set(PUGS_PREFIX_PATH "" CACHE STRING "pugs intall dir")
+
+string(COMPARE EQUAL "${PUGS_PREFIX_PATH}" "" pugs_prefix_undefined)
+
+if (pugs_prefix_undefined)
+ message(FATAL_ERROR "PUGS_PREFIX_PATH must be defined")
+endif()
+
+# CMake utils
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
+
+# Forbids in-source builds
+include(CheckNotInSources)
+
+# use PkgConfig to find packages
+find_package(PkgConfig REQUIRED)
+find_package(Pugs REQUIRED)
+
+list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/Kokkos")
+include(KokkosConfig)
+
+set(HDF5_PREFER_PARALLEL TRUE)
+list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/HighFive")
+include(HighFiveConfig)
+
+list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/pugs")
+include(PugsTargets)
+include(PugsCompileFlags)
+
+#------------------------------------------------------
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${PUGS_CMAKE_CXX_FLAGS}")
+set(CMAKE_CXX_STANDARD "${PUGS_CMAKE_CXX_STANDARD}")
+
+# -----------------------------------------------------
+# dynamic libraries
+
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+link_libraries("-rdynamic")
+set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+
+#------------------------------------------------------
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+#------------------------------------------------------
+
+set(_PLUGIN_UP__SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}")
+set(_PLUGIN_UP__BINARY_DIR "${CMAKE_CURRENT_BINARY_DIR}")
+
+#------------------------------------------------------
+# use pugs compilation settings
+
+set(CMAKE_BUILD_TYPE "${PUGS_CMAKE_BUILD_TYPE}" CACHE STRING "" FORCE)
+set(CMAKE_CXX_COMPILER "${PUGS_CMAKE_CXX_COMPILER}" CACHE STRING "" FORCE)
+set(CMAKE_C_COMPILER "${PUGS_CMAKE_C_COMPILER}" CACHE STRING "" FORCE)
+
+mark_as_advanced(CMAKE_BUILD_TYPE CMAKE_CXX_COMPILER CMAKE_C_COMPILER)
+
+#------------------------------------------------------
+# default build shared libraries
+set(BUILD_SHARED_LIBS ON CACHE STRING "" FORCE)
+
+#------------------------------------------------------
+
+# Checks if compiler version is compatible with Pugs sources
+set(GNU_CXX_MIN_VERSION "10.0.0")
+set(CLANG_CXX_MIN_VERSION "11.0.0")
+
+#------------------------------------------------------
+
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}")
+include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include")
+include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include/kokkos")
+include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include/tao/")
+include_directories(SYSTEM "${MPI_CXX_INCLUDE_DIRS}")
+
+get_target_property(_prop Pugs::PugsAlgebra INTERFACE_INCLUDE_DIRECTORIES)
+set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
+get_target_property(_prop Pugs::PugsUtils INTERFACE_INCLUDE_DIRECTORIES)
+set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
+get_target_property(_prop Pugs::pugs INTERFACE_INCLUDE_DIRECTORIES)
+set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
+
+include_directories(SYSTEM ${PUGS_INC_DIR})
+link_directories(${PUGS_PREFIX_PATH}/lib)
+
+#------------------------------------------------------
+
+if(${PUGS_HAS_MPI})
+ set(MPIEXEC_OPTION_FLAGS --oversubscribe)
+ if (NOT "$ENV{GITLAB_CI}" STREQUAL "")
+ set(MPIEXEC_OPTION_FLAGS ${MPIEXEC_OPTION_FLAGS} --allow-run-as-root)
+ endif()
+ set(MPIEXEC_NUMPROC 3)
+ set(MPIEXEC_PATH_FLAG --path)
+ set(MPI_LAUNCHER ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_NUMPROC} ${MPIEXEC_OPTION_FLAGS})
+endif()
+
+add_custom_target(all_unit_tests
+ DEPENDS unit_tests mpi_unit_tests
+)
+
+add_custom_target(check
+ DEPENDS test
+ )
+
+add_custom_target(test
+ DEPENDS run_all_unit_tests
+ )
+
+add_custom_target(run_all_unit_tests
+ DEPENDS run_mpi_unit_tests
+ )
+
+if(PUGS_HAS_MPI)
+ set(RUN_MPI_UNIT_TESTS_COMMENT "Running mpi_unit_tests [using ${MPIEXEC_NUMPROC} procs]")
+else()
+ set(RUN_MPI_UNIT_TESTS_COMMENT "Running mpi_unit_tests [sequentially]")
+endif()
+
+add_custom_target(run_mpi_unit_tests
+ COMMAND ${MPI_LAUNCHER} "${_PLUGIN_UP__BINARY_DIR}/mpi_unit_tests" --allow-running-no-tests
+ WORKING_DIRECTORY ${_PLUGIN_UP__BINARY_DIR}
+ DEPENDS run_unit_tests
+ COMMENT ${RUN_MPI_UNIT_TESTS_COMMENT}
+ )
+
+add_custom_target(run_unit_tests
+ COMMAND "${_PLUGIN_UP__BINARY_DIR}/unit_tests" --allow-running-no-tests
+ DEPENDS all_unit_tests
+ COMMENT "Running unit_tests"
+ )
+
+#------------------------------------------------------
+
+add_library(_PLUGIN_NAME_
+ _PLUGIN_NAME_Module.cpp
+ # add cpp sources files here
+)
+
+#------------------------------------------------------
+
+add_subdirectory(tests)
+
+#------------------------------------------------------
+
+install(TARGETS _PLUGIN_NAME_)
diff --git a/tools/plugin-template/FindPugs.cmake b/tools/plugin-template/FindPugs.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..6269bed8a9055339894cdd9fa9e63b2b479c9b70
--- /dev/null
+++ b/tools/plugin-template/FindPugs.cmake
@@ -0,0 +1,10 @@
+# Finds for the pugs installation directory
+
+find_path(PUGS_PREFIX_PATH include/utils/pugs_version.hpp
+ HINTS
+ $ENV{PUGS_INSTALL_DIR}
+ /usr/local/pugs
+ NO_DEFAULT_PATH
+)
+
+find_package_handle_standard_args(Pugs REQUIRED_VARS PUGS_PREFIX_PATH )
diff --git a/tools/plugin-template/Module.cpp-template b/tools/plugin-template/Module.cpp-template
new file mode 100644
index 0000000000000000000000000000000000000000..dd264e03631590e3cf0db152779e2118df86d474
--- /dev/null
+++ b/tools/plugin-template/Module.cpp-template
@@ -0,0 +1,28 @@
+#include <_PLUGIN_NAME_Module.hpp>
+
+#include <language/modules/ModuleRepository.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+
+_PLUGIN_NAME_Module::_PLUGIN_NAME_Module() : BuiltinModule(false)
+{
+ // Simple hello world example
+ this->_addBuiltinFunction("_PLUGIN_LOW__hello", std::function(
+
+ []() -> void { std::cout << "_PLUGIN_NAME_: Hello world\n"; }
+
+ ));
+}
+
+void
+_PLUGIN_NAME_Module::registerOperators() const
+{
+ // Kept empty for basic use
+}
+
+void
+_PLUGIN_NAME_Module::registerCheckpointResume() const
+{
+ // kept empty for basic use
+}
+
+ModuleRepository::Subscribe<_PLUGIN_NAME_Module> _PLUGIN_LOW__module;
diff --git a/tools/plugin-template/Module.hpp-template b/tools/plugin-template/Module.hpp-template
new file mode 100644
index 0000000000000000000000000000000000000000..757e3b5e650390251d15515cfdcf26a9c970a102
--- /dev/null
+++ b/tools/plugin-template/Module.hpp-template
@@ -0,0 +1,22 @@
+#ifndef _PLUGIN_UP__MODULE_HPP
+#define _PLUGIN_UP__MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+
+class _PLUGIN_NAME_Module : public BuiltinModule
+{
+ public:
+ std::string_view
+ name() const final
+ {
+ return "_PLUGIN_LOW_";
+ }
+
+ void registerOperators() const final;
+ void registerCheckpointResume() const final;
+
+ _PLUGIN_NAME_Module();
+ ~_PLUGIN_NAME_Module() = default;
+};
+
+#endif // _PLUGIN_UP__MODULE_HPP
diff --git a/tools/plugin-template/README.md-template b/tools/plugin-template/README.md-template
new file mode 100644
index 0000000000000000000000000000000000000000..4db1071714e56642f3114dad4bc3451bbb274b56
--- /dev/null
+++ b/tools/plugin-template/README.md-template
@@ -0,0 +1,67 @@
+`pugs`'s plugin `_PLUGIN_NAME_`
+===============================
+
+# Building `_PLUGIN_NAME_`
+
+## `pugs` installation
+
+Building this plugin requires an **installed** version of `pugs`.
+`pugs` follows standard `cmake` installation recipes.
+
+Before building `pugs` one should define its installation directory.
+In the `pugs` compilation directory one should execute
+```shell
+cmake -DCMAKE_INSTALL_PREFIX=pugs_install_dir ...
+```
+where `pugs_install_dir` is the chosen installation directory.
+
+Then one simply runs
+```shell
+make install
+```
+
+## building the plugin `_PLUGIN_NAME_`
+
+> **Warning**:<br>
+> Building `_PLUGIN_NAME_` in its source directory is
+> **forbidden**. Trying to do so will result in a failure. However it
+> generally leaves some garbage files in the source directory, namely
+> the `CMakeCache.txt` and the `CMakeFiles` directory. `CMake` itself
+> is not able to remove them, to avoid the risk of compilation issues,
+> one has to dot it manually...
+
+In the build directory one runs
+```shell
+cmake -DPUGS_PREFIX_PATH=pugs_install_dir _PLUGIN_LOW__dir
+```
+where `pugs_install_dir` has the same value as above and `_PLUGIN_LOW__dir`
+is the directory that contains this `README.md` file.
+
+Then to build the plugin, one runs
+```shell
+make
+```
+
+If anything runs fine, the dynamic library `lib_PLUGIN_NAME_.so` is
+built.
+
+# Using `_PLUGIN_NAME_`
+
+In order to use the created plugin, one simply has to give the
+location of `lib_PLUGIN_NAME_.so` to `pugs`. This is done by means of
+environment variables. There are two possibilities:
+- `PUGS_PLUGIN` contains a semicolumn separated list of plugin
+ libraries,
+- `PUGS_PLUGIN_DIR` contains a semicolumn separated list of path to
+ plugin libraries.
+
+Example
+```shell
+export PUGS_PLUGIN="/pathtoplugin/lib_PLUGIN_NAME_.so"
+```
+or
+```shell
+export PUGS_PLUGIN_DIR="/pathtoplugin1;/pathtoplugin2"
+```
+
+Then one launches `pugs` classically.
diff --git a/tools/plugin-template/tests-CMakeLists.txt-template b/tools/plugin-template/tests-CMakeLists.txt-template
new file mode 100644
index 0000000000000000000000000000000000000000..026ff2c2962ed988d35e2680b8b97919801c7a4b
--- /dev/null
+++ b/tools/plugin-template/tests-CMakeLists.txt-template
@@ -0,0 +1,95 @@
+set(EXECUTABLE_OUTPUT_PATH ${_PLUGIN_UP__BINARY_DIR})
+
+include_directories(${_PLUGIN_UP__SOURCE_DIR}/tests)
+
+add_executable (unit_tests
+ test_main.cpp
+
+ # add unit tests here
+ )
+
+ set(_PLUGIN_UP__checkpointing_TESTS
+ )
+
+if(PUGS_HAS_HDF5)
+ list(APPEND _PLUGIN_UP__checkpointing_TESTS
+ )
+endif(PUGS_HAS_HDF5)
+
+add_executable (mpi_unit_tests
+ mpi_test_main.cpp
+ ${_PLUGIN_UP__checkpointing_TESTS}
+
+ # add mpi unit tests here
+)
+
+add_library(test_MeshDataBase
+ MeshDataBaseForTests.cpp)
+
+add_library(test_ParallelCheckerTester
+ ParallelCheckerTester.cpp)
+
+target_link_libraries (test_ParallelCheckerTester
+ ${HIGHFIVE_TARGET})
+
+target_link_libraries (unit_tests
+ test_MeshDataBase
+ test_ParallelCheckerTester
+ PugsAlgebra
+ PugsAnalysis
+ PugsUtils
+ PugsLanguage
+ PugsLanguageAST
+ PugsLanguageModules
+ PugsMesh
+ PugsAlgebra
+ PugsUtils
+ PugsLanguageUtils
+ PugsScheme
+ PugsOutput
+ PugsUtils
+ PugsCheckpointing
+ PugsDev
+ PugsAlgebra
+ PugsMesh
+ Kokkos::kokkos
+ ${PARMETIS_LIBRARIES}
+ ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+ ${PETSC_LIBRARIES}
+ Catch2
+ ${PUGS_STD_LINK_FLAGS}
+ ${HIGHFIVE_TARGET}
+ ${SLURM_LIBRARY}
+ stdc++fs
+ )
+
+target_link_libraries (mpi_unit_tests
+ test_MeshDataBase
+ test_ParallelCheckerTester
+ PugsAlgebra
+ PugsAnalysis
+ PugsUtils
+ PugsLanguage
+ PugsLanguageAST
+ PugsLanguageModules
+ PugsMesh
+ PugsAlgebra
+ PugsUtils
+ PugsLanguageUtils
+ PugsScheme
+ PugsOutput
+ PugsUtils
+ PugsCheckpointing
+ PugsDev
+ PugsAlgebra
+ PugsMesh
+ Kokkos::kokkos
+ ${PARMETIS_LIBRARIES}
+ ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+ ${PETSC_LIBRARIES}
+ Catch2
+ ${PUGS_STD_LINK_FLAGS}
+ ${HIGHFIVE_TARGET}
+ ${SLURM_LIBRARY}
+ stdc++fs
+ )