[PL] SmallDeformation process/fem implementation.

Documentation/ProjectFile/process/SMALL_DEFORMATION/c_SMALL_DEFORMATION.md

+Deformation process with linear kinematics.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/c_LinearElasticIsotropic.md

+Isotropic linear elastic material model.

Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/t_process_variable.md

+The displacement vector field.

Documentation/ProjectFile/process/SMALL_DEFORMATION/t_dimension.md

+Displacement vector dimension. The displacement dimension must be equal to the mesh dimension.

MaterialLib/SolidModels/MechanicsBase.h

@@ -22,6 +22,11 @@ class Element;
 
 namespace Solids
 {
+/// Interface for mechanical solid material models. Provides updates of the
+/// stress for a given current state and also a tangent at that position. If the
+/// implemented material model stores an internal state, the nested
+/// MaterialStateVariables class should be used; it's only responsibility is to
+/// provide state's push back possibility.
 template <int DisplacementDim>
 struct MechanicsBase
 {
@@ -36,9 +41,16 @@ struct MechanicsBase
         virtual void pushBackState() = 0;
     };
 
+    /// Polymorphic creator for MaterialStateVariables objects specific for a
+    /// material model.
+    virtual std::unique_ptr<MaterialStateVariables>
+    createMaterialStateVariables() = 0;
+
     using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
     using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
 
+    /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
+    /// constitutive relation compute function.
     void computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x,
@@ -76,6 +88,10 @@ struct MechanicsBase
         C = C_;
     }
 
+    /// Computation of the constitutive relation for specific material model.
+    /// This should be implemented in the derived model. Fixed Kelvin vector and
+    /// matrix size version; for dynamic size arguments there is an overloaded
+    /// wrapper function.
     virtual void computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x,
@@ -87,9 +103,6 @@ struct MechanicsBase
         KelvinMatrix& C,
         MaterialStateVariables& material_state_variables) = 0;
 
-    virtual std::unique_ptr<MaterialStateVariables>
-    createMaterialStateVariables() = 0;
-
     virtual ~MechanicsBase() = default;
 };
 }

ProcessLib/CMakeLists.txt

@@ -12,6 +12,8 @@ APPEND_SOURCE_FILES(SOURCES Parameter)
 add_subdirectory(GroundwaterFlow)
 APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
 
+APPEND_SOURCE_FILES(SOURCES SmallDeformation)
+
 add_subdirectory(TES)
 APPEND_SOURCE_FILES(SOURCES TES)
 

