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 <ostream> 18 #include <algorithm> 19 20 #include "TFEL/Raise.hxx" 21 #include "TFEL/Math/tmatrix.hxx" 22 #include "TFEL/Math/stensor.hxx" 23 #include "TFEL/Math/tensor.hxx" 24 #include "TFEL/Math/t2tost2.hxx" 25 #include "TFEL/Math/st2tost2.hxx" 26 #include "TFEL/Math/Tensor/TensorView.hxx" 27 #include "TFEL/Math/Stensor/StensorView.hxx" 28 #include "TFEL/Math/ST2toST2/ST2toST2View.hxx" 29 #include "TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx" 30 #include "TFEL/System/ExternalLibraryManager.hxx" 31 #include "MFront/MFrontLogStream.hxx" 32 #include "MTest/Evolution.hxx" 33 #include "MTest/CurrentState.hxx" 34 #include "MTest/BehaviourWorkSpace.hxx" 35 #include "MFront/GenericBehaviour/BehaviourData.hxx" 36 #include "MTest/GenericBehaviour.hxx" 37 38 namespace mtest { 39 40 template <unsigned short N> convertToSecondPiolaKirchhoffStress(real * const Sv,const real * const sv,const real * const Fv)41 void convertToSecondPiolaKirchhoffStress(real* const Sv, 42 const real* const sv, 43 const real* const Fv) { 44 tfel::math::StensorView<N, real> S(Sv); 45 tfel::math::tensor<N, real> F; 46 tfel::math::stensor<N, real> s; 47 std::copy(Fv, Fv + F.size(), F.begin()); 48 std::copy(sv, sv + s.size(), s.begin()); 49 S = tfel::math::convertCauchyStressToSecondPiolaKirchhoffStress(s, F); 50 } // end of convertToSecondPiolaKirchhoffStress 51 52 template <unsigned short N> convertFromSecondPiolaKirchhoffStress(real * const sv,const real * const Sv,const real * const Fv)53 void convertFromSecondPiolaKirchhoffStress(real* const sv, 54 const real* const Sv, 55 const real* const Fv) { 56 tfel::math::StensorView<N, real> s(sv); 57 tfel::math::tensor<N, real> F; 58 tfel::math::stensor<N, real> S; 59 std::copy(Fv, Fv + F.size(), F.begin()); 60 std::copy(Sv, Sv + S.size(), S.begin()); 61 s = tfel::math::convertSecondPiolaKirchhoffStressToCauchyStress(S, F); 62 } // end of convertFromSecondPiolaKirchhoffStress 63 64 template <unsigned short N> convertToFirstPiolaKirchhoffStress(real * const pkv,const real * const sv,const real * const Fv)65 void convertToFirstPiolaKirchhoffStress(real* const pkv, 66 const real* const sv, 67 const real* const Fv) { 68 tfel::math::TensorView<N, real> pk(pkv); 69 tfel::math::tensor<N, real> F; 70 tfel::math::stensor<N, real> s; 71 std::copy(Fv, Fv + F.size(), F.begin()); 72 std::copy(sv, sv + s.size(), s.begin()); 73 pk = tfel::math::convertCauchyStressToFirstPiolaKirchhoffStress(s, F); 74 } // end of convertToFirstPiolaKirchhoffStress 75 76 template <unsigned short N> convertFromFirstPiolaKirchhoffStress(real * const sv,const real * const pkv,const real * const Fv)77 void convertFromFirstPiolaKirchhoffStress(real* const sv, 78 const real* const pkv, 79 const real* const Fv) { 80 tfel::math::StensorView<N, real> s(sv); 81 tfel::math::tensor<N, real> F; 82 tfel::math::tensor<N, real> pk; 83 std::copy(Fv, Fv + F.size(), F.begin()); 84 std::copy(pkv, pkv + pk.size(), pk.begin()); 85 s = tfel::math::convertFirstPiolaKirchhoffStressToCauchyStress(pk, F); 86 } // end of convertFromFirstPiolaKirchhoffStress 87 88 template <unsigned short N> convertFromDSDEGL(real * const K,const real * const F0,const real * const F1,const real * const s)89 void convertFromDSDEGL(real* const K, 90 const real* const F0, 91 const real* const F1, 92 const real* const s) { 93 using FSTOBase = tfel::material::FiniteStrainBehaviourTangentOperatorBase; 94 tfel::math::st2tost2<N, real> Cse; 95 std::copy(K, K + Cse.size(), Cse.begin()); 96 const auto Kv = 97 tfel::material::convert<FSTOBase::DSIG_DF, FSTOBase::DS_DEGL>( 98 Cse, tfel::math::tensor<N, real>(F0), 99 tfel::math::tensor<N, real>(F1), tfel::math::stensor<N, real>(s)); 100 std::copy(Kv.begin(), Kv.end(), K); 101 } // end of convertFromDSDEGL 102 103 template <unsigned short N> convertFromDPK1DF(real * const K,const real * const F0,const real * const F1,const real * const s)104 void convertFromDPK1DF(real* const K, 105 const real* const F0, 106 const real* const F1, 107 const real* const s) { 108 using FSTOBase = tfel::material::FiniteStrainBehaviourTangentOperatorBase; 109 tfel::math::t2tot2<N, real> dP; 110 std::copy(K, K + dP.size(), dP.begin()); 111 const auto Kv = 112 tfel::material::convert<FSTOBase::DSIG_DF, FSTOBase::DPK1_DF>( 113 dP, tfel::math::tensor<N, real>(F0), 114 tfel::math::tensor<N, real>(F1), tfel::math::stensor<N, real>(s)); 115 std::copy(Kv.begin(), Kv.end(), K); 116 } // end of convertFromDPK1DF 117 GenericBehaviour(const Hypothesis h,const std::string & l,const std::string & b)118 GenericBehaviour::GenericBehaviour(const Hypothesis h, 119 const std::string& l, 120 const std::string& b) 121 : StandardBehaviourBase(h, l, b) { 122 using tfel::system::ExternalLibraryManager; 123 auto& elm = ExternalLibraryManager::getExternalLibraryManager(); 124 const auto f = b + "_" + ModellingHypothesis::toString(h); 125 this->fct = elm.getGenericBehaviourFunction(l, f); 126 if (this->stype == 1u) { 127 // load the rotation functions 128 this->rg_fct = elm.getGenericBehaviourRotateGradientsFunction( 129 l, f + "_rotateGradients"); 130 if (this->btype == 2u) { 131 // finite strain behaviour 132 this->rtf_fct = elm.getGenericBehaviourRotateThermodynamicForcesFunction( 133 l, f + "_rotateThermodynamicForces_CauchyStress"); 134 this->rto_fct = 135 elm.getGenericBehaviourRotateTangentOperatorBlocksFunction( 136 l, f + "_rotateTangentOperatorBlocks_dsig_dF"); 137 } else { 138 this->rtf_fct = elm.getGenericBehaviourRotateThermodynamicForcesFunction( 139 l, f + "_rotateThermodynamicForces"); 140 this->rto_fct = 141 elm.getGenericBehaviourRotateTangentOperatorBlocksFunction( 142 l, f + "_rotateTangentOperatorBlocks"); 143 } 144 } 145 // additional material properties to compute the elastic stiffness 146 auto mps = std::vector<std::string>{}; 147 if (this->requiresStiffnessTensor) { 148 if (this->etype == 0u) { 149 mps.insert(mps.end(), {"YoungModulus", "PoissonRatio"}); 150 } else if (this->etype == 1u) { 151 if ((h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN) || 152 (h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS)) { 153 mps.insert(mps.end(), 154 {"YoungModulus1", "YoungModulus2", "YoungModulus3", 155 "PoissonRatio12", "PoissonRatio23", "PoissonRatio13"}); 156 } else if ((h == ModellingHypothesis::PLANESTRESS) || 157 (h == ModellingHypothesis::PLANESTRAIN) || 158 (h == ModellingHypothesis::AXISYMMETRICAL) || 159 (h == ModellingHypothesis::GENERALISEDPLANESTRAIN)) { 160 mps.insert(mps.end(), 161 {"YoungModulus1", "YoungModulus2", "YoungModulus3", 162 "PoissonRatio12", "PoissonRatio23", "PoissonRatio13", 163 "ShearModulus12"}); 164 } else if (h == ModellingHypothesis::TRIDIMENSIONAL) { 165 mps.insert(mps.end(), 166 {"YoungModulus1", "YoungModulus2", "YoungModulus3", 167 "PoissonRatio12", "PoissonRatio23", "PoissonRatio13", 168 "ShearModulus12", "ShearModulus23", "ShearModulus13"}); 169 } else { 170 tfel::raise( 171 "GenericBehaviour::GenericBehaviour : " 172 "unsupported modelling hypothesis"); 173 } 174 } else { 175 tfel::raise( 176 "GenericBehaviour::GenericBehaviour : " 177 "unsupported symmetry type " 178 "(neither isotropic nor orthotropic)"); 179 } 180 } 181 if (this->requiresThermalExpansionCoefficientTensor) { 182 if (this->stype == 0u) { 183 mps.push_back("ThermalExpansion"); 184 } else if (this->stype == 1u) { 185 mps.insert(mps.end(), {"ThermalExpansion1", "ThermalExpansion2", 186 "ThermalExpansion3"}); 187 } else { 188 tfel::raise( 189 "GenericBehaviour::GenericBehaviour : " 190 "unsupported symmetry type " 191 "(neither isotropic nor orthotropic)"); 192 } 193 } 194 this->mpnames.insert(this->mpnames.begin(), mps.begin(), mps.end()); 195 } // end of GenericBehaviour::GenericBehaviour 196 GenericBehaviour(const Hypothesis h,const std::string & l,const std::string & b,const std::map<std::string,Parameters> & params)197 GenericBehaviour::GenericBehaviour( 198 const Hypothesis h, 199 const std::string& l, 200 const std::string& b, 201 const std::map<std::string, Parameters>& params) 202 : GenericBehaviour(h, l, b) { 203 if (params.empty()) { 204 return; 205 } 206 tfel::raise_if(this->btype != 2u, 207 "GenericBehaviour::GenericBehaviour: " 208 "no parameter expected"); 209 for (const auto& p : params) { 210 if ((p.first != "stress_measure") && (p.first != "tangent_operator")) { 211 tfel::raise( 212 "GenericBehaviour::GenericBehaviour: " 213 "unexpected parameter '" + 214 p.first + "'"); 215 } 216 tfel::raise_if(!p.second.is<std::string>(), 217 "GenericBehaviour::GenericBehaviour: " 218 "unexpected type for parameter '" + 219 p.first + "'"); 220 if (p.first == "stress_measure") { 221 const auto& ss = p.second.get<std::string>(); 222 if ((ss == "CAUCHY") || (ss == "CauchyStress")) { 223 } else if ((ss == "PK2") || (ss == "SecondPiolaKirchoffStress")) { 224 this->stress_measure = PK2; 225 } else if ((ss == "PK1") || (ss == "FirstPiolaKirchoffStress")) { 226 this->stress_measure = PK1; 227 } else { 228 tfel::raise( 229 "GenericBehaviour::GenericBehaviour: " 230 "unsupported stress measure'" + 231 ss + "'"); 232 } 233 } else { 234 const auto& to = p.second.get<std::string>(); 235 if ((to == "DSIG_DF") || (to == "DCAUCHY_DF")) { 236 } else if (to == "DS_DEGL") { 237 this->fsto = DS_DEGL; 238 } else if (to == "DPK1_DF") { 239 this->fsto = DPK1_DF; 240 } else { 241 tfel::raise( 242 "GenericBehaviour::GenericBehaviour: " 243 "unsupported tangent operator '" + 244 to + "'"); 245 } 246 } 247 } 248 } // end of GenericBehaviour::GenericBehaviour 249 allocate(BehaviourWorkSpace & wk) const250 void GenericBehaviour::allocate(BehaviourWorkSpace& wk) const { 251 const auto ndv = this->getGradientsSize(); 252 const auto nth = this->getThermodynamicForcesSize(); 253 const auto nstatev = this->getInternalStateVariablesSize(); 254 const auto nevs = this->getExternalStateVariablesSize(); 255 if (this->btype == 2u) { 256 // allocating an array large enough to store 257 // the tangent operator and its conversion to 258 // DSIG_DF 259 if (this->fsto == DSIG_DF) { 260 wk.D.resize(nth, ndv); 261 wk.k.resize(nth, ndv); 262 wk.kt.resize(nth, ndv); 263 } else if (this->fsto == DS_DEGL) { 264 wk.D.resize(nth, ndv); 265 wk.k.resize(nth, ndv); 266 wk.kt.resize(nth, ndv); 267 } else if (this->fsto == DPK1_DF) { 268 wk.D.resize(ndv, ndv); 269 wk.k.resize(ndv, ndv); 270 wk.kt.resize(ndv, ndv); 271 } else { 272 tfel::raise( 273 "GenericBehaviour::allocate: " 274 "unsupported tangent operator type"); 275 } 276 } else { 277 wk.D.resize(nth, ndv); 278 wk.k.resize(nth, ndv); 279 wk.kt.resize(nth, ndv); 280 } 281 wk.nk.resize(nth, ndv); 282 wk.ne.resize(ndv); 283 wk.ns.resize(nth); 284 wk.mps.resize(this->getMaterialPropertiesSize()); 285 wk.ivs0.resize(nstatev); 286 wk.ivs.resize(nstatev); 287 wk.nivs.resize(nstatev); 288 wk.evs0.resize(nevs); 289 wk.evs.resize(nevs); 290 if ((this->stype == 1u) || (this->btype == 1u)) { 291 wk.e0.resize(ndv); 292 wk.e1.resize(ndv); 293 wk.s0.resize(nth); 294 } 295 mtest::allocate(wk.cs, this->shared_from_this()); 296 if (this->btype == 2u) { 297 if (this->stress_measure == PK1) { 298 wk.pk0.resize(ndv); 299 wk.pk1.resize(ndv); 300 } else if (this->stress_measure == PK2) { 301 wk.S0.resize(ndv); 302 wk.S1.resize(ndv); 303 } 304 } 305 } // end f GenericBehaviour::allocate 306 getGradientsDefaultInitialValues(tfel::math::vector<real> & v) const307 void GenericBehaviour::getGradientsDefaultInitialValues( 308 tfel::math::vector<real>& v) const { 309 std::fill(v.begin(), v.end(), real(0)); 310 if (this->btype == 2u) { 311 // finite strain behaviour 312 v[0] = v[1] = v[2] = real(1); 313 } 314 } // end of GenericBehaviour::setGradientsDefaultInitialValue 315 computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const316 std::pair<bool, real> GenericBehaviour::computePredictionOperator( 317 BehaviourWorkSpace& wk, 318 const CurrentState& s, 319 const StiffnessMatrixType ktype) const { 320 wk.cs = s; 321 return this->call_behaviour(wk.kt, wk.cs, wk, real(1), ktype, false); 322 } // end of GenericBehaviour::computePredictionOperator 323 integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const324 std::pair<bool, real> GenericBehaviour::integrate( 325 CurrentState& s, 326 BehaviourWorkSpace& wk, 327 const real dt, 328 const StiffnessMatrixType ktype) const { 329 return this->call_behaviour(wk.k, s, wk, dt, ktype, true); 330 } // end of GenericBehaviour::integrate 331 call_behaviour(tfel::math::matrix<real> & Kt,CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype,const bool b) const332 std::pair<bool, real> GenericBehaviour::call_behaviour( 333 tfel::math::matrix<real>& Kt, 334 CurrentState& s, 335 BehaviourWorkSpace& wk, 336 const real dt, 337 const StiffnessMatrixType ktype, 338 const bool b) const { 339 using namespace std; 340 using namespace tfel::math; 341 using tfel::math::vector; 342 using size_type = tfel::math::matrix<real>::size_type; 343 auto throw_if = [](const bool c, const std::string& m) { 344 tfel::raise_if(c, "GenericBehaviour::call_behaviour: " + m); 345 }; 346 auto init_ptr = [](vector<real>& t, const vector<real>& v) -> real* { 347 if (v.empty()) { 348 return nullptr; 349 } 350 t = v; 351 return &t[0]; 352 }; 353 throw_if(wk.mps.size() != s.mprops1.size(), 354 "temporary material properties vector was not allocated properly"); 355 throw_if(wk.ivs0.size() != s.iv0.size(), 356 "temporary internal state variable vector was not allocated " 357 "properly"); 358 throw_if(wk.ivs.size() != s.iv1.size(), 359 "temporary internal state variable vector was not allocated " 360 "properly"); 361 throw_if(wk.evs0.size() != s.esv0.size(), 362 "temporary external state variable vector was not allocated " 363 "properly"); 364 throw_if(wk.evs.size() != s.esv0.size(), 365 "temporary external state variable vector was not allocated " 366 "properly"); 367 const auto ir = tfel::math::transpose(s.r); 368 if (this->btype == 2u) { 369 if (this->fsto == DSIG_DF) { 370 throw_if((wk.D.getNbRows() != Kt.getNbRows()) || 371 (wk.D.getNbCols() != Kt.getNbCols()), 372 "the memory has not been allocated correctly"); 373 } else if (this->fsto == DS_DEGL) { 374 const auto ndv = this->getGradientsSize(); 375 const auto nth = this->getThermodynamicForcesSize(); 376 throw_if((wk.D.getNbRows() != nth) || (wk.D.getNbCols() != ndv), 377 "the memory has not been allocated correctly"); 378 } else if (this->fsto == DPK1_DF) { 379 const auto ndv = this->getGradientsSize(); 380 throw_if((wk.D.getNbRows() != ndv) || (wk.D.getNbCols() != ndv), 381 "the memory has not been allocated correctly"); 382 } else { 383 throw_if(true, "unsupported tangent operator"); 384 } 385 } else { 386 throw_if((wk.D.getNbRows() != Kt.getNbRows()) || 387 (wk.D.getNbCols() != Kt.getNbCols()), 388 "the memory has not been allocated correctly"); 389 } 390 std::fill(wk.D.begin(), wk.D.end(), 0.); 391 mfront::gb::BehaviourData d; 392 if (this->stype == 1u) { 393 // orthotropic behaviour 394 std::copy(s.e0.begin(), s.e0.end(), wk.e0.begin()); 395 std::copy(s.e1.begin(), s.e1.end(), wk.e1.begin()); 396 std::copy(s.s0.begin(), s.s0.end(), wk.s0.begin()); 397 if (this->btype == 1u) { 398 // small strain behaviour 399 for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) { 400 wk.e0(i) -= s.e_th0(i); 401 wk.e1(i) -= s.e_th1(i); 402 } 403 } 404 this->rg_fct(&(wk.e0[0]), &(wk.e0[0]), s.r.begin()); 405 this->rg_fct(&(wk.e1[0]), &(wk.e1[0]), s.r.begin()); 406 this->rtf_fct(&(wk.s0[0]), &(wk.s0[0]), ir.begin()); 407 d.s0.gradients = &(wk.e0[0]); 408 d.s1.gradients = &(wk.e1[0]); 409 d.s0.thermodynamic_forces = &(wk.s0[0]); 410 d.s1.thermodynamic_forces = &s.s1[0]; 411 } else { 412 if (this->btype == 1u) { 413 // small strain behaviour 414 std::copy(s.e0.begin(), s.e0.end(), wk.e0.begin()); 415 std::copy(s.e1.begin(), s.e1.end(), wk.e1.begin()); 416 for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) { 417 wk.e0(i) -= s.e_th0(i); 418 wk.e1(i) -= s.e_th1(i); 419 } 420 d.s0.gradients = &(wk.e0[0]); 421 d.s1.gradients = &(wk.e1[0]); 422 } else { 423 d.s0.gradients = &(s.e0[0]); 424 d.s1.gradients = &(s.e1[0]); 425 } 426 d.s0.thermodynamic_forces = &s.s0[0]; 427 d.s1.thermodynamic_forces = &s.s1[0]; 428 } 429 d.s0.material_properties = init_ptr(wk.mps, s.mprops1); 430 d.s1.material_properties = d.s0.material_properties; 431 d.s0.internal_state_variables = init_ptr(wk.ivs0, s.iv0); 432 d.s1.internal_state_variables = init_ptr(wk.ivs, s.iv1); 433 d.s0.stored_energy = &s.se0; 434 d.s1.stored_energy = &s.se1; 435 d.s0.dissipated_energy = &s.de0; 436 d.s1.dissipated_energy = &s.de1; 437 d.s0.external_state_variables = init_ptr(wk.evs0, s.esv0); 438 if (!s.esv0.empty()) { 439 for (decltype(s.esv0.size()) i = 0; i != s.esv0.size(); ++i) { 440 wk.evs[i] = s.esv0[i] + s.desv[i]; 441 } 442 d.s1.external_state_variables = &wk.evs[0]; 443 } else { 444 d.s1.external_state_variables = nullptr; 445 } 446 d.dt = dt; 447 d.rdt = 1; 448 // type of integration to be performed 449 StandardBehaviourBase::initializeTangentOperator(wk.D, ktype, b); 450 d.K = &(wk.D(0, 0)); 451 // saving the point the thermodynamic forces in the material frame 452 auto* const s0 = d.s0.thermodynamic_forces; 453 auto* const s1 = d.s1.thermodynamic_forces; 454 if (this->btype == 2u) { 455 // selecting the stress measure 456 this->executeFiniteStrainBehaviourStressPreProcessing(wk, d); 457 // choosing the type of stiffness matrix 458 this->executeFiniteStrainBehaviourTangentOperatorPreProcessing(d, ktype); 459 } 460 // calling the behaviour 461 const auto r = (this->fct)(&d); 462 if (r != 1) { 463 return {false, d.rdt}; 464 } 465 if (mfront::getVerboseMode() >= mfront::VERBOSE_DEBUG) { 466 auto& log = mfront::getLogStream(); 467 log << "Consistent tangent operator returned by the behaviour:\n"; 468 const auto ndv = this->getGradientsSize(); 469 const auto nth = this->getThermodynamicForcesSize(); 470 for (size_type i = 0; i != ndv; ++i) { 471 for (size_type j = 0; j != nth;) { 472 log << *(&(wk.D(0, 0)) + i * nth + j); 473 if (++j != nth) { 474 log << " "; 475 } 476 } 477 log << '\n'; 478 } 479 } 480 // turn back the gradients in the global frame 481 if (stype == 1u) { 482 // here we use the fact that d.s1.gradients is a pointer to &(wk.e0[0]) 483 this->rg_fct(&(wk.e0[0]), &(wk.e0[0]), ir.begin()); 484 this->rg_fct(&(wk.e1[0]), &(wk.e1[0]), ir.begin()); 485 } 486 if (b) { 487 std::copy(wk.ivs.begin(), wk.ivs.end(), s.iv1.begin()); 488 if (this->btype == 2u) { 489 d.s0.thermodynamic_forces = s0; 490 d.s1.thermodynamic_forces = s1; 491 this->executeFiniteStrainBehaviourStressPostProcessing(wk, d); 492 } 493 } 494 // tangent operator (...) 495 if (ktype != StiffnessMatrixType::NOSTIFFNESS) { 496 const auto ndv = this->getGradientsSize(); 497 const auto nth = this->getThermodynamicForcesSize(); 498 if (this->btype == 2u) { 499 this->executeFiniteStrainBehaviourTangentOperatorPostProcessing(d); 500 } 501 if (this->stype == 1u) { 502 this->rto_fct(&(wk.D(0, 0)), &(wk.D(0, 0)), s.r.begin()); 503 } 504 if ((this->gtypes.size() == 1u) && (this->thtypes.size() == 1u)) { 505 for (unsigned short i = 0; i != nth; ++i) { 506 for (unsigned short j = 0; j != ndv; ++j) { 507 Kt(i, j) = wk.D(i, j); 508 } 509 } 510 } else { 511 const auto h = this->getHypothesis(); 512 auto to_offset = size_type{}; 513 for (const auto& to : this->getTangentOperatorBlocks()) { 514 auto getVariableOffset = [&h](const size_type pos, 515 const std::vector<int>& types) { 516 auto o = size_type{}; 517 for (size_type p = 0; p != pos; ++p) { 518 o += mtest::getVariableSize(types[p], h); 519 } 520 return o; 521 }; 522 const auto ptf = 523 std::find(this->thnames.begin(), this->thnames.end(), to.first); 524 const auto piv = 525 std::find(this->ivnames.begin(), this->ivnames.end(), to.first); 526 const auto pg = 527 std::find(this->gnames.begin(), this->gnames.end(), to.second); 528 const auto pev = 529 std::find(this->evnames.begin(), this->evnames.end(), to.second); 530 const auto to_ro = [&, this]() -> size_type { 531 if (pev != this->evnames.end()) { 532 return 1; 533 } 534 if (pg == this->gnames.end()) { 535 tfel::raise( 536 "GenericBehaviour::call_behaviour(1): " 537 "invalid tangent operator block ('" + 538 to.first + "'" + " vs '" + to.second + "')"); 539 } 540 return mtest::getVariableSize( 541 this->gtypes[pg - this->gnames.begin()], h); 542 }(); 543 const auto to_co = [&, this]() -> size_type { 544 if (piv != this->ivnames.end()) { 545 return mtest::getVariableSize( 546 this->ivtypes[piv - this->ivnames.begin()], h); 547 } 548 if (ptf == this->thnames.end()) { 549 tfel::raise( 550 "GenericBehaviour::call_behaviour(2): " 551 "invalid tangent operator block ('" + 552 to.first + "'" + " vs '" + to.second + "')"); 553 } 554 return mtest::getVariableSize( 555 this->thtypes[ptf - this->thnames.begin()], h); 556 }(); 557 if ((ptf == this->thnames.end()) || (pg == this->gnames.end())) { 558 to_offset += to_ro * to_co; 559 continue; 560 } 561 const auto tfpos = ptf - this->thnames.begin(); 562 const auto gpos = pg - this->gnames.begin(); 563 const auto og = getVariableOffset(gpos, this->gtypes); 564 const auto otf = getVariableOffset(tfpos, this->thtypes); 565 const auto g_size = mtest::getVariableSize(this->gtypes[gpos], h); 566 const auto th_size = mtest::getVariableSize(this->thtypes[tfpos], h); 567 const auto p = wk.D.begin(); 568 for (size_type i = 0; i != g_size; ++i) { 569 for (size_type j = 0; j != th_size; ++j) { 570 Kt(og + i, otf + j) = p[to_offset + i * th_size + j]; 571 } 572 } 573 to_offset += to_ro * to_co; 574 } 575 } 576 if (b && (this->stype == 1u)) { 577 this->rtf_fct(&s.s1[0], &s.s1[0], s.r.begin()); 578 } 579 } 580 return {true, d.rdt}; 581 } // end of GenericBehaviour::call_behaviour 582 executeFiniteStrainBehaviourStressPreProcessing(BehaviourWorkSpace & wk,mfront::gb::BehaviourData & d) const583 void GenericBehaviour::executeFiniteStrainBehaviourStressPreProcessing( 584 BehaviourWorkSpace& wk, mfront::gb::BehaviourData& d) const { 585 auto throw_if = [](const bool c, const std::string& m) { 586 tfel::raise_if(c, 587 "GenericBehaviour::" 588 "executeFiniteStrainBehaviourStressPreProcessing: " + 589 m); 590 }; 591 if (this->stress_measure == CAUCHY) { 592 d.K[1] = real(0); 593 } else if (this->stress_measure == PK2) { 594 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 595 d.K[1] = real(1); 596 if (n == 1u) { 597 if (this->getHypothesis() == 598 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 599 tfel::raise( 600 "GenericBehaviour::call_behaviour: " 601 "using DS_DEGL in plane strain is not supported yet"); 602 } 603 convertToSecondPiolaKirchhoffStress<1u>( 604 &wk.S0[0], d.s0.thermodynamic_forces, d.s0.gradients); 605 } else if (n == 2u) { 606 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 607 throw_if(true, "using DS_DEGL in plane strain is not supported yet"); 608 } 609 convertToSecondPiolaKirchhoffStress<2u>( 610 &wk.S0[0], d.s0.thermodynamic_forces, d.s0.gradients); 611 } else if (n == 3u) { 612 convertToSecondPiolaKirchhoffStress<3u>( 613 &wk.S0[0], d.s0.thermodynamic_forces, d.s0.gradients); 614 } else { 615 throw_if(true, "internal error, unsupported space dimension"); 616 } 617 d.s0.thermodynamic_forces = &wk.S0[0]; 618 d.s1.thermodynamic_forces = &wk.S1[0]; 619 } else if (this->stress_measure == PK1) { 620 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 621 d.K[1] = real(2); 622 if (n == 1u) { 623 if (this->getHypothesis() == 624 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 625 tfel::raise( 626 "GenericBehaviour::call_behaviour: " 627 "using DS_DEGL in plane strain is not supported yet"); 628 } 629 convertToFirstPiolaKirchhoffStress<1u>( 630 &wk.pk0[0], d.s0.thermodynamic_forces, d.s0.gradients); 631 } else if (n == 2u) { 632 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 633 throw_if(true, "using DS_DEGL in plane strain is not supported yet"); 634 } 635 convertToFirstPiolaKirchhoffStress<2u>( 636 &wk.pk0[0], d.s0.thermodynamic_forces, d.s0.gradients); 637 } else if (n == 3u) { 638 convertToFirstPiolaKirchhoffStress<3u>( 639 &wk.pk0[0], d.s0.thermodynamic_forces, d.s0.gradients); 640 } else { 641 throw_if(true, "internal error, unsupported space dimension"); 642 } 643 d.s0.thermodynamic_forces = &wk.pk0[0]; 644 d.s1.thermodynamic_forces = &wk.pk1[0]; 645 } else { 646 throw_if(true, "internal error, unexpected stress measure"); 647 } 648 } // end of 649 // GenericBehaviour::executeFiniteStrainBehaviourStressPreProcessing 650 651 void executeFiniteStrainBehaviourTangentOperatorPreProcessing(mfront::gb::BehaviourData & d,const StiffnessMatrixType ktype) const652 GenericBehaviour::executeFiniteStrainBehaviourTangentOperatorPreProcessing( 653 mfront::gb::BehaviourData& d, const StiffnessMatrixType ktype) const { 654 d.K[2] = real(0); 655 if (ktype != StiffnessMatrixType::NOSTIFFNESS) { 656 if (this->fsto == DSIG_DF) { 657 d.K[2] = real(0); 658 } else if (this->fsto == DS_DEGL) { 659 d.K[2] = real(1); 660 } else if (this->fsto == DPK1_DF) { 661 d.K[2] = real(2); 662 } else { 663 tfel::raise( 664 "GenericBehaviour::" 665 "executeFiniteStrainBehaviourTangentOperatorPreProcessing: " 666 "internal error, unexpected tangent operator type"); 667 } 668 } 669 } // end of 670 // GenericBehaviour::executeFiniteStrainBehaviourTangentOperatorPreProcessing 671 executeFiniteStrainBehaviourStressPostProcessing(BehaviourWorkSpace & wk,mfront::gb::BehaviourData & d) const672 void GenericBehaviour::executeFiniteStrainBehaviourStressPostProcessing( 673 BehaviourWorkSpace& wk, mfront::gb::BehaviourData& d) const { 674 auto throw_if = [](const bool c, const std::string& m) { 675 tfel::raise_if(c, 676 "GenericBehaviour::" 677 "executeFiniteStrainBehaviourStressPostProcessing: " + 678 m); 679 }; 680 if (this->stress_measure == CAUCHY) { 681 } else if (this->stress_measure == PK2) { 682 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 683 if (n == 1u) { 684 if (this->getHypothesis() == 685 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 686 tfel::raise( 687 "GenericBehaviour::call_behaviour: " 688 "using DS_DEGL in plane strain is not supported yet"); 689 } 690 convertFromSecondPiolaKirchhoffStress<1u>(d.s1.thermodynamic_forces, 691 &wk.S1[0], d.s1.gradients); 692 } else if (n == 2u) { 693 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 694 throw_if(true, "using DS_DEGL in plane strain is not supported yet"); 695 } 696 convertFromSecondPiolaKirchhoffStress<2u>(d.s1.thermodynamic_forces, 697 &wk.S1[0], d.s1.gradients); 698 } else if (n == 3u) { 699 convertFromSecondPiolaKirchhoffStress<3u>(d.s1.thermodynamic_forces, 700 &wk.S1[0], d.s1.gradients); 701 } else { 702 throw_if(true, "internal error, unsupported space dimension"); 703 } 704 } else if (this->stress_measure == PK1) { 705 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 706 if (n == 1u) { 707 if (this->getHypothesis() == 708 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 709 tfel::raise( 710 "GenericBehaviour::call_behaviour: " 711 "using DS_DEGL in plane strain is not supported yet"); 712 } 713 convertFromFirstPiolaKirchhoffStress<1u>(d.s1.thermodynamic_forces, 714 &wk.pk1[0], d.s1.gradients); 715 } else if (n == 2u) { 716 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 717 throw_if(true, "using DS_DEGL in plane strain is not supported yet"); 718 } 719 convertFromFirstPiolaKirchhoffStress<2u>(d.s1.thermodynamic_forces, 720 &wk.pk1[0], d.s1.gradients); 721 } else if (n == 3u) { 722 convertFromFirstPiolaKirchhoffStress<3u>(d.s1.thermodynamic_forces, 723 &wk.pk1[0], d.s1.gradients); 724 } else { 725 throw_if(true, "internal error, unsupported space dimension"); 726 } 727 } else { 728 throw_if(true, "internal error, unexpected stress measure"); 729 } 730 } // end of 731 // GenericBehaviour::executeFiniteStrainBehaviourStressPostProcessing 732 733 void executeFiniteStrainBehaviourTangentOperatorPostProcessing(mfront::gb::BehaviourData & d) const734 GenericBehaviour::executeFiniteStrainBehaviourTangentOperatorPostProcessing( 735 mfront::gb::BehaviourData& d) const { 736 if (this->fsto == DSIG_DF) { 737 // nothing to be done 738 } else if (this->fsto == DS_DEGL) { 739 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 740 if (n == 1u) { 741 if (this->getHypothesis() == 742 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 743 tfel::raise( 744 "GenericBehaviour::call_behaviour: " 745 "using DS_DEGL in plane strain is not supported yet"); 746 } else { 747 convertFromDSDEGL<1u>(d.K, d.s0.gradients, d.s1.gradients, 748 d.s1.thermodynamic_forces); 749 } 750 } else if (n == 2u) { 751 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 752 tfel::raise( 753 "GenericBehaviour::call_behaviour: " 754 "using DS_DEGL in plane strain is not supported yet"); 755 } else { 756 convertFromDSDEGL<2u>(d.K, d.s0.gradients, d.s1.gradients, 757 d.s1.thermodynamic_forces); 758 } 759 } else if (n == 3u) { 760 convertFromDSDEGL<3u>(d.K, d.s0.gradients, d.s1.gradients, 761 d.s1.thermodynamic_forces); 762 } else { 763 tfel::raise( 764 "GenericBehaviour::call_behaviour: " 765 "invalid space dimensions"); 766 } 767 } else if (this->fsto == DPK1_DF) { 768 const auto n = tfel::material::getSpaceDimension(this->getHypothesis()); 769 if (n == 1u) { 770 if (this->getHypothesis() == 771 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) { 772 tfel::raise( 773 "GenericBehaviour::call_behaviour: " 774 "using DS_DEGL in plane strain is not supported yet"); 775 } else { 776 convertFromDPK1DF<1u>(d.K, d.s0.gradients, d.s1.gradients, 777 d.s1.thermodynamic_forces); 778 } 779 } else if (n == 2u) { 780 if (this->getHypothesis() == ModellingHypothesis::PLANESTRESS) { 781 tfel::raise( 782 "GenericBehaviour::call_behaviour: " 783 "using DS_DEGL in plane strain is not supported yet"); 784 } else { 785 convertFromDPK1DF<2u>(d.K, d.s0.gradients, d.s1.gradients, 786 d.s1.thermodynamic_forces); 787 } 788 } else if (n == 3u) { 789 convertFromDPK1DF<3u>(d.K, d.s0.gradients, d.s1.gradients, 790 d.s1.thermodynamic_forces); 791 } else { 792 tfel::raise( 793 "GenericBehaviour::call_behaviour: " 794 "invalid space dimensions"); 795 } 796 } else { 797 tfel::raise( 798 "GenericBehaviour::call_behaviour: " 799 "internal error, unexpected tangent operator type"); 800 } 801 } // end of 802 // GenericBehaviour::executeFiniteStrainBehaviourTangentOperatorPostProcessing 803 getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> & r) const804 tfel::math::tmatrix<3u, 3u, real> GenericBehaviour::getRotationMatrix( 805 const tfel::math::vector<real>&, 806 const tfel::math::tmatrix<3u, 3u, real>& r) const { 807 return r; 808 } // end of GenericBehaviour::getRotationMatrix 809 getOptionalMaterialProperties() const810 std::vector<std::string> GenericBehaviour::getOptionalMaterialProperties() 811 const { 812 return {}; 813 } // end of GenericBehaviour::getOptionalMaterialProperties 814 setOptionalMaterialPropertiesDefaultValues(EvolutionManager &,const EvolutionManager &) const815 void GenericBehaviour::setOptionalMaterialPropertiesDefaultValues( 816 EvolutionManager&, const EvolutionManager&) const { 817 } // end of GenericBehaviour::setOptionalMaterialPropertiesDefaultValues 818 getDefaultStiffnessMatrixType() const819 StiffnessMatrixType GenericBehaviour::getDefaultStiffnessMatrixType() const { 820 return StiffnessMatrixType::CONSISTENTTANGENTOPERATOR; 821 } // end of GenericBehaviour::getDefaultStiffnessMatrixType 822 823 GenericBehaviour::~GenericBehaviour() = default; 824 825 } // end of namespace mtest 826