Merge pull request #1452 from norihiro-w/sd-lie

[SD-LIE] Small deformation process with lower-dimensional interface elements

Applications/ApplicationsLib/ProjectData.cpp

@@ -36,6 +36,7 @@
 #include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
 #include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
+#include "ProcessLib/SmallDeformationWithLIE/CreateSmallDeformationProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
 
 namespace detail
@@ -330,6 +331,23 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                         "given dimension");
             }
         }
+        else if (type == "SMALL_DEFORMATION_WITH_LIE")
+        {
+            switch (process_config.getConfigParameter<int>("dimension"))
+            {
+                case 2:
+                    process = ProcessLib::SmallDeformationWithLIE::
+                        createSmallDeformationProcess<2>(
+                            *_mesh_vec[0], std::move(jacobian_assembler),
+                            _process_variables, _parameters, integration_order,
+                            process_config);
+                    break;
+                default:
+                    OGS_FATAL(
+                        "SMALL_DEFORMATION_WITH_LIE process does not support "
+                        "given dimension");
+            }
+        }
         else if (type == "RICHARDS_FLOW")
         {
             process = ProcessLib::RichardsFlow::createRichardsFlowProcess(

Applications/CLI/Tests.cmake

@@ -503,6 +503,20 @@ if(NOT OGS_USE_MPI)
     #     EXECUTABLE_ARGS tes-inert-wedge.prj
     # )
 
+    # LIE; Small deformation
+    AddTest(
+        NAME LIE_M_single_joint
+        PATH LIE/Mechanics
+        EXECUTABLE ogs
+        EXECUTABLE_ARGS single_joint.prj
+        WRAPPER time
+        TESTER vtkdiff
+        ABSTOL 1e-16 RELTOL 1e-16
+        DIFF_DATA
+        single_joint_expected_pcs_0_ts_1_t_1.000000.vtu single_joint_pcs_0_ts_1_t_1.000000.vtu displacement displacement
+        single_joint_expected_pcs_0_ts_1_t_1.000000.vtu single_joint_pcs_0_ts_1_t_1.000000.vtu displacement_jump1 displacement_jump1
+    )
+
 else()
     # MPI groundwater flow tests
     AddTest(

NumLib/DOF/LocalToGlobalIndexMap.cpp

@@ -313,6 +313,7 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
         << " rows\n";
     for (std::size_t e=0; e<map.size(); ++e)
     {
+        os << "== e " << e << " ==\n";
         for (std::size_t c=0; c<map.getNumberOfComponents(); ++c)
         {
             auto const& line = map._rows(e, c);

ProcessLib/Process.h

@@ -154,8 +154,10 @@ private:
         return NumLib::IterationResult::SUCCESS;
     }
 
-    void constructDofTable();
+protected:
+    virtual void constructDofTable();
 
+private:
     void initializeExtrapolator();
 
     /// Finishes the \c _named_function_caller and adds a secondary variable for

ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h

 /**
  * \copyright
- * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * 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

ProcessLib/SmallDeformationWithLIE/Common/FractureProperty.h

@@ -17,7 +17,6 @@
 
 #include "ProcessLib/Parameter/Parameter.h"
 
-#include "FractureProperty.h"
 #include "Utils.h"
 
 

ProcessLib/SmallDeformationWithLIE/CreateSmallDeformationProcess.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 "CreateSmallDeformationProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/FractureModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/FractureModels/CreateMohrCoulomb.h"
+
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"  // required for findParameter
+
+#include "SmallDeformationProcess.h"
+#include "SmallDeformationProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+template <int DisplacementDim>
+class SmallDeformationProcess;
+
+extern template class SmallDeformationProcess<2>;
+
+template <int DisplacementDim>
+std::unique_ptr<Process>
+createSmallDeformationProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{process__type}
+    config.checkConfigParameter("type", "SMALL_DEFORMATION_WITH_LIE");
+    DBUG("Create SmallDeformationProcess with LIE.");
+
+    // Process variables
+    auto const pv_conf = config.getConfigSubtree("process_variables");
+    auto range = pv_conf.getConfigParameterList<std::string>("process_variable");
+    std::vector<std::reference_wrapper<ProcessVariable>> process_variables;
+    for (std::string const& pv_name : range)
+    {
+        if (pv_name != "displacement"
+            && pv_name.find("displacement_jump")==std::string::npos)
+            OGS_FATAL("Found a process variable name '%s'. It should be 'displacement' or 'displacement_jumpN'");
+        auto variable = std::find_if(
+            variables.cbegin(), variables.cend(),
+            [&pv_name](ProcessVariable const& v) { return v.getName() == pv_name; });
+
+        if (variable == variables.end())
+        {
+            OGS_FATAL(
+                "Could not find process variable '%s' in the provided variables "
+                "list for config tag <%s>.",
+                pv_name.c_str(), "process_variable");
+        }
+        DBUG("Found process variable \'%s\' for config tag <%s>.",
+             variable->getName().c_str(), "process_variable");
+
+        process_variables.emplace_back(const_cast<ProcessVariable&>(*variable));
+    }
+    if (process_variables.size() > 2)
+        OGS_FATAL("Currently only one displacement jump is supported");
+
+    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_WITH_LIE__constitutive_relation}
+        config.getConfigSubtree("constitutive_relation");
+
+    auto const type =
+        constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>> material = nullptr;
+    if (type == "LinearElasticIsotropic")
+    {
+        material = MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
+            parameters, constitutive_relation_config);
+    }
+    else
+    {
+        OGS_FATAL(
+            "Cannot construct constitutive relation of given type \'%s\'.",
+            type.c_str());
+    }
+
+    // Fracture constitutive relation.
+    // read type;
+    auto const fracture_constitutive_relation_config =
+        //! \ogs_file_param{process__SMALL_DEFORMATION_WITH_LIE__constitutive_relation}
+        config.getConfigSubtree("fracture_constitutive_relation");
+
+    auto const frac_type =
+        fracture_constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+    std::unique_ptr<MaterialLib::Fracture::FractureModelBase<DisplacementDim>> fracture_model = nullptr;
+    if (frac_type == "LinearElasticIsotropic")
+    {
+        fracture_model = MaterialLib::Fracture::createLinearElasticIsotropic<DisplacementDim>(
+            parameters, fracture_constitutive_relation_config);
+    }
+    else if (frac_type == "MohrCoulomb")
+    {
+        fracture_model = MaterialLib::Fracture::createMohrCoulomb<DisplacementDim>(
+            parameters, fracture_constitutive_relation_config);
+    }
+    else
+    {
+        OGS_FATAL(
+            "Cannot construct fracture constitutive relation of given type \'%s\'.",
+            frac_type.c_str());
+    }
+
+    // Fracture properties
+    //! \ogs_file_param{process__SMALL_DEFORMATION_WITH_LIE__fracture_properties}
+    auto fracture_properties_config = config.getConfigSubtree("fracture_properties");
+    auto &para_b0 = ProcessLib::findParameter<double>(fracture_properties_config, "initial_aperture", parameters, 1);
+    std::unique_ptr<FractureProperty> frac_prop(new FractureProperty());
+    frac_prop->mat_id = fracture_properties_config.getConfigParameter<int>("material_id");
+    frac_prop->aperture0 = &para_b0;
+
+
+    SmallDeformationProcessData<DisplacementDim> process_data(
+        std::move(material), std::move(fracture_model), std::move(frac_prop));
+
+    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, std::move(jacobian_assembler), parameters, integration_order,
+            std::move(process_variables), std::move(process_data),
+            std::move(secondary_variables), std::move(named_function_caller)}};
+}
+
+template
+std::unique_ptr<Process>
+createSmallDeformationProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+

ProcessLib/SmallDeformationWithLIE/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_SMALLDEFORMATION_WITH_LIE_CREATESMALLDEFORMATIONPROCESS_H_
+#define PROCESS_LIB_SMALLDEFORMATION_WITH_LIE_CREATESMALLDEFORMATIONPROCESS_H_
+
+#include <memory>
+#include <vector>
+
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createSmallDeformationProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_SMALLDEFORMATION_WITH_LIE_CREATESMALLDEFORMATIONPROCESS_H_

ProcessLib/SmallDeformationWithLIE/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h

@@ -147,7 +147,10 @@ private:
         cache.reserve(_ip_data.size());
 
         for (auto const& ip_data : _ip_data) {
-            cache.push_back(ip_data._sigma[component]);
+            if (component < 3)  // xx, yy, zz components
+                cache.push_back(ip_data._sigma[component]);
+            else    // mixed xy, yz, xz components
+                cache.push_back(ip_data._sigma[component] / std::sqrt(2));
         }
 
         return cache;

ProcessLib/SmallDeformationWithLIE/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h

@@ -148,7 +148,10 @@ private:
         cache.reserve(_ip_data.size());
 
         for (auto const& ip_data : _ip_data) {
-            cache.push_back(ip_data._sigma[component]);
+            if (component < 3)  // xx, yy, zz components
+                cache.push_back(ip_data._sigma[component]);
+            else    // mixed xy, yz, xz components
+                cache.push_back(ip_data._sigma[component] / std::sqrt(2));
         }
 
         return cache;

ProcessLib/SmallDeformationWithLIE/SmallDeformationProcess-fwd.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_WITH_LIE_PROCESS_FWD_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_PROCESS_FWD_H_
+
+#include "SmallDeformationProcess.h"
+
+extern template class ProcessLib::SmallDeformationWithLIE::SmallDeformationProcess<2>;
+
+#endif  // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_PROCESS_FWD_H_

ProcessLib/SmallDeformationWithLIE/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"
+
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <vector>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/ElementCoordinatesMappingLocal.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+#include "MeshLib/Properties.h"
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h"
+#include "ProcessLib/SmallDeformationWithLIE/Common/LevelSetFunction.h"
+#include "ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h"
+#include "ProcessLib/SmallDeformationWithLIE/Common/Utils.h"
+#include "ProcessLib/SmallDeformationWithLIE/LocalAssembler/LocalAssemblerDataMatrix.h"
+#include "ProcessLib/SmallDeformationWithLIE/LocalAssembler/LocalAssemblerDataMatrixNearFracture.h"
+#include "ProcessLib/SmallDeformationWithLIE/LocalAssembler/LocalAssemblerDataFracture.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+template <int DisplacementDim>
+SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+        jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::reference_wrapper<ProcessVariable>>&&
+        process_variables,
+    SmallDeformationProcessData<DisplacementDim>&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller)
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              integration_order,
+              std::move(process_variables),
+              std::move(secondary_variables),
+              std::move(named_function_caller)),
+      _process_data(std::move(process_data))
+{
+    getFractureMatrixDataInMesh(mesh,
+                                _vec_matrix_elements,
+                                _vec_fracture_elements,
+                                _vec_fracture_matrix_elements,
+                                _vec_fracture_nodes);
+
+    // set fracture property assuming a fracture forms a straight line
+    setFractureProperty(DisplacementDim,
+                        *_vec_fracture_elements[0],
+                        *_process_data._fracture_property.get());
+
+    // need to use a custom Neumann BC assembler for displacement jumps
+    for (ProcessVariable& pv : getProcessVariables())
+    {
+        if (pv.getName().find("displacement_jump") == std::string::npos)
+            continue;
+        pv.setBoundaryConditionBuilder(
+                    std::unique_ptr<ProcessLib::BoundaryConditionBuilder>(
+                        new BoundaryConditionBuilder(*_process_data._fracture_property.get())));
+    }
+
+}
+
+
+template <int DisplacementDim>
+void SmallDeformationProcess<DisplacementDim>::constructDofTable()
+{
+    //------------------------------------------------------------
+    // prepare mesh subsets to define DoFs
+    //------------------------------------------------------------
+    // for extrapolation
+    _mesh_subset_all_nodes.reset(new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
+    // regular u
+    _mesh_subset_matrix_nodes.reset(new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
+    // u jump
+    _mesh_subset_fracture_nodes.reset(new MeshLib::MeshSubset(_mesh, &_vec_fracture_nodes));
+
+    // Collect the mesh subsets in a vector.
+    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::generate_n(
+        std::back_inserter(all_mesh_subsets),
+        DisplacementDim,
+        [&]() {
+            return std::unique_ptr<MeshLib::MeshSubsets>{
+                new MeshLib::MeshSubsets{_mesh_subset_matrix_nodes.get()}};
+        });
+    std::generate_n(
+        std::back_inserter(all_mesh_subsets),
+        DisplacementDim,
+        [&]() {
+            return std::unique_ptr<MeshLib::MeshSubsets>{
+                new MeshLib::MeshSubsets{_mesh_subset_fracture_nodes.get()}};
+        });
+
+    std::vector<unsigned> const vec_n_components(2, DisplacementDim);
+
+    std::vector<std::vector<MeshLib::Element*>const*> vec_var_elements;
+    vec_var_elements.push_back(&_vec_matrix_elements);
+    vec_var_elements.push_back(&_vec_fracture_matrix_elements);
+
+    _local_to_global_index_map.reset(
+        new NumLib::LocalToGlobalIndexMap(
+            std::move(all_mesh_subsets),
+            vec_n_components,
+            vec_var_elements,
+            NumLib::ComponentOrder::BY_COMPONENT));
+}
+
+
+template <int DisplacementDim>
+void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    ProcessLib::SmallDeformationWithLIE::createLocalAssemblers
+            <DisplacementDim,
+             LocalAssemblerDataMatrix,
+             LocalAssemblerDataMatrixNearFracture,
+             LocalAssemblerDataFracture>(
+        mesh.getDimension(), mesh.getElements(), dof_table,
+        _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+        _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_zz", 1,
+        makeExtrapolator(
+            getExtrapolator(), _local_assemblers,
+            &SmallDeformationLocalAssemblerInterface::getIntPtSigmaZZ));
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma_xy", 1,
+        makeExtrapolator(
+            getExtrapolator(), _local_assemblers,
+            &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY));
+
+    auto mesh_prop_sigma_xx = const_cast<MeshLib::Mesh&>(mesh)
+                                  .getProperties()
+                                  .template createNewPropertyVector<double>(
+                                      "stress_xx", MeshLib::MeshItemType::Cell);
+    mesh_prop_sigma_xx->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_stress_xx = mesh_prop_sigma_xx;
+
+    auto mesh_prop_sigma_yy = const_cast<MeshLib::Mesh&>(mesh)
+                                  .getProperties()
+                                  .template createNewPropertyVector<double>(
+                                      "stress_yy", MeshLib::MeshItemType::Cell);
+    mesh_prop_sigma_yy->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_stress_yy = mesh_prop_sigma_yy;
+
+    auto mesh_prop_sigma_xy = const_cast<MeshLib::Mesh&>(mesh)
+                                  .getProperties()
+                                  .template createNewPropertyVector<double>(
+                                      "stress_xy", MeshLib::MeshItemType::Cell);
+    mesh_prop_sigma_xy->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_stress_xy = mesh_prop_sigma_xy;
+
+    auto mesh_prop_epsilon_xx =
+        const_cast<MeshLib::Mesh&>(mesh)
+            .getProperties()
+            .template createNewPropertyVector<double>(
+                "strain_xx", MeshLib::MeshItemType::Cell);
+    mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_strain_xx = mesh_prop_epsilon_xx;
+
+    auto mesh_prop_epsilon_yy =
+        const_cast<MeshLib::Mesh&>(mesh)
+            .getProperties()
+            .template createNewPropertyVector<double>(
+                "strain_yy", MeshLib::MeshItemType::Cell);
+    mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_strain_yy = mesh_prop_epsilon_yy;
+
+    auto mesh_prop_epsilon_xy =
+        const_cast<MeshLib::Mesh&>(mesh)
+            .getProperties()
+            .template createNewPropertyVector<double>(
+                "strain_xy", MeshLib::MeshItemType::Cell);
+    mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements());
+    _process_data._mesh_prop_strain_xy = mesh_prop_epsilon_xy;
+
+    auto mesh_prop_levelset = const_cast<MeshLib::Mesh&>(mesh)
+                                  .getProperties()
+                                  .template createNewPropertyVector<double>(
+                                      "levelset1", MeshLib::MeshItemType::Cell);
+    mesh_prop_levelset->resize(mesh.getNumberOfElements());
+    for (MeshLib::Element const* e : _mesh.getElements())
+    {
+        if (e->getDimension() < DisplacementDim)
+            continue;
+
+        double const levelsets =
+            calculateLevelSetFunction(*_process_data._fracture_property,
+                                      e->getCenterOfGravity().getCoords());
+        (*mesh_prop_levelset)[e->getID()] = levelsets;
+    }
+
+    auto mesh_prop_b = const_cast<MeshLib::Mesh&>(mesh)
+                           .getProperties()
+                           .template createNewPropertyVector<double>(
+                               "aperture", MeshLib::MeshItemType::Cell);
+    mesh_prop_b->resize(mesh.getNumberOfElements());
+    auto mesh_prop_matid = mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+    auto frac = _process_data._fracture_property.get();
+    for (MeshLib::Element const* e : _mesh.getElements())
+    {
+        if (e->getDimension() == DisplacementDim)
+            continue;
+        if ((*mesh_prop_matid)[e->getID()] != frac->mat_id)
+            continue;
+        ProcessLib::SpatialPosition x;
+        x.setElementID(e->getID());
+        (*mesh_prop_b)[e->getID()] = (*frac->aperture0)(0, x)[0];
+    }
+    _process_data._mesh_prop_b = mesh_prop_b;
+}
+
+
+template <int DisplacementDim>
+void SmallDeformationProcess<DisplacementDim>::postTimestepConcreteProcess(GlobalVector const& x)
+{
+    DBUG("PostTimestep SmallDeformationProcess.");
+
+    GlobalExecutor::executeMemberOnDereferenced(
+        &SmallDeformationLocalAssemblerInterface::postTimestep,
+        _local_assemblers, *_local_to_global_index_map, x);
+
+}
+
+
+// ------------------------------------------------------------------------------------
+// template instantiation
+// ------------------------------------------------------------------------------------
+template class SmallDeformationProcess<2>;
+
+}   // namespace SmallDeformationWithLIE
+}   // namespace ProcessLib

ProcessLib/SmallDeformationWithLIE/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 PROCESSLIB_SMALLDEFORMATION_WITH_LIE_PROCESS_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_PROCESS_H_
+
+#include <cassert>
+
+#include "ProcessLib/Process.h"
+
+#include "ProcessLib/SmallDeformationWithLIE/LocalAssembler/CreateLocalAssemblers.h"
+#include "ProcessLib/SmallDeformationWithLIE/LocalAssembler/SmallDeformationLocalAssemblerInterface.h"
+
+#include "SmallDeformationProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+template <int DisplacementDim>
+class SmallDeformationProcess final : public Process
+{
+    using Base = Process;
+
+    static_assert(DisplacementDim==2,
+                  "Currently SmallDeformationWithLIE::SmallDeformationProcess "
+                  "supports only 2D.");
+
+public:
+    SmallDeformationProcess(
+        MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        unsigned const integration_order,
+        std::vector<std::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        SmallDeformationProcessData<DisplacementDim>&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller);
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override { return false; }
+    //! @}
+
+private:
+    using LocalAssemblerInterface = SmallDeformationLocalAssemblerInterface;
+
+    void constructDofTable() override;
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    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::executeMemberDereferenced(
+            _global_assembler, &VectorMatrixAssembler::assemble,
+            _local_assemblers, *_local_to_global_index_map, t, x, M, K, b);
+    }
+
+    void assembleWithJacobianConcreteProcess(
+        const double t, GlobalVector const& x, GlobalVector const& xdot,
+        const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override
+    {
+        DBUG("AssembleWithJacobian SmallDeformationProcess.");
+
+        // Call global assembler for each local assembly item.
+        GlobalExecutor::executeMemberDereferenced(
+            _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+            _local_assemblers, *_local_to_global_index_map, t, x, xdot,
+            dxdot_dx, dx_dx, M, K, b, Jac);
+    }
+
+    void preTimestepConcreteProcess(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 postTimestepConcreteProcess(GlobalVector const& x) override;
+
+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;
+
+    std::vector<MeshLib::Element*> _vec_matrix_elements;
+    std::vector<MeshLib::Element*> _vec_fracture_elements;
+    std::vector<MeshLib::Element*> _vec_fracture_matrix_elements;
+    std::vector<MeshLib::Node*> _vec_fracture_nodes;
+
+    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_fracture_nodes;
+    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_matrix_nodes;
+};
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_PROCESS_H_