ProcessLib/SmallDeformation/CreateLocalAssemblers.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_
+#define PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_
+
+#include <vector>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "LocalDataInitializer.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+namespace detail
+{
+template <unsigned GlobalDim, int DisplacementDim,
+          template <typename, typename, unsigned, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    unsigned const integration_order,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    // Shape matrices initializer
+    using LocalDataInitializer =
+        LocalDataInitializer<LocalAssemblerInterface,
+                             LocalAssemblerImplementation, GlobalDim,
+                             DisplacementDim, ExtraCtorArgs...>;
+
+    DBUG("Create local assemblers.");
+    // Populate the vector of local assemblers.
+    local_assemblers.resize(mesh_elements.size());
+
+    LocalDataInitializer initializer(dof_table);
+
+    DBUG("Calling local assembler builder for all mesh elements.");
+    GlobalExecutor::transformDereferenced(
+        initializer,
+        mesh_elements,
+        local_assemblers,
+        integration_order,
+        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 <int DisplacementDim,
+          template <typename, typename, unsigned, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    const unsigned dimension,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    unsigned const integration_order,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    DBUG("Create local assemblers.");
+
+    switch (dimension)
+    {
+        case 2:
+            detail::createLocalAssemblers<2, DisplacementDim,
+                                          LocalAssemblerImplementation>(
+                dof_table, mesh_elements, integration_order, local_assemblers,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
+        case 3:
+            detail::createLocalAssemblers<3, DisplacementDim,
+                                          LocalAssemblerImplementation>(
+                dof_table, mesh_elements, integration_order, local_assemblers,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
+        default:
+            OGS_FATAL(
+                "Meshes with dimension different than two and three are not "
+                "supported.");
+    }
+}
+}  // SmallDeformation
+
+}  // ProcessLib
+
+#endif  // PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_

ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_
+#define PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+
+#include "SmallDeformationProcess.h"
+#include "SmallDeformationProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+template <int DisplacementDim>
+class SmallDeformationProcess;
+
+extern template class SmallDeformationProcess<2>;
+extern template class SmallDeformationProcess<3>;
+
+template <int DisplacementDim>
+std::unique_ptr<SmallDeformationProcess<DisplacementDim>>
+createSmallDeformationProcess(
+    MeshLib::Mesh& mesh,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{process__type}
+    config.checkConfigParameter("type", "SMALL_DEFORMATION");
+    DBUG("Create SmallDeformationProcess.");
+
+    // Process variable.
+    auto process_variables = findProcessVariables(
+        variables, config,
+        {//! \ogs_file_param_special{process__SMALL_DEFORMATION__process_variables__process_variable}
+         "process_variable"});
+
+    DBUG("Associate displacement with process variable \'%s\'.",
+         process_variables.back().get().getName().c_str());
+
+    if (process_variables.back().get().getNumberOfComponents() !=
+        DisplacementDim)
+    {
+        OGS_FATAL(
+            "Number of components of the process variable '%s' is different "
+            "from the displacement dimension: got %d, expected %d",
+            process_variables.back().get().getName().c_str(),
+            process_variables.back().get().getNumberOfComponents(),
+            DisplacementDim);
+    }
+
+    // Constitutive relation.
+    // read type;
+    auto const constitutive_relation_config =
+        //! \ogs_file_param{process__SMALL_DEFORMATION__constitutive_relation}
+        config.getConfigSubtree("constitutive_relation");
+
+    auto const type =
+        constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+    std::unique_ptr<Solids::MechanicsBase<DisplacementDim>> material = nullptr;
+    if (type == "LinearElasticIsotropic")
+    {
+        material = Solids::createLinearElasticIsotropic<DisplacementDim>(
+            parameters, constitutive_relation_config);
+    }
+    else
+    {
+        OGS_FATAL(
+            "Cannot construct constitutive relation of given type \'%s\'.",
+            type.c_str());
+    }
+
+    SmallDeformationProcessData<DisplacementDim> process_data{
+        std::move(material)};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"SmallDeformation_displacement"});
+
+    ProcessLib::parseSecondaryVariables(config, secondary_variables,
+                                        named_function_caller);
+
+    return std::unique_ptr<SmallDeformationProcess<DisplacementDim>>{
+        new SmallDeformationProcess<DisplacementDim>{
+            mesh, parameters, std::move(process_variables),
+            std::move(process_data), std::move(secondary_variables),
+            std::move(named_function_caller)}};
+}
+}  // namespace SmallDeformation
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_

