diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index faa16291ef1d36d57afa619b3ee2ae668a389fe3..e4ff5913c6a6eee165483282a8887265761645d3 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -81,6 +81,9 @@
#ifdef OGS_BUILD_PROCESS_SMALLDEFORMATION
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#endif
+#ifdef OGS_BUILD_PROCESS_SMALLDEFORMATIONNONLOCAL
+#include "ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.h"
+#endif
#ifdef OGS_BUILD_PROCESS_TES
#include "ProcessLib/TES/CreateTESProcess.h"
#endif
@@ -530,6 +533,34 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
}
else
#endif
+#ifdef OGS_BUILD_PROCESS_SMALLDEFORMATIONNONLOCAL
+ if (type == "SMALL_DEFORMATION_NONLOCAL")
+ {
+ switch (_mesh_vec[0]->getDimension())
+ {
+ case 2:
+ process = ProcessLib::SmallDeformationNonlocal::
+ createSmallDeformationNonlocalProcess<2>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ break;
+ case 3:
+ process = ProcessLib::SmallDeformationNonlocal::
+ createSmallDeformationNonlocalProcess<3>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ break;
+ default:
+ OGS_FATAL(
+ "SMALL_DEFORMATION_NONLOCAL process does not support "
+ "given dimension %d",
+ _mesh_vec[0]->getDimension());
+ }
+ }
+ else
+#endif
#ifdef OGS_BUILD_PROCESS_LIE
if (type == "SMALL_DEFORMATION_WITH_LIE")
{
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a72111375c3854b6b44a890dbaa341152b597c74..6cee71c71463810957f7ffcdee16b727c1f3e357 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -98,6 +98,7 @@ set(ProcessesList
RichardsFlow
RichardsMechanics
SmallDeformation
+ SmallDeformationNonlocal
TES
ThermalTwoPhaseFlowWithPP
ThermoMechanicalPhaseField
diff --git a/ProcessLib/Deformation/Divergence.h b/ProcessLib/Deformation/Divergence.h
new file mode 100644
index 0000000000000000000000000000000000000000..adbc0ca42db01bbbfa8a25d587aa96400326ad44
--- /dev/null
+++ b/ProcessLib/Deformation/Divergence.h
@@ -0,0 +1,34 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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
+
+namespace ProcessLib
+{
+namespace Deformation
+{
+/// Divergence of displacement, the volumetric strain.
+template <int DisplacementDim, int NPOINTS, typename DNDX_Type>
+double divergence(
+ const Eigen::Ref<Eigen::Matrix<double, NPOINTS * DisplacementDim, 1> const>&
+ u,
+ DNDX_Type const& dNdx)
+{
+ double divergence = 0;
+ for (int i = 0; i < DisplacementDim; ++i)
+ {
+ divergence += dNdx.template block<1, NPOINTS>(i, 0) *
+ u.template segment<NPOINTS>(i * NPOINTS);
+ }
+ return divergence;
+}
+} // namespace Deformation
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/CMakeLists.txt b/ProcessLib/SmallDeformationNonlocal/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f99ef53c547de078cdd9a0d49275a1a88629a72c
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/CMakeLists.txt
@@ -0,0 +1,6 @@
+APPEND_SOURCE_FILES(SOURCES)
+
+add_library(SmallDeformationNonlocal ${SOURCES})
+target_link_libraries(SmallDeformationNonlocal PUBLIC ProcessLib)
+
+include(Tests.cmake)
diff --git a/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.cpp b/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..775ed844a65d9bbb2b14b89df73f35a99f669114
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.cpp
@@ -0,0 +1,141 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "CreateSmallDeformationNonlocalProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "SmallDeformationNonlocalProcess.h"
+#include "SmallDeformationNonlocalProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createSmallDeformationNonlocalProcess(
+ 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{prj__processes__process__type}
+ config.checkConfigParameter("type", "SMALL_DEFORMATION_NONLOCAL");
+ DBUG("Create SmallDeformationNonlocalProcess.");
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ auto per_process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__process_variables__process_variable}
+ "process_variable"});
+
+ DBUG("Associate displacement with process variable '%s'.",
+ per_process_variables.back().get().getName().c_str());
+
+ if (per_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",
+ per_process_variables.back().get().getName().c_str(),
+ per_process_variables.back().get().getNumberOfComponents(),
+ DisplacementDim);
+ }
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+ process_variables;
+ process_variables.push_back(std::move(per_process_variables));
+
+ auto solid_constitutive_relations =
+ MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
+ parameters, config);
+
+ // Solid density
+ auto& solid_density = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__solid_density}
+ "solid_density", parameters, 1);
+ DBUG("Use '%s' as solid density parameter.", solid_density.name.c_str());
+
+ // Specific body force
+ Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+ {
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__specific_body_force}
+ config.getConfigParameter<std::vector<double>>(
+ "specific_body_force");
+ if (b.size() != DisplacementDim)
+ OGS_FATAL(
+ "The size of the specific body force vector does not match the "
+ "displacement dimension. Vector size is %d, displacement "
+ "dimension is %d",
+ b.size(), DisplacementDim);
+
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+ }
+
+ // Reference temperature
+ const auto& reference_temperature =
+ //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__reference_temperature}
+ config.getConfigParameter<double>(
+ "reference_temperature", std::numeric_limits<double>::quiet_NaN());
+
+ auto const internal_length =
+ //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_NONLOCAL__internal_length}
+ config.getConfigParameter<double>("internal_length");
+
+ SmallDeformationNonlocalProcessData<DisplacementDim> process_data{
+ materialIDs(mesh), std::move(solid_constitutive_relations),
+ solid_density, specific_body_force,
+ reference_temperature, internal_length};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"SmallDeformationNonlocal_displacement"});
+
+ ProcessLib::createSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::make_unique<SmallDeformationNonlocalProcess<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> createSmallDeformationNonlocalProcess<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);
+
+template std::unique_ptr<Process> createSmallDeformationNonlocalProcess<3>(
+ 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 SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.h b/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..06e5723be1249911950b28f0ce1f82d8c0454ecc
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/CreateSmallDeformationNonlocalProcess.h
@@ -0,0 +1,65 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+struct ParameterBase;
+class Process;
+class ProcessVariable;
+} // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createSmallDeformationNonlocalProcess(
+ 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);
+
+extern template std::unique_ptr<Process>
+createSmallDeformationNonlocalProcess<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);
+
+extern template std::unique_ptr<Process>
+createSmallDeformationNonlocalProcess<3>(
+ 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 SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/IntegrationPointData.h b/ProcessLib/SmallDeformationNonlocal/IntegrationPointData.h
new file mode 100644
index 0000000000000000000000000000000000000000..cd83bb19da69ccbcfe34c929044a6de5067034f7
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/IntegrationPointData.h
@@ -0,0 +1,73 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "IntegrationPointDataNonlocalInterface.h"
+#include "MaterialLib/SolidModels/Ehlers.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <typename BMatricesType, typename ShapeMatricesType,
+ int DisplacementDim>
+struct IntegrationPointData final : public IntegrationPointDataNonlocalInterface
+{
+ explicit IntegrationPointData(
+ MaterialLib::Solids::Ehlers::SolidEhlers<DisplacementDim>&
+ solid_material)
+ : solid_material(solid_material),
+ material_state_variables(
+ solid_material.createMaterialStateVariables())
+ {
+ auto const msv =
+ static_cast<typename MaterialLib::Solids::Ehlers::StateVariables<
+ DisplacementDim>&>(*material_state_variables);
+
+ eps_p_V = &msv.eps_p.V;
+ eps_p_D_xx = &msv.eps_p.D[0];
+ }
+
+ typename BMatricesType::BMatrixType b_matrices;
+ typename BMatricesType::KelvinVectorType sigma, sigma_prev;
+ typename BMatricesType::KelvinVectorType eps, eps_prev;
+ double free_energy_density = 0;
+ double damage = 0; ///< isotropic damage
+ double damage_prev = 0; ///< \copydoc damage
+ double kappa_d_prev = 0; ///< \copydoc kappa_d
+
+ MaterialLib::Solids::Ehlers::SolidEhlers<DisplacementDim>& solid_material;
+ std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>
+ material_state_variables;
+
+ typename BMatricesType::KelvinMatrixType C;
+ typename ShapeMatricesType::NodalRowVectorType N;
+ typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
+
+ double const* eps_p_V; // Used for the secondary variables output.
+ double const* eps_p_D_xx; // Used for the secondary variables output.
+
+ void pushBackState()
+ {
+ eps_prev = eps;
+ sigma_prev = sigma;
+ damage_prev = damage;
+ kappa_d_prev = kappa_d;
+ material_state_variables->pushBackState();
+ }
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/IntegrationPointDataNonlocalInterface.h b/ProcessLib/SmallDeformationNonlocal/IntegrationPointDataNonlocalInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..98a874bdda3015a33c7aeb1fc6eefc09db4da740
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/IntegrationPointDataNonlocalInterface.h
@@ -0,0 +1,43 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+struct IntegrationPointDataNonlocalInterface;
+
+struct NonlocalIP final
+{
+ IntegrationPointDataNonlocalInterface* const ip_l_pointer;
+ double alpha_kl_times_w_l;
+ double distance2; ///< Squared distance to current integration point.
+};
+
+struct IntegrationPointDataNonlocalInterface
+{
+ virtual ~IntegrationPointDataNonlocalInterface() = default;
+
+ std::vector<NonlocalIP> non_local_assemblers;
+
+ double kappa_d = 0; ///< damage driving variable.
+ double integration_weight;
+ double nonlocal_internal_length;
+ Eigen::Vector3d coordinates;
+ bool active_self = false;
+ bool activated = false;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/IntegrationPointWriter.h b/ProcessLib/SmallDeformationNonlocal/IntegrationPointWriter.h
new file mode 100644
index 0000000000000000000000000000000000000000..18c911cb3f8ba96b93cc91d1a1aa7b20024759f8
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/IntegrationPointWriter.h
@@ -0,0 +1,87 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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>
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+struct SigmaIntegrationPointWriter final : public IntegrationPointWriter
+{
+ SigmaIntegrationPointWriter(
+ int const n_components,
+ int const integration_order,
+ std::function<std::vector<std::vector<double>>()>
+ callback)
+ : _callback(callback),
+ _integration_order(integration_order),
+ _n_components(n_components)
+ {
+ }
+
+ int numberOfComponents() const override { return _n_components; }
+ int integrationOrder() const override { return _integration_order; }
+
+ std::string name() const override
+ {
+ // TODO (naumov) remove ip suffix. Probably needs modification of the
+ // mesh properties, s.t. there is no "overlapping" with cell/point data.
+ // See getOrCreateMeshProperty.
+ return "sigma_ip";
+ }
+
+ std::vector<std::vector<double>> values() const override
+ {
+ return _callback();
+ }
+
+private:
+ std::function<std::vector<std::vector<double>>()> _callback;
+ int const _integration_order;
+ int const _n_components;
+};
+
+struct KappaDIntegrationPointWriter final : public IntegrationPointWriter
+{
+ KappaDIntegrationPointWriter(
+ int const integration_order,
+ std::function<std::vector<std::vector<double>>()>
+ callback)
+ : _callback(callback), _integration_order(integration_order)
+ {
+ }
+
+ int numberOfComponents() const override { return 1; }
+ int integrationOrder() const override { return _integration_order; }
+
+ std::string name() const override
+ {
+ // TODO (naumov) remove ip suffix. Probably needs modification of the
+ // mesh properties, s.t. there is no "overlapping" with cell/point data.
+ // See getOrCreateMeshProperty.
+ return "kappa_d_ip";
+ }
+
+ std::vector<std::vector<double>> values() const override
+ {
+ return _callback();
+ }
+
+private:
+ std::function<std::vector<std::vector<double>>()> _callback;
+ int const _integration_order;
+};
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/LocalAssemblerInterface.h b/ProcessLib/SmallDeformationNonlocal/LocalAssemblerInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..6cf49d7917cef2eeb60ea5be3a84b70fdb61c05c
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/LocalAssemblerInterface.h
@@ -0,0 +1,104 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <memory>
+#include <vector>
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+#include "IntegrationPointDataNonlocalInterface.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+struct SmallDeformationNonlocalLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+ virtual std::size_t setIPDataInitialConditions(
+ std::string const& name, double const* values,
+ int const integration_order) = 0;
+
+ virtual void setIPDataInitialConditionsFromCellData(
+ std::string const& name, std::vector<double> const& value) = 0;
+
+ virtual void computeCrackIntegral(
+ std::size_t mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table, GlobalVector const& x,
+ double& crack_volume) = 0;
+
+ virtual std::vector<double> const& getIntPtEpsPV(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+ virtual std::vector<double> const& getIntPtEpsPDXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+ virtual std::vector<double> getKappaD() const = 0;
+ virtual std::vector<double> const& getIntPtDamage(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtFreeEnergyDensity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> getSigma() const = 0;
+ virtual std::vector<double> const& getIntPtSigma(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilon(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const = 0;
+
+ // TODO move to NumLib::ExtrapolatableElement
+ virtual unsigned getNumberOfIntegrationPoints() const = 0;
+
+ virtual typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const&
+ getMaterialStateVariablesAt(int integration_point) const = 0;
+
+ virtual std::vector<double> const& getNodalValues(
+ std::vector<double>& nodal_values) const = 0;
+
+ virtual void nonlocal(
+ std::size_t const mesh_item_id,
+ std::vector<std::unique_ptr<
+ SmallDeformationNonlocalLocalAssemblerInterface>> const&
+ local_assemblers) = 0;
+
+ virtual void getIntegrationPointCoordinates(
+ Eigen::Vector3d const& coords,
+ std::vector<double>& distances) const = 0;
+
+ virtual IntegrationPointDataNonlocalInterface* getIPDataPtr(
+ int const ip) = 0;
+};
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..5fa904f23a5006af1de5d93f604bc756fb521883
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
@@ -0,0 +1,867 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <algorithm>
+#include <limits>
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/Ehlers.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "MeshLib/findElementsWithinRadius.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/Divergence.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "Damage.h"
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "SmallDeformationNonlocalProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+/// Used for the extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+ std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+class SmallDeformationNonlocalLocalAssembler
+ : public SmallDeformationNonlocalLocalAssemblerInterface<DisplacementDim>
+{
+public:
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ 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;
+
+ SmallDeformationNonlocalLocalAssembler(
+ SmallDeformationNonlocalLocalAssembler const&) = delete;
+ SmallDeformationNonlocalLocalAssembler(
+ SmallDeformationNonlocalLocalAssembler&&) = delete;
+
+ SmallDeformationNonlocalLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ SmallDeformationNonlocalProcessData<DisplacementDim>& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, DisplacementDim>(
+ e, is_axially_symmetric, _integration_method);
+
+ auto& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials,
+ _process_data.material_ids,
+ e.getID());
+ auto* ehlers_solid_material = dynamic_cast<
+ MaterialLib::Solids::Ehlers::SolidEhlers<DisplacementDim>*>(
+ &solid_material);
+ if (ehlers_solid_material == nullptr)
+ {
+ OGS_FATAL(
+ "The SmallDeformationNonlocalLocal process supports only "
+ "Ehlers material at the moment. For other materials the "
+ "interface must be extended first.");
+ }
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data.emplace_back(*ehlers_solid_material);
+ auto& ip_data = _ip_data[ip];
+ auto const& sm = shape_matrices[ip];
+ _ip_data[ip].integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ;
+
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
+
+ // Initialize current time step values
+ ip_data.sigma.setZero(MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value);
+ ip_data.eps.setZero(MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value);
+
+ // Previous time step values are not initialized and are set later.
+ ip_data.sigma_prev.resize(
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value);
+ ip_data.eps_prev.resize(
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value);
+
+ _secondary_data.N[ip] = shape_matrices[ip].N;
+
+ ip_data.coordinates = getSingleIntegrationPointCoordinates(ip);
+ }
+ }
+
+ std::size_t setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order) override
+ {
+ if (integration_order !=
+ static_cast<int>(_integration_method.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element %d is different from "
+ "the integration order in the initial condition.",
+ _element.getID());
+ }
+
+ if (name == "sigma_ip")
+ {
+ return setSigma(values);
+ }
+
+ if (name == "kappa_d_ip")
+ {
+ return setKappaD(values);
+ }
+
+ return 0;
+ }
+
+ void setIPDataInitialConditionsFromCellData(
+ std::string const& name, std::vector<double> const& value) override
+ {
+ if (name == "kappa_d_ip")
+ {
+ if (value.size() != 1)
+ {
+ OGS_FATAL(
+ "CellData for kappa_d initial conditions has wrong number "
+ "of components. 1 expected, got %d.",
+ value.size());
+ }
+ setKappaD(value[0]);
+ }
+ }
+
+ double alpha_0(double const distance2) const
+ {
+ double const internal_length2 = _process_data.internal_length_squared;
+ return (distance2 > internal_length2)
+ ? 0
+ : (1 - distance2 / (internal_length2)) *
+ (1 - distance2 / (internal_length2));
+ }
+
+ void nonlocal(
+ std::size_t const /*mesh_item_id*/,
+ std::vector<
+ std::unique_ptr<SmallDeformationNonlocalLocalAssemblerInterface<
+ DisplacementDim>>> const& local_assemblers) override
+ {
+ auto const search_element_ids = MeshLib::findElementsWithinRadius(
+ _element, _process_data.internal_length_squared);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ std::vector<double> distances; // Cache for ip-ip distances.
+ //
+ // For every integration point in this element collect the neighbouring
+ // integration points falling in given radius (internal length) and
+ // compute the alpha_kl weight.
+ //
+ for (unsigned k = 0; k < n_integration_points; k++)
+ {
+ //
+ // Collect the integration points.
+ //
+
+ auto const& xyz = _ip_data[k].coordinates;
+
+ // For all neighbors of element
+ for (auto const search_element_id : search_element_ids)
+ {
+ auto const& la = local_assemblers[search_element_id];
+ la->getIntegrationPointCoordinates(xyz, distances);
+ for (int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
+ {
+ if (distances[ip] >= _process_data.internal_length_squared)
+ continue;
+ // save into current ip_k
+ _ip_data[k].non_local_assemblers.push_back(
+ {la->getIPDataPtr(ip),
+ std::numeric_limits<double>::quiet_NaN(),
+ distances[ip]});
+ }
+ }
+ if (_ip_data[k].non_local_assemblers.size() == 0)
+ {
+ OGS_FATAL("no neighbours found!");
+ }
+
+ double a_k_sum_m = 0;
+ for (auto const& tuple : _ip_data[k].non_local_assemblers)
+ {
+ double const distance2_m = tuple.distance2;
+
+ auto const& w_m = tuple.ip_l_pointer->integration_weight;
+
+ a_k_sum_m += w_m * alpha_0(distance2_m);
+ }
+
+ //
+ // Calculate alpha_kl =
+ // alpha_0(|x_k - x_l|) / int_{m \in ip} alpha_0(|x_k - x_m|)
+ //
+ for (auto& tuple : _ip_data[k].non_local_assemblers)
+ {
+ double const distance2_l = tuple.distance2;
+ double const a_kl = alpha_0(distance2_l) / a_k_sum_m;
+
+ // Store the a_kl already multiplied with the integration
+ // weight of that l integration point.
+ auto const w_l = tuple.ip_l_pointer->integration_weight;
+ tuple.alpha_kl_times_w_l = a_kl * w_l;
+ }
+ }
+ }
+
+ Eigen::Vector3d getSingleIntegrationPointCoordinates(
+ int integration_point) const
+ {
+ auto const& N = _secondary_data.N[integration_point];
+
+ Eigen::Vector3d xyz = Eigen::Vector3d::Zero(); // Resulting coordinates
+ auto* nodes = _element.getNodes();
+ for (int i = 0; i < N.size(); ++i)
+ {
+ Eigen::Map<Eigen::Vector3d const> node_coordinates{
+ nodes[i]->getCoords(), 3};
+ xyz += node_coordinates * N[i];
+ }
+ return xyz;
+ }
+
+ /// For each of the current element's integration points the squared
+ /// distance from the current integration point is computed and stored in
+ /// the given distances cache.
+ void getIntegrationPointCoordinates(
+ Eigen::Vector3d const& coords,
+ std::vector<double>& distances) const override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ distances.resize(n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& xyz = _ip_data[ip].coordinates;
+ distances[ip] = (xyz - coords).squaredNorm();
+ }
+ }
+
+ void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& /*local_b_data*/) override
+ {
+ OGS_FATAL(
+ "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
+ "is not "
+ "implemented.");
+ }
+
+ void preAssemble(double const t,
+ std::vector<double> const& local_x) override
+ {
+ auto const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
+ 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& state = _ip_data[ip].material_state_variables;
+ double const& damage_prev = _ip_data[ip].damage_prev;
+
+ eps.noalias() =
+ B *
+ Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
+ local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
+
+ // sigma is for plastic part only.
+ std::unique_ptr<
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>>
+ new_C;
+ std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>
+ new_state;
+
+ // Compute sigma_eff from damage total stress sigma
+ using KelvinVectorType = typename BMatricesType::KelvinVectorType;
+ KelvinVectorType const sigma_eff_prev =
+ sigma_prev / (1. - damage_prev);
+
+ auto&& solution = _ip_data[ip].solid_material.integrateStress(
+ t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
+ *state, _process_data.reference_temperature);
+
+ if (!solution)
+ OGS_FATAL("Computation of local constitutive relation failed.");
+
+ std::tie(sigma, state, C) = std::move(*solution);
+
+ /// Compute only the local kappa_d.
+ {
+ auto const& ehlers_material =
+ static_cast<MaterialLib::Solids::Ehlers::SolidEhlers<
+ DisplacementDim> const&>(_ip_data[ip].solid_material);
+ auto const damage_properties =
+ ehlers_material.evaluatedDamageProperties(t, x_position);
+ auto const material_properties =
+ ehlers_material.evaluatedMaterialProperties(t, x_position);
+
+ // Ehlers material state variables
+ auto& state_vars =
+ static_cast<MaterialLib::Solids::Ehlers::StateVariables<
+ DisplacementDim>&>(
+ *_ip_data[ip].material_state_variables);
+
+ double const eps_p_eff_diff =
+ state_vars.eps_p.eff - state_vars.eps_p_prev.eff;
+
+ _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
+ eps_p_eff_diff, sigma, _ip_data[ip].kappa_d_prev,
+ damage_properties.h_d, material_properties);
+
+ if (!_ip_data[ip].active_self)
+ {
+ _ip_data[ip].active_self |= _ip_data[ip].kappa_d > 0;
+ if (_ip_data[ip].active_self)
+ {
+ for (auto const& tuple :
+ _ip_data[ip].non_local_assemblers)
+ {
+ // Activate the integration point.
+ tuple.ip_l_pointer->activated = true;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void assembleWithJacobian(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double> const& /*local_xdot*/,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_b_data,
+ std::vector<double>& local_Jac_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+
+ auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
+ local_b_data, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ // Non-local integration.
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
+
+ auto& sigma = _ip_data[ip].sigma;
+ auto& C = _ip_data[ip].C;
+ double& damage = _ip_data[ip].damage;
+
+ {
+ double nonlocal_kappa_d = 0;
+
+ if (_ip_data[ip].active_self || _ip_data[ip].activated)
+ {
+ for (auto const& tuple : _ip_data[ip].non_local_assemblers)
+ {
+ // Get local variable for the integration point l.
+ double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
+ double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
+ nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
+ }
+ }
+
+ auto const& ehlers_material =
+ static_cast<MaterialLib::Solids::Ehlers::SolidEhlers<
+ DisplacementDim> const&>(_ip_data[ip].solid_material);
+
+ //
+ // Overnonlocal formulation
+ //
+ // See (Di Luzio & Bazant 2005, IJSS) for details.
+ // The implementation would go here and would be for a given
+ // gamma_nonlocal:
+ //
+ // Update nonlocal damage with local damage (scaled with 1 -
+ // \gamma_{nonlocal}) for the current integration point and the
+ // nonlocal integral part.
+ // nonlocal_kappa_d = (1. - gamma_nonlocal) * kappa_d +
+ // gamma_nonlocal * nonlocal_kappa_d;
+
+ nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
+
+ // Update damage based on nonlocal kappa_d
+ {
+ auto const damage_properties =
+ ehlers_material.evaluatedDamageProperties(t,
+ x_position);
+ damage = calculateDamage(nonlocal_kappa_d,
+ damage_properties.alpha_d,
+ damage_properties.beta_d);
+ damage = std::max(0., damage);
+ }
+ sigma = sigma * (1. - damage);
+ }
+
+ local_b.noalias() -= B.transpose() * sigma * w;
+ local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
+ }
+ }
+
+ 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();
+ }
+ }
+
+ void computeCrackIntegral(std::size_t mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ GlobalVector const& x,
+ double& crack_volume) override
+ {
+ auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+ auto local_x = x.get(indices);
+
+ auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
+ local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
+
+ int const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& dNdx = _ip_data[ip].dNdx;
+ auto const& d = _ip_data[ip].damage;
+ auto const& w = _ip_data[ip].integration_weight;
+
+ double const div_u =
+ Deformation::divergence<DisplacementDim,
+ ShapeFunction::NPOINTS>(u, dNdx);
+ crack_volume += div_u * d * w;
+ }
+ }
+
+ 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& getNodalValues(
+ std::vector<double>& nodal_values) const override
+ {
+ nodal_values.clear();
+ auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
+ nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
+ auto& sigma = _ip_data[ip].sigma;
+
+ local_b.noalias() += B.transpose() * sigma * w;
+ }
+
+ return nodal_values;
+ }
+
+ std::vector<double> const& getIntPtFreeEnergyDensity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ cache.push_back(ip_data.free_energy_density);
+ }
+
+ return cache;
+ }
+
+ std::vector<double> const& getIntPtEpsPV(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ cache.push_back(*ip_data.eps_p_V);
+ }
+
+ return cache;
+ }
+
+ std::vector<double> const& getIntPtEpsPDXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ cache.push_back(*ip_data.eps_p_D_xx);
+ }
+
+ return cache;
+ }
+
+ std::vector<double> const& getIntPtSigma(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ static const int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
+
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ auto const& sigma = _ip_data[ip].sigma;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+ }
+
+ return cache;
+ }
+
+ virtual std::vector<double> const& getIntPtEpsilon(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
+
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ auto const& eps = _ip_data[ip].eps;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+ }
+
+ return cache;
+ }
+
+ std::size_t setSigma(double const* values)
+ {
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const n_integration_points = _ip_data.size();
+
+ std::vector<double> ip_sigma_values;
+ auto sigma_values =
+ Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
+ Eigen::ColMajor> const>(
+ values, kelvin_vector_size, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ _ip_data[ip].sigma =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector(
+ sigma_values.col(ip));
+ }
+
+ return n_integration_points;
+ }
+
+ // TODO (naumov) This method is same as getIntPtSigma but for arguments and
+ // the ordering of the cache_mat.
+ // There should be only one.
+ std::vector<double> getSigma() const override
+ {
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const n_integration_points = _ip_data.size();
+
+ std::vector<double> ip_sigma_values;
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+ ip_sigma_values, n_integration_points, kelvin_vector_size);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& sigma = _ip_data[ip].sigma;
+ cache_mat.row(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+ }
+
+ return ip_sigma_values;
+ }
+
+ std::size_t setKappaD(double const* values)
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ _ip_data[ip].kappa_d = values[ip];
+ }
+ return n_integration_points;
+ }
+
+ void setKappaD(double value)
+ {
+ for (auto& ip_data : _ip_data)
+ {
+ ip_data.kappa_d = value;
+ }
+ }
+ std::vector<double> getKappaD() const override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ std::vector<double> result_values;
+ result_values.resize(n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ result_values[ip] = _ip_data[ip].kappa_d;
+ }
+
+ return result_values;
+ }
+
+ std::vector<double> const& getIntPtDamage(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const override
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ cache.push_back(ip_data.damage);
+ }
+
+ return cache;
+ }
+
+ unsigned getNumberOfIntegrationPoints() const override
+ {
+ return _integration_method.getNumberOfPoints();
+ }
+
+ typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const&
+ getMaterialStateVariablesAt(int const integration_point) const override
+ {
+ return *_ip_data[integration_point].material_state_variables;
+ }
+
+private:
+ std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
+ std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ 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;
+ }
+
+ std::vector<double> const& getIntPtEpsilon(
+ std::vector<double>& cache, std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ if (component < 3) // xx, yy, zz components
+ cache.push_back(ip_data.eps[component]);
+ else // mixed xy, yz, xz components
+ cache.push_back(ip_data.eps[component] / std::sqrt(2));
+ }
+
+ return cache;
+ }
+
+ IntegrationPointDataNonlocalInterface*
+ getIPDataPtr(int const ip) override
+ {
+ return &_ip_data[ip];
+ }
+
+private:
+ SmallDeformationNonlocalProcessData<DisplacementDim>& _process_data;
+
+ std::vector<
+ IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ BMatricesType, ShapeMatricesType, DisplacementDim>>>
+ _ip_data;
+
+ IntegrationMethod const _integration_method;
+ MeshLib::Element const& _element;
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+ bool const _is_axially_symmetric;
+
+ static const int displacement_size =
+ ShapeFunction::NPOINTS * DisplacementDim;
+};
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c65ae6dacd3954ee7f7e002809061a623692dee4
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
@@ -0,0 +1,316 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "SmallDeformationNonlocalProcess.h"
+
+#include <nlohmann/json.hpp>
+#include <iostream>
+
+// Reusing local assembler creation code.
+#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+SmallDeformationNonlocalProcess<DisplacementDim>::
+ SmallDeformationNonlocalProcess(
+ 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::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ SmallDeformationNonlocalProcessData<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))
+{
+ _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<SigmaIntegrationPointWriter>(
+ static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+ integration_order, [this]() {
+ // Result containing integration point data for each local
+ // assembler.
+ std::vector<std::vector<double>> result;
+ result.resize(_local_assemblers.size());
+
+ for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+ {
+ auto const& local_asm = *_local_assemblers[i];
+
+ result[i] = local_asm.getSigma();
+ }
+
+ return result;
+ }));
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<KappaDIntegrationPointWriter>(
+ integration_order, [this]() {
+ // Result containing integration point data for each local
+ // assembler.
+ std::vector<std::vector<double>> result;
+ result.resize(_local_assemblers.size());
+
+ for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+ {
+ auto const& local_asm = *_local_assemblers[i];
+
+ result[i] = local_asm.getKappaD();
+ }
+
+ return result;
+ }));
+}
+
+template <int DisplacementDim>
+bool SmallDeformationNonlocalProcess<DisplacementDim>::isLinear() const
+{
+ return false;
+}
+
+template <int DisplacementDim>
+void SmallDeformationNonlocalProcess<DisplacementDim>::
+ initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ // Reusing local assembler creation code.
+ ProcessLib::SmallDeformation::createLocalAssemblers<
+ DisplacementDim, SmallDeformationNonlocalLocalAssembler>(
+ 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<MeshLib::MeshSubset> all_mesh_subsets_single_component{
+ *_mesh_subset_all_nodes};
+ _local_to_global_index_map_single_component =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets_single_component),
+ // by location order is needed for output
+ NumLib::ComponentOrder::BY_LOCATION);
+
+ Process::_secondary_variables.addSecondaryVariable(
+ "sigma",
+ makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+ DisplacementDim>::RowsAtCompileTime,
+ getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigma));
+
+ Process::_secondary_variables.addSecondaryVariable(
+ "epsilon",
+ makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+ DisplacementDim>::RowsAtCompileTime,
+ getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilon));
+
+ Process::_secondary_variables.addSecondaryVariable(
+ "eps_p_V",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsPV));
+ Process::_secondary_variables.addSecondaryVariable(
+ "eps_p_D_xx",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsPDXX));
+
+ Process::_secondary_variables.addSecondaryVariable(
+ "damage",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtDamage));
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::nonlocal, _local_assemblers,
+ _local_assemblers);
+
+
+ // Set initial conditions for integration point data.
+ for (auto const& ip_writer : _integration_point_writer)
+ {
+ auto const& name = ip_writer->name();
+ // First check the field data, which is used for restart.
+ if (mesh.getProperties().existsPropertyVector<double>(name))
+ {
+ auto const& mesh_property =
+ *mesh.getProperties().template getPropertyVector<double>(name);
+
+ // The mesh property must be defined on integration points.
+ if (mesh_property.getMeshItemType() !=
+ MeshLib::MeshItemType::IntegrationPoint)
+ {
+ continue;
+ }
+
+ auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
+
+ // Check the number of components.
+ if (ip_meta_data.n_components !=
+ mesh_property.getNumberOfComponents())
+ {
+ OGS_FATAL(
+ "Different number of components in meta data (%d) than in "
+ "the integration point field data for '%s': %d.",
+ ip_meta_data.n_components, name.c_str(),
+ mesh_property.getNumberOfComponents());
+ }
+
+ // Now we have a properly named vtk's field data array and the
+ // corresponding meta data.
+ std::size_t position = 0;
+ for (auto& local_asm : _local_assemblers)
+ {
+ std::size_t const integration_points_read =
+ local_asm->setIPDataInitialConditions(
+ name, &mesh_property[position],
+ ip_meta_data.integration_order);
+ if (integration_points_read == 0)
+ {
+ OGS_FATAL(
+ "No integration points read in the integration point "
+ "initial conditions set function.");
+ }
+ position += integration_points_read * ip_meta_data.n_components;
+ }
+ }
+ else if (mesh.getProperties().existsPropertyVector<double>(name +
+ "_ic"))
+ { // Try to find cell data with '_ic' suffix
+ auto const& mesh_property =
+ *mesh.getProperties().template getPropertyVector<double>(name +
+ "_ic");
+ if (mesh_property.getMeshItemType() != MeshLib::MeshItemType::Cell)
+ {
+ continue;
+ }
+
+ // Now we have a vtk's cell data array containing the initial
+ // conditions for the corresponding integration point writer.
+
+ // For each assembler use the single cell value for all
+ // integration points.
+ for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+ {
+ auto& local_asm = _local_assemblers[i];
+
+ std::vector<double> value(
+ &mesh_property[i],
+ &mesh_property[i] + mesh_property.getNumberOfComponents());
+ // TODO (naumov) Check sizes / read size / etc.
+ // OR reconstruct dimensions from size / component =
+ // ip_points
+ local_asm->setIPDataInitialConditionsFromCellData(name, value);
+ }
+ }
+ }
+}
+
+template <int DisplacementDim>
+void SmallDeformationNonlocalProcess<DisplacementDim>::assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble SmallDeformationNonlocalProcess.");
+
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_table = {std::ref(*_local_to_global_index_map)};
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ dof_table, t, x, M, K, b, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+void SmallDeformationNonlocalProcess<
+ DisplacementDim>::preAssembleConcreteProcess(const double t,
+ GlobalVector const& x)
+{
+ DBUG("preAssemble SmallDeformationNonlocalProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::preAssemble,
+ _local_assemblers, *_local_to_global_index_map, t, x);
+}
+
+template <int DisplacementDim>
+void SmallDeformationNonlocalProcess<DisplacementDim>::
+ 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)
+{
+ DBUG("AssembleWithJacobian SmallDeformationNonlocalProcess.");
+
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_table = {std::ref(*_local_to_global_index_map)};
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, dof_table, t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac,
+ _coupled_solutions);
+
+ b.copyValues(*_nodal_forces);
+ std::transform(_nodal_forces->begin(), _nodal_forces->end(),
+ _nodal_forces->begin(), [](double val) { return -val; });
+}
+
+template <int DisplacementDim>
+void SmallDeformationNonlocalProcess<
+ DisplacementDim>::preTimestepConcreteProcess(GlobalVector const& x,
+ double const t,
+ double const dt,
+ int const /*process_id*/)
+{
+ DBUG("PreTimestep SmallDeformationNonlocalProcess.");
+
+ _process_data.dt = dt;
+ _process_data.t = t;
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::preTimestep, _local_assemblers,
+ *_local_to_global_index_map, x, t, dt);
+}
+
+template <int DisplacementDim>
+NumLib::IterationResult
+SmallDeformationNonlocalProcess<DisplacementDim>::postIterationConcreteProcess(
+ GlobalVector const& x)
+{
+ _process_data.crack_volume_old = _process_data.crack_volume;
+ _process_data.crack_volume = 0.0;
+
+ DBUG("PostNonLinearSolver crack volume computation.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::computeCrackIntegral, _local_assemblers,
+ *_local_to_global_index_map, x, _process_data.crack_volume);
+
+ INFO("Integral of crack: %g", _process_data.crack_volume);
+
+ return NumLib::IterationResult::SUCCESS;
+}
+
+template class SmallDeformationNonlocalProcess<2>;
+template class SmallDeformationNonlocalProcess<3>;
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..438f9e9379ec1196a9e2ea92f50c8b825720fb54
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
@@ -0,0 +1,94 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <cassert>
+
+#include "BaseLib/Functional.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Process.h"
+
+#include "IntegrationPointWriter.h"
+#include "SmallDeformationNonlocalFEM.h"
+#include "SmallDeformationNonlocalProcessData.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+class SmallDeformationNonlocalProcess final : public Process
+{
+public:
+ SmallDeformationNonlocalProcess(
+ 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::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ SmallDeformationNonlocalProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller);
+
+ //! \name ODESystem interface
+ //! @{
+
+ bool isLinear() const override;
+ //! @}
+
+private:
+ 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;
+
+ void preAssembleConcreteProcess(const double t,
+ GlobalVector const& x) override;
+
+ 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;
+
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt,
+ int const /*process_id*/) override;
+
+ NumLib::IterationResult postIterationConcreteProcess(
+ GlobalVector const& x) override;
+
+private:
+ SmallDeformationNonlocalProcessData<DisplacementDim> _process_data;
+
+ using LocalAssemblerInterface =
+ SmallDeformationNonlocalLocalAssemblerInterface<DisplacementDim>;
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+ _local_to_global_index_map_single_component;
+
+ MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+};
+
+extern template class ProcessLib::SmallDeformationNonlocal::
+ SmallDeformationNonlocalProcess<2>;
+extern template class ProcessLib::SmallDeformationNonlocal::
+ SmallDeformationNonlocalProcess<3>;
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcessData.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..716eaaea2338289c8507caf9b5803bc631b32774
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcessData.h
@@ -0,0 +1,92 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <memory>
+#include <utility>
+
+#include <Eigen/Eigen>
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+}
+} // namespace MaterialLib
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+template <int DisplacementDim>
+struct SmallDeformationNonlocalProcessData
+{
+ SmallDeformationNonlocalProcessData(
+ MeshLib::PropertyVector<int> const* const material_ids_,
+ std::map<int,
+ std::unique_ptr<
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
+ solid_materials_,
+ Parameter<double> const& solid_density_,
+ Eigen::Matrix<double, DisplacementDim, 1>
+ specific_body_force_,
+ double const reference_temperature_,
+ double const internal_length_)
+ : material_ids(material_ids_),
+ solid_materials{std::move(solid_materials_)},
+ solid_density(solid_density_),
+ specific_body_force(std::move(specific_body_force_)),
+ reference_temperature(reference_temperature_),
+ internal_length_squared(internal_length_ * internal_length_)
+ {
+ }
+
+ SmallDeformationNonlocalProcessData(
+ SmallDeformationNonlocalProcessData&& other) = default;
+
+ //! Copies are forbidden.
+ SmallDeformationNonlocalProcessData(
+ SmallDeformationNonlocalProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(SmallDeformationNonlocalProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(SmallDeformationNonlocalProcessData&&) = delete;
+
+ MeshLib::PropertyVector<int> const* const material_ids;
+
+ std::map<
+ int,
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
+ solid_materials;
+ /// Solid's density. A scalar quantity, Parameter<double>.
+ Parameter<double> const& solid_density;
+ /// Specific body forces applied to the solid.
+ /// It is usually used to apply gravitational forces.
+ /// A vector of displacement dimension's length.
+ Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+
+ double crack_volume_old = 0.0;
+ double crack_volume = 0.0;
+
+ double dt = 0;
+ double t = 0;
+ double const reference_temperature;
+ double const internal_length_squared;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib