1 /*! 2 * \file mfront/src/InelasticFlowBase.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 15/03/2018 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 "TFEL/Raise.hxx" 15 #include "TFEL/Glossary/Glossary.hxx" 16 #include "TFEL/Glossary/GlossaryEntry.hxx" 17 #include "TFEL/Utilities/Data.hxx" 18 #include "MFront/ImplicitDSLBase.hxx" 19 #include "MFront/NonLinearSystemSolver.hxx" 20 #include "MFront/BehaviourBrick/BrickUtilities.hxx" 21 #include "MFront/BehaviourBrick/StressPotential.hxx" 22 #include "MFront/BehaviourBrick/StressCriterion.hxx" 23 #include "MFront/BehaviourBrick/OptionDescription.hxx" 24 #include "MFront/BehaviourBrick/StressCriterionFactory.hxx" 25 #include "MFront/BehaviourBrick/IsotropicHardeningRule.hxx" 26 #include "MFront/BehaviourBrick/IsotropicHardeningRuleFactory.hxx" 27 #include "MFront/BehaviourBrick/KinematicHardeningRule.hxx" 28 #include "MFront/BehaviourBrick/KinematicHardeningRuleFactory.hxx" 29 #include "MFront/BehaviourBrick/InelasticFlowBase.hxx" 30 31 namespace mfront { 32 33 namespace bbrick { 34 getOptions() const35 std::vector<OptionDescription> InelasticFlowBase::getOptions() const { 36 std::vector<OptionDescription> opts; 37 opts.emplace_back( 38 "criterion", 39 "stress criterion to be used (Mises, Hill, Hosford, Barlat)", 40 OptionDescription::DATASTRUCTURE); 41 opts.emplace_back("flow_criterion", 42 "stress criterion used to build the plastic potential " 43 "(Mises, Hill, Hosford, Barlat)", 44 OptionDescription::DATASTRUCTURE); 45 opts.emplace_back("isotropic_hardening", 46 "choice of an isotropic hardening rule", 47 OptionDescription::DATASTRUCTURES); 48 opts.emplace_back("kinematic_hardening", 49 "description of an hardening rule", 50 OptionDescription::DATASTRUCTURES); 51 opts.emplace_back( 52 "save_porosity_increase", 53 "if appropriate, save the porosity increase induced " 54 "by this inelastic flow in a dedicated auxiliary state variable", 55 OptionDescription::BOOLEAN); 56 opts.emplace_back( 57 "porosity_effect_on_equivalent_plastic_strain", 58 "specify the effect of the porosity of the flow rule. " 59 "Valid strings are 'StandardPorosityEffect' (or equivalently " 60 "'standard_porosity_effect') or 'None' (or equivalently 'none' " 61 "or 'false')'", 62 OptionDescription::STRING); 63 opts.emplace_back("porosity_evolution_algorithm", 64 "reserved for internal use", OptionDescription::STRING); 65 opts.emplace_back("cosine_threshold", 66 "minimum value of the cosine of the angle between two " 67 "successive estimates of the flow direction. If the " 68 "computed angle is lower than this threshold, the " 69 "Newton step is rejected. This parameter must be in " 70 "the range [-1:1].", 71 OptionDescription::REAL); 72 opts.emplace_back( 73 "cosine_check_minimum_iteration_factor", 74 "a factor $\\alpha$ which gives the threshold below which the " 75 "check on the cosine of the angle between two successive estimates " 76 "of the flow direction is performed, i.e. the test is performed if " 77 "the iteration counter is greater than $\\alpha \\cdot i_{\\max{}}$ " 78 "where $i_{\\max{}}$ is the maximum number of iterations. This " 79 "parameter must be in the range [0:1].", 80 OptionDescription::REAL); 81 opts.emplace_back( 82 "cosine_check_maximum_iteration_factor", 83 "a factor $\\alpha$ which gives the threshold upper which the " 84 "check on the cosine of the angle between two successive estimates " 85 "of the flow direction is performed, i.e. the test is performed if " 86 "the iteration counter is lower than $\\alpha \\cdot i_{\\max{}}$ " 87 "where $i_{\\max{}}$ is the maximum number of iterations. This " 88 "parameter must be in the range [0:1].", 89 OptionDescription::REAL); 90 return opts; 91 } // end of InelasticFlowBase::getOptions() 92 initialize(BehaviourDescription & bd,AbstractBehaviourDSL & dsl,const std::string & id,const DataMap & d)93 void InelasticFlowBase::initialize(BehaviourDescription& bd, 94 AbstractBehaviourDSL& dsl, 95 const std::string& id, 96 const DataMap& d) { 97 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 98 auto raise = [](const std::string& m) { 99 tfel::raise("InelasticFlowBase::initialize: " + m); 100 }; // end of raise 101 auto getDataStructure = [&raise](const std::string& n, const Data& ds) { 102 if (ds.is<std::string>()) { 103 tfel::utilities::DataStructure nds; 104 nds.name = ds.get<std::string>(); 105 return nds; 106 } 107 if (!ds.is<tfel::utilities::DataStructure>()) { 108 raise("invalid data type for entry '" + n + "'"); 109 } 110 return ds.get<tfel::utilities::DataStructure>(); 111 }; // end of getDataStructure 112 // checking options 113 mfront::bbrick::check(d, this->getOptions()); 114 // parsing options 115 for (const auto& e : d) { 116 if (e.first == "criterion") { 117 if (this->sc != nullptr) { 118 raise("criterion has already been defined"); 119 } 120 const auto ds = getDataStructure(e.first, e.second); 121 auto& cf = StressCriterionFactory::getFactory(); 122 this->sc = cf.generate(ds.name); 123 if (d.count("flow_criterion") != 0) { 124 this->sc->initialize(bd, dsl, id, ds.data, 125 StressCriterion::STRESSCRITERION); 126 } else { 127 this->sc->initialize(bd, dsl, id, ds.data, 128 StressCriterion::STRESSANDFLOWCRITERION); 129 } 130 } else if (e.first == "flow_criterion") { 131 if (this->fc != nullptr) { 132 raise("criterion has already been defined"); 133 } 134 const auto ds = getDataStructure(e.first, e.second); 135 auto& cf = StressCriterionFactory::getFactory(); 136 this->fc = cf.generate(ds.name); 137 this->fc->initialize(bd, dsl, id, ds.data, 138 StressCriterion::FLOWCRITERION); 139 } else if (e.first == "isotropic_hardening") { 140 auto add_isotropic_hardening_rule = 141 [this, &bd, &dsl, &id](const tfel::utilities::DataStructure& ds) { 142 auto& rf = IsotropicHardeningRuleFactory::getFactory(); 143 const auto ihr = rf.generate(ds.name); 144 ihr->initialize(bd, dsl, id, std::to_string(this->ihrs.size()), 145 ds.data); 146 this->ihrs.push_back(ihr); 147 }; 148 if (e.second.is<std::vector<Data>>()) { 149 for (const auto& ird : e.second.get<std::vector<Data>>()) { 150 add_isotropic_hardening_rule(getDataStructure(e.first, ird)); 151 } 152 } else { 153 auto& rf = IsotropicHardeningRuleFactory::getFactory(); 154 const auto ds = getDataStructure(e.first, e.second); 155 const auto ihr = rf.generate(ds.name); 156 ihr->initialize(bd, dsl, id, "", ds.data); 157 this->ihrs.push_back(ihr); 158 } 159 } else if (e.first == "kinematic_hardening") { 160 auto add_kinematic_hardening_rule = 161 [this, &bd, &dsl, &id](const tfel::utilities::DataStructure& ds) { 162 auto& kf = KinematicHardeningRuleFactory::getFactory(); 163 const auto k = kf.generate(ds.name); 164 k->initialize(bd, dsl, id, std::to_string(this->khrs.size()), 165 ds.data); 166 this->khrs.push_back(k); 167 }; 168 if (e.second.is<std::vector<Data>>()) { 169 for (const auto& krd : e.second.get<std::vector<Data>>()) { 170 add_kinematic_hardening_rule(getDataStructure(e.first, krd)); 171 } 172 } else { 173 add_kinematic_hardening_rule(getDataStructure(e.first, e.second)); 174 } 175 } else if (e.first == "save_porosity_increase") { 176 if (!e.second.is<bool>()) { 177 raise("'save_porosity_increase' is not a boolean"); 178 } 179 this->save_porosity_increase = e.second.get<bool>(); 180 } else if (e.first == "porosity_evolution_algorithm") { 181 if (!e.second.is<std::string>()) { 182 raise("'porosity_evolution_algorithm' is not a boolean"); 183 } 184 const auto& a = e.second.get<std::string>(); 185 if (a == "standard_implicit_scheme") { 186 this->porosity_evolution_algorithm = 187 PorosityEvolutionAlgorithm::STANDARD_IMPLICIT_SCHEME; 188 } else if (a == "staggered_scheme") { 189 this->porosity_evolution_algorithm = 190 PorosityEvolutionAlgorithm::STAGGERED_SCHEME; 191 } else { 192 raise("internal error: unsupported porosity evolution algorithm"); 193 } 194 } else if (e.first == "porosity_effect_on_equivalent_plastic_strain") { 195 if (!e.second.is<std::string>()) { 196 raise("'porosity_effect_on_equivalent_plastic_strain' is not a string"); 197 } 198 const auto& pe = e.second.get<std::string>(); 199 if ((pe == "StandardPorosityEffect") || 200 (pe == "standard_porosity_effect")) { 201 this->porosity_effect_on_equivalent_plastic_strain = 202 STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN; 203 } else if ((pe == "None") || (pe == "none") || (pe == "false")) { 204 this->porosity_effect_on_equivalent_plastic_strain = 205 NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN; 206 } else if (pe != "Undefined") { 207 raise( 208 "invalid value of 'porosity_effect_on_equivalent_plastic_strain'. Expected " 209 "'StandardPorosityEffect' (or equivalently " 210 "'standard_porosity_effect') or 'None' (or equivalently 'none' " 211 "or 'false'), but got '" + 212 pe + "'"); 213 } 214 } else if (e.first == "cosine_threshold") { 215 this->cosine_threshold = tfel::utilities::convert<double>(e.second); 216 if ((this->cosine_threshold < -1) || (this->cosine_threshold > 1)) { 217 raise( 218 "invalid value for the option `cosine_threshold` (value must " 219 "be in range [-1:1]"); 220 } 221 } else if (e.first == "cosine_check_minimum_iteration_factor") { 222 if (d.count("cosine_threshold") == 0) { 223 raise( 224 "option `cosine_check_minimum_iteration_factor` is only " 225 "meaningful if the option `cosine_threshold` is also " 226 "provided"); 227 } 228 this->cosine_check_minimum_iteration_factor = 229 tfel::utilities::convert<double>(e.second); 230 if ((this->cosine_check_minimum_iteration_factor < -1) || 231 (this->cosine_check_minimum_iteration_factor > 1)) { 232 raise( 233 "invalid value for the option " 234 "`cosine_check_minimum_iteration_factor` (value must be in " 235 "range [-1:1]"); 236 } 237 if (this->cosine_check_maximum_iteration_factor < 238 this->cosine_check_minimum_iteration_factor) { 239 raise( 240 "the value of the option " 241 "`cosine_check_maximum_iteration_factor` must be greater than " 242 "the vaue of `cosine_check_minmum_iteration_factor`"); 243 } 244 } else if (e.first == "cosine_check_maximum_iteration_factor") { 245 if (d.count("cosine_threshold") == 0) { 246 raise( 247 "option `cosine_check_maximum_iteration_factor` is only " 248 "meaningful if the option `cosine_threshold` is also " 249 "provided"); 250 } 251 this->cosine_check_maximum_iteration_factor = 252 tfel::utilities::convert<double>(e.second); 253 if ((this->cosine_check_maximum_iteration_factor < -1) || 254 (this->cosine_check_maximum_iteration_factor > 1)) { 255 raise( 256 "invalid value for the option " 257 "`cosine_check_maximum_iteration_factor` (value must be in " 258 "range [-1:1]"); 259 } 260 if (this->cosine_check_maximum_iteration_factor < 261 this->cosine_check_minimum_iteration_factor) { 262 raise( 263 "the value of the option " 264 "`cosine_check_maximum_iteration_factor` must be greater than " 265 "the vaue of `cosine_check_minmum_iteration_factor`"); 266 } 267 } // other options must be treated in child classes 268 } 269 if (this->cosine_threshold < 1.5) { 270 addLocalVariable(bd, "Stensor", "np" + id); 271 CodeBlock init; 272 init.code = "this->np" + id + 273 " = Stensor(real(0));\n"; 274 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, init, 275 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 276 } 277 if (this->sc == nullptr) { 278 raise("criterion has not been defined"); 279 } 280 if (!this->ihrs.empty()) { 281 addLocalVariable(bd, "bool", "bpl" + id); 282 if (this->ihrs.size() != 1) { 283 const auto Rel = "Rel" + id; 284 const auto R = "R" + id; 285 const auto dR = "d" + R + "_ddp" + id; 286 bd.reserveName(uh, Rel); 287 bd.reserveName(uh, R); 288 bd.reserveName(uh, dR); 289 } 290 } 291 if (this->isCoupledWithPorosityEvolution()) { 292 if (this->fc != nullptr) { 293 if (this->fc->isNormalDeviatoric()) { 294 raise( 295 "InelasticFlowBase::initialize: " 296 "the flow is coupled with the porosity evolution, but " 297 "the flow direction is deviatoric"); 298 } 299 } 300 if (this->sc->isNormalDeviatoric()) { 301 raise( 302 "InelasticFlowBase::initialize: " 303 "the flow is coupled with the porosity evolution by " 304 "the flow direction is deviatoric"); 305 } 306 if (!this->khrs.empty()) { 307 raise( 308 "InelasticFlowBase::initialize: " 309 "kinematic hardening rules are not supported " 310 "when coupled with a porosity evolution"); 311 } 312 } 313 if (this->contributesToPorosityGrowth()) { 314 addLocalVariable(bd, "real", "trace_n" + id); 315 if (this->save_porosity_increase) { 316 addLocalVariable(bd, "real", "dfg" + id); 317 VariableDescription fg("real", "fg" + id, 1u, 0u); 318 const auto& g = 319 tfel::glossary::Glossary::PorosityIncreaseDueToInelasticFlow; 320 if (id.empty()) { 321 fg.setGlossaryName(g); 322 } else { 323 fg.setEntryName(g.getKey() + id); 324 } 325 bd.addAuxiliaryStateVariable(uh, fg); 326 } 327 } 328 } // end of InelasticFlowBase::initialize 329 setPorosityEvolutionHandled(const bool b)330 void InelasticFlowBase::setPorosityEvolutionHandled(const bool b) { 331 if ((!b) && (this->isCoupledWithPorosityEvolution())) { 332 tfel::raise( 333 "InelasticFlowBase::setPorosityEvolutionHandled: " 334 "internal error (consistent situation)"); 335 } 336 this->porosity_evolution_explicitely_handled = b; 337 } // end of InelasticFlowBase::setPorosityEvolutionHandled 338 contributesToPorosityGrowth() const339 bool InelasticFlowBase::contributesToPorosityGrowth() const { 340 if ((this->porosity_evolution_explicitely_handled == true) || 341 (this->isCoupledWithPorosityEvolution())) { 342 if (this->sc == nullptr) { 343 tfel::raise( 344 "InelasticFlowBase::contributesToPorosityGrowth: " 345 "unitialised object"); 346 } 347 const bool bn = [this] { 348 if (this->fc != nullptr) { 349 return !this->fc->isNormalDeviatoric(); 350 } 351 return !this->sc->isNormalDeviatoric(); 352 }(); 353 return bn; 354 } 355 return false; 356 } // end of InelasticFlowBase::contributesToPorosityGrowth 357 completeVariableDeclaration(BehaviourDescription &,const AbstractBehaviourDSL &,const std::string &) const358 void InelasticFlowBase::completeVariableDeclaration( 359 BehaviourDescription&, 360 const AbstractBehaviourDSL&, 361 const std::string&) const { 362 } // end of InelasticFlowBase::completeVariableDeclaration 363 requiresActivationState() const364 bool InelasticFlowBase::requiresActivationState() const { 365 return !this->ihrs.empty(); 366 } // end of InelasticFlowBase::requiresActivationState 367 computeInitialActivationState(BehaviourDescription & bd,const StressPotential & sp,const std::string & id) const368 void InelasticFlowBase::computeInitialActivationState( 369 BehaviourDescription& bd, 370 const StressPotential& sp, 371 const std::string& id) const { 372 if (!this->ihrs.empty()) { 373 CodeBlock i; 374 if (this->khrs.empty()) { 375 i.code += "const auto& sel" + id + " = sigel;\n"; 376 } else { 377 auto kid = decltype(khrs.size()){}; 378 for (const auto& khr : khrs) { 379 i.code += khr->computeKinematicHardeningsInitialValues( 380 id, std::to_string(kid)); 381 ++kid; 382 } 383 // effective stress 384 i.code += "const auto sel" + id + " = eval(sigel"; 385 kid = decltype(khrs.size()){}; 386 for (const auto& khr : khrs) { 387 for (const auto& X : khr->getKinematicHardeningsVariables( 388 id, std::to_string(kid))) { 389 i.code += "-" + X; 390 ++kid; 391 } 392 } 393 i.code += ");\n"; 394 } 395 i.code += this->sc->computeElasticPrediction(id, bd, sp); 396 i.code += computeElasticLimitInitialValue(this->ihrs, id); 397 i.code += "this->bpl" + id + " = seqel" + id + " > Rel" + id + ";\n"; 398 bd.setCode(ModellingHypothesis::UNDEFINEDHYPOTHESIS, 399 BehaviourData::BeforeInitializeLocalVariables, i, 400 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 401 } 402 } // end of InelasticFlowBase::computeInitialActivationState 403 computeEffectiveStress(const std::string & id) const404 std::string InelasticFlowBase::computeEffectiveStress( 405 const std::string& id) const { 406 if (this->khrs.empty()) { 407 return "const auto& s" + id + " = this->sig;\n"; 408 } 409 auto c = std::string{}; 410 auto kid = decltype(khrs.size()){}; 411 for (const auto& khr : khrs) { 412 c += khr->computeKinematicHardenings(id, std::to_string(kid)); 413 ++kid; 414 } 415 // effective stress 416 c += "const auto s" + id + " = eval(this->sig"; 417 kid = decltype(khrs.size()){}; 418 for (const auto& khr : khrs) { 419 for (const auto& X : 420 khr->getKinematicHardeningsVariables(id, std::to_string(kid))) { 421 c += "-" + X + "_"; 422 } 423 ++kid; 424 } 425 c += ");\n"; 426 return c; 427 } // end of InelasticFlowBase::computeEffectiveStress 428 endTreatment(BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const StressPotential & sp,const std::string & id) const429 void InelasticFlowBase::endTreatment(BehaviourDescription& bd, 430 const AbstractBehaviourDSL& dsl, 431 const StressPotential& sp, 432 const std::string& id) const { 433 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 434 const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(dsl); 435 if (this->fc != nullptr) { 436 this->sc->endTreatment(bd, dsl, id, StressCriterion::STRESSCRITERION); 437 this->fc->endTreatment(bd, dsl, id, StressCriterion::FLOWCRITERION); 438 } else { 439 this->sc->endTreatment(bd, dsl, id, 440 StressCriterion::STRESSANDFLOWCRITERION); 441 } 442 const bool is_deviatoric = (this->fc != nullptr) 443 ? this->fc->isNormalDeviatoric() 444 : this->sc->isNormalDeviatoric(); 445 auto iid = decltype(ihrs.size()){}; 446 if (this->ihrs.size() == 1) { 447 ihrs[0]->endTreatment(bd, dsl, id, ""); 448 } else { 449 for (const auto& ihr : this->ihrs) { 450 ihr->endTreatment(bd, dsl, id, std::to_string(iid)); 451 ++iid; 452 } 453 } 454 auto kid = decltype(khrs.size()){}; 455 for (const auto& khr : this->khrs) { 456 khr->endTreatment(bd, dsl, id, std::to_string(kid)); 457 ++kid; 458 } 459 const auto requiresAnalyticalJacobian = 460 ((idsl.getSolver().usesJacobian()) && 461 (!idsl.getSolver().requiresNumericalJacobian())); 462 // implicit equation associated with the elastic strain 463 CodeBlock ib; 464 if (!this->ihrs.empty()) { 465 ib.code += "if(this->bpl" + id + "){\n"; 466 } 467 if (idsl.getSolver().usesJacobian()) { 468 ib.code += "if(!perturbatedSystemEvaluation){\n"; 469 } 470 // ib.code += "this->dp" + id + " = this->dp" + id + ", 471 // strain(0));\n"; 472 if (idsl.getSolver().usesJacobian()) { 473 ib.code += "}\n"; 474 } 475 ib.code += this->computeEffectiveStress(id); 476 if (requiresAnalyticalJacobian) { 477 if (this->fc == nullptr) { 478 ib.code += this->sc->computeNormalDerivative( 479 id, bd, sp, StressCriterion::STRESSANDFLOWCRITERION); 480 } else { 481 ib.code += this->sc->computeNormal(id, bd, sp, 482 StressCriterion::STRESSCRITERION); 483 ib.code += this->fc->computeNormalDerivative( 484 id, bd, sp, StressCriterion::FLOWCRITERION); 485 } 486 } else { 487 if (this->fc == nullptr) { 488 ib.code += this->sc->computeNormal( 489 id, bd, sp, StressCriterion::STRESSANDFLOWCRITERION); 490 } else { 491 ib.code += this->sc->computeCriterion(id, bd, sp); 492 ib.code += this->fc->computeNormal(id, bd, sp, 493 StressCriterion::FLOWCRITERION); 494 } 495 } 496 // check on the flow direction 497 if (this->cosine_threshold < 1.5) { 498 const auto imin = 499 std::to_string(this->cosine_check_minimum_iteration_factor) + 500 " * this->iterMax"; 501 const auto imax = 502 std::to_string(this->cosine_check_maximum_iteration_factor) + 503 " * this->iterMax"; 504 const auto seps = sp.getEquivalentStressLowerBound(bd); 505 if (this->fc == nullptr) { 506 ib.code += "if((seq" + id + ">" + seps + ")&&"; 507 } else { 508 ib.code += "if((seqf" + id + ">" + seps + ")&&"; 509 } 510 ib.code += "(this->iter > " + imin + ")&&"; 511 ib.code += "(this->iter < " + imax + ")){\n"; 512 ib.code += "if (std::abs(n" + id + " | this->np" + id + ") < "; 513 ib.code += "std::sqrt((n" + id + ") | (n" + id + ")) * "; 514 ib.code += "std::sqrt((this->np" + id + ") | (this->np" + id + ")) * " + 515 std::to_string(this->cosine_threshold) + ") {\n"; 516 ib.code += "return false;\n"; 517 ib.code += "}\n"; 518 ib.code += "}\n"; 519 ib.code += "this->np" + id + " = n" + id + ";\n"; 520 } 521 if (this->isCoupledWithPorosityEvolution()) { 522 ib.code += "this->trace_n" + id + " = trace(n" + id + ");\n"; 523 } 524 // elasticity 525 if (this->getPorosityEffectOnEquivalentPlasticStrain() == 526 STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN) { 527 const auto& f = 528 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 529 tfel::glossary::Glossary::Porosity); 530 const auto f_ = f.name + "_"; 531 ib.code += 532 "feel += (1 - " + f_ + ") * this->dp" + id + "* n" + id + ";\n"; 533 if (requiresAnalyticalJacobian) { 534 // jacobian terms 535 ib.code += "dfeel_ddp" + id + " = (1 - " + f_ + ") * n" + id + ";\n"; 536 ib.code += sp.generateImplicitEquationDerivatives( 537 bd, "StrainStensor", "eel", 538 "(1 - " + f_ + ") * (this->dp" + id + ") * dn" + id + "_ds" + id, 539 is_deviatoric); 540 kid = decltype(khrs.size()){}; 541 for (const auto& khr : khrs) { 542 ib.code += khr->generateImplicitEquationDerivatives( 543 "eel", "- (1 - " + f_ + ") * (this->dp" + id + ") * dn" + id + 544 "_ds" + id, 545 id, std::to_string(kid)); 546 ++kid; 547 } 548 ib.code += "dfeel_dd" + f.name + " += "; 549 ib.code += "- theta * this->dp" + id + "* n" + id; 550 ib.code += "+ theta * (1 - " + f_ + ") * this->dp" + id + "* dn" + 551 id + "_d" + f.name + ";\n"; 552 } 553 } else { 554 ib.code += "feel += this->dp" + id + "* n" + id + ";\n"; 555 if (requiresAnalyticalJacobian) { 556 // jacobian terms 557 ib.code += "dfeel_ddp" + id + " = n" + id + ";\n"; 558 ib.code += sp.generateImplicitEquationDerivatives( 559 bd, "StrainStensor", "eel", 560 "(this->dp" + id + ") * dn" + id + "_ds" + id, is_deviatoric); 561 kid = decltype(khrs.size()){}; 562 for (const auto& khr : khrs) { 563 ib.code += khr->generateImplicitEquationDerivatives( 564 "eel", "-(this->dp" + id + ") * dn" + id + "_ds" + id, id, 565 std::to_string(kid)); 566 ++kid; 567 } 568 if (this->isCoupledWithPorosityEvolution()) { 569 const auto& f = bd.getBehaviourData(uh) 570 .getStateVariableDescriptionByExternalName( 571 tfel::glossary::Glossary::Porosity); 572 ib.code += "dfeel_dd" + f.name + " += "; 573 ib.code += 574 "theta * this->dp" + id + " * dn" + id + "_d" + f.name + ";\n"; 575 } 576 } 577 } 578 // inelastic flow 579 ib.code += this->buildFlowImplicitEquations(bd, sp, id, 580 requiresAnalyticalJacobian); 581 // hardening rules 582 kid = decltype(khrs.size()){}; 583 for (const auto& khr : khrs) { 584 if (this->fc != nullptr) { 585 ib.code += khr->buildBackStrainImplicitEquations( 586 bd, sp, *(this->fc), this->khrs, id, std::to_string(kid), 587 requiresAnalyticalJacobian); 588 } else { 589 ib.code += khr->buildBackStrainImplicitEquations( 590 bd, sp, *(this->sc), this->khrs, id, std::to_string(kid), 591 requiresAnalyticalJacobian); 592 } 593 ++kid; 594 } 595 // porosity evolution 596 if (this->contributesToPorosityGrowth()) { 597 if (this->porosity_evolution_algorithm == 598 PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 599 ib.code += "if("; 600 ib.code += StandardElastoViscoPlasticityBrick:: 601 computeStandardSystemOfImplicitEquations; 602 ib.code += "){\n"; 603 this->addFlowContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 604 ib, bd, dsl, sp, id); 605 ib.code += "}\n"; 606 } else if (this->porosity_evolution_algorithm == 607 PorosityEvolutionAlgorithm::STANDARD_IMPLICIT_SCHEME) { 608 this->addFlowContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 609 ib, bd, dsl, sp, id); 610 } else { 611 tfel::raise( 612 "InelasticFlowBase::endTreatment: internal error " 613 "(unsupported porosity evolution algorithm)"); 614 } 615 } 616 if (!this->ihrs.empty()) { 617 ib.code += "} // end if(this->bpl" + id + ")\n"; 618 } 619 bd.setCode(uh, BehaviourData::Integrator, ib, 620 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 621 // additional checks 622 if (!this->ihrs.empty()) { 623 CodeBlock acc; 624 acc.code += "if (converged) {\n"; 625 acc.code += "// checking status consistency\n"; 626 acc.code += "if (this->bpl" + id + ") {\n"; 627 acc.code += "if (this->dp" + id + " < -2*this->epsilon) {\n"; 628 acc.code += "// desactivating this system\n"; 629 acc.code += "converged = this->bpl" + id + " = false;\n"; 630 acc.code += "}\n"; 631 acc.code += "} else {\n"; 632 if (this->isCoupledWithPorosityEvolution()) { 633 const auto& f = 634 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 635 tfel::glossary::Glossary::Porosity); 636 acc.code += "const auto " + f.name + "_ = "; 637 acc.code += "this->" + f.name + " + (" + bd.getClassName() + 638 "::theta) * (this->d" + f.name + ");\n"; 639 } 640 acc.code += this->computeEffectiveStress(id); 641 acc.code += this->sc->computeCriterion(id, bd, sp); 642 acc.code += computeElasticLimit(this->ihrs, id); 643 acc.code += "if(seq" + id + " > R" + id + ") {\n"; 644 acc.code += "converged = false;\n"; 645 acc.code += "this->bpl" + id + " = true;\n"; 646 acc.code += "}\n"; 647 acc.code += "}\n"; 648 acc.code += "} // end of if(converged)\n"; 649 bd.setCode(uh, BehaviourData::AdditionalConvergenceChecks, acc, 650 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 651 } 652 // update auxiliary state variables 653 if ((this->contributesToPorosityGrowth()) && 654 (this->save_porosity_increase)) { 655 CodeBlock uav; 656 uav.code += "fg" + id + " += this->dfg" + id + ";\n"; 657 bd.setCode(uh, BehaviourData::UpdateAuxiliaryStateVariables, uav, 658 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 659 } 660 } // end of InelasticFlowBase::endTreatment 661 updateNextEstimateOfThePorosityIncrement(const BehaviourDescription & bd,const std::string & id) const662 std::string InelasticFlowBase::updateNextEstimateOfThePorosityIncrement( 663 const BehaviourDescription& bd, const std::string& id) const { 664 if (!this->contributesToPorosityGrowth()) { 665 return ""; 666 } 667 const auto df = [this, &bd, &id] { 668 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 669 const auto& f = 670 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 671 tfel::glossary::Glossary::Porosity); 672 const auto theta = bd.getClassName() + "::theta"; 673 const auto f_ = "((this->" + f.name + ") + (" + theta + ") * (this->" + 674 std::string(StandardElastoViscoPlasticityBrick:: 675 currentEstimateOfThePorosityIncrement) + 676 "))"; 677 if (this->getPorosityEffectOnEquivalentPlasticStrain() == 678 STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN) { 679 return "power<2>(1 - " + f_ + ") * " + // 680 "(this->dp" + id + ") * (this->trace_n" + id + ")"; 681 } 682 return "(1 - " + f_ + ") * " + // 683 "(this->dp" + id + ") * (this->trace_n" + id + ")"; 684 }(); 685 auto c = std::string{}; 686 c = "if (this->bpl" + id + "){\n"; 687 if (this->save_porosity_increase) { 688 c += "this->dfg" + id + " = " + df + ";\n"; 689 c += StandardElastoViscoPlasticityBrick:: 690 nextEstimateOfThePorosityIncrement; 691 c += " += this->dfg" + id + ";\n"; 692 } else { 693 c += StandardElastoViscoPlasticityBrick:: 694 nextEstimateOfThePorosityIncrement; 695 c += " += " + df + ";\n"; 696 } 697 if (this->save_porosity_increase) { 698 c += "} else {\n"; 699 c += "this->dfg" + id + " = real{};\n"; 700 } 701 c += "}\n"; 702 return c; 703 } // end of updateNextEstimateOfThePorosityIncrement 704 705 void InelasticFlowBase:: addFlowContributionToTheImplicitEquationAssociatedWithPorosityEvolution(CodeBlock & ib,const BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const StressPotential & sp,const std::string & id) const706 addFlowContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 707 CodeBlock& ib, 708 const BehaviourDescription& bd, 709 const AbstractBehaviourDSL& dsl, 710 const StressPotential& sp, 711 const std::string& id) const { 712 if (!this->contributesToPorosityGrowth()) { 713 return; 714 } 715 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 716 const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(dsl); 717 const auto requiresAnalyticalJacobian = 718 ((idsl.getSolver().usesJacobian()) && 719 (!idsl.getSolver().requiresNumericalJacobian())); 720 const auto& f = 721 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 722 tfel::glossary::Glossary::Porosity); 723 const auto& broken = 724 bd.getBehaviourData(uh) 725 .getAuxiliaryStateVariableDescriptionByExternalName( 726 tfel::glossary::Glossary::Broken); 727 const auto f_ = f.name + "_"; 728 ib.code += "if(2 * (this->" + broken.name + ") > 1){\n"; 729 if (this->save_porosity_increase) { 730 ib.code += "this->dfg" + id + " = real{};\n"; 731 ib.code += "f" + f.name + " -= this->dfg" + id + ";\n"; 732 } else { 733 ib.code += "f" + f.name + " -= real{};\n"; 734 } 735 ib.code += "} else {\n"; 736 if (this->getPorosityEffectOnEquivalentPlasticStrain() == 737 STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN) { 738 if (this->save_porosity_increase) { 739 ib.code += "this->dfg" + id + " = power<2>(1 - " + f_ + 740 ") * (this->dp" + id + ") * (this->trace_n" + id + ");\n"; 741 ib.code += "f" + f.name + " -= this->dfg" + id + ";\n"; 742 } else { 743 ib.code += "f" + f.name; 744 ib.code += " -= power<2>(1 - " + f_ + ") * (this->dp" + id + 745 ") * (this->trace_n" + id + ");\n"; 746 } 747 } else { 748 if (this->save_porosity_increase) { 749 ib.code += "this->dfg" + id + " = (1 - " + f_ + ") * (this->dp" + id + 750 ") * (this->trace_n" + id + ");\n"; 751 ib.code += "f" + f.name + " -= this->dfg" + id + ";\n"; 752 } else { 753 ib.code += "f" + f.name; 754 ib.code += " -= (1 - " + f_ + ") * (this->dp" + id + 755 ") * (this->trace_n" + id + ");\n"; 756 } 757 } 758 if (requiresAnalyticalJacobian) { 759 if (this->getPorosityEffectOnEquivalentPlasticStrain() == 760 STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN) { 761 ib.code += "df" + f.name + "_dd" + f.name; 762 ib.code += " += 2 * theta * (1 - " + f_ + ") * (this->dp" + id + 763 ") * (this->trace_n" + id + ")"; 764 ib.code += " - theta * power<2>(1 - " + f_ + ")* (this->dp" + 765 id + ") * (Stensor::Id() | dn_d" + f.name + ");\n"; 766 // derivatives with respect to stress 767 for (const auto& dsig : sp.getStressDerivatives(bd)) { 768 const auto& dsig_ddv = std::get<0>(dsig); 769 const auto& v = std::get<1>(dsig); 770 const auto& t = std::get<2>(dsig); 771 if ((t != SupportedTypes::SCALAR) && 772 (t != SupportedTypes::STENSOR)) { 773 tfel::raise( 774 "InelasticFlowBase::endTreatment: " 775 "stress dependency on variable '" + 776 v + "' is not unsupported "); 777 } 778 ib.code += "df" + f.name + "_dd" + v + " -= "; 779 ib.code += "power<2>(1 - " + f_ + ") * (this->dp" + id + ") * "; 780 ib.code += "(Stensor::Id() | (dn" + id + "_ds" + id + " * (" + 781 dsig_ddv + ")));\n"; 782 } 783 auto kid = decltype(khrs.size()){}; 784 for (const auto& khr : this->khrs) { 785 const auto an = khr->getBackStrainVariable(id, std::to_string(kid)); 786 const auto dX = 787 khr->getBackStressDerivative(id, std::to_string(kid)); 788 ib.code += "df" + f.name + "_dd" + an + " += "; 789 ib.code += "power<2>(1 - " + f_ + ") * (this->dp" + id + ") * "; 790 ib.code += "(Stensor::Id() | " + dX + ");\n"; 791 ++kid; 792 } 793 ib.code += "df" + f.name + "_ddp" + id; 794 ib.code += 795 " = - power<2>(1 - " + f_ + ") * (this->trace_n" + id + ");\n"; 796 } else { 797 ib.code += "df" + f.name + "_dd" + f.name; 798 ib.code += 799 " += theta *(this->dp" + id + ") * (this->trace_n" + id + ")"; 800 ib.code += " - theta * (1 - " + f_ + ")* (this->dp" + id + 801 ") * (Stensor::Id() | dn_d" + f.name + ");\n"; 802 // derivatives with respect to stress 803 for (const auto& dsig : sp.getStressDerivatives(bd)) { 804 const auto& dsig_ddv = std::get<0>(dsig); 805 const auto& v = std::get<1>(dsig); 806 const auto& t = std::get<2>(dsig); 807 if ((t != SupportedTypes::SCALAR) && 808 (t != SupportedTypes::STENSOR)) { 809 tfel::raise( 810 "InelasticFlowBase::endTreatment: " 811 "stress dependency on variable '" + 812 v + "' is not unsupported "); 813 } 814 ib.code += "df" + f.name + "_dd" + v + " -= "; 815 ib.code += "(1 - " + f_ + ") * (this->dp" + id + ") * "; 816 ib.code += "(Stensor::Id() | (dn" + id + "_ds" + id + " * (" + 817 dsig_ddv + ")));\n"; 818 } 819 auto kid = decltype(khrs.size()){}; 820 for (const auto& khr : this->khrs) { 821 const auto an = khr->getBackStrainVariable(id, std::to_string(kid)); 822 const auto dX = 823 khr->getBackStressDerivative(id, std::to_string(kid)); 824 ib.code += "df" + f.name + "_dd" + an + " += "; 825 ib.code += "(1 - " + f_ + ") * (this->dp" + id + ") * "; 826 ib.code += "(Stensor::Id() | " + dX + ");\n"; 827 ++kid; 828 } 829 ib.code += "df" + f.name + "_ddp" + id; 830 ib.code += " = - (1 - " + f_ + ") * (this->trace_n" + id + ");\n"; 831 } 832 } 833 ib.code += "}\n"; 834 } // end of 835 // InelasticFlowBase::addFlowContributionToTheImplicitEquationAssociatedWithPorosityEvolution 836 isCoupledWithPorosityEvolution() const837 bool InelasticFlowBase::isCoupledWithPorosityEvolution() const { 838 if (this->sc == nullptr) { 839 tfel::raise( 840 "InelasticFlowBase::isCoupledWithPorosityEvolution: " 841 "uninitialized flow"); 842 } 843 if (this->sc->isCoupledWithPorosityEvolution()) { 844 return true; 845 } 846 if (this->fc != nullptr) { 847 return this->fc->isCoupledWithPorosityEvolution(); 848 } 849 return false; 850 } // end of InelasticFlowBase::isCoupledWithPorosityEvolution 851 852 InelasticFlowBase::PorosityEffectOnFlowRule getPorosityEffectOnEquivalentPlasticStrain() const853 InelasticFlowBase::getPorosityEffectOnEquivalentPlasticStrain() const { 854 if (this->porosity_effect_on_equivalent_plastic_strain == 855 UNDEFINED_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN) { 856 if (this->sc == nullptr) { 857 tfel::raise( 858 "InelasticFlowBase::getPorosityEffectOnEquivalentPlasticStrain:" 859 "uninitialised flow"); 860 } 861 const auto scpe = this->sc->getPorosityEffectOnEquivalentPlasticStrain(); 862 if ((scpe != StressCriterion::NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN) && 863 (scpe != 864 StressCriterion::STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN)) { 865 tfel::raise( 866 "InelasticFlowBase::getPorosityEffectOnEquivalentPlasticStrain:" 867 "unsupported porosity effect defined by the stress criterion"); 868 } 869 if (this->fc != nullptr) { 870 const auto fcpe = this->fc->getPorosityEffectOnEquivalentPlasticStrain(); 871 if ((fcpe != StressCriterion::NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN) && 872 (fcpe != 873 StressCriterion::STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN)) { 874 tfel::raise( 875 "InelasticFlowBase::getPorosityEffectOnEquivalentPlasticStrain:" 876 "unsupported porosity effect defined by the flow criterion"); 877 } 878 if ((scpe == 879 StressCriterion::STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN) || 880 (fcpe == 881 StressCriterion::STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN)) { 882 return STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN; 883 } 884 return NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN; 885 } else { 886 if (scpe == StressCriterion::NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN) { 887 return NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN; 888 } 889 return STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN; 890 } 891 } 892 return this->porosity_effect_on_equivalent_plastic_strain; 893 } // end of InelasticFlowBase::getPorosityEffectOnEquivalentPlasticStrain 894 updatePorosityUpperBound(const BehaviourDescription & bd,const std::string & id) const895 std::string InelasticFlowBase::updatePorosityUpperBound( 896 const BehaviourDescription& bd, const std::string& id) const { 897 if (this->sc == nullptr) { 898 tfel::raise( 899 "InelasticFlowBase::updatePorosityUpperBound: " 900 "uninitialized flow"); 901 } 902 return this->sc->updatePorosityUpperBound( 903 bd, id, StressCriterion::STRESSCRITERION); 904 } // end of InelasticFlowBase::updatePorosityUpperBound 905 906 InelasticFlowBase::~InelasticFlowBase() = default; 907 908 } // end of namespace bbrick 909 910 } // end of namespace mfront 911