diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
index 406c38c7e3abf8ab6ce3eafdec4ac14cc97aa132..8f18986a8a928cdd538eb546dba30cc2b78fe1e0 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -17,7 +17,7 @@
#include "MeshLib/MeshSearch/NodeSearch.h" // for getUniqueNodes
#include "MeshLib/Node.h"
#include "ParameterLib/Utils.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
namespace ProcessLib
{
@@ -100,7 +100,7 @@ ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
const int shape_function_order = 1;
- ProcessLib::createLocalAssemblers<
+ BoundaryConditionOrSourceTerm::createLocalAssemblers<
ConstraintDirichletBoundaryConditionLocalAssembler>(
_bulk_mesh.getDimension(), _bc_mesh.getElements(), *_dof_table_boundary,
shape_function_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 0c91154d7c2f9a6d7de8e3e456c52c474642dfaf..dfd07225c6f4d23179df06a81cdd489a44fb48d9 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -9,7 +9,7 @@
*/
#include "GenericNaturalBoundaryConditionLocalAssembler.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
namespace ProcessLib
{
@@ -65,7 +65,8 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
_dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
variable_id, {component_id}, std::move(bc_mesh_subset)));
- createLocalAssemblers<LocalAssemblerImplementation>(
+ BoundaryConditionOrSourceTerm::createLocalAssemblers<
+ LocalAssemblerImplementation>(
global_dim, _bc_mesh.getElements(), *_dof_table_boundary,
shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
integration_order, _data);
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
index 28706ab75f2bb66b4c686ca9134849bd4bdd2d83..460fdc8d2f70957e62504e8c28ca08d7431d6665 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
@@ -13,10 +13,9 @@
#include <numeric>
#include "MeshLib/MeshSearch/NodeSearch.h"
-#include "ParameterLib/Utils.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
-
#include "NormalTractionBoundaryConditionLocalAssembler.h"
+#include "ParameterLib/Utils.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
namespace ProcessLib
{
@@ -52,7 +51,8 @@ NormalTractionBoundaryCondition<GlobalDim, LocalAssemblerImplementation>::
_dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
variable_id, component_ids, std::move(bc_mesh_subset)));
- createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
+ BoundaryConditionOrSourceTerm::detail::createLocalAssemblers<
+ GlobalDim, LocalAssemblerImplementation>(
*_dof_table_boundary, shapefunction_order, _bc_mesh.getElements(),
_local_assemblers, _bc_mesh.isAxiallySymmetric(), _integration_order,
_pressure);
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
index 1ba9c7ecad38e2ce295978192dcc288fe9d47f0c..deed7068845ef7cdd7c238123797804d96588748 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
@@ -16,7 +16,7 @@
#include "BaseLib/ConfigTree.h"
#include "MeshLib/MeshSearch/NodeSearch.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
#include "PythonBoundaryConditionLocalAssembler.h"
namespace
@@ -75,7 +75,8 @@ PythonBoundaryCondition::PythonBoundaryCondition(
_dof_table_boundary = _bc_data.dof_table_bulk.deriveBoundaryConstrainedMap(
std::move(bc_mesh_subset));
- createLocalAssemblers<PythonBoundaryConditionLocalAssembler>(
+ BoundaryConditionOrSourceTerm::createLocalAssemblers<
+ PythonBoundaryConditionLocalAssembler>(
global_dim, _bc_data.boundary_mesh.getElements(), *_dof_table_boundary,
shapefunction_order, _local_assemblers,
_bc_data.boundary_mesh.isAxiallySymmetric(), integration_order,
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 2c417a1d6d8e1aaa246deb9b984e299074daec8a..f1b2e308d5da780a973698c7bfa3bc245e2caea9 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -17,6 +17,7 @@ append_source_files(SOURCES SurfaceFlux)
append_source_files(SOURCES Output)
append_source_files(SOURCES SourceTerms)
append_source_files(SOURCES Utils)
+append_source_files(SOURCES Utils/ForBoundaryConditionOrSourceTerm)
ogs_add_library(ProcessLib ${SOURCES})
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
index 35b4bbc35fa7ba89d60e52854a07c2f097a108fb..12cb687bd2087dda6850bbbdff03c77f8964c5ef 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
@@ -16,7 +16,7 @@
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
#include "PythonSourceTermLocalAssembler.h"
namespace
@@ -71,7 +71,8 @@ PythonSourceTerm::PythonSourceTerm(
_source_term_data(std::move(source_term_data)),
_flush_stdout(flush_stdout)
{
- createLocalAssemblers<PythonSourceTermLocalAssembler>(
+ BoundaryConditionOrSourceTerm::createLocalAssemblers<
+ PythonSourceTermLocalAssembler>(
global_dim, _source_term_data.source_term_mesh.getElements(),
*_source_term_dof_table, shapefunction_order, _local_assemblers,
_source_term_data.source_term_mesh.isAxiallySymmetric(),
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
index 6901d7537c80e74beacb3b324772d3c2e847b254..62be8d56432bdb5a0d6684e8067e50cd633a9348 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
@@ -10,7 +10,7 @@
#include "VolumetricSourceTerm.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h"
namespace ProcessLib
{
@@ -22,7 +22,8 @@ VolumetricSourceTerm::VolumetricSourceTerm(
: SourceTerm(std::move(source_term_dof_table)),
_source_term_parameter(source_term_parameter)
{
- ProcessLib::createLocalAssemblers<VolumetricSourceTermLocalAssembler>(
+ BoundaryConditionOrSourceTerm::createLocalAssemblers<
+ VolumetricSourceTermLocalAssembler>(
bulk_mesh_dimension, source_term_mesh.getElements(),
*_source_term_dof_table, shapefunction_order, _local_assemblers,
source_term_mesh.isAxiallySymmetric(), integration_order,
diff --git a/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h b/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h
new file mode 100644
index 0000000000000000000000000000000000000000..b12f7808662f9126d5b55ec30c3c6e4b4abe150d
--- /dev/null
+++ b/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h
@@ -0,0 +1,101 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "BaseLib/Logging.h"
+#include "LocalDataInitializer.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace ProcessLib
+{
+namespace BoundaryConditionOrSourceTerm
+{
+namespace detail
+{
+template <int GlobalDim,
+ template <typename, typename, int> class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ static_assert(
+ GlobalDim == 1 || GlobalDim == 2 || GlobalDim == 3,
+ "Meshes with dimension greater than three are not supported.");
+ // Shape matrices initializer
+ using LocalDataInitializer =
+ LocalDataInitializer<LocalAssemblerInterface,
+ LocalAssemblerImplementation, GlobalDim,
+ ExtraCtorArgs...>;
+
+ DBUG("Create local assemblers.");
+ // Populate the vector of local assemblers.
+ local_assemblers.resize(mesh_elements.size());
+
+ LocalDataInitializer initializer(dof_table, shapefunction_order);
+
+ DBUG("Calling local assembler builder for all mesh elements.");
+ GlobalExecutor::transformDereferenced(
+ initializer, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+} // namespace detail
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ * Those arguments will be passed to the constructor of
+ * \c LocalAssemblerImplementation.
+ *
+ * The first two template parameters cannot be deduced from the arguments.
+ * Therefore they always have to be provided manually.
+ */
+template <template <typename, typename, int> class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ const unsigned dimension,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ DBUG("Create local assemblers.");
+
+ switch (dimension)
+ {
+ case 1:
+ detail::createLocalAssemblers<1, LocalAssemblerImplementation>(
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ break;
+ case 2:
+ detail::createLocalAssemblers<2, LocalAssemblerImplementation>(
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ break;
+ case 3:
+ detail::createLocalAssemblers<3, LocalAssemblerImplementation>(
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ break;
+ default:
+ OGS_FATAL(
+ "Meshes with dimension greater than three are not supported.");
+ }
+}
+} // namespace BoundaryConditionOrSourceTerm
+} // namespace ProcessLib
diff --git a/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/LocalDataInitializer.h b/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/LocalDataInitializer.h
new file mode 100644
index 0000000000000000000000000000000000000000..032af2890b6b911177875cd52b2d48169a7bd6d5
--- /dev/null
+++ b/ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/LocalDataInitializer.h
@@ -0,0 +1,401 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+
+#ifndef OGS_MAX_ELEMENT_DIM
+static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
+#endif
+
+#ifndef OGS_MAX_ELEMENT_ORDER
+static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
+#endif
+
+// The following macros decide which element types will be compiled, i.e.
+// which element types will be available for use in simulations.
+
+#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
+#else
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_CUBOID
+#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
+#else
+#define ENABLED_ELEMENT_TYPE_CUBOID 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PRISM
+#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
+#else
+#define ENABLED_ELEMENT_TYPE_PRISM 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PYRAMID
+#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
+#else
+#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
+#endif
+
+// Dependent element types.
+// Faces of tets, pyramids and prisms are triangles
+#define ENABLED_ELEMENT_TYPE_TRI \
+ ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+ (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of hexes, pyramids and prisms are quads
+#define ENABLED_ELEMENT_TYPE_QUAD \
+ ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+ (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of triangles and quads are lines
+#define ENABLED_ELEMENT_TYPE_LINE \
+ ((ENABLED_ELEMENT_TYPE_TRI) | (ENABLED_ELEMENT_TYPE_QUAD))
+
+// All enabled element types
+#define OGS_ENABLED_ELEMENTS \
+ ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
+ (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
+
+// Include only what is needed (Well, the conditions are not sharp).
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
+#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
+#endif
+
+namespace ProcessLib
+{
+namespace BoundaryConditionOrSourceTerm
+{
+/// The LocalDataInitializer is a functor creating a local assembler data with
+/// corresponding to the mesh element type shape functions and calling
+/// initialization of the new local assembler data.
+/// For example for MeshLib::Quad a local assembler data with template argument
+/// NumLib::ShapeQuad4 is created.
+template <typename LocalAssemblerInterface,
+ template <typename, typename, int> class LocalAssemblerData,
+ int GlobalDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+ using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+ LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order)
+ : _dof_table(dof_table)
+ {
+ if (shapefunction_order < 1 || 2 < shapefunction_order)
+ OGS_FATAL("The given shape function order {:d} is not supported",
+ shapefunction_order);
+
+ if (shapefunction_order == 1)
+ {
+ // /// Lines and points ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 0 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Point))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePoint1>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Line))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine2>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Line3))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine2>();
+#endif
+
+ // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Quad))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Hex))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex8>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Quad8))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
+ _builder[std::type_index(typeid(MeshLib::Quad9))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Hex20))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex8>();
+#endif
+
+ // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Tri))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri3>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Tet))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet4>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tri6))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri3>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tet10))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet4>();
+#endif
+
+ // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Prism))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Prism15))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism6>();
+#endif
+
+ // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Pyramid))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra5>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra5>();
+#endif
+ }
+ else if (shapefunction_order == 2)
+ {
+// /// Lines and points ///////////////////////////////////
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 0 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Point))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePoint1>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Line3))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
+#endif
+
+ // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Quad8))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+ _builder[std::type_index(typeid(MeshLib::Quad9))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Hex20))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+ // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tri6))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tet10))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+ // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Prism15))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+ // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
+#endif
+ }
+ }
+
+ /// Returns data pointer to the newly created local assembler data.
+ ///
+ /// \attention
+ /// The index \c id is not necessarily the mesh item's id. Especially when
+ /// having multiple meshes it will differ from the latter.
+ LADataIntfPtr operator()(std::size_t const id,
+ MeshLib::Element const& mesh_item,
+ ConstructorArgs&&... args) const
+ {
+ auto const type_idx = std::type_index(typeid(mesh_item));
+ auto const it = _builder.find(type_idx);
+
+ if (it != _builder.end())
+ {
+ auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
+ return it->second(mesh_item, num_local_dof,
+ std::forward<ConstructorArgs>(args)...);
+ }
+ OGS_FATAL(
+ "You are trying to build a local assembler for an unknown mesh "
+ "element type ({:s})."
+ " Maybe you have disabled this mesh element type in your build "
+ "configuration.",
+ type_idx.name());
+ }
+
+private:
+ using LADataBuilder =
+ std::function<LADataIntfPtr(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ ConstructorArgs&&...)>;
+
+ template <typename ShapeFunction>
+ using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+ typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+ template <typename ShapeFunction>
+ using LAData =
+ LocalAssemblerData<ShapeFunction, IntegrationMethod<ShapeFunction>,
+ GlobalDim>;
+
+ /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+ /// depending whether the global dimension is less than the shape function's
+ /// dimension or not.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder()
+ {
+ return makeLocalAssemblerBuilder<ShapeFunction>(
+ static_cast<std::integral_constant<
+ bool, (GlobalDim >= ShapeFunction::DIM)>*>(nullptr));
+ }
+
+ /// Mapping of element types to local assembler constructors.
+ std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+ NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+ // local assembler builder implementations.
+private:
+ /// Generates a function that creates a new LocalAssembler of type
+ /// LAData<ShapeFunction>. Only functions with shape function's dimension
+ /// less or equal to the global dimension are instantiated, e.g. following
+ /// combinations of shape functions and global dimensions: (Line2, 1),
+ /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder(std::true_type* /*unused*/)
+ {
+ return [](MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ ConstructorArgs&&... args)
+ {
+ return LADataIntfPtr{new LAData<ShapeFunction>{
+ e, local_matrix_size, std::forward<ConstructorArgs>(args)...}};
+ };
+ }
+
+ /// Returns nullptr for shape functions whose dimensions are less than the
+ /// global dimension.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder(std::false_type* /*unused*/)
+ {
+ return nullptr;
+ }
+};
+
+} // namespace BoundaryConditionOrSourceTerm
+} // namespace ProcessLib
+#undef ENABLED_ELEMENT_TYPE_SIMPLEX
+#undef ENABLED_ELEMENT_TYPE_CUBOID
+#undef ENABLED_ELEMENT_TYPE_PYRAMID
+#undef ENABLED_ELEMENT_TYPE_PRISM
+#undef ENABLED_ELEMENT_TYPE_LINE
+#undef ENABLED_ELEMENT_TYPE_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
index 52a3cf317a93d8763ab9ad224b5d8214cbbd90e1..df1733e6dd2359191a0cc68b79502b2ff15f5d25 100644
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ b/ProcessLib/Utils/LocalDataInitializer.h
@@ -78,7 +78,6 @@ static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0
#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
-#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
@@ -139,11 +138,6 @@ public:
{
// /// Lines and points ///////////////////////////////////
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 0 && OGS_MAX_ELEMENT_ORDER >= 1
- _builder[std::type_index(typeid(MeshLib::Point))] =
- makeLocalAssemblerBuilder<NumLib::ShapePoint1>();
-#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 1
@@ -241,12 +235,6 @@ public:
}
else if (shapefunction_order == 2)
{
-// /// Lines and points ///////////////////////////////////
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 0 && OGS_MAX_ELEMENT_ORDER >= 1
- _builder[std::type_index(typeid(MeshLib::Point))] =
- makeLocalAssemblerBuilder<NumLib::ShapePoint1>();
-#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 && \
OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2