ProcessLib/SmallDeformation/LocalDataInitializer.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_
+#define PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_
+
+#include <functional>
+#include <memory>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/Integration/GaussIntegrationPolicy.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))
+
+// 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_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
+{
+
+/// 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, unsigned, int> class LocalAssemblerData,
+    unsigned GlobalDim, int DisplacementDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+    explicit LocalDataInitializer(
+        NumLib::LocalToGlobalIndexMap const& dof_table)
+        : _dof_table(dof_table)
+    {
+        // /// 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::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 >= 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::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 >= 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::ShapePrism15>();
+#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::ShapePyra13>();
+#endif
+    }
+
+    /// Sets the provided \c data_ptr 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.
+    void operator()(
+            std::size_t const id,
+            MeshLib::Element const& mesh_item,
+            LADataIntfPtr& data_ptr,
+            unsigned const integration_order,
+            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);
+            data_ptr = it->second(
+                           mesh_item, num_local_dof, integration_order,
+                           std::forward<ConstructorArgs>(args)...);
+        } else {
+            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,
+            unsigned const integration_order,
+            ConstructorArgs&&...
+        )>;
+
+    template <typename ShapeFunction>
+    using IntegrationMethod = typename NumLib::GaussIntegrationPolicy<
+                typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+    template <typename ShapeFunction>
+    using LAData =
+        LocalAssemblerData<ShapeFunction, IntegrationMethod<ShapeFunction>,
+                           GlobalDim, DisplacementDim>;
+
+    /// Generates a function that creates a new LocalAssembler of type LAData<SHAPE_FCT>
+    template<typename ShapeFct>
+    static LADataBuilder makeLocalAssemblerBuilder()
+    {
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  unsigned const integration_order,
+                  ConstructorArgs&&... args)
+        {
+            return LADataIntfPtr{
+                new LAData<ShapeFct>{
+                    e, local_matrix_size, integration_order,
+                    std::forward<ConstructorArgs>(args)...
+                }};
+        };
+    }
+
+    /// Mapping of element types to local assembler constructors.
+    std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+};
+
+}   // 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_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS
+
+#endif  // PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_

