12#include <MGIS/Behaviour/Behaviour.hxx>
13#include <MGIS/Behaviour/BehaviourData.hxx>
14#include <MGIS/Behaviour/Integrate.hxx>
15#include <boost/mp11.hpp>
35template <
typename Derived>
39 Eigen::Matrix<
typename Derived::Scalar, Derived::RowsAtCompileTime,
40 Derived::ColsAtCompileTime>;
42 return Matrix::NullaryExpr(
43 matrix.rows(), matrix.cols(),
44 [&m = matrix.derived()](Eigen::Index
const row, Eigen::Index
const col)
46 constexpr std::ptrdiff_t result[6] = {0, 1, 2, 3, 5, 4};
47 return m(result[row], result[col]);
58template <
int DisplacementDim,
typename Derived>
62 Eigen::Matrix<
typename Derived::Scalar, Derived::RowsAtCompileTime,
63 Derived::ColsAtCompileTime>;
65 if constexpr (DisplacementDim == 2)
67 return Matrix::NullaryExpr(
68 matrix.rows(), matrix.cols(),
69 [&m = matrix.derived()](Eigen::Index
const row,
70 Eigen::Index
const col)
72 constexpr std::ptrdiff_t result[5] = {0, 3, 4, 1, 2};
73 return m(result[row], result[col]);
76 if constexpr (DisplacementDim == 3)
78 return Matrix::NullaryExpr(
79 matrix.rows(), matrix.cols(),
80 [&m = matrix.derived()](Eigen::Index
const row,
81 Eigen::Index
const col)
83 constexpr std::ptrdiff_t result[9] = {0, 4, 8, 1, 3,
85 return m(result[row], result[col]);
99template <
int DisplacementDim>
101 std::vector<double>
const& mfront_data,
104 mgis::behaviour::Behaviour
const& behaviour);
106extern template OGSMFrontTangentOperatorData tangentOperatorDataMFrontToOGS<2>(
107 std::vector<double>
const& mfront_data,
109 mgis::behaviour::Behaviour
const& behaviour);
110extern template OGSMFrontTangentOperatorData tangentOperatorDataMFrontToOGS<3>(
111 std::vector<double>
const& mfront_data,
113 mgis::behaviour::Behaviour
const& behaviour);
116template <
int DisplacementDim>
121 int const equivalent_plastic_strain_offset,
122 mgis::behaviour::Behaviour
const& b)
123 : equivalent_plastic_strain_offset_(equivalent_plastic_strain_offset),
140 if (equivalent_plastic_strain_offset_ >= 0)
142 return _behaviour_data.s1
143 .internal_state_variables[
static_cast<mgis::size_type
>(
144 equivalent_plastic_strain_offset_)];
154template <
int DisplacementDim, mgis::behaviour::Variable::Type MFrontType>
157template <
int DisplacementDim>
163template <
int DisplacementDim>
169template <
int DisplacementDim>
175template <
int DisplacementDim, mgis::behaviour::Variable::Type MFrontType>
180template <
int DisplacementDim>
188 template <
typename Grad>
191 auto constexpr num_comp = Grad::template size<DisplacementDim>();
193 if constexpr (Grad::type == mgis::behaviour::Variable::Type::SCALAR)
195 *target = variable_array.*Grad::mpl_var;
197 else if constexpr (Grad::type ==
198 mgis::behaviour::Variable::Type::STENSOR)
201 auto const& grad_ogs =
202 std::get<MPLType>(variable_array.*Grad::mpl_var);
204 auto const grad_mfront =
207 std::copy_n(grad_mfront.data(), num_comp, target);
209 else if constexpr (Grad::type ==
210 mgis::behaviour::Variable::Type::TENSOR)
213 auto const& grad_ogs =
214 std::get<MPLType>(variable_array.*Grad::mpl_var);
218 OGS_FATAL(
"Rotations of tensors are not implemented.");
221 Eigen::Map<Eigen::Vector<double, num_comp>>{target} =
222 ogsTensorToMFrontTensor<DisplacementDim>(grad_ogs);
226 OGS_FATAL(
"Unsupported gradient type {}.",
237template <
int DisplacementDim,
240 typename ExtStateVars>
243 static_assert(boost::mp11::mp_is_set<Gradients>::value);
244 static_assert(boost::mp11::mp_is_set<TDynForces>::value);
245 static_assert(boost::mp11::mp_is_set<ExtStateVars>::value);
248 boost::mp11::mp_append<Gradients, ExtStateVars>;
250 static_assert(boost::mp11::mp_is_set<GradientsAndExtStateVars>::value);
261 mgis::behaviour::Behaviour&& behaviour,
265 state_variables_initial_properties,
266 std::optional<ParameterLib::CoordinateSystem>
const&
267 local_coordinate_system)
268 : _behaviour(std::move(behaviour)),
269 equivalent_plastic_strain_offset_(
271 _material_properties(std::move(material_properties)),
272 _state_variables_initial_properties(
273 std::move(state_variables_initial_properties)),
274 _local_coordinate_system(local_coordinate_system
275 ? &local_coordinate_system.value()
279 auto check_gradient = [&grads = _behaviour.gradients,
280 hyp = _behaviour.hypothesis,
281 i = 0]<
typename Grad>(Grad)
mutable
284 if (grads[i].name != Grad::name)
287 "OGS expects the {}th gradient to be {} but MFront "
289 i, Grad::name, grads[i].name);
292 if (grads[i].type != Grad::type)
295 "The behaviour's {}th driver ({}) must be of type {}.",
299 if (mgis::behaviour::getVariableSize(grads[i], hyp) !=
300 Grad::template size<DisplacementDim>())
303 "The behaviour's {}th driver's ({}) size in OGS is {} "
306 Grad::template size<DisplacementDim>(),
307 mgis::behaviour::getVariableSize(grads[i], hyp));
313 if (_behaviour.gradients.size() !=
314 boost::mp11::mp_size<Gradients>::value)
316 "The behaviour must have exactly {} gradients as input.",
317 boost::mp11::mp_size<Gradients>::value);
319 boost::mp11::mp_for_each<Gradients>(check_gradient);
323 auto check_tdyn_force = [&tdfs = _behaviour.thermodynamic_forces,
324 hyp = _behaviour.hypothesis,
325 i = 0]<
typename TDF>(TDF)
mutable
327 if (tdfs[i].name != TDF::name)
330 "OGS expects the {}th thermodynamic force to be {} but "
331 "MFront provides {}.",
332 i, TDF::name, tdfs[i].name);
335 if (tdfs[i].type != TDF::type)
338 "The behaviour's {}th thermodynamic force ({}) must be "
343 if (mgis::behaviour::getVariableSize(tdfs[i], hyp) !=
344 TDF::template size<DisplacementDim>())
347 "The behaviour's {}th thermodynamic force ({}) must "
348 "have size {} instead of {}.",
349 i, tdfs[i].name, TDF::template size<DisplacementDim>(),
350 mgis::behaviour::getVariableSize(tdfs[i], hyp));
356 if (_behaviour.thermodynamic_forces.size() !=
357 boost::mp11::mp_size<TDynForces>::value)
359 "The behaviour must compute exactly {} thermodynamic "
361 boost::mp11::mp_size<TDynForces>::value);
363 boost::mp11::mp_for_each<TDynForces>(check_tdyn_force);
366 auto const hypothesis = _behaviour.hypothesis;
369 std::is_same_v<ExtStateVars, boost::mp11::mp_list<Temperature>>,
370 "Temperature is the only allowed external state variable.");
372 if (!_behaviour.esvs.empty())
374 if (_behaviour.esvs[0].name !=
"Temperature")
377 "Only temperature is supported as external state "
381 if (mgis::behaviour::getVariableSize(_behaviour.esvs[0],
384 "Temperature must be a scalar instead of having {:d} "
386 mgis::behaviour::getVariableSize(
387 _behaviour.thermodynamic_forces[0], hypothesis));
390 if (_behaviour.mps.size() != _material_properties.size())
392 ERR(
"There are {:d} material properties in the loaded behaviour:",
393 _behaviour.mps.size());
394 for (auto const& mp : _behaviour.mps)
396 ERR(
"\t{:s}", mp.name);
398 OGS_FATAL(
"But the number of passed material properties is {:d}.",
399 _material_properties.size());
404 typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
407 return std::make_unique<MaterialStateVariablesMFront<DisplacementDim>>(
408 equivalent_plastic_strain_offset_, _behaviour);
415 material_state_variables)
const
418 &material_state_variables));
422 material_state_variables);
424 auto const& ivs = getInternalVariables();
426 for (
auto& [name, parameter] : _state_variables_initial_properties)
432 {
return iv.name == name; },
434 {
OGS_FATAL(
"Internal variable `{:s}' not found.", name); });
437 std::vector<double> values = (*parameter)(t, x);
440 auto const values_span = iv.reference(state);
441 assert(values.size() == values_span.size());
442 std::copy_n(begin(values), values_span.size(), values_span.begin());
445 auto const& s1 = state._behaviour_data.s1.internal_state_variables;
446 auto& s0 = state._behaviour_data.s0.internal_state_variables;
447 std::copy(begin(s1), end(s1), begin(s0));
452 DisplacementDim>::MaterialStateVariables>,
461 material_state_variables)
const
467 &material_state_variables));
469 auto state = std::make_unique<
472 material_state_variables));
475 behaviour_data.dt = dt;
476 behaviour_data.rdt = 1.0;
477 behaviour_data.K[0] =
480 if (behaviour_data.s0.b.btype ==
481 mgis::behaviour::Behaviour::STANDARDFINITESTRAINBEHAVIOUR)
483 behaviour_data.K[1] = 1.0;
485 behaviour_data.K[2] =
494 auto out = behaviour_data.s0.material_properties.begin();
495 for (
auto* param : _material_properties)
497 auto const& vals = (*param)(t - dt, x);
498 out = std::copy(vals.begin(), vals.end(), out);
500 assert(out == behaviour_data.s0.material_properties.end());
504 auto out = behaviour_data.s1.material_properties.begin();
505 for (
auto* param : _material_properties)
507 auto const& vals = (*param)(t, x);
508 out = std::copy(vals.begin(), vals.end(), out);
510 assert(out == behaviour_data.s1.material_properties.end());
516 if (!behaviour_data.s0.external_state_variables.empty())
521 behaviour_data.s0.external_state_variables[0] =
525 if (!behaviour_data.s1.external_state_variables.empty())
530 behaviour_data.s1.external_state_variables[0] =
535 std::optional<KelvinMatrixType<DisplacementDim>> Q;
536 if (_local_coordinate_system)
538 Q = fourthOrderRotationMatrix(
539 _local_coordinate_system->transformation<DisplacementDim>(x));
542 boost::mp11::mp_for_each<Gradients>(
544 variable_array_prev, Q, behaviour_data.s0.gradients.data()});
545 boost::mp11::mp_for_each<Gradients>(
547 variable_array, Q, behaviour_data.s1.gradients.data()});
552 boost::mp11::mp_for_each<TDynForces>(
554 variable_array_prev, Q,
555 behaviour_data.s0.thermodynamic_forces.data()});
556 boost::mp11::mp_for_each<TDynForces>(
558 variable_array_prev, Q,
559 behaviour_data.s1.thermodynamic_forces.data()});
561 auto v = mgis::behaviour::make_view(behaviour_data);
562 auto const status = mgis::behaviour::integrate(v, _behaviour);
566 "MFront: integration failed with status " +
567 std::to_string(status) +
".");
571 tdyn_forces_data.
data.resize(
573 behaviour_data.s1.thermodynamic_forces.size());
575 boost::mp11::mp_for_each<TDynForces>(
576 [&out_data = tdyn_forces_data,
577 &in_data = behaviour_data.s1.thermodynamic_forces,
578 &Q]<
typename TDF>(TDF tdf)
580 OGSMFrontThermodynamicForcesView<DisplacementDim, TDynForces>
583 if constexpr (TDF::type ==
584 mgis::behaviour::Variable::Type::STENSOR)
588 view.block(tdf, out_data) =
589 *Q * eigenSwap45View(view.block(tdf, in_data));
593 view.block(tdf, out_data) =
594 eigenSwap45View(view.block(tdf, in_data));
597 else if constexpr (TDF::type ==
598 mgis::behaviour::Variable::Type::SCALAR)
600 view.block(tdf, out_data) = view.block(tdf, in_data);
608 return std::make_optional(
611 DisplacementDim>::MaterialStateVariables>,
613 std::move(tdyn_forces_data),
615 tangentOperatorDataMFrontToOGS<DisplacementDim>(
616 behaviour_data.K, Q, _behaviour)));
621 std::vector<InternalVariable> internal_variables;
623 for (
auto const& iv : _behaviour.isvs)
625 auto const name = iv.name;
626 auto const offset = mgis::behaviour::getVariableOffset(
627 _behaviour.isvs, name, _behaviour.hypothesis);
629 mgis::behaviour::getVariableSize(iv, _behaviour.hypothesis);
637 name,
static_cast<int>(size),
640 DisplacementDim>::MaterialStateVariables
const& state,
641 std::vector<double>& cache) -> std::vector<double>
const&
644 DisplacementDim
> const*>(&state) !=
nullptr);
645 auto const& internal_state_variables =
647 DisplacementDim
> const&>(state)
648 ._behaviour_data.s1.internal_state_variables;
651 std::copy_n(internal_state_variables.data() + offset,
657 DisplacementDim>::MaterialStateVariables& state)
661 DisplacementDim
> const*>(&state) !=
nullptr);
662 auto& internal_state_variables =
666 ._behaviour_data.s1.internal_state_variables;
668 return {internal_state_variables.data() + offset, size};
670 internal_variables.push_back(new_variable);
673 return internal_variables;
684 _behaviour.to_blocks};
700 "MFront::getBulkModulus() requires the tangent stiffness C "
702 "argument to be valid.");
706 DisplacementDim)>::identity2;
707 return 1. / 9. * identity2.transpose() * *C * identity2;
720 return std::numeric_limits<double>::quiet_NaN();
727 std::map<std::string, ParameterLib::Parameter<double>
const*>
std::optional< std::tuple< OGSMFrontThermodynamicForcesData, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, OGSMFrontTangentOperatorData > > integrateStress(MaterialPropertyLib::VariableArray const &variable_array_prev, MaterialPropertyLib::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x, double const dt, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const
int const equivalent_plastic_strain_offset_
mgis::behaviour::Behaviour _behaviour
OGSMFrontTangentOperatorBlocksView< DisplacementDim, ForcesGradsCombinations > createTangentOperatorBlocksView() const
std::vector< InternalVariable > getInternalVariables() const
std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables > createMaterialStateVariables() const
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
double computeFreeEnergyDensity(double const, ParameterLib::SpatialPosition const &, double const, KelvinVector const &, KelvinVector const &, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &) const
OGSMFrontThermodynamicForcesView< DisplacementDim, TDynForces > createThermodynamicForcesView() const
std::vector< ParameterLib::Parameter< double > const * > _material_properties
void initializeInternalStateVariables(double const t, ParameterLib::SpatialPosition const &x, typename MechanicsBase< DisplacementDim >::MaterialStateVariables &material_state_variables) const
typename MechanicsBase< DisplacementDim >::InternalVariable InternalVariable
std::map< std::string, ParameterLib::Parameter< double > const * > _state_variables_initial_properties
boost::mp11::mp_append< Gradients, ExtStateVars > GradientsAndExtStateVars
MFrontGeneric(mgis::behaviour::Behaviour &&behaviour, std::vector< ParameterLib::Parameter< double > const * > &&material_properties, std::map< std::string, ParameterLib::Parameter< double > const * > &&state_variables_initial_properties, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system)
double getBulkModulus(double const, ParameterLib::SpatialPosition const &, KelvinMatrix const *const C) const
ParameterLib::CoordinateSystem const *const _local_coordinate_system
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
ranges::range_reference_t< Range > findElementOrError(Range &range, std::predicate< ranges::range_reference_t< Range > > auto &&predicate, std::invocable auto error_callback)
typename MapToMPLType< DisplacementDim, MFrontType >::type MapToMPLType_t
OGSMFrontTangentOperatorData tangentOperatorDataMFrontToOGS(std::vector< double > const &mfront_data, std::optional< MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > const &Q, mgis::behaviour::Behaviour const &behaviour)
constexpr auto ogsTensorToMFrontTensor(Eigen::MatrixBase< Derived > const &matrix)
constexpr auto eigenSwap45View(Eigen::MatrixBase< Derived > const &matrix)
const char * varTypeToString(int v)
int getEquivalentPlasticStrainOffset(mgis::behaviour::Behaviour const &b)
Eigen::Matrix< double, tensorSize(Dim), 1 > Tensor
Eigen::Matrix< double, symmetric_tensor_size< GlobalDim >, 1 > SymmetricTensor
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
int const equivalent_plastic_strain_offset_
MaterialStateVariablesMFront(int const equivalent_plastic_strain_offset, mgis::behaviour::Behaviour const &b)
MaterialStateVariablesMFront(MaterialStateVariablesMFront< DisplacementDim > &&)=delete
void pushBackState() override
mgis::behaviour::BehaviourData _behaviour_data
MaterialStateVariablesMFront(MaterialStateVariablesMFront< DisplacementDim > const &)=default
double getEquivalentPlasticStrain() const override
Used for disambiguation with MFront's thermodynamic forces data.
Used for disambiguation with MFront's tangent operator blocks data.
std::vector< double > data
MaterialPropertyLib::Tensor< DisplacementDim > type
MaterialPropertyLib::SymmetricTensor< DisplacementDim > type
A struct mapping a MGIS variable type to its corresponding OGS type.
std::optional< MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > const & Q
MaterialPropertyLib::VariableArray const & variable_array
Helper type for providing access to internal variables.
A local coordinate system used for tensor transformations.