1 /*! 2 * \file mfront/src/StandardElastoViscoPlasticityBrick.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 17/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/MFrontDebugMode.hxx" 19 #include "MFront/ImplicitDSLBase.hxx" 20 #include "MFront/NonLinearSystemSolver.hxx" 21 #include "MFront/BehaviourBrick/BrickUtilities.hxx" 22 #include "MFront/BehaviourBrick/StressPotential.hxx" 23 #include "MFront/BehaviourBrick/StressPotentialFactory.hxx" 24 #include "MFront/BehaviourBrick/OptionDescription.hxx" 25 #include "MFront/BehaviourBrick/InelasticFlow.hxx" 26 #include "MFront/BehaviourBrick/InelasticFlowFactory.hxx" 27 #include "MFront/BehaviourBrick/PorosityNucleationModel.hxx" 28 #include "MFront/BehaviourBrick/PorosityNucleationModelFactory.hxx" 29 #include "MFront/AbstractBehaviourDSL.hxx" 30 #include "MFront/StandardElastoViscoPlasticityBrick.hxx" 31 32 namespace mfront { 33 getId(const size_t i,const size_t m)34 static std::string getId(const size_t i, const size_t m) { 35 if (m == 1) { 36 return ""; 37 } 38 return std::to_string(i); 39 } // end of getId 40 declarePorosityUpperBoundSafetyFactorParameter(BehaviourDescription & bd,const double v)41 static void declarePorosityUpperBoundSafetyFactorParameter( 42 BehaviourDescription& bd, const double v) { 43 constexpr const auto uh = 44 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 45 const auto n = 46 StandardElastoViscoPlasticityBrick::porosityUpperBoundSafetyFactor; 47 VariableDescription p("real", n, 1u, 0u); 48 p.description = "a safety factor for the porosity upper bound"; 49 bd.addParameter(uh, p, BehaviourData::UNREGISTRED); 50 bd.setParameterDefaultValue(uh, n, v); 51 } // end of declarePorosityUpperBoundSafetyFactor 52 53 static void declarePorosityUpperBoundSafetyFactorForFractureDetectionParameter(BehaviourDescription & bd,const double v)54 declarePorosityUpperBoundSafetyFactorForFractureDetectionParameter( 55 BehaviourDescription& bd, const double v) { 56 constexpr const auto uh = 57 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 58 const auto n = StandardElastoViscoPlasticityBrick:: 59 porosityUpperBoundSafetyFactorForFractureDetection; 60 VariableDescription p("real", n, 1u, 0u); 61 p.description = "a safety factor for the porosity upper bound"; 62 bd.addParameter(uh, p, BehaviourData::UNREGISTRED); 63 bd.setParameterDefaultValue(uh, n, v); 64 } // end of declarePorosityUpperBoundSafetyFactorForFractureDetection 65 66 const char* const StandardElastoViscoPlasticityBrick:: 67 computeStandardSystemOfImplicitEquations = 68 "compute_standard_system_of_implicit_equations"; 69 const char* const StandardElastoViscoPlasticityBrick::brokenVariable = 70 "broken"; 71 const char* const StandardElastoViscoPlasticityBrick:: 72 currentEstimateOfThePorosityIncrement = 73 "current_estimate_of_the_porosity_increment"; 74 const char* const 75 StandardElastoViscoPlasticityBrick::nextEstimateOfThePorosityIncrement = 76 "next_estimate_of_the_porosity_increment"; 77 const char* const StandardElastoViscoPlasticityBrick::fixedPointConverged = 78 "fixed_point_converged"; 79 const char* const StandardElastoViscoPlasticityBrick:: 80 nextEstimateOfThePorosityAtTheEndOfTheTimeStep = 81 "next_estimate_of_the_porosity_at_the_end_of_the_time_step"; 82 const char* const StandardElastoViscoPlasticityBrick::porosityUpperBound = 83 "upper_bound_of_the_porosity"; 84 const char* const 85 StandardElastoViscoPlasticityBrick::porosityUpperBoundSafetyFactor = 86 "safety_factor_for_the_upper_bound_of_the_porosity"; 87 const char* const StandardElastoViscoPlasticityBrick:: 88 porosityUpperBoundSafetyFactorForFractureDetection = 89 "safety_factor_for_the_upper_bound_of_the_porosity_for_fracture_" 90 "detection"; 91 const char* const StandardElastoViscoPlasticityBrick:: 92 differenceBetweenSuccessiveEstimatesOfThePorosityIncrement = 93 "difference_between_successive_estimates_of_the_porosity_increment"; 94 const char* const 95 StandardElastoViscoPlasticityBrick::staggeredSchemeConvergenceCriterion = 96 "staggered_scheme_porosity_criterion"; 97 const char* const 98 StandardElastoViscoPlasticityBrick::staggeredSchemeIterationCounter = 99 "staggered_scheme_iteration_counter"; 100 const char* const StandardElastoViscoPlasticityBrick:: 101 maximumNumberOfIterationsOfTheStaggeredScheme = 102 "staggered_scheme_maximum_number_of_iterations"; 103 const char* const 104 StandardElastoViscoPlasticityBrick::staggeredSchemeBissectionAlgorithm = 105 "staggered_scheme_bissection_algorithm"; 106 const char* const StandardElastoViscoPlasticityBrick:: 107 staggeredSchemeAitkenAccelerationAlgorithm = 108 "staggered_scheme_aitken_acceleration_algorithm"; 109 StandardElastoViscoPlasticityBrick(AbstractBehaviourDSL & dsl_,BehaviourDescription & mb_)110 StandardElastoViscoPlasticityBrick::StandardElastoViscoPlasticityBrick( 111 AbstractBehaviourDSL& dsl_, BehaviourDescription& mb_) 112 : BehaviourBrickBase(dsl_, mb_) {} // end of 113 // StandardElastoViscoPlasticityBrick::StandardElastoViscoPlasticityBrick 114 getName() const115 std::string StandardElastoViscoPlasticityBrick::getName() const { 116 return "ElastoViscoPlasticity"; 117 } // end of StandardElastoViscoPlasticityBrick::getName 118 getDescription() const119 BehaviourBrickDescription StandardElastoViscoPlasticityBrick::getDescription() 120 const { 121 auto d = BehaviourBrickDescription{}; 122 d.behaviourType = 123 tfel::material::MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR; 124 d.integrationScheme = IntegrationScheme::IMPLICITSCHEME; 125 d.supportedModellingHypotheses = 126 ModellingHypothesis::getModellingHypotheses(); 127 d.supportedBehaviourSymmetries = {mfront::ISOTROPIC, mfront::ORTHOTROPIC}; 128 return d; 129 } // end of StandardElastoViscoPlasticityBrick::getDescription 130 131 std::vector<bbrick::OptionDescription> getOptions(const bool) const132 StandardElastoViscoPlasticityBrick::getOptions(const bool) const { 133 auto opts = std::vector<bbrick::OptionDescription>{}; 134 opts.emplace_back("stress_potential", 135 "Decare the stress potential (required)", 136 bbrick::OptionDescription::STRING); 137 opts.emplace_back("inelastic_flow", "Declare another inelastic flow", 138 bbrick::OptionDescription::STRING); 139 opts.emplace_back( 140 "porosity_evolution", 141 "state if the porosity evolution must be taken into account", 142 bbrick::OptionDescription::STRING); 143 return opts; 144 } // end of StandardElastoViscoPlasticityBrick::getOptions 145 initialize(const Parameters &,const DataMap & d)146 void StandardElastoViscoPlasticityBrick::initialize(const Parameters&, 147 const DataMap& d) { 148 auto raise = [](const std::string& m) { 149 tfel::raise("StandardElastoViscoPlasticityBrick::initialize: " + m); 150 }; // end of raise 151 auto& iff = bbrick::InelasticFlowFactory::getFactory(); 152 auto getDataStructure = [&raise](const std::string& n, const Data& ds) { 153 if (ds.is<std::string>()) { 154 tfel::utilities::DataStructure nds; 155 nds.name = ds.get<std::string>(); 156 return nds; 157 } 158 if (!ds.is<tfel::utilities::DataStructure>()) { 159 raise("invalid data type for entry '" + n + "'"); 160 } 161 return ds.get<tfel::utilities::DataStructure>(); 162 }; // end of getDataStructure 163 auto getStressPotential = [&d, &getDataStructure, &raise, 164 this](const char* const n) { 165 if (d.count(n) != 0) { 166 const auto ds = getDataStructure(n, d.at(n)); 167 if (this->stress_potential != nullptr) { 168 raise("the stress potential has already been defined"); 169 } 170 auto& spf = bbrick::StressPotentialFactory::getFactory(); 171 this->stress_potential = spf.generate(ds.name); 172 this->stress_potential->initialize(this->bd, this->dsl, ds.data); 173 } 174 }; 175 // 176 this->checkOptionsNames(d); 177 // 178 getStressPotential("elastic_potential"); 179 getStressPotential("stress_potential"); 180 if (this->stress_potential == nullptr) { 181 raise("no stress potential defined"); 182 } 183 auto save_individual_porosity_increase = false; 184 if (d.count("porosity_evolution") != 0) { 185 save_individual_porosity_increase = 186 this->treatPorosityEvolutionSection(d.at("porosity_evolution")); 187 } 188 for (const auto& e : d) { 189 if ((e.first == "elastic_potential") || (e.first == "stress_potential")) { 190 // already treated 191 } else if (e.first == "inelastic_flow") { 192 auto append_flow = [this, &iff, getDataStructure, raise, 193 &save_individual_porosity_increase]( 194 const Data& ifd, const size_t msize) { 195 const auto ds = getDataStructure("inelatic_flow", ifd); 196 auto iflow = iff.generate(ds.name); 197 auto data = ds.data; 198 if (data.count("save_porosity_increase") == 0) { 199 data["save_porosity_increase"] = save_individual_porosity_increase; 200 } 201 if (data.count("porosity_evolution_algorithm") != 0) { 202 raise( 203 "the `porosity_evolution_algorithm` entry is reserved and " 204 "shall not be set by the user"); 205 } 206 if (this->porosity_evolution_algorithm == 207 mfront::bbrick::PorosityEvolutionAlgorithm:: 208 STANDARD_IMPLICIT_SCHEME) { 209 data["porosity_evolution_algorithm"] = 210 std::string("standard_implicit_scheme"); 211 } else if (this->porosity_evolution_algorithm == 212 mfront::bbrick::PorosityEvolutionAlgorithm:: 213 STAGGERED_SCHEME) { 214 data["porosity_evolution_algorithm"] = 215 std::string("staggered_scheme"); 216 } else { 217 raise("internal error (unsupported porosity algorithm)"); 218 } 219 iflow->initialize(this->bd, this->dsl, 220 getId(this->flows.size(), msize), data); 221 this->flows.push_back(iflow); 222 }; 223 if (e.second.is<std::vector<Data>>()) { 224 // multiple inelastic flows are defined 225 const auto& ifs = e.second.get<std::vector<Data>>(); 226 for (const auto& iflow : ifs) { 227 append_flow(iflow, ifs.size()); 228 } 229 } else { 230 append_flow(e.second, 1u); 231 } 232 } else if (e.first == "porosity_evolution") { 233 } else { 234 raise("unsupported entry '" + e.first + "'"); 235 } 236 } 237 // say to every flows if the porosity is handled by the brick 238 const auto pe = this->isCoupledWithPorosityEvolution(); 239 for (auto& f : this->flows) { 240 f->setPorosityEvolutionHandled(pe); 241 } 242 // 243 if (this->porosity_evolution_algorithm == 244 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 245 this->bd.appendToIncludes( 246 "#include \"TFEL/Math/RootFinding/BissectionAlgorithmBase.hxx\"\n"); 247 mfront::bbrick::addLocalVariable( 248 bd, "tfel::math::BissectionAlgorithmBase<real>", 249 staggeredSchemeBissectionAlgorithm); 250 if (this->staggered_scheme_parameters.acceleration_algorithm == 251 StaggeredSchemeParameters::AITKEN) { 252 this->bd.appendToIncludes( 253 "#include " 254 "\"TFEL/Math/AccelerationAlgorithms/" 255 "AitkenAccelerationAlgorithm.hxx\"\n"); 256 mfront::bbrick::addLocalVariable( 257 bd, "tfel::math::AitkenAccelerationAlgorithm<real>", 258 staggeredSchemeAitkenAccelerationAlgorithm); 259 } 260 } 261 } // end of StandardElastoViscoPlasticityBrick::initialize 262 treatPorosityEvolutionSection(const Data & e)263 bool StandardElastoViscoPlasticityBrick::treatPorosityEvolutionSection( 264 const Data& e) { 265 auto save_individual_porosity_increase = false; 266 auto& nmf = bbrick::PorosityNucleationModelFactory::getFactory(); 267 auto raise = [](const std::string& m) { 268 tfel::raise( 269 "StandardElastoViscoPlasticityBrick::" 270 "treatPorosityEvolutionSection: " + 271 m); 272 }; // end of raise 273 auto getDataStructure = [&raise](const std::string& n, const Data& ds) { 274 if (ds.is<std::string>()) { 275 tfel::utilities::DataStructure nds; 276 nds.name = ds.get<std::string>(); 277 return nds; 278 } 279 if (!ds.is<tfel::utilities::DataStructure>()) { 280 raise("invalid data type for entry '" + n + "'"); 281 } 282 return ds.get<tfel::utilities::DataStructure>(); 283 }; // end of getDataStructure 284 if (!e.is<DataMap>()) { 285 raise("invalid data type for entry 'porosity_evolution'"); 286 } 287 const auto ed = e.get<DataMap>(); 288 if (ed.count("save_individual_porosity_increase") != 0) { 289 const auto b = ed.at("save_individual_porosity_increase"); 290 if (!b.is<bool>()) { 291 raise( 292 "the 'save_individual_porosity_increase' option is not a boolean " 293 "value"); 294 } 295 save_individual_porosity_increase = b.get<bool>(); 296 } 297 for (const auto& ped : ed) { 298 if (ped.first == "save_individual_porosity_increase") { 299 // already treated 300 } else if (ped.first == "elastic_contribution") { 301 const auto b = ed.at("elastic_contribution"); 302 if (!b.is<bool>()) { 303 raise("the 'elastic_contribution' option is not a boolean value"); 304 } 305 this->elastic_contribution = b.get<bool>(); 306 } else if (ped.first == "nucleation_model") { 307 auto append_nucleation_model = [this, &nmf, getDataStructure, raise, 308 &save_individual_porosity_increase]( 309 const Data& nmd, const size_t msize) { 310 const auto ds = getDataStructure("nucleation_model", nmd); 311 auto nm = nmf.generate(ds.name); 312 auto data = ds.data; 313 if (data.count("save_individual_porosity_increase") == 0) { 314 data["save_individual_porosity_increase"] = 315 save_individual_porosity_increase; 316 } 317 if (data.count("porosity_evolution_algorithm") != 0) { 318 raise( 319 "the `porosity_evolution_algorithm` entry is reserved and " 320 "shall not be set by the user"); 321 } 322 if (porosity_evolution_algorithm == 323 mfront::bbrick::PorosityEvolutionAlgorithm:: 324 STANDARD_IMPLICIT_SCHEME) { 325 data["porosity_evolution_algorithm"] = 326 std::string("standard_implicit_scheme"); 327 } else if (porosity_evolution_algorithm == 328 mfront::bbrick::PorosityEvolutionAlgorithm:: 329 STAGGERED_SCHEME) { 330 data["porosity_evolution_algorithm"] = 331 std::string("staggered_scheme"); 332 } else { 333 raise("internal error (unsupported porosity algorithm)"); 334 } 335 nm->initialize(this->bd, this->dsl, 336 getId(this->nucleation_models.size(), msize), data); 337 this->nucleation_models.push_back(nm); 338 }; 339 if (ped.second.is<std::vector<Data>>()) { 340 // multiple inelastic nucleation_models are defined 341 const auto& nms = ped.second.get<std::vector<Data>>(); 342 for (const auto& nm : nms) { 343 append_nucleation_model(nm, nms.size()); 344 } 345 } else { 346 append_nucleation_model(ped.second, 1u); 347 } 348 } else if (ped.first == "algorithm") { 349 this->treatPorosityEvolutionAlgorithmSection(ped.second); 350 } else if (ped.first == porosityUpperBoundSafetyFactor) { 351 const auto value = tfel::utilities::convert<double>(ped.second); 352 declarePorosityUpperBoundSafetyFactorParameter(bd, value); 353 } else if (ped.first == 354 porosityUpperBoundSafetyFactorForFractureDetection) { 355 const auto value = tfel::utilities::convert<double>(ped.second); 356 declarePorosityUpperBoundSafetyFactorForFractureDetectionParameter( 357 bd, value); 358 } else { 359 raise("invalid entry '" + ped.first + "' in 'porosity_evolution'"); 360 } 361 } // end of the loop over the entry of the porosity_evolution options 362 // Since a porosity_evolution section has been defined, 363 // we consider that all the non purely deviatoric inelastic flow 364 // contributes to the porosity evolution 365 if (this->porosity_growth_policy == UNDEFINEDPOROSITYGROWTHPOLICY) { 366 this->porosity_growth_policy = STANDARDVISCOPLASTICPOROSITYGROWTHPOLICY; 367 } 368 return save_individual_porosity_increase; 369 } // end of 370 // StandardElastoViscoPlasticityBrick::treatPorosityEvolutionSection 371 372 void treatPorosityEvolutionAlgorithmSection(const Data & d)373 StandardElastoViscoPlasticityBrick::treatPorosityEvolutionAlgorithmSection( 374 const Data& d) { 375 auto raise = [](const std::string& m) { 376 tfel::raise( 377 "StandardElastoViscoPlasticityBrick::" 378 "treatPorosityEvolutionAlgorithmSection: " + 379 m); 380 }; // end of raise 381 auto get_algorithm = [&raise](const std::string& a) { 382 if ((a == "standard implicit scheme") || 383 (a == "standard_implicit_scheme") || 384 (a == "StandardImplicitScheme")) { 385 return mfront::bbrick::PorosityEvolutionAlgorithm:: 386 STANDARD_IMPLICIT_SCHEME; 387 } else if (!((a == "staggered scheme") || (a == "staggered_scheme") || 388 (a == "StaggeredScheme"))) { 389 raise( 390 "error while treating entry 'algorithm' in " 391 "'porosity_evolution': unhandled algorithm '" + 392 a + 393 "'. Known algorithms are:\n" 394 "- 'standard implicit scheme' (or equivalently " 395 "'standard_implicit_scheme' or 'StandardImplicitScheme')\n" 396 "- 'staggered scheme' (or equivalently 'staggered_scheme' or " 397 "'StaggeredScheme'"); 398 } 399 return mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME; 400 }; 401 if (d.is<std::string>()) { 402 this->porosity_evolution_algorithm = get_algorithm(d.get<std::string>()); 403 } else if (d.is<tfel::utilities::DataStructure>()) { 404 const auto& a = d.get<tfel::utilities::DataStructure>(); 405 const auto& algorithm_parameters = a.data; 406 this->porosity_evolution_algorithm = get_algorithm(a.name); 407 if (this->porosity_evolution_algorithm == 408 mfront::bbrick::PorosityEvolutionAlgorithm:: 409 STANDARD_IMPLICIT_SCHEME) { 410 if (!algorithm_parameters.empty()) { 411 raise( 412 "no algorithm parameter for the 'standard implicit " 413 "scheme'"); 414 } 415 } else if (this->porosity_evolution_algorithm == 416 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 417 this->treatStaggeredPorosityEvolutionAlgorithmParameters( 418 algorithm_parameters); 419 } else { 420 raise("internal error (unsupported porosity evolution algorithm)"); 421 } 422 } else { 423 raise("invalid type for entry 'algorithm' in 'porosity_evolution'"); 424 } 425 } // end of 426 // StandardElastoViscoPlasticityBrick::treatPorosityEvolutionAlgorithmSection 427 428 void StandardElastoViscoPlasticityBrick:: treatStaggeredPorosityEvolutionAlgorithmParameters(const DataMap & algorithm_parameters)429 treatStaggeredPorosityEvolutionAlgorithmParameters( 430 const DataMap& algorithm_parameters) { 431 auto raise = [](const std::string& m) { 432 tfel::raise( 433 "StandardElastoViscoPlasticityBrick::" 434 "treatStaggeredPorosityEvolutionAlgorithmParameters: " + 435 m); 436 }; // end of raise 437 for (const auto& p : algorithm_parameters) { 438 if (p.first == "convergence_criterion") { 439 if (!p.second.is<double>()) { 440 raise("invalid type for the `convergence_criterion` parameter"); 441 } 442 const auto& c = p.second.get<double>(); 443 if (c <= 0) { 444 raise("invalid value of the `convergence_criterion` parameter"); 445 } 446 this->staggered_scheme_parameters.convergence_criterion = c; 447 } else if (p.first == "maximum_number_of_iterations") { 448 if (!p.second.is<int>()) { 449 raise( 450 "invalid type for the `maximum_number_of_iterations` " 451 "parameter"); 452 } 453 const auto& i = p.second.get<int>(); 454 if (i < 1) { 455 raise( 456 "invalid value of the `maximum_number_of_iterations` " 457 "parameter"); 458 } 459 this->staggered_scheme_parameters.maximum_number_of_iterations = 460 static_cast<unsigned int>(i); 461 } else if (p.first == "acceleration_algorithm"){ 462 if (!p.second.is<std::string>()) { 463 raise( 464 "invalid type for the `acceleration_algorithm`" 465 "parameter"); 466 } 467 const auto& a = p.second.get<std::string>(); 468 if (a == "Aitken") { 469 this->staggered_scheme_parameters.acceleration_algorithm = 470 StaggeredSchemeParameters::AITKEN; 471 } else if ((a == "Relaxation") || (a == "relaxation")) { 472 this->staggered_scheme_parameters.acceleration_algorithm = 473 StaggeredSchemeParameters::RELAXATION; 474 } else { 475 raise("unsupported acceleration algorithm '" + a + 476 "' (expected 'Aitken' or 'Relaxation')"); 477 } 478 } else { 479 raise("invalid parameter '" + p.first + "'"); 480 } 481 } 482 } // end of treatStaggeredPorosityEvolutionAlgorithmParameters 483 isCoupledWithPorosityEvolution() const484 bool StandardElastoViscoPlasticityBrick::isCoupledWithPorosityEvolution() 485 const { 486 for (const auto& f : this->flows) { 487 if (f->isCoupledWithPorosityEvolution()) { 488 return true; 489 } 490 } 491 if (!this->nucleation_models.empty()) { 492 return true; 493 } 494 // no flow coupled with porosity evolution. However, if 495 // `porosity_growth_policy` is not equal to 496 // `UNDEFINEDPOROSITYGROWTHPOLICY`, 497 // the user has defined an `porosity_evolution` section, which means that 498 // he wants that any non purely deviatoric inelastic flow contributes to 499 // the 500 // porosity growth 501 return this->porosity_growth_policy != UNDEFINEDPOROSITYGROWTHPOLICY; 502 } // end if 503 // StandardElastoViscoPlasticityBrick::isCoupledWithPorosityEvolution 504 505 std::vector<StandardElastoViscoPlasticityBrick::Hypothesis> getSupportedModellingHypotheses() const506 StandardElastoViscoPlasticityBrick::getSupportedModellingHypotheses() const { 507 return this->stress_potential->getSupportedModellingHypotheses(this->bd, 508 this->dsl); 509 } 510 511 std::map<std::string, std::shared_ptr<mfront::bbrick::InelasticFlow>> buildInelasticFlowsMap() const512 StandardElastoViscoPlasticityBrick::buildInelasticFlowsMap() const { 513 auto i = size_t{}; 514 auto m = std::map<std::string, std::shared_ptr<bbrick::InelasticFlow>>{}; 515 for (const auto& f : this->flows) { 516 m.insert({getId(i, this->flows.size()), f}); 517 } 518 return m; 519 } // end of StandardElastoViscoPlasticityBrick::buildInelasticFlowsMap() 520 completeVariableDeclaration() const521 void StandardElastoViscoPlasticityBrick::completeVariableDeclaration() const { 522 constexpr const auto uh = 523 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 524 if (this->isCoupledWithPorosityEvolution()) { 525 mfront::bbrick::addStateVariableIfNotDefined( 526 bd, "real", "f", tfel::glossary::Glossary::Porosity, 1u, true); 527 const auto& f = 528 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 529 tfel::glossary::Glossary::Porosity); 530 if (!f.hasPhysicalBounds()) { 531 VariableBoundsDescription fbounds; 532 fbounds.boundsType = VariableBoundsDescription::LOWERANDUPPER; 533 fbounds.lowerBound = 0; 534 fbounds.upperBound = 1; 535 this->bd.setPhysicalBounds(uh, f.name, fbounds); 536 } 537 this->bd.reserveName(uh, f.name + "_"); 538 // 539 mfront::bbrick::addAuxiliaryStateVariableIfNotDefined( 540 bd, "real", brokenVariable, tfel::glossary::Glossary::Broken, 1u, 541 true); 542 } 543 this->stress_potential->completeVariableDeclaration(this->bd, this->dsl); 544 auto i = size_t{}; 545 for (const auto& f : this->flows) { 546 f->completeVariableDeclaration(this->bd, this->dsl, 547 getId(i, this->flows.size())); 548 ++i; 549 } 550 i = size_t{}; 551 for (const auto& nm : this->nucleation_models) { 552 nm->completeVariableDeclaration(this->bd, this->dsl, 553 this->buildInelasticFlowsMap(), 554 getId(i, this->nucleation_models.size())); 555 ++i; 556 } 557 // automatic registration of the porosity for porous flows 558 if (this->isCoupledWithPorosityEvolution()) { 559 // 560 if (!bd.hasParameter(uh, porosityUpperBoundSafetyFactor)) { 561 declarePorosityUpperBoundSafetyFactorParameter(bd, 0.985); 562 } 563 // 564 if (!bd.hasParameter( 565 uh, porosityUpperBoundSafetyFactorForFractureDetection)) { 566 declarePorosityUpperBoundSafetyFactorForFractureDetectionParameter( 567 bd, 0.984); 568 } 569 mfront::bbrick::addLocalVariable(bd, "real", porosityUpperBound); 570 // 571 if (this->porosity_evolution_algorithm == 572 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 573 mfront::bbrick::addLocalVariable( 574 bd, "bool", computeStandardSystemOfImplicitEquations); 575 mfront::bbrick::addLocalVariable(bd, "real", 576 currentEstimateOfThePorosityIncrement); 577 mfront::bbrick::addLocalVariable(bd, "ushort", 578 staggeredSchemeIterationCounter); 579 // convergence criterion of the staggered scheme 580 VariableDescription eps("real", staggeredSchemeConvergenceCriterion, 1u, 581 0u); 582 eps.description = "stopping criterion value of the staggered scheme"; 583 bd.addParameter(uh, eps, BehaviourData::UNREGISTRED); 584 bd.setParameterDefaultValue( 585 uh, staggeredSchemeConvergenceCriterion, 586 this->staggered_scheme_parameters.convergence_criterion); 587 // maximum number of iterations of the staggered scheme 588 VariableDescription iterMax( 589 "ushort", maximumNumberOfIterationsOfTheStaggeredScheme, 1u, 0u); 590 iterMax.description = 591 "maximum number of iterations of the staggered scheme allowed"; 592 bd.addParameter(uh, iterMax, BehaviourData::UNREGISTRED); 593 bd.setParameterDefaultValue( 594 uh, maximumNumberOfIterationsOfTheStaggeredScheme, 595 static_cast<unsigned short>(this->staggered_scheme_parameters 596 .maximum_number_of_iterations)); 597 } 598 } 599 } // end of StandardElastoViscoPlasticityBrick::completeVariableDeclaration 600 endTreatment() const601 void StandardElastoViscoPlasticityBrick::endTreatment() const { 602 constexpr const auto uh = 603 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 604 // 605 if (this->isCoupledWithPorosityEvolution()) { 606 if (this->porosity_evolution_algorithm == 607 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 608 const auto& broken = 609 bd.getBehaviourData(uh) 610 .getAuxiliaryStateVariableDescriptionByExternalName( 611 tfel::glossary::Glossary::Broken); 612 CodeBlock init; 613 // initialize the porosity status 614 init.code += "this->"; 615 init.code += StandardElastoViscoPlasticityBrick:: 616 computeStandardSystemOfImplicitEquations; 617 init.code += " = 2 * (this->" + broken.name + ") > 1;\n"; 618 // initialize the porosity increment 619 init.code += StandardElastoViscoPlasticityBrick:: 620 currentEstimateOfThePorosityIncrement; 621 init.code += " = real{};\n"; 622 // initialize the iteration counter of the staggered scheme 623 init.code += "this->"; 624 init.code += 625 StandardElastoViscoPlasticityBrick::staggeredSchemeIterationCounter; 626 init.code += " = static_cast<unsigned short>(0);\n"; 627 if ((this->porosity_evolution_algorithm == 628 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) && 629 (this->staggered_scheme_parameters.acceleration_algorithm == 630 StaggeredSchemeParameters::AITKEN)) { 631 init.code += staggeredSchemeAitkenAccelerationAlgorithm; 632 init.code += ".initialize(real(0));\n"; 633 } 634 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, init, 635 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 636 } 637 const auto& f = 638 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 639 tfel::glossary::Glossary::Porosity); 640 const auto f_ = f.name + "_"; 641 // value of the porosity at t+theta*dt 642 CodeBlock ib; 643 ib.code += "const auto " + f_ + " = "; 644 ib.code += "std::max(std::min(this->" + f.name + " + (" + 645 this->bd.getClassName() + "::theta) * (this->d" + f.name + 646 "), (this->" + std::string(porosityUpperBoundSafetyFactor) + 647 ") * " + std::string(porosityUpperBound) + "), real(0));\n "; 648 ib.code += "static_cast<void>(" + f_ + ");\n"; 649 if (this->elastic_contribution) { 650 if (this->porosity_evolution_algorithm == 651 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 652 ib.code += "if("; 653 ib.code += StandardElastoViscoPlasticityBrick:: 654 computeStandardSystemOfImplicitEquations; 655 ib.code += "){\n"; 656 this->addElasticContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 657 ib); 658 ib.code += "}\n"; 659 } else if (this->porosity_evolution_algorithm == 660 mfront::bbrick::PorosityEvolutionAlgorithm:: 661 STANDARD_IMPLICIT_SCHEME) { 662 this->addElasticContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 663 ib); 664 } else { 665 tfel::raise( 666 "StandardElastoViscoPlasticityBrick::endTreatment: internal " 667 "error " 668 "(unsupported porosity evolution algorithm)"); 669 } 670 } 671 bd.setCode(uh, BehaviourData::Integrator, ib, 672 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 673 // implicit equation associated with the porosity 674 if (this->porosity_evolution_algorithm == 675 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 676 ib.code = "if(!this->"; 677 ib.code += StandardElastoViscoPlasticityBrick:: 678 computeStandardSystemOfImplicitEquations; 679 ib.code += "){\n"; 680 ib.code += "f" + f.name + " -= " + 681 StandardElastoViscoPlasticityBrick:: 682 currentEstimateOfThePorosityIncrement + 683 ";\n"; 684 ib.code += "}\n"; 685 bd.setCode(uh, BehaviourData::Integrator, ib, 686 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 687 } 688 // modification of the implicit equation to impose the bounds 689 if (this->porosity_evolution_algorithm == 690 mfront::bbrick::PorosityEvolutionAlgorithm:: 691 STANDARD_IMPLICIT_SCHEME) { 692 const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(this->dsl); 693 const auto requiresAnalyticalJacobian = 694 ((idsl.getSolver().usesJacobian()) && 695 (!idsl.getSolver().requiresNumericalJacobian())); 696 ib.code = "if(this->" + f.name + " + this->d" + f.name + // 697 " - f" + f.name + " < 0){\n"; 698 ib.code += "f" + f.name + " = " + // 699 "this->d" + f.name + " - this->" + f.name + ";\n"; 700 if (requiresAnalyticalJacobian) { 701 ib.code += 702 "for(unsigned short idx = 0; idx!= " // 703 "this->jacobian.getNbCols(); ++idx){\n"; 704 ib.code += "this->jacobian(" + f.name + "_offset, idx) = 0;\n"; 705 ib.code += "}\n"; 706 ib.code += "this->jacobian(" + f.name + "_offset, " + f.name + 707 "_offset) = 1;\n"; 708 } 709 ib.code += "}\n"; 710 ib.code += "if(this->" + f.name + " + this->d" + f.name + // 711 " - f" + f.name + " > " + porosityUpperBound + "){\n"; 712 ib.code += "f" + f.name + " = this->d" + f.name + // 713 " - (" + porosityUpperBound + " - this->" + f.name + ");\n"; 714 if (requiresAnalyticalJacobian) { 715 ib.code += 716 "for(unsigned short idx = 0; idx!= " // 717 "this->jacobian.getNbCols(); ++idx){\n"; 718 ib.code += "this->jacobian(" + f.name + "_offset, idx) = 0;\n"; 719 ib.code += "}\n"; 720 ib.code += "this->jacobian(" + f.name + "_offset, " + f.name + 721 "_offset) = 1;\n"; 722 } 723 ib.code += "}\n"; 724 bd.setCode(uh, BehaviourData::Integrator, ib, 725 BehaviourData::CREATEORAPPEND, BehaviourData::AT_END); 726 } 727 // additional convergence checks 728 if (this->porosity_evolution_algorithm == 729 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 730 const auto& broken = 731 bd.getBehaviourData(uh) 732 .getAuxiliaryStateVariableDescriptionByExternalName( 733 tfel::glossary::Glossary::Broken); 734 // additional convergence checks 735 CodeBlock acc; 736 std::ostringstream os; 737 // We check if we did not built the implicit system only to have the 738 // jacobian of the full implicit system after the convergence of the 739 // staggered scheme. If this is the case, the converged flag is set 740 // to true and the the resolution is stopped 741 os << "if (("; 742 os << computeStandardSystemOfImplicitEquations; 743 os << ") || (2 * (this->" + broken.name + ") > 1)) {\n"; 744 os << "converged = true;\n"; 745 os << "return;\n"; 746 os << "}\n"; 747 acc.code = os.str(); 748 bd.setCode(uh, BehaviourData::AdditionalConvergenceChecks, acc, 749 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 750 } 751 } 752 this->stress_potential->endTreatment(this->bd, this->dsl); 753 auto i = size_t{}; 754 for (const auto& f : this->flows) { 755 f->endTreatment(this->bd, this->dsl, *(this->stress_potential), 756 getId(i, this->flows.size())); 757 ++i; 758 } 759 i = size_t{}; 760 for (const auto& nm : this->nucleation_models) { 761 nm->endTreatment(this->bd, this->dsl, *(this->stress_potential), 762 this->buildInelasticFlowsMap(), 763 getId(i, this->nucleation_models.size())); 764 ++i; 765 } 766 // At this stage, one expects to be able to compute the upper porosity bound 767 if (this->isCoupledWithPorosityEvolution()) { 768 CodeBlock init; 769 init.code += "this->"; 770 init.code += porosityUpperBound; 771 init.code += " = real{1};\n"; 772 auto flow_id = size_t{}; 773 for (const auto& pf : this->flows) { 774 const auto id = getId(flow_id, this->flows.size()); 775 init.code += pf->updatePorosityUpperBound(bd, id); 776 ++flow_id; 777 } 778 // init.code += "if(2 * (this->" + broken.name + ") > 1){\n"; 779 // init.code += "this->" + f.name + " = (this->" + 780 // std::string(porosityUpperBoundSafetyFactor) + ") * 781 // (this->" + 782 // std::string(porosityUpperBound) + ");\n"; 783 // init.code += "}\n"; 784 // if (this->porosity_evolution_algorithm == 785 // mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME) 786 // { 787 // init.code += porosityIncrementLowerBound + " = -this->f;\n"; 788 // init.code += porosityIncrementUpperBound + " = "; 789 // init.code += porosityUpperBound + "- (this->f);\n"; 790 // } 791 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, init, 792 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 793 } 794 // at this stage, one assumes that the various components of the inelastic 795 // flow (stress_potential, isotropic hardening rule) have added the 796 // initialization of their material properties the 797 // `BeforeInitializeLocalVariables`. We then ask the inelastic flows if 798 // they 799 // require an activation state (in practice, it mean that an isotropic 800 // hardening rule has been defined). If so, the initialization of the 801 // activation states requires the the computation of an elastic prediction 802 // of the stress. The brik asks the stress potential to compute it in a 803 // variable called sigel and the inelastic flows shall use it to compute 804 // their initial state. All thoses steps must be added to the 805 // `BeforeInitializeLocalVariables` code block. 806 const bool bep = [this] { 807 for (const auto& pf : this->flows) { 808 if (pf->requiresActivationState()) { 809 return true; 810 } 811 } 812 return false; 813 }(); 814 if (bep) { 815 // compute the elastic prediction 816 this->stress_potential->computeElasticPrediction(bd); 817 i = size_t{}; 818 for (const auto& pf : this->flows) { 819 if (pf->requiresActivationState()) { 820 pf->computeInitialActivationState(bd, *(this->stress_potential), 821 getId(i, this->flows.size())); 822 } 823 ++i; 824 } 825 } 826 if ((this->isCoupledWithPorosityEvolution()) && 827 (this->porosity_evolution_algorithm == 828 mfront::bbrick::PorosityEvolutionAlgorithm::STAGGERED_SCHEME)) { 829 const auto& f = 830 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 831 tfel::glossary::Glossary::Porosity); 832 const auto& broken = 833 bd.getBehaviourData(uh) 834 .getAuxiliaryStateVariableDescriptionByExternalName( 835 tfel::glossary::Glossary::Broken); 836 const auto df = 837 differenceBetweenSuccessiveEstimatesOfThePorosityIncrement; 838 const auto f_ets = nextEstimateOfThePorosityAtTheEndOfTheTimeStep; 839 // additional convergence checks 840 CodeBlock acc; 841 // defining the `nextEstimateOfThePorosityIncrement` variable before the 842 // user code, so that the user can update the porosity with its own 843 // contribution 844 acc.code = "auto " + std::string(fixedPointConverged) + " = true;\n"; 845 acc.code += "auto " + std::string(nextEstimateOfThePorosityIncrement) + 846 " = real{};\n"; 847 bd.setCode(uh, BehaviourData::AdditionalConvergenceChecks, acc, 848 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 849 acc.code.clear(); 850 std::ostringstream os; 851 auto generate_debug_message = [&os, &f, &df](const bool b) { 852 if (!getDebugMode()) { 853 return; 854 } 855 if (b) { 856 os << "std::cout " 857 << "<< \"convergence of the staggered algorithm, iteration \""; 858 } else { 859 os << "std::cout " 860 << "<< \"non convergence of the staggered algorithm:, iteration " 861 "\""; 862 } 863 os << " << this->" << staggeredSchemeIterationCounter << " << '\\n'\n" 864 << "<< \"- current estimate of the porosity at the end of the time " 865 << "step: \"" 866 << " << ((this->" << f.name << ")+(this->" 867 << currentEstimateOfThePorosityIncrement << ")) << '\\n'\n" 868 << "<< \"- current estimate of the porosity increment: \" << " 869 "(this->" 870 << currentEstimateOfThePorosityIncrement << ") << '\\n'\n" 871 << "<< \"- next estimate of the porosity increment: \" << " 872 << nextEstimateOfThePorosityIncrement << " << '\\n'\n" 873 << "<< \"- difference: \" << " 874 << "std::abs(" << df << ") << '\\n';\n"; 875 }; 876 // now we check that the staggered scheme has convergenced 877 os << "if (converged && (!"; 878 os << computeStandardSystemOfImplicitEquations; 879 os << ")) {\n"; 880 if (this->elastic_contribution) { 881 const auto f_ = "((this->" + f.name + ") + (" + 882 this->bd.getClassName() + "::theta) * (this->" + 883 std::string(StandardElastoViscoPlasticityBrick:: 884 currentEstimateOfThePorosityIncrement) + 885 "))"; 886 os << nextEstimateOfThePorosityIncrement 887 << " += (1-" + f_ + ") * trace(this->deel);\n"; 888 } 889 // 1. compute a new estimate of the porosity increment and get the maximum 890 // allowed porosity (onset of fracture) 891 auto flow_id = size_t{}; 892 for (const auto& pf : this->flows) { 893 const auto id = getId(flow_id, this->flows.size()); 894 os << pf->updateNextEstimateOfThePorosityIncrement(bd, id); 895 ++flow_id; 896 } 897 auto nucleation_id = size_t{}; 898 for (const auto& nm : this->nucleation_models) { 899 const auto id = getId(nucleation_id, this->nucleation_models.size()); 900 os << nm->updateNextEstimateOfThePorosityIncrement( 901 bd, this->buildInelasticFlowsMap(), id); 902 ++nucleation_id; 903 } 904 auto apply_next_estimate_bounds = [&os, f_ets] { 905 // next estimate of the porosity at the end of the time step 906 os << "{\n" 907 << "const auto " << f_ets << " = this->f + " 908 << nextEstimateOfThePorosityIncrement << ";\n"; 909 // Treat the case when the newly computed porosity exceeds the 910 // porosity upper bounds. 911 // 912 // The maximum allowed increment is the current value of the upper of 913 // the porosity minus the value at the beginning of the time step. We 914 // then use a dichotomic approach by choosing an increment which is an 915 // average of the previous increment and this maximum allowed increment 916 os << "if(" << f_ets << " > (this->" << porosityUpperBoundSafetyFactor 917 << " ) * ( this->" << porosityUpperBound << ")){\n" 918 << nextEstimateOfThePorosityIncrement << " = " 919 << "(this->" << currentEstimateOfThePorosityIncrement << " + (this->" 920 << porosityUpperBound << " - this->f))/2;\n" 921 // Treat the case when the newly computed porosity is negative 922 // 923 // The maximum allowed increment is the opposite of the initial 924 // porosity. We then use a dichotomic approach by choosing an 925 // increment 926 // which is an average of the previous increment and this maximum 927 // allowed increment 928 << "} else if(" << f_ets << " < 0){\n" 929 << nextEstimateOfThePorosityIncrement << " = " 930 << "(" << currentEstimateOfThePorosityIncrement 931 << " - (this->f))/2;\n" 932 << "}\n" 933 << "}\n"; 934 }; 935 // 2. apply the bounds on the next estimate 936 apply_next_estimate_bounds(); 937 // 3. check if the staggered algorithm has converged 938 os << "const auto " << df << " = " << nextEstimateOfThePorosityIncrement 939 << " - (this->" << currentEstimateOfThePorosityIncrement << ");\n"; 940 os << std::string(fixedPointConverged) << " = (" 941 << std::string(fixedPointConverged) << ") && " 942 << "(std::abs(" << df << ") < this->" 943 << staggeredSchemeConvergenceCriterion << ");\n"; 944 os << "if(" + std::string(fixedPointConverged) + "){\n"; 945 generate_debug_message(true); 946 // try to dectect failure 947 os << "if (this->f + this->" << currentEstimateOfThePorosityIncrement 948 << " > (this->" << porosityUpperBoundSafetyFactorForFractureDetection 949 << ") * (this->" << porosityUpperBound << ")){\n" 950 << "this->" << broken.name << "= 1;\n" 951 << "}\n" 952 // makes one more iteration to compute the exact consistent tangent 953 // operator if required 954 << "converged = smt != CONSISTENTTANGENTOPERATOR;\n" 955 // the staggered scheme has converged 956 << computeStandardSystemOfImplicitEquations << " = true;\n" 957 << "} else {\n"; 958 generate_debug_message(false); 959 // the staggered scheme has not converged 960 os << "if(this->" << staggeredSchemeIterationCounter << " > this->" 961 << maximumNumberOfIterationsOfTheStaggeredScheme << "){\n" 962 << "tfel::raise<DivergenceException>(\"" 963 << "maximum number of iterations of " 964 << "the staggered algorithm reached\");\n" 965 << "}\n"; 966 // update the current estimate 967 os << staggeredSchemeBissectionAlgorithm << ".updateBounds(this->" 968 << currentEstimateOfThePorosityIncrement << "," 969 << "difference_between_successive_estimates_of_the_porosity_increment" 970 << ");\n"; 971 if (this->staggered_scheme_parameters.acceleration_algorithm == 972 StaggeredSchemeParameters::RELAXATION) { 973 os << staggeredSchemeBissectionAlgorithm << ".iterate(" 974 << nextEstimateOfThePorosityIncrement << ");\n"; 975 os << "this->" << currentEstimateOfThePorosityIncrement << " = " 976 << "(this->" << currentEstimateOfThePorosityIncrement << " + " 977 << nextEstimateOfThePorosityIncrement << ")/2;\n"; 978 } else if (this->staggered_scheme_parameters.acceleration_algorithm == 979 StaggeredSchemeParameters::AITKEN) { 980 os << staggeredSchemeAitkenAccelerationAlgorithm << ".accelerate(" 981 << nextEstimateOfThePorosityIncrement << ");\n"; 982 os << staggeredSchemeBissectionAlgorithm << ".iterate(" 983 << nextEstimateOfThePorosityIncrement << ");\n"; 984 apply_next_estimate_bounds(); 985 os << "this->" << currentEstimateOfThePorosityIncrement << " = " 986 << nextEstimateOfThePorosityIncrement << ";\n"; 987 } 988 // update the staggered iteration number 989 os << "++(this->" << staggeredSchemeIterationCounter << ");\n" 990 << "this->iter = 0;\n" 991 << "converged = false;\n" 992 << "}\n" 993 << "} // end of if(converged&&...)\n"; 994 acc.code = os.str(); 995 bd.setCode(uh, BehaviourData::AdditionalConvergenceChecks, acc, 996 BehaviourData::CREATEORAPPEND, BehaviourData::AT_END); 997 } 998 // 999 if ((this->isCoupledWithPorosityEvolution()) && 1000 (this->porosity_evolution_algorithm == 1001 mfront::bbrick::PorosityEvolutionAlgorithm:: 1002 STANDARD_IMPLICIT_SCHEME)) { 1003 const auto& broken = 1004 bd.getBehaviourData(uh) 1005 .getAuxiliaryStateVariableDescriptionByExternalName( 1006 tfel::glossary::Glossary::Broken); 1007 auto uasv = CodeBlock{}; 1008 const auto fmax = 1009 "(this->" + 1010 std::string(porosityUpperBoundSafetyFactorForFractureDetection) + 1011 ") * " + std::string(porosityUpperBound); 1012 uasv.code += "if(this->f > " + fmax + "){\n"; 1013 uasv.code += "this->" + broken.name + " = true;\n"; 1014 uasv.code += "}\n"; 1015 bd.setCode(uh, BehaviourData::UpdateAuxiliaryStateVariables, uasv, 1016 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 1017 } 1018 } // end of StandardElastoViscoPlasticityBrick::endTreatment 1019 1020 void StandardElastoViscoPlasticityBrick:: addElasticContributionToTheImplicitEquationAssociatedWithPorosityEvolution(CodeBlock & ib) const1021 addElasticContributionToTheImplicitEquationAssociatedWithPorosityEvolution( 1022 CodeBlock& ib) const { 1023 constexpr const auto uh = 1024 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 1025 const auto& f = 1026 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 1027 tfel::glossary::Glossary::Porosity); 1028 const auto f_ = f.name + "_"; 1029 const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(this->dsl); 1030 const auto requiresAnalyticalJacobian = 1031 ((idsl.getSolver().usesJacobian()) && 1032 (!idsl.getSolver().requiresNumericalJacobian())); 1033 const auto& broken = 1034 bd.getBehaviourData(uh) 1035 .getAuxiliaryStateVariableDescriptionByExternalName( 1036 tfel::glossary::Glossary::Broken); 1037 ib.code += "if(2 * (this->" + broken.name + ") < 1){\n"; 1038 ib.code += "f" + f.name + " -= (1 - " + f_ + ") * trace(this->deel);\n"; 1039 if (requiresAnalyticalJacobian) { 1040 ib.code += "df" + f.name + "_ddeel -= " + // 1041 "(1 - " + f_ + ") * Stensor::Id();\n"; 1042 ib.code += "df" + f.name + "_dd" + f.name + " += " + // 1043 "(" + this->bd.getClassName() + 1044 "::theta) * trace(this->deel);\n"; 1045 } 1046 ib.code += "}\n"; 1047 } // end of 1048 // addElasticContributionToTheImplicitEquationAssociatedWithPorosityEvolution 1049 1050 StandardElastoViscoPlasticityBrick::~StandardElastoViscoPlasticityBrick() = 1051 default; 1052 1053 } // end of namespace mfront 1054