1 /*!
2  * \file   mtest/src/GenericBehaviour.cxx
3  * \brief
4  * \author Thomas Helfer
5  * \brief 07 avril 2013
6  * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights
7  * reserved.
8  * This project is publicly released under either the GNU GPL Licence
9  * or the CECILL-A licence. A copy of thoses licences are delivered
10  * with the sources of TFEL. CEA or EDF may also distribute this
11  * project under specific licensing conditions.
12  */
13 
14 #include <cmath>
15 #include <limits>
16 #include <cstring>
17 #include <algorithm>
18 
19 #include "TFEL/Raise.hxx"
20 #include "TFEL/Math/tmatrix.hxx"
21 #include "TFEL/Math/stensor.hxx"
22 #include "TFEL/Math/tensor.hxx"
23 #include "TFEL/Math/t2tost2.hxx"
24 #include "TFEL/Math/st2tost2.hxx"
25 #include "TFEL/Math/ST2toST2/ST2toST2View.hxx"
26 #include "TFEL/System/ExternalLibraryManager.hxx"
27 #include "MFront/MFrontLogStream.hxx"
28 #include "MTest/Evolution.hxx"
29 #include "MTest/CurrentState.hxx"
30 #include "MTest/BehaviourWorkSpace.hxx"
31 #include "MFront/GenericBehaviour/BehaviourData.hxx"
32 #include "MTest/GenericBehaviour.hxx"
33 
34 namespace mtest {
35 
applyRotation(real * const v,const std::vector<int> & types,const GenericBehaviour::Hypothesis h,const tfel::math::tmatrix<3u,3u,real> & r)36   static void applyRotation(real* const v,
37                             const std::vector<int>& types,
38                             const GenericBehaviour::Hypothesis h,
39                             const tfel::math::tmatrix<3u, 3u, real>& r) {
40     auto o = size_t{};
41     const auto n = tfel::material::getSpaceDimension(h);
42     for (const auto& type : types) {
43       if (type == 0) {
44         o += 1;
45       } else if (type == 1) {
46         if (n == 2u) {
47           tfel::math::stensor<2u, real> s(v + o);
48           const auto rs = change_basis(s, r);
49           tfel::fsalgo::copy<4u>::exe(rs.begin(), v + o);
50         } else if (n == 3u) {
51           tfel::math::stensor<3u, real> s(v + o);
52           const auto rs = change_basis(s, r);
53           tfel::fsalgo::copy<6u>::exe(rs.begin(), v + o);
54         }
55         o += tfel::material::getStensorSize(h);
56       } else if (type == 2) {
57         tfel::raise("applyRotation: vector are not supported yet");
58         o += tfel::material::getSpaceDimension(h);
59       } else if (type == 3) {
60         if (n == 2u) {
61           tfel::math::tensor<2u, real> t(v + o);
62           const auto rt = change_basis(t, r);
63           tfel::fsalgo::copy<5u>::exe(rt.begin(), v + o);
64         } else if (n == 3u) {
65           tfel::math::tensor<3u, real> t(v + o);
66           const auto rt = change_basis(t, r);
67           tfel::fsalgo::copy<9u>::exe(rt.begin(), v + o);
68         }
69         o += tfel::material::getTensorSize(h);
70       }
71     }
72   }  // end of applyRotation
73 
applyRotation(real * const v,const std::vector<int> & dvtypes,const std::vector<int> & thtypes,const GenericBehaviour::Hypothesis h,const tfel::math::tmatrix<3u,3u,real> & r)74   static void applyRotation(real* const v,
75                             const std::vector<int>& dvtypes,
76                             const std::vector<int>& thtypes,
77                             const GenericBehaviour::Hypothesis h,
78                             const tfel::math::tmatrix<3u, 3u, real>& r) {
79     auto o = size_t{};
80     const auto n = tfel::material::getSpaceDimension(h);
81     tfel::raise_if(dvtypes.size() != thtypes.size(),
82                    "applyRotation: the number of driving variables "
83                    "does not match the number of thermodynamic fores");
84     tfel::raise_if(dvtypes.size() != 1u, "applyRotation: unsupported case");
85     for (decltype(dvtypes.size()) i = 0; i != dvtypes.size(); ++i) {
86       if (dvtypes[i] == 1) {
87         // symmetric tensors
88         tfel::raise_if(thtypes[i] != 1,
89                        "applyRotation: "
90                        "unsupported case");
91         if (n == 2u) {
92           tfel::math::st2tost2<2u, real> k;
93           std::copy(v + o, v + o + k.size(), k.begin());
94           const auto rk = change_basis(k, r);
95           std::copy(rk.begin(), rk.end(), v + o);
96         } else if (n == 3u) {
97           tfel::math::st2tost2<3u, real> k;
98           std::copy(v + o, v + o + k.size(), k.begin());
99           const auto rk = change_basis(k, r);
100           std::copy(rk.begin(), rk.end(), v + o);
101         }
102       } else if (dvtypes[i] == 2) {
103         tfel::raise_if(dvtypes[i] != 1,
104                        "applyRotation: "
105                        "unsupported case");
106         if (n == 2u) {
107           tfel::math::t2tost2<2u, real> k;
108           std::copy(v + o, v + o + k.size(), k.begin());
109           const auto rk = change_basis(k, r);
110           std::copy(rk.begin(), rk.end(), v + o);
111         } else if (n == 3u) {
112           tfel::math::t2tost2<3u, real> k;
113           std::copy(v + o, v + o + k.size(), k.begin());
114           const auto rk = change_basis(k, r);
115           std::copy(rk.begin(), rk.end(), v + o);
116         }
117       } else {
118         tfel::raise(
119             "applyRotation: "
120             "unsupported driving variable type");
121       }
122     }
123   }  // end of applyRotation
124 
GenericBehaviour(const Hypothesis h,const std::string & l,const std::string & b)125   GenericBehaviour::GenericBehaviour(const Hypothesis h,
126                                      const std::string& l,
127                                      const std::string& b)
128       : StandardBehaviourBase(h, l, b) {
129     using tfel::system::ExternalLibraryManager;
130     auto& elm = ExternalLibraryManager::getExternalLibraryManager();
131     const auto f = b + "_" + ModellingHypothesis::toString(h);
132     this->fct = elm.getGenericBehaviourFunction(l, f);
133     // additional material properties to compute the elastic stiffness
134     auto mps = std::vector<std::string>{};
135     if (this->requiresStiffnessTensor) {
136       if (this->etype == 0u) {
137         mps.insert(mps.end(), {"YoungModulus", "PoissonRatio"});
138       } else if (this->etype == 1u) {
139         if ((h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN) ||
140             (h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS)) {
141           mps.insert(mps.end(),
142                      {"YoungModulus1", "YoungModulus2", "YoungModulus3",
143                       "PoissonRatio12", "PoissonRatio23", "PoissonRatio13"});
144         } else if ((h == ModellingHypothesis::PLANESTRESS) ||
145                    (h == ModellingHypothesis::PLANESTRAIN) ||
146                    (h == ModellingHypothesis::AXISYMMETRICAL) ||
147                    (h == ModellingHypothesis::GENERALISEDPLANESTRAIN)) {
148           mps.insert(mps.end(),
149                      {"YoungModulus1", "YoungModulus2", "YoungModulus3",
150                       "PoissonRatio12", "PoissonRatio23", "PoissonRatio13",
151                       "ShearModulus12"});
152         } else if (h == ModellingHypothesis::TRIDIMENSIONAL) {
153           mps.insert(mps.end(),
154                      {"YoungModulus1", "YoungModulus2", "YoungModulus3",
155                       "PoissonRatio12", "PoissonRatio23", "PoissonRatio13",
156                       "ShearModulus12", "ShearModulus23", "ShearModulus13"});
157         } else {
158           tfel::raise(
159               "GenericBehaviour::GenericBehaviour : "
160               "unsupported modelling hypothesis");
161         }
162       } else {
163         tfel::raise(
164             "GenericBehaviour::GenericBehaviour : "
165             "unsupported symmetry type "
166             "(neither isotropic nor orthotropic)");
167       }
168     }
169     if (this->requiresThermalExpansionCoefficientTensor) {
170       if (this->stype == 0u) {
171         mps.push_back("ThermalExpansion");
172       } else if (this->stype == 1u) {
173         mps.insert(mps.end(), {"ThermalExpansion1", "ThermalExpansion2",
174                                "ThermalExpansion3"});
175       } else {
176         tfel::raise(
177             "GenericBehaviour::GenericBehaviour : "
178             "unsupported symmetry type "
179             "(neither isotropic nor orthotropic)");
180       }
181     }
182     this->mpnames.insert(this->mpnames.begin(), mps.begin(), mps.end());
183   }  // end of GenericBehaviour::GenericBehaviour
184 
allocate(BehaviourWorkSpace & wk) const185   void GenericBehaviour::allocate(BehaviourWorkSpace& wk) const {
186     const auto ndv = this->getGradientsSize();
187     const auto nth = this->getThermodynamicForcesSize();
188     const auto nstatev = this->getInternalStateVariablesSize();
189     const auto nevs = this->getExternalStateVariablesSize();
190     wk.D.resize(nth, ndv);
191     wk.k.resize(nth, ndv);
192     wk.kt.resize(nth, ndv);
193     wk.nk.resize(nth, ndv);
194     wk.ne.resize(ndv);
195     wk.ns.resize(nth);
196     wk.mps.resize(this->getMaterialPropertiesSize());
197     wk.ivs0.resize(nstatev);
198     wk.ivs.resize(nstatev);
199     wk.nivs.resize(nstatev);
200     wk.evs0.resize(nevs);
201     wk.evs.resize(nevs);
202     if ((this->stype == 1u) || (this->btype == 1u)) {
203       wk.e0.resize(ndv);
204       wk.e1.resize(ndv);
205       wk.s0.resize(nth);
206     }
207     mtest::allocate(wk.cs, this->shared_from_this());
208   }  // end f GenericBehaviour::allocate
209 
getGradientsDefaultInitialValues(tfel::math::vector<real> & v) const210   void GenericBehaviour::getGradientsDefaultInitialValues(
211       tfel::math::vector<real>& v) const {
212     std::fill(v.begin(), v.end(), real(0));
213     if (this->btype == 2u) {
214       // finite strain behaviour
215       v[0] = v[1] = v[2] = real(1);
216     }
217   }  // end of GenericBehaviour::setGradientsDefaultInitialValue
218 
computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const219   std::pair<bool, real> GenericBehaviour::computePredictionOperator(
220       BehaviourWorkSpace& wk,
221       const CurrentState& s,
222       const StiffnessMatrixType ktype) const {
223     wk.cs = s;
224     return this->call_behaviour(wk.kt, wk.cs, wk, real(1), ktype, false);
225   }  // end of GenericBehaviour::computePredictionOperator
226 
integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const227   std::pair<bool, real> GenericBehaviour::integrate(
228       CurrentState& s,
229       BehaviourWorkSpace& wk,
230       const real dt,
231       const StiffnessMatrixType ktype) const {
232     return this->call_behaviour(wk.k, s, wk, dt, ktype, true);
233   }  // end of GenericBehaviour::integrate
234 
call_behaviour(tfel::math::matrix<real> & Kt,CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype,const bool b) const235   std::pair<bool, real> GenericBehaviour::call_behaviour(
236       tfel::math::matrix<real>& Kt,
237       CurrentState& s,
238       BehaviourWorkSpace& wk,
239       const real dt,
240       const StiffnessMatrixType ktype,
241       const bool b) const {
242     using namespace std;
243     using namespace tfel::math;
244     using tfel::math::vector;
245     auto throw_if = [](const bool c, const std::string& m) {
246       tfel::raise_if(c, "GenericBehaviour::call_behaviour: " + m);
247     };
248     auto init_ptr = [](vector<real>& t, const vector<real>& v) -> real* {
249       if (v.empty()) {
250         return nullptr;
251       }
252       t = v;
253       return &t[0];
254     };
255     throw_if(wk.mps.size() != s.mprops1.size(),
256              "temporary material properties vector was not allocated properly");
257     throw_if(
258         wk.ivs0.size() != s.iv0.size(),
259         "temporary internal state variable vector was not allocated properly");
260     throw_if(
261         wk.ivs.size() != s.iv1.size(),
262         "temporary internal state variable vector was not allocated properly");
263     throw_if(
264         wk.evs0.size() != s.esv0.size(),
265         "temporary external state variable vector was not allocated properly");
266     throw_if(
267         wk.evs.size() != s.esv0.size(),
268         "temporary external state variable vector was not allocated properly");
269     throw_if((wk.D.getNbRows() != Kt.getNbRows()) ||
270                  (wk.D.getNbCols() != Kt.getNbCols()),
271              "the memory has not been allocated correctly");
272     std::fill(wk.D.begin(), wk.D.end(), 0.);
273     mfront::gb::BehaviourData d;
274     if (this->stype == 1u) {
275       // orthotropic behaviour
276       std::copy(s.e0.begin(), s.e0.end(), wk.e0.begin());
277       std::copy(s.e1.begin(), s.e1.end(), wk.e1.begin());
278       std::copy(s.s0.begin(), s.s0.end(), wk.s0.begin());
279       if (this->btype == 1u) {
280         // small strain behaviour
281         for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) {
282           wk.e0(i) -= s.e_th0(i);
283           wk.e1(i) -= s.e_th1(i);
284         }
285       }
286       d.s0.gradients = &(wk.e0[0]);
287       d.s1.gradients = &(wk.e1[0]);
288       d.s0.thermodynamic_forces = &(wk.s0[0]);
289       d.s1.thermodynamic_forces = &s.s1[0];
290       applyRotation(d.s0.gradients, this->dvtypes, this->getHypothesis(), s.r);
291       applyRotation(d.s1.gradients, this->dvtypes, this->getHypothesis(), s.r);
292       applyRotation(d.s0.thermodynamic_forces, this->thtypes,
293                     this->getHypothesis(), s.r);
294     } else {
295       if (this->btype == 1u) {
296         // small strain behaviour
297         std::copy(s.e0.begin(), s.e0.end(), wk.e0.begin());
298         std::copy(s.e1.begin(), s.e1.end(), wk.e1.begin());
299         for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) {
300           wk.e0(i) -= s.e_th0(i);
301           wk.e1(i) -= s.e_th1(i);
302         }
303         d.s0.gradients = &(wk.e0[0]);
304         d.s1.gradients = &(wk.e1[0]);
305       } else {
306         d.s0.gradients = &(s.e0[0]);
307         d.s1.gradients = &(s.e1[0]);
308       }
309       d.s0.thermodynamic_forces = &s.s0[0];
310       d.s1.thermodynamic_forces = &s.s1[0];
311     }
312     d.s0.material_properties = init_ptr(wk.mps, s.mprops1);
313     d.s1.material_properties = d.s0.material_properties;
314     d.s0.internal_state_variables = init_ptr(wk.ivs0, s.iv0);
315     d.s1.internal_state_variables = init_ptr(wk.ivs, s.iv1);
316     d.s0.stored_energy = &s.se0;
317     d.s1.stored_energy = &s.se1;
318     d.s0.dissipated_energy = &s.de0;
319     d.s1.dissipated_energy = &s.de1;
320     d.s0.external_state_variables = init_ptr(wk.evs0, s.esv0);
321     if (!s.esv0.empty()) {
322       for (decltype(s.esv0.size()) i = 0; i != s.esv0.size(); ++i) {
323         wk.evs[i] = s.esv0[i] + s.desv[i];
324       }
325       d.s1.external_state_variables = &wk.evs[0];
326     } else {
327       d.s1.external_state_variables = nullptr;
328     }
329     d.dt = dt;
330     d.rdt = 1;
331     // choosing the type of stiffness matrix
332     StandardBehaviourBase::initializeTangentOperator(wk.D, ktype, b);
333     d.K = &(wk.D(0, 0));
334     // calling the behaviour
335     const auto r = (this->fct)(&d);
336     if (r != 1) {
337       return {false, d.rdt};
338     }
339     // tangent operator (...)
340     if (ktype != StiffnessMatrixType::NOSTIFFNESS) {
341       const auto ndv = this->getGradientsSize();
342       const auto nth = this->getThermodynamicForcesSize();
343       if (this->stype == 1u) {
344         if ((this->btype == 1u) || (this->btype == 2u)) {
345           applyRotation(&(wk.D(0, 0)), this->dvtypes, this->thtypes,
346                         this->getHypothesis(), transpose(s.r));
347         } else {
348           tfel::raise(
349               "GenericBehaviour::call_behaviour: "
350               "orthotropic behaviours are only "
351               "supported for small or finite strain behaviours");
352         }
353       }
354       for (unsigned short i = 0; i != nth; ++i) {
355         for (unsigned short j = 0; j != ndv; ++j) {
356           Kt(i, j) = wk.D(i, j);
357         }
358       }
359     }
360     if (b) {
361       if (stype == 1u) {
362         applyRotation(d.s1.thermodynamic_forces, this->thtypes,
363                       this->getHypothesis(), transpose(s.r));
364       }
365       std::copy(wk.ivs.begin(), wk.ivs.end(), s.iv1.begin());
366     }
367     return {true, d.rdt};
368   }  // end of GenericBehaviour::call_behaviour
369 
getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> & r) const370   tfel::math::tmatrix<3u, 3u, real> GenericBehaviour::getRotationMatrix(
371       const tfel::math::vector<real>&,
372       const tfel::math::tmatrix<3u, 3u, real>& r) const {
373     return r;
374   }  // end of GenericBehaviour::getRotationMatrix
375 
setOptionalMaterialPropertiesDefaultValues(EvolutionManager &,const EvolutionManager &) const376   void GenericBehaviour::setOptionalMaterialPropertiesDefaultValues(
377       EvolutionManager&, const EvolutionManager&) const {
378   }  // end of GenericBehaviour::setOptionalMaterialPropertiesDefaultValues
379 
getDefaultStiffnessMatrixType() const380   StiffnessMatrixType GenericBehaviour::getDefaultStiffnessMatrixType() const {
381     return StiffnessMatrixType::CONSISTENTTANGENTOPERATOR;
382   }  // end of GenericBehaviour::getDefaultStiffnessMatrixType
383 
384   GenericBehaviour::~GenericBehaviour() = default;
385 
386 }  // end of namespace mtest
387