OGS 6.1.0-1721-g6382411ad
MFront.cpp
Go to the documentation of this file.
1 
10 #include "MFront.h"
11 
12 #include <MGIS/Behaviour/Integrate.hxx>
13 
14 namespace
15 {
17 constexpr std::ptrdiff_t OGSToMFront(std::ptrdiff_t i)
18 {
19  // MFront: 11 22 33 12 13 23
20  // OGS: 11 22 33 12 23 13
21  if (i < 4)
22  return i;
23  if (i == 4)
24  return 5;
25  return 4;
26 }
28 constexpr std::ptrdiff_t MFrontToOGS(std::ptrdiff_t i)
29 {
30  // MFront: 11 22 33 12 13 23
31  // OGS: 11 22 33 12 23 13
32  return OGSToMFront(i); // Same algorithm: indices 4 and 5 swapped.
33 }
34 
36 template <typename Derived>
37 typename Derived::PlainObject OGSToMFront(Eigen::DenseBase<Derived> const& m)
38 {
39  static_assert(Derived::RowsAtCompileTime != Eigen::Dynamic, "Error");
40  static_assert(Derived::ColsAtCompileTime != Eigen::Dynamic, "Error");
41 
42  typename Derived::PlainObject n;
43 
44  // optimal for row-major storage order
45  for (std::ptrdiff_t r = 0; r < Eigen::DenseBase<Derived>::RowsAtCompileTime;
46  ++r)
47  {
48  auto const R = OGSToMFront(r);
49  for (std::ptrdiff_t c = 0;
50  c < Eigen::DenseBase<Derived>::ColsAtCompileTime;
51  ++c)
52  {
53  auto const C = OGSToMFront(c);
54  n(R, C) = m(r, c);
55  }
56  }
57 
58  return n;
59 }
60 
62 template <typename Derived>
63 typename Derived::PlainObject MFrontToOGS(Eigen::DenseBase<Derived> const& m)
64 {
65  static_assert(Derived::RowsAtCompileTime != Eigen::Dynamic, "Error");
66  static_assert(Derived::ColsAtCompileTime != Eigen::Dynamic, "Error");
67 
68  typename Derived::PlainObject n;
69 
70  // optimal for row-major storage order
71  for (std::ptrdiff_t r = 0; r < Eigen::DenseBase<Derived>::RowsAtCompileTime;
72  ++r)
73  {
74  auto const R = MFrontToOGS(r);
75  for (std::ptrdiff_t c = 0;
76  c < Eigen::DenseBase<Derived>::ColsAtCompileTime;
77  ++c)
78  {
79  auto const C = MFrontToOGS(c);
80  n(R, C) = m(r, c);
81  }
82  }
83 
84  return n;
85 }
86 
87 } // namespace
88 
89 namespace MaterialLib
90 {
91 namespace Solids
92 {
93 namespace MFront
94 {
95 const char* toString(mgis::behaviour::Behaviour::Kinematic kin)
96 {
97  using K = mgis::behaviour::Behaviour::Kinematic;
98  switch (kin)
99  {
100  case K::UNDEFINEDKINEMATIC:
101  return "UNDEFINEDKINEMATIC";
102  case K::SMALLSTRAINKINEMATIC:
103  return "SMALLSTRAINKINEMATIC";
104  case K::COHESIVEZONEKINEMATIC:
105  return "COHESIVEZONEKINEMATIC";
106  case K::FINITESTRAINKINEMATIC_F_CAUCHY:
107  return "FINITESTRAINKINEMATIC_F_CAUCHY";
108  case K::FINITESTRAINKINEMATIC_ETO_PK1:
109  return "FINITESTRAINKINEMATIC_ETO_PK1";
110  }
111 
112  OGS_FATAL("Unknown kinematic %d.", kin);
113 }
114 
115 const char* toString(mgis::behaviour::Behaviour::Symmetry sym)
116 {
117  using S = mgis::behaviour::Behaviour::Symmetry;
118  switch (sym)
119  {
120  case S::ISOTROPIC:
121  return "ISOTROPIC";
122  case S::ORTHOTROPIC:
123  return "ORTHOTROPIC";
124  }
125 
126  OGS_FATAL("Unknown symmetry %d.", sym);
127 }
128 const char* btypeToString(int btype)
129 {
130  using B = mgis::behaviour::Behaviour;
131  if (btype == B::GENERALBEHAVIOUR)
132  return "GENERALBEHAVIOUR";
133  if (btype == B::STANDARDSTRAINBASEDBEHAVIOUR)
134  return "STANDARDSTRAINBASEDBEHAVIOUR";
135  if (btype == B::STANDARDFINITESTRAINBEHAVIOUR)
136  return "STANDARDFINITESTRAINBEHAVIOUR";
137  if (btype == B::COHESIVEZONEMODEL)
138  return "COHESIVEZONEMODEL";
139 
140  OGS_FATAL("Unknown behaviour type %d.", btype);
141 }
142 const char* varTypeToString(int v)
143 {
144  using V = mgis::behaviour::Variable;
145  if (v == V::SCALAR)
146  return "SCALAR";
147  if (v == V::VECTOR)
148  return "VECTOR";
149  if (v == V::STENSOR)
150  return "STENSOR";
151  if (v == V::TENSOR)
152  return "TENSOR";
153 
154  OGS_FATAL("Unknown variable type %d.", v);
155 }
156 
157 template <int DisplacementDim>
159  mgis::behaviour::Behaviour&& behaviour,
160  std::vector<ProcessLib::Parameter<double> const*>&& material_properties)
161  : _behaviour(std::move(behaviour)),
162  _material_properties(std::move(material_properties))
163 {
164  auto const hypothesis = behaviour.hypothesis;
165 
166  if (_behaviour.symmetry != mgis::behaviour::Behaviour::Symmetry::ISOTROPIC)
167  OGS_FATAL(
168  "The storage order of the stiffness matrix is not tested, yet. "
169  "Thus, we cannot be sure if we compute the behaviour of "
170  "anisotropic materials correctly. Therefore, currently only "
171  "isotropic materials are allowed.");
172 
173  if (_behaviour.gradients.size() != 1)
174  OGS_FATAL(
175  "The behaviour must have exactly a single gradient as input.");
176 
177  if (_behaviour.gradients[0].name != "Strain")
178  OGS_FATAL("The behaviour must be driven by strain.");
179 
180  if (_behaviour.gradients[0].type != mgis::behaviour::Variable::STENSOR)
181  OGS_FATAL("Strain must be a symmetric tensor.");
182 
183  if (mgis::behaviour::getVariableSize(_behaviour.gradients[0], hypothesis) !=
185  OGS_FATAL("Strain must have %ld components instead of %lu.",
187  mgis::behaviour::getVariableSize(_behaviour.gradients[0],
188  hypothesis));
189 
190  if (_behaviour.thermodynamic_forces.size() != 1)
191  OGS_FATAL(
192  "The behaviour must compute exactly one thermodynamic force.");
193 
194  if (_behaviour.thermodynamic_forces[0].name != "Stress")
195  OGS_FATAL("The behaviour must compute stress.");
196 
197  if (_behaviour.thermodynamic_forces[0].type !=
198  mgis::behaviour::Variable::STENSOR)
199  OGS_FATAL("Stress must be a symmetric tensor.");
200 
201  if (mgis::behaviour::getVariableSize(_behaviour.thermodynamic_forces[0],
202  hypothesis) !=
204  OGS_FATAL("Stress must have %ld components instead of %lu.",
206  mgis::behaviour::getVariableSize(
207  _behaviour.thermodynamic_forces[0], hypothesis));
208 
209  if (!_behaviour.esvs.empty())
210  {
211  if (_behaviour.esvs[0].name != "Temperature")
212  {
213  OGS_FATAL(
214  "Only temperature is supported as external state variable.");
215  }
216 
217  if (mgis::behaviour::getVariableSize(_behaviour.esvs[0], hypothesis) !=
218  1)
219  OGS_FATAL(
220  "Temperature must be a scalar instead of having %lu "
221  "components.",
222  mgis::behaviour::getVariableSize(
223  _behaviour.thermodynamic_forces[0], hypothesis));
224  }
225 }
226 
227 template <int DisplacementDim>
228 std::unique_ptr<typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
230 {
231  return std::make_unique<MaterialStateVariables>(_behaviour);
232 }
233 
234 template <int DisplacementDim>
235 boost::optional<std::tuple<typename MFront<DisplacementDim>::KelvinVector,
236  std::unique_ptr<typename MechanicsBase<
237  DisplacementDim>::MaterialStateVariables>,
240  double const t,
242  double const dt,
243  KelvinVector const& /*eps_prev*/,
244  KelvinVector const& eps,
245  KelvinVector const& /*sigma_prev*/,
247  material_state_variables,
248  double const T) const
249 {
250  assert(
251  dynamic_cast<MaterialStateVariables const*>(&material_state_variables));
252  auto& d =
253  static_cast<MaterialStateVariables const&>(material_state_variables)
254  ._data;
255 
256  // TODO add a test of material behaviour where the value of dt matters.
257  d.dt = dt;
258  d.rdt = 1.0;
259  d.K[0] = 4.0; // if K[0] is greater than 3.5, the consistent tangent
260  // operator must be computed.
261 
262  // evaluate parameters at (t, x)
263  {
264  auto out = d.s1.material_properties.begin();
265  for (auto* param : _material_properties)
266  {
267  auto const& vals = (*param)(t, x);
268  out = std::copy(vals.begin(), vals.end(), out);
269  }
270  }
271 
272  if (!d.s1.external_state_variables.empty())
273  {
274  // assuming that there is only temperature
275  d.s1.external_state_variables[0] = T;
276  }
277 
278  auto v = mgis::behaviour::make_view(d);
279 
280  auto const eps_MFront = OGSToMFront(eps);
281  for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
282  {
283  v.s1.gradients[i] = eps_MFront[i];
284  }
285 
286  auto const status = mgis::behaviour::integrate(v, _behaviour);
287  if (status != 1)
288  {
289  OGS_FATAL("Integration failed with status %i.", status);
290  }
291 
292  KelvinVector sigma;
293  for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
294  {
295  sigma[i] = d.s1.thermodynamic_forces[i];
296  }
297  sigma = MFrontToOGS(sigma);
298 
299  // TODO row- vs. column-major storage order. This should only matter for
300  // anisotropic materials.
301  if (d.K.size() !=
302  KelvinMatrix::RowsAtCompileTime * KelvinMatrix::ColsAtCompileTime)
303  OGS_FATAL("Stiffness matrix has wrong size.");
304 
305  KelvinMatrix C = MFrontToOGS(Eigen::Map<KelvinMatrix>(d.K.data()));
306 
307  // TODO avoid copying the state
308  auto state_copy = std::make_unique<MaterialStateVariables>(
309  static_cast<MaterialStateVariables const&>(material_state_variables));
310  std::unique_ptr<
312  state_upcast(state_copy.release());
313 
314  return {std::make_tuple(std::move(sigma), std::move(state_upcast), std::move(C))};
315 }
316 
317 template <int DisplacementDim>
319  double const /*t*/,
320  ProcessLib::SpatialPosition const& /*x*/,
321  double const /*dt*/,
322  KelvinVector const& /*eps*/,
323  KelvinVector const& /*sigma*/,
325  /*material_state_variables*/) const
326 {
327  // TODO implement
328  return std::numeric_limits<double>::quiet_NaN();
329 }
330 
331 template class MFront<2>;
332 template class MFront<3>;
333 
334 } // namespace MFront
335 } // namespace Solids
336 } // namespace MaterialLib
static const double r
mgis::behaviour::Behaviour _behaviour
Definition: MFront.h:107
const char * varTypeToString(int v)
Converts MGIS variable type to a string representation.
Definition: MFront.cpp:142
std::vector< ProcessLib::Parameter< double > const * > _material_properties
Definition: MFront.h:108
Derived::PlainObject OGSToMFront(Eigen::DenseBase< Derived > const &m)
Converts between OGSes and MFront&#39;s Kelvin vectors and matrices.
Definition: MFront.cpp:37
std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables > createMaterialStateVariables() const override
Definition: MFront.cpp:229
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
Definition: MFront.cpp:318
const char * btypeToString(int btype)
Converts MGIS behaviour type to a string representation.
Definition: MFront.cpp:128
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition: MechanicsBase.h:74
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition: MFront.h:72
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
const char * toString(mgis::behaviour::Behaviour::Kinematic kin)
Converts MGIS kinematic to a string representation.
Definition: MFront.cpp:95
Derived::PlainObject MFrontToOGS(Eigen::DenseBase< Derived > const &m)
Converts between OGSes and MFront&#39;s Kelvin vectors and matrices.
Definition: MFront.cpp:63
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
Definition: MechanicsBase.h:72
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
Definition: MFront.cpp:239