ProcessLib/SmallDeformation/SmallDeformationFEM.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_SMALLDEFORMATION_FEM_H_
+#define PROCESS_LIB_SMALLDEFORMATION_FEM_H_
+
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "SmallDeformationProcessData.h"
+
+template <typename BMatricesType, int DisplacementDim>
+struct IntegrationPointData final
+{
+    explicit IntegrationPointData(
+        Solids::MechanicsBase<DisplacementDim>& solid_material)
+        : _solid_material(solid_material),
+          _material_state_variables(
+              _solid_material.createMaterialStateVariables())
+    {
+    }
+
+#if defined(_MSC_VER) && _MSC_VER < 1900
+    // The default generated move-ctor is correctly generated for other
+    // compilers.
+    explicit IntegrationPointData(IntegrationPointData&& other)
+        : _b_matrices(std::move(other._b_matrices)),
+          _sigma(std::move(other._sigma)),
+          _sigma_prev(std::move(other._sigma_prev)),
+          _eps(std::move(other._eps)),
+          _eps_prev(std::move(other._eps_prev)),
+          _solid_material(other._solid_material),
+          _material_state_variables(std::move(other._material_state_variables)),
+          _C(std::move(other._C)),
+          _detJ(std::move(other._detJ))
+    {
+    }
+#endif  // _MSC_VER
+
+    typename BMatricesType::BMatrixType _b_matrices;
+    typename BMatricesType::KelvinVectorType _sigma, _sigma_prev;
+    typename BMatricesType::KelvinVectorType _eps, _eps_prev;
+
+    Solids::MechanicsBase<DisplacementDim>& _solid_material;
+    std::unique_ptr<
+        typename Solids::MechanicsBase<DisplacementDim>::MaterialStateVariables>
+        _material_state_variables;
+
+    typename BMatricesType::KelvinMatrixType _C;
+    double _detJ;
+
+    void pushBackState()
+    {
+        _eps_prev = _eps;
+        _sigma_prev = _sigma;
+        _material_state_variables->pushBackState();
+    }
+};
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+    std::vector<ShapeMatrixType> N;
+    std::vector<double> _sigmaXX;
+    std::vector<double> _sigmaYY;
+    std::vector<double> _sigmaXY;
+};
+
+struct SmallDeformationLocalAssemblerInterface
+    : public ProcessLib::LocalAssemblerInterface,
+      public NumLib::ExtrapolatableElement
+{
+    virtual std::vector<double> const& getIntPtSigmaXX(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaYY(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtSigmaXY(
+        std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          int DisplacementDim>
+class SmallDeformationLocalAssembler
+    : public SmallDeformationLocalAssemblerInterface
+{
+public:
+    using ShapeMatricesType =
+        ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
+    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+    using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+    using BMatrixType = typename BMatricesType::BMatrixType;
+    using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType;
+    using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
+    using NodalDisplacementVectorType =
+        typename BMatricesType::NodalForceVectorType;
+
+    SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) =
+        delete;
+    SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete;
+
+    SmallDeformationLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const local_matrix_size,
+        unsigned const integration_order,
+        SmallDeformationProcessData<DisplacementDim>& process_data)
+        : _process_data(process_data),
+          _localA(local_matrix_size, local_matrix_size),
+          _localRhs(local_matrix_size),
+          _integration_method(integration_order),
+          _element(e)
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        _ip_data.reserve(n_integration_points);
+
+        _secondary_data.N.resize(n_integration_points);
+        _secondary_data._sigmaXX.resize(n_integration_points);
+        _secondary_data._sigmaYY.resize(n_integration_points);
+        _secondary_data._sigmaXY.resize(n_integration_points);
+
+        auto const shape_matrices =
+            initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                              IntegrationMethod, DisplacementDim>(
+                e, _integration_method);
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data.emplace_back(*_process_data._material);
+            auto& ip_data = _ip_data[ip];
+            ip_data._detJ = shape_matrices[ip].detJ;
+            ip_data._b_matrices.resize(
+                KelvinVectorDimensions<DisplacementDim>::value,
+                ShapeFunction::NPOINTS * DisplacementDim);
+            LinearBMatrix::computeBMatrix<
+                DisplacementDim, ShapeFunction::NPOINTS,
+                typename ShapeMatricesType::GlobalDimNodalMatrixType,
+                BMatrixType>(shape_matrices[ip].dNdx, ip_data._b_matrices);
+
+            ip_data._sigma.resize(
+                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data._sigma_prev.resize(
+                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data._eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data._eps_prev.resize(
+                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data._C.resize(KelvinVectorDimensions<DisplacementDim>::value,
+                              KelvinVectorDimensions<DisplacementDim>::value);
+
+            _secondary_data.N[ip] = shape_matrices[ip].N;
+        }
+    }
+
+    void assembleConcrete(
+        double const t, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& /*M*/, GlobalMatrix& /*K*/, GlobalVector& b) override
+    {
+        _localRhs.setZero();
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        NodalDisplacementVectorType local_displacement(ShapeFunction::NPOINTS *
+                                                       DisplacementDim);
+
+        SpatialPosition x_position;
+        x_position.setElementID(_element.getID());
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            x_position.setIntegrationPoint(ip);
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            auto const& detJ = _ip_data[ip]._detJ;
+
+            auto const& B = _ip_data[ip]._b_matrices;
+            auto const& eps_prev = _ip_data[ip]._eps_prev;
+            auto const& sigma_prev = _ip_data[ip]._sigma_prev;
+
+            auto& eps = _ip_data[ip]._eps;
+            auto& sigma = _ip_data[ip]._sigma;
+            auto& C = _ip_data[ip]._C;
+            auto& material_state_variables =
+                *_ip_data[ip]._material_state_variables;
+
+            eps.noalias() =
+                B *
+                Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
+                    local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
+
+            _ip_data[ip]._solid_material.computeConstitutiveRelation(
+                t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
+                sigma, C, material_state_variables);
+
+            _localRhs.noalias() -=
+                B.transpose() * sigma * detJ * wp.getWeight();
+
+            // TODO: Reuse _ip_data[ip]._sigma
+            _secondary_data._sigmaXX[ip] = sigma[0];
+            _secondary_data._sigmaYY[ip] = sigma[1];
+            _secondary_data._sigmaXY[ip] = sigma[3];
+        }
+
+        b.add(indices.rows, _localRhs);
+    }
+
+    void assembleJacobianConcrete(
+        double const /*t*/,
+        std::vector<double> const& /*local_x*/,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& Jac) override
+    {
+        _localA.setZero();
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            auto const& B = _ip_data[ip]._b_matrices;
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            auto& C = _ip_data[ip]._C;
+            auto const& detJ = _ip_data[ip]._detJ;
+
+            _localA.noalias() += B.transpose() * C * B * detJ * wp.getWeight();
+        }
+
+        Jac.add(indices, _localA);
+    }
+
+    void preTimestepConcrete(std::vector<double> const& /*local_x*/,
+                             double const /*t*/,
+                             double const /*delta_t*/) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data[ip].pushBackState();
+        }
+    }
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N = _secondary_data.N[integration_point];
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+    }
+
+    std::vector<double> const& getIntPtSigmaXX(
+        std::vector<double>& /*cache*/) const override
+    {
+        return _secondary_data._sigmaXX;
+    }
+
+    std::vector<double> const& getIntPtSigmaYY(
+        std::vector<double>& /*cache*/) const override
+    {
+        return _secondary_data._sigmaYY;
+    }
+
+    std::vector<double> const& getIntPtSigmaXY(
+        std::vector<double>& /*cache*/) const override
+    {
+        return _secondary_data._sigmaXY;
+    }
+
+private:
+    SmallDeformationProcessData<DisplacementDim>& _process_data;
+
+    std::vector<IntegrationPointData<BMatricesType, DisplacementDim>> _ip_data;
+
+    StiffnessMatrixType _localA;
+    NodalForceVectorType _localRhs;
+
+    IntegrationMethod _integration_method;
+    MeshLib::Element const& _element;
+    SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim, int DisplacementDim>
+class LocalAssemblerData final
+    : public SmallDeformationLocalAssembler<ShapeFunction, IntegrationMethod,
+                                            DisplacementDim>
+{
+public:
+    LocalAssemblerData(LocalAssemblerData const&) = delete;
+    LocalAssemblerData(LocalAssemblerData&&) = delete;
+
+    LocalAssemblerData(
+        MeshLib::Element const& e,
+        std::size_t const local_matrix_size,
+        unsigned const integration_order,
+        SmallDeformationProcessData<DisplacementDim>& process_data)
+        : SmallDeformationLocalAssembler<ShapeFunction, IntegrationMethod,
+                                         DisplacementDim>(
+              e, local_matrix_size, integration_order, process_data)
+    {
+    }
+};
+
+}  // namespace SmallDeformation
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_SMALLDEFORMATION_FEM_H_

ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_
+#define PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_
+
+#include "SmallDeformationProcess.h"
+
+extern template class ProcessLib::SmallDeformation::SmallDeformationProcess<2>;
+extern template class ProcessLib::SmallDeformation::SmallDeformationProcess<3>;
+
+#endif  // PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_

