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