[Process] Separate createLocalAssemblers for IC/BC and the equation, respectively

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(),

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);

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);

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,

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})
 

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(),

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,

ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/CreateLocalAssemblers.h

+/**
+ * \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

ProcessLib/Utils/ForBoundaryConditionOrSourceTerm/LocalDataInitializer.h

+/**
+ * \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

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