ProcessLib/SmallDeformation/SmallDeformationProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "SmallDeformationProcess-fwd.h"
+#include "SmallDeformationProcess.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+
+template class SmallDeformationProcess<2>;
+template class SmallDeformationProcess<3>;
+
+}   // namespace SmallDeformation
+}   // namespace ProcessLib

ProcessLib/SmallDeformation/SmallDeformationProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_
+#define PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_
+
+#include <cassert>
+
+#include "ProcessLib/Process.h"
+#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+
+#include "SmallDeformationFEM.h"
+#include "SmallDeformationProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+template <int DisplacementDim>
+class SmallDeformationProcess final : public Process
+{
+    using Base = Process;
+
+public:
+    SmallDeformationProcess(
+        MeshLib::Mesh& mesh,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        std::vector<std::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        SmallDeformationProcessData<DisplacementDim>&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller)
+        : Process(mesh, parameters, std::move(process_variables),
+                  std::move(secondary_variables),
+                  std::move(named_function_caller)),
+          _process_data(std::move(process_data))
+    {
+    }
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override { return false; }
+    //! @}
+
+private:
+    using LocalAssemblerInterface = SmallDeformationLocalAssemblerInterface;
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override
+    {
+        ProcessLib::SmallDeformation::createLocalAssemblers<DisplacementDim,
+                                                            LocalAssemblerData>(
+            mesh.getDimension(), mesh.getElements(), dof_table,
+            integration_order, _local_assemblers, _process_data);
+
+        // TODO move the two data members somewhere else.
+        // for extrapolation of secondary variables
+        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
+            all_mesh_subsets_single_component;
+        all_mesh_subsets_single_component.emplace_back(
+            new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+        _local_to_global_index_map_single_component.reset(
+            new NumLib::LocalToGlobalIndexMap(
+                std::move(all_mesh_subsets_single_component),
+                // by location order is needed for output
+                NumLib::ComponentOrder::BY_LOCATION));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "sigma_xx", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXX));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "sigma_yy", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtSigmaYY));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "sigma_xy", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY));
+    }
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override
+    {
+        DBUG("Assemble SmallDeformationProcess.");
+
+        // Call global assembler for each local assembly item.
+        GlobalExecutor::executeMemberOnDereferenced(
+            &SmallDeformationLocalAssemblerInterface::assemble,
+            _local_assemblers, *_local_to_global_index_map, t, x, M, K, b);
+    }
+
+    void assembleJacobianConcreteProcess(const double t, GlobalVector const& x,
+                                         GlobalVector const& /*xdot*/,
+                                         const double /*dxdot_dx*/,
+                                         GlobalMatrix const& /*M*/,
+                                         const double /*dx_dx*/,
+                                         GlobalMatrix const& /*K*/,
+                                         GlobalMatrix& Jac) override
+    {
+        DBUG("AssembleJacobian SmallDeformationProcess.");
+
+        // Call global assembler for each local assembly item.
+        GlobalExecutor::executeMemberOnDereferenced(
+            &SmallDeformationLocalAssemblerInterface::assembleJacobian,
+            _local_assemblers, *_local_to_global_index_map, t, x, Jac);
+    }
+
+    void preTimestep(GlobalVector const& x, double const t,
+                     double const dt) override
+    {
+        DBUG("PreTimestep SmallDeformationProcess.");
+
+        _process_data.dt = dt;
+        _process_data.t = t;
+
+        GlobalExecutor::executeMemberOnDereferenced(
+            &SmallDeformationLocalAssemblerInterface::preTimestep,
+            _local_assemblers, *_local_to_global_index_map, x, t, dt);
+    }
+
+    void postTimestep(GlobalVector const& x) override
+    {
+        DBUG("PostTimestep SmallDeformationProcess.");
+
+        GlobalExecutor::executeMemberOnDereferenced(
+            &SmallDeformationLocalAssemblerInterface::postTimestep,
+            _local_assemblers, *_local_to_global_index_map, x);
+    }
+
+private:
+    SmallDeformationProcessData<DisplacementDim> _process_data;
+
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_single_component;
+};
+
+}  // namespace SmallDeformation
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_

ProcessLib/SmallDeformation/SmallDeformationProcessData.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_
+#define PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+template <int DisplacementDim>
+struct SmallDeformationProcessData
+{
+    SmallDeformationProcessData(
+        std::unique_ptr<Solids::MechanicsBase<DisplacementDim>>&& material)
+        : _material{std::move(material)}
+    {
+    }
+
+    SmallDeformationProcessData(SmallDeformationProcessData&& other)
+        : _material{std::move(other._material)}
+    {
+    }
+
+    //! Copies are forbidden.
+    SmallDeformationProcessData(SmallDeformationProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(SmallDeformationProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(SmallDeformationProcessData&&) = delete;
+
+    std::unique_ptr<Solids::MechanicsBase<DisplacementDim>> _material;
+    double dt;
+    double t;
+};
+
+}  // namespace SmallDeformation
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_