From 573e2b426b0e61366252e3488d790e0124446fb6 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Wed, 7 Nov 2018 10:44:09 +0100
Subject: [PATCH] [Mat] added a solid model interfacing with MFront
---
MaterialLib/SolidModels/MFront/CMakeLists.txt | 26 ++
.../SolidModels/MFront/CreateMFront.cpp | 222 ++++++++++++
MaterialLib/SolidModels/MFront/CreateMFront.h | 38 ++
MaterialLib/SolidModels/MFront/MFront.cpp | 336 ++++++++++++++++++
MaterialLib/SolidModels/MFront/MFront.h | 116 ++++++
5 files changed, 738 insertions(+)
create mode 100644 MaterialLib/SolidModels/MFront/CMakeLists.txt
create mode 100644 MaterialLib/SolidModels/MFront/CreateMFront.cpp
create mode 100644 MaterialLib/SolidModels/MFront/CreateMFront.h
create mode 100644 MaterialLib/SolidModels/MFront/MFront.cpp
create mode 100644 MaterialLib/SolidModels/MFront/MFront.h
diff --git a/MaterialLib/SolidModels/MFront/CMakeLists.txt b/MaterialLib/SolidModels/MFront/CMakeLists.txt
new file mode 100644
index 00000000000..a693ff0a812
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/CMakeLists.txt
@@ -0,0 +1,26 @@
+set(SOURCES CreateMFront.cpp CreateMFront.h)
+
+if(OGS_USE_MFRONT)
+ list(APPEND SOURCES MFront.cpp MFront.h)
+endif()
+
+add_library(MaterialLib_SolidModels_MFront ${SOURCES})
+
+if(BUILD_SHARED_LIBS)
+ install(TARGETS MaterialLib_SolidModels_MFront LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+
+target_link_libraries(MaterialLib_SolidModels_MFront
+ PUBLIC BaseLib NumLib logog
+ PRIVATE MathLib MeshLib
+)
+
+if (OGS_USE_MFRONT)
+ target_include_directories(MaterialLib_SolidModels_MFront
+ PUBLIC ${MGIS_INCLUDE_DIR}
+ )
+ target_link_libraries(MaterialLib_SolidModels_MFront
+ PUBLIC ${MGIS_LIBRARY}
+ )
+ target_compile_definitions(MaterialLib_SolidModels_MFront PRIVATE OGS_USE_MFRONT)
+endif()
diff --git a/MaterialLib/SolidModels/MFront/CreateMFront.cpp b/MaterialLib/SolidModels/MFront/CreateMFront.cpp
new file mode 100644
index 00000000000..0cc87c43b46
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/CreateMFront.cpp
@@ -0,0 +1,222 @@
+/**
+ * \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 "CreateMFront.h"
+
+#ifdef OGS_USE_MFRONT
+
+#include "MFront.h"
+
+#include "BaseLib/FileTools.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace
+{
+/// Prints info about MFront variables.
+void varInfo(std::string const& msg,
+ std::vector<mgis::behaviour::Variable> const& vars,
+ mgis::behaviour::Hypothesis hypothesis)
+{
+ INFO("#%s: %lu (array size %lu).",
+ msg.c_str(),
+ vars.size(),
+ mgis::behaviour::getArraySize(vars, hypothesis));
+ for (auto& var : vars)
+ {
+ INFO(" --> type `%s' with name `%s', size %lu, offset %lu.",
+ MaterialLib::Solids::MFront::varTypeToString(var.type),
+ var.name.c_str(),
+ mgis::behaviour::getVariableSize(var, hypothesis),
+ mgis::behaviour::getVariableOffset(vars, var.name, hypothesis));
+ }
+}
+} // anonymous namespace
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ INFO("### MFRONT ########################################################");
+
+ //! \ogs_file_param{material__solid__constitutive_relation__type}
+ config.checkConfigParameter("type", "MFront");
+
+ auto const lib_path = BaseLib::joinPaths(
+ BaseLib::getProjectDirectory(),
+ //! \ogs_file_param{material__solid__constitutive_relation__MFront__library}
+ config.getConfigParameter<std::string>("library"));
+ auto const behaviour_name =
+ //! \ogs_file_param{material__solid__constitutive_relation__MFront__behaviour}
+ config.getConfigParameter<std::string>("behaviour");
+
+ static_assert(DisplacementDim == 2 || DisplacementDim == 3,
+ "Given DisplacementDim not supported.");
+
+ mgis::behaviour::Hypothesis hypothesis;
+ if (DisplacementDim == 2)
+ {
+ // TODO support the axial symmetry modelling hypothesis.
+ WARN(
+ "The model is defined in 2D. On the material level currently a "
+ "plane strain setting is used. In particular it is not checked if "
+ "axial symmetry or plane stress are assumed. Special material "
+ "behaviour for these settings is currently not supported.");
+ hypothesis = mgis::behaviour::Hypothesis::PLANESTRAIN;
+ }
+ else if (DisplacementDim == 3)
+ {
+ hypothesis = mgis::behaviour::Hypothesis::TRIDIMENSIONAL;
+ }
+
+ auto behaviour =
+ mgis::behaviour::load(lib_path, behaviour_name, hypothesis);
+
+ INFO("Behaviour: `%s'.", behaviour.behaviour.c_str());
+ INFO("Hypothesis: `%s'.", mgis::behaviour::toString(hypothesis));
+ INFO("Source: `%s'.", behaviour.source.c_str());
+ INFO("TFEL version: `%s'.", behaviour.tfel_version.c_str());
+ INFO("Behaviour type: `%s'.", btypeToString(behaviour.btype));
+ INFO("Kinematic: `%s'.", toString(behaviour.kinematic));
+ INFO("Symmetry: `%s'.", toString(behaviour.symmetry));
+
+ varInfo("Mat. props.", behaviour.mps, hypothesis);
+ varInfo("Gradients", behaviour.gradients, hypothesis);
+ varInfo("Thdyn. forces", behaviour.thermodynamic_forces, hypothesis);
+ varInfo("Int. StVars.", behaviour.isvs, hypothesis);
+ varInfo("Ext. StVars.", behaviour.esvs, hypothesis);
+
+ // TODO read parameters from prj file, not yet (2018-11-05) supported by
+ // MGIS library.
+ varInfo("Parameters", behaviour.parameters, hypothesis);
+
+ std::vector<ProcessLib::Parameter<double> const*> material_properties;
+
+ if (!behaviour.mps.empty())
+ {
+ std::map<std::string, std::string> map_name_to_param;
+
+ // gather material properties from the prj file
+ //! \ogs_file_param{material__solid__constitutive_relation__MFront__material_properties}
+ auto const mps_config = config.getConfigSubtree("material_properties");
+ for (
+ auto const mp_config :
+ //! \ogs_file_param{material__solid__constitutive_relation__MFront__material_properties__material_property}
+ mps_config.getConfigParameterList("material_property"))
+ {
+ //! \ogs_file_attr{material__solid__constitutive_relation__MFront__material_properties__material_property__name}
+ auto name = mp_config.getConfigAttribute<std::string>("name");
+ auto const param_name =
+ //! \ogs_file_attr{material__solid__constitutive_relation__MFront__material_properties__material_property__parameter}
+ mp_config.getConfigAttribute<std::string>("parameter");
+
+ map_name_to_param.emplace(std::move(name), std::move(param_name));
+ }
+
+ for (auto& mp : behaviour.mps)
+ {
+ auto const it = map_name_to_param.find(mp.name);
+ if (it == map_name_to_param.end())
+ OGS_FATAL(
+ "Material Property `%s' has not been configured in the "
+ "project file.",
+ mp.name.c_str());
+
+ auto const param_name = it->second;
+ auto const num_comp =
+ mgis::behaviour::getVariableSize(mp, hypothesis);
+ auto const* param = &ProcessLib::findParameter<double>(
+ param_name, parameters, num_comp);
+
+ INFO("Using OGS parameter `%s' for material property `%s'.",
+ param_name.c_str(), mp.name.c_str());
+
+ using V = mgis::behaviour::Variable;
+ if (mp.type == V::STENSOR || mp.type == V::TENSOR)
+ {
+ WARN(
+ "Material property `%s' is a tensorial quantity. You, the "
+ "user, have to make sure that the component order of "
+ "parameter `%s' matches the one required by MFront!",
+ mp.name.c_str(), param_name.c_str());
+ }
+
+ material_properties.push_back(param);
+ map_name_to_param.erase(it);
+ }
+
+ if (!map_name_to_param.empty())
+ {
+ ERR("Some material parameters that were configured are not used by "
+ "the material model.");
+ ERR("These parameters are:");
+
+ for (auto& e : map_name_to_param)
+ {
+ ERR(" name: `%s', parameter: `%s'.", e.first.c_str(),
+ e.second.c_str());
+ }
+
+ OGS_FATAL(
+ "Configuration errors occured. Please fix the project file.");
+ }
+ }
+
+ INFO("### MFRONT END ####################################################");
+
+ return std::make_unique<MFront<DisplacementDim>>(
+ std::move(behaviour), std::move(material_properties));
+}
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
+
+#else // OGS_USE_MFRONT
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
+ std::vector<
+ std::unique_ptr<ProcessLib::ParameterBase>> const& /*parameters*/,
+ BaseLib::ConfigTree const& /*config*/)
+{
+ OGS_FATAL("OpenGeoSys has not been build with MFront support.");
+}
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
+
+#endif // OGS_USE_MFRONT
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+template std::unique_ptr<MechanicsBase<2>> createMFront<2>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+template std::unique_ptr<MechanicsBase<3>> createMFront<3>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/MFront/CreateMFront.h b/MaterialLib/SolidModels/MFront/CreateMFront.h
new file mode 100644
index 00000000000..4df0a9a8236
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/CreateMFront.h
@@ -0,0 +1,38 @@
+/**
+ * \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 "BaseLib/ConfigTree.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MechanicsBase<2>> createMFront<2>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+extern template std::unique_ptr<MechanicsBase<3>> createMFront<3>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/MFront/MFront.cpp b/MaterialLib/SolidModels/MFront/MFront.cpp
new file mode 100644
index 00000000000..a6c15cdbcf8
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MFront.cpp
@@ -0,0 +1,336 @@
+/**
+ * \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 "MFront.h"
+
+#include <MGIS/Behaviour/Integrate.hxx>
+
+namespace
+{
+/// Converts between OGSes and MFront's Kelvin vector indices.
+constexpr std::ptrdiff_t OGSToMFront(std::ptrdiff_t i)
+{
+ // MFront: 11 22 33 12 13 23
+ // OGS: 11 22 33 12 23 13
+ if (i < 4)
+ return i;
+ if (i == 4)
+ return 5;
+ return 4;
+}
+/// Converts between OGSes and MFront's Kelvin vector indices.
+constexpr std::ptrdiff_t MFrontToOGS(std::ptrdiff_t i)
+{
+ // MFront: 11 22 33 12 13 23
+ // OGS: 11 22 33 12 23 13
+ return OGSToMFront(i); // Same algorithm: indices 4 and 5 swapped.
+}
+
+/// Converts between OGSes and MFront's Kelvin vectors and matrices.
+template <typename Derived>
+typename Derived::PlainObject OGSToMFront(Eigen::DenseBase<Derived> const& m)
+{
+ static_assert(Derived::RowsAtCompileTime != Eigen::Dynamic, "Error");
+ static_assert(Derived::ColsAtCompileTime != Eigen::Dynamic, "Error");
+
+ typename Derived::PlainObject n;
+
+ // optimal for row-major storage order
+ for (std::ptrdiff_t r = 0; r < Eigen::DenseBase<Derived>::RowsAtCompileTime;
+ ++r)
+ {
+ auto const R = OGSToMFront(r);
+ for (std::ptrdiff_t c = 0;
+ c < Eigen::DenseBase<Derived>::ColsAtCompileTime;
+ ++c)
+ {
+ auto const C = OGSToMFront(c);
+ n(R, C) = m(r, c);
+ }
+ }
+
+ return n;
+}
+
+/// Converts between OGSes and MFront's Kelvin vectors and matrices.
+template <typename Derived>
+typename Derived::PlainObject MFrontToOGS(Eigen::DenseBase<Derived> const& m)
+{
+ static_assert(Derived::RowsAtCompileTime != Eigen::Dynamic, "Error");
+ static_assert(Derived::ColsAtCompileTime != Eigen::Dynamic, "Error");
+
+ typename Derived::PlainObject n;
+
+ // optimal for row-major storage order
+ for (std::ptrdiff_t r = 0; r < Eigen::DenseBase<Derived>::RowsAtCompileTime;
+ ++r)
+ {
+ auto const R = MFrontToOGS(r);
+ for (std::ptrdiff_t c = 0;
+ c < Eigen::DenseBase<Derived>::ColsAtCompileTime;
+ ++c)
+ {
+ auto const C = MFrontToOGS(c);
+ n(R, C) = m(r, c);
+ }
+ }
+
+ return n;
+}
+
+} // namespace
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+const char* toString(mgis::behaviour::Behaviour::Kinematic kin)
+{
+ using K = mgis::behaviour::Behaviour::Kinematic;
+ switch (kin)
+ {
+ case K::UNDEFINEDKINEMATIC:
+ return "UNDEFINEDKINEMATIC";
+ case K::SMALLSTRAINKINEMATIC:
+ return "SMALLSTRAINKINEMATIC";
+ case K::COHESIVEZONEKINEMATIC:
+ return "COHESIVEZONEKINEMATIC";
+ case K::FINITESTRAINKINEMATIC_F_CAUCHY:
+ return "FINITESTRAINKINEMATIC_F_CAUCHY";
+ case K::FINITESTRAINKINEMATIC_ETO_PK1:
+ return "FINITESTRAINKINEMATIC_ETO_PK1";
+ }
+
+ OGS_FATAL("Unknown kinematic %d.", kin);
+}
+
+const char* toString(mgis::behaviour::Behaviour::Symmetry sym)
+{
+ using S = mgis::behaviour::Behaviour::Symmetry;
+ switch (sym)
+ {
+ case S::ISOTROPIC:
+ return "ISOTROPIC";
+ case S::ORTHOTROPIC:
+ return "ORTHOTROPIC";
+ }
+
+ OGS_FATAL("Unknown symmetry %d.", sym);
+}
+const char* btypeToString(int btype)
+{
+ using B = mgis::behaviour::Behaviour;
+ if (btype == B::GENERALBEHAVIOUR)
+ return "GENERALBEHAVIOUR";
+ if (btype == B::STANDARDSTRAINBASEDBEHAVIOUR)
+ return "STANDARDSTRAINBASEDBEHAVIOUR";
+ if (btype == B::STANDARDFINITESTRAINBEHAVIOUR)
+ return "STANDARDFINITESTRAINBEHAVIOUR";
+ if (btype == B::COHESIVEZONEMODEL)
+ return "COHESIVEZONEMODEL";
+
+ OGS_FATAL("Unknown behaviour type %d.", btype);
+}
+const char* varTypeToString(int v)
+{
+ using V = mgis::behaviour::Variable;
+ if (v == V::SCALAR)
+ return "SCALAR";
+ if (v == V::VECTOR)
+ return "VECTOR";
+ if (v == V::STENSOR)
+ return "STENSOR";
+ if (v == V::TENSOR)
+ return "TENSOR";
+
+ OGS_FATAL("Unknown variable type %d.", v);
+}
+
+template <int DisplacementDim>
+MFront<DisplacementDim>::MFront(
+ mgis::behaviour::Behaviour&& behaviour,
+ std::vector<ProcessLib::Parameter<double> const*>&& material_properties)
+ : _behaviour(std::move(behaviour)),
+ _material_properties(std::move(material_properties))
+{
+ auto const hypothesis = behaviour.hypothesis;
+
+ if (_behaviour.symmetry != mgis::behaviour::Behaviour::Symmetry::ISOTROPIC)
+ OGS_FATAL(
+ "The storage order of the stiffness matrix is not tested, yet. "
+ "Thus, we cannot be sure if we compute the behaviour of "
+ "anisotropic materials correctly. Therefore, currently only "
+ "isotropic materials are allowed.");
+
+ if (_behaviour.gradients.size() != 1)
+ OGS_FATAL(
+ "The behaviour must have exactly a single gradient as input.");
+
+ if (_behaviour.gradients[0].name != "Strain")
+ OGS_FATAL("The behaviour must be driven by strain.");
+
+ if (_behaviour.gradients[0].type != mgis::behaviour::Variable::STENSOR)
+ OGS_FATAL("Strain must be a symmetric tensor.");
+
+ if (mgis::behaviour::getVariableSize(_behaviour.gradients[0], hypothesis) !=
+ MFront<DisplacementDim>::KelvinVector::SizeAtCompileTime)
+ OGS_FATAL("Strain must have %ld components instead of %lu.",
+ MFront<DisplacementDim>::KelvinVector::SizeAtCompileTime,
+ mgis::behaviour::getVariableSize(_behaviour.gradients[0],
+ hypothesis));
+
+ if (_behaviour.thermodynamic_forces.size() != 1)
+ OGS_FATAL(
+ "The behaviour must compute exactly one thermodynamic force.");
+
+ if (_behaviour.thermodynamic_forces[0].name != "Stress")
+ OGS_FATAL("The behaviour must compute stress.");
+
+ if (_behaviour.thermodynamic_forces[0].type !=
+ mgis::behaviour::Variable::STENSOR)
+ OGS_FATAL("Stress must be a symmetric tensor.");
+
+ if (mgis::behaviour::getVariableSize(_behaviour.thermodynamic_forces[0],
+ hypothesis) !=
+ MFront<DisplacementDim>::KelvinVector::SizeAtCompileTime)
+ OGS_FATAL("Stress must have %ld components instead of %lu.",
+ MFront<DisplacementDim>::KelvinVector::SizeAtCompileTime,
+ mgis::behaviour::getVariableSize(
+ _behaviour.thermodynamic_forces[0], hypothesis));
+
+ if (!_behaviour.esvs.empty())
+ {
+ if (_behaviour.esvs[0].name != "Temperature")
+ {
+ OGS_FATAL(
+ "Only temperature is supported as external state variable.");
+ }
+
+ if (mgis::behaviour::getVariableSize(_behaviour.esvs[0], hypothesis) !=
+ 1)
+ OGS_FATAL(
+ "Temperature must be a scalar instead of having %lu "
+ "components.",
+ mgis::behaviour::getVariableSize(
+ _behaviour.thermodynamic_forces[0], hypothesis));
+ }
+}
+
+template <int DisplacementDim>
+std::unique_ptr<typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+MFront<DisplacementDim>::createMaterialStateVariables() const
+{
+ return std::make_unique<MaterialStateVariables>(_behaviour);
+}
+
+template <int DisplacementDim>
+boost::optional<std::tuple<typename MFront<DisplacementDim>::KelvinVector,
+ std::unique_ptr<typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>,
+ typename MFront<DisplacementDim>::KelvinMatrix>>
+MFront<DisplacementDim>::integrateStress(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& /*eps_prev*/,
+ KelvinVector const& eps,
+ KelvinVector const& /*sigma_prev*/,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables,
+ double const T) const
+{
+ assert(
+ dynamic_cast<MaterialStateVariables const*>(&material_state_variables));
+ auto& d =
+ static_cast<MaterialStateVariables const&>(material_state_variables)
+ ._data;
+
+ // TODO add a test of material behaviour where the value of dt matters.
+ d.dt = dt;
+ d.rdt = 1.0;
+ d.K[0] = 4.0; // if K[0] is greater than 3.5, the consistent tangent
+ // operator must be computed.
+
+ // evaluate parameters at (t, x)
+ {
+ auto out = d.s1.material_properties.begin();
+ for (auto* param : _material_properties)
+ {
+ auto const& vals = (*param)(t, x);
+ out = std::copy(vals.begin(), vals.end(), out);
+ }
+ }
+
+ if (!d.s1.external_state_variables.empty())
+ {
+ // assuming that there is only temperature
+ d.s1.external_state_variables[0] = T;
+ }
+
+ auto v = mgis::behaviour::make_view(d);
+
+ auto const eps_MFront = OGSToMFront(eps);
+ for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
+ {
+ v.s1.gradients[i] = eps_MFront[i];
+ }
+
+ auto const status = mgis::behaviour::integrate(v, _behaviour);
+ if (status != 1)
+ {
+ OGS_FATAL("Integration failed with status %i.", status);
+ }
+
+ KelvinVector sigma;
+ for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
+ {
+ sigma[i] = d.s1.thermodynamic_forces[i];
+ }
+ sigma = MFrontToOGS(sigma);
+
+ // TODO row- vs. column-major storage order. This should only matter for
+ // anisotropic materials.
+ if (d.K.size() !=
+ KelvinMatrix::RowsAtCompileTime * KelvinMatrix::ColsAtCompileTime)
+ OGS_FATAL("Stiffness matrix has wrong size.");
+
+ KelvinMatrix C = MFrontToOGS(Eigen::Map<KelvinMatrix>(d.K.data()));
+
+ // TODO avoid copying the state
+ auto state_copy = std::make_unique<MaterialStateVariables>(
+ static_cast<MaterialStateVariables const&>(material_state_variables));
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+ state_upcast(state_copy.release());
+
+ return {{std::move(sigma), std::move(state_upcast), std::move(C)}};
+}
+
+template <int DisplacementDim>
+double MFront<DisplacementDim>::computeFreeEnergyDensity(
+ double const /*t*/,
+ ProcessLib::SpatialPosition const& /*x*/,
+ double const /*dt*/,
+ KelvinVector const& /*eps*/,
+ KelvinVector const& /*sigma*/,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ /*material_state_variables*/) const
+{
+ // TODO implement
+ return std::numeric_limits<double>::quiet_NaN();
+}
+
+template class MFront<2>;
+template class MFront<3>;
+
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/MFront/MFront.h b/MaterialLib/SolidModels/MFront/MFront.h
new file mode 100644
index 00000000000..85d26a25c42
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MFront.h
@@ -0,0 +1,116 @@
+/**
+ * \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 "MaterialLib/SolidModels/MechanicsBase.h"
+
+#include <MGIS/Behaviour/Behaviour.hxx>
+#include <MGIS/Behaviour/BehaviourData.hxx>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+namespace MFront
+{
+/// Converts MGIS kinematic to a string representation.
+const char* toString(mgis::behaviour::Behaviour::Kinematic);
+
+/// Converts MGIS symmetry to a string representation.
+const char* toString(mgis::behaviour::Behaviour::Symmetry);
+
+/// Converts MGIS behaviour type to a string representation.
+const char* btypeToString(int);
+
+/// Converts MGIS variable type to a string representation.
+const char* varTypeToString(int);
+
+/// Uses a material model provided by MFront (via MFront's generic interface and
+/// the MGIS library).
+template <int DisplacementDim>
+class MFront : public MechanicsBase<DisplacementDim>
+{
+public:
+ struct MaterialStateVariables
+ : public MechanicsBase<DisplacementDim>::MaterialStateVariables
+ {
+ explicit MaterialStateVariables(mgis::behaviour::Behaviour const& b)
+ : _data{b}
+ {
+ }
+
+ MaterialStateVariables(MaterialStateVariables const&) = default;
+
+ void pushBackState() override { mgis::behaviour::update(_data); }
+
+ MaterialStateVariables& operator=(MaterialStateVariables const&) =
+ default;
+
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ operator=(typename MechanicsBase<DisplacementDim>::
+ MaterialStateVariables const& state) noexcept override
+ {
+ return operator=(static_cast<MaterialStateVariables const&>(state));
+ }
+
+ // TODO fix: this has to be mutable.
+ mutable mgis::behaviour::BehaviourData _data;
+ };
+
+ using KelvinVector =
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix =
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+
+ MFront(mgis::behaviour::Behaviour&& behaviour,
+ std::vector<ProcessLib::Parameter<double> const*>&&
+ material_properties);
+
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() const override;
+
+ boost::optional<std::tuple<KelvinVector,
+ std::unique_ptr<typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>,
+ KelvinMatrix>>
+ integrateStress(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables,
+ double const T) const override;
+
+ double computeFreeEnergyDensity(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& eps,
+ KelvinVector const& sigma,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables) const override;
+
+private:
+ mgis::behaviour::Behaviour _behaviour;
+ std::vector<ProcessLib::Parameter<double> const*> _material_properties;
+};
+
+extern template class MFront<2>;
+extern template class MFront<3>;
+
+} // namespace MFront
+} // namespace Solids
+} // namespace MaterialLib
--
GitLab