1 /*! 2 * \file mfront/src/StrainBasedPorosityNucleationModelBase.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 05/04/2020 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 "MFront/ImplicitDSLBase.hxx" 18 #include "MFront/NonLinearSystemSolver.hxx" 19 #include "MFront/StandardElastoViscoPlasticityBrick.hxx" 20 #include "MFront/BehaviourBrick/BrickUtilities.hxx" 21 #include "MFront/BehaviourBrick/OptionDescription.hxx" 22 #include "MFront/BehaviourBrick/StrainBasedPorosityNucleationModelBase.hxx" 23 24 namespace mfront { 25 26 namespace bbrick { 27 28 bool StrainBasedPorosityNucleationModelBase:: requiresSavingNucleatedPorosity() const29 requiresSavingNucleatedPorosity() const { 30 return this->requiresLimitOnNucleationPorosity(); 31 } // end of 32 // StrainBasedPorosityNucleationModelBase::requiresSavingNucleatedPorosity() 33 34 std::vector<OptionDescription> getOptions() const35 StrainBasedPorosityNucleationModelBase::getOptions() const { 36 auto opts = PorosityNucleationModelBase::getOptions(); 37 if (this->requiresLimitOnNucleationPorosity()) { 38 opts.emplace_back("fmax", "maximum value of the nucleated porosity", 39 OptionDescription::MATERIALPROPERTY); 40 } 41 return opts; 42 } // end of StrainBasedPorosityNucleationModelBase::getOptions() 43 initialize(BehaviourDescription & bd,AbstractBehaviourDSL & dsl,const std::string & id,const DataMap & d)44 void StrainBasedPorosityNucleationModelBase::initialize( 45 BehaviourDescription& bd, 46 AbstractBehaviourDSL& dsl, 47 const std::string& id, 48 const DataMap& d) { 49 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 50 const auto mn = this->getModelName(); 51 PorosityNucleationModelBase::initialize(bd, dsl, id, d); 52 bd.appendToIncludes("#include \"TFEL/Material/" + mn + ".hxx\""); 53 const auto parameters = 54 PorosityNucleationModel::getVariableId("parameters", id); 55 addLocalVariable(bd, mn + "Parameters<real>", parameters); 56 if (this->requiresLimitOnNucleationPorosity()) { 57 const auto fmax_n = PorosityNucleationModel::getVariableId("fmax", id); 58 tfel::raise_if( 59 d.count("fmax") == 0, 60 "ChuNeedleman1980StressBasedPorosityNucleationModel::initialize: " 61 "material property 'fmax' is not defined"); 62 this->fmax = 63 getBehaviourDescriptionMaterialProperty(dsl, "fmax", d.at("fmax")); 64 declareParameterOrLocalVariable(bd, this->fmax, "real", fmax_n); 65 } 66 // porosity factor 67 const auto fn_n = "fn" + id; 68 bd.reserveName(uh, "d" + fn_n + "_ddp"); 69 } // end of StrainBasedPorosityNucleationModelBase::initialize 70 endTreatment(BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const StressPotential & sp,const std::map<std::string,std::shared_ptr<bbrick::InelasticFlow>> & iflows,const std::string & id) const71 void StrainBasedPorosityNucleationModelBase::endTreatment( 72 BehaviourDescription& bd, 73 const AbstractBehaviourDSL& dsl, 74 const StressPotential& sp, 75 const std::map<std::string, std::shared_ptr<bbrick::InelasticFlow>>& 76 iflows, 77 const std::string& id) const { 78 PorosityNucleationModelBase::endTreatment(bd, dsl, sp, iflows, id); 79 // 80 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 81 const auto mn = this->getModelName(); 82 const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(dsl); 83 const auto requiresAnalyticalJacobian = 84 ((idsl.getSolver().usesJacobian()) && 85 (!idsl.getSolver().requiresNumericalJacobian())); 86 // material coefficients initialisation 87 const auto parameters = 88 PorosityNucleationModel::getVariableId("parameters", id); 89 CodeBlock init; 90 for (const auto& mc : this->getMaterialCoefficientDescriptions()) { 91 const auto mc_n = PorosityNucleationModel::getVariableId(mc.name, id); 92 init.code += "this->" + parameters + "." + mc.name + " = this->" + mc_n + ";\n"; 93 } 94 if (this->requiresLimitOnNucleationPorosity()) { 95 const auto fmax_n = PorosityNucleationModel::getVariableId("fmax", id); 96 init.code += generateMaterialPropertyInitializationCode(dsl, bd, fmax_n, 97 this->fmax); 98 } 99 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, init, 100 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 101 // 102 const auto dfn_n = PorosityNucleationModel::getVariableId("dfn", id); 103 const auto dfn_ddp = 104 PorosityNucleationModel::getVariableId("dfn_ddp", id); 105 // sum of the equivalent plastic strains 106 const auto p = [&iflows, &bd] { 107 auto first = true; 108 auto v = std::string{}; 109 for (const auto& flow : iflows) { 110 if (!first) { 111 v += " + "; 112 } 113 v += "this->p" + flow.first; 114 first = false; 115 } 116 return v; 117 }(); 118 // sum of the increments of the equivalent plastic strains 119 const auto dp = [&iflows, &bd] { 120 auto first = true; 121 auto v = std::string{}; 122 for (const auto& flow : iflows) { 123 if (!first) { 124 v += " + "; 125 } 126 v += "std::max(this->dp" + flow.first + ", real(0))"; 127 first = false; 128 } 129 return v; 130 }(); 131 // 132 // porosity evolution 133 const auto& f = 134 bd.getBehaviourData(uh).getStateVariableDescriptionByExternalName( 135 tfel::glossary::Glossary::Porosity); 136 CodeBlock i; 137 if (this->porosity_evolution_algorithm == 138 PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 139 i.code += "if("; 140 i.code += StandardElastoViscoPlasticityBrick:: 141 computeStandardSystemOfImplicitEquations; 142 i.code += "){\n"; 143 } 144 auto call = "compute" + mn + "PorosityIncrementAndDerivative" + // 145 "(this->" + parameters + "," + p + ", " + dp + ", " + 146 bd.getClassName() + "::theta)"; 147 i.code += "#if __cplusplus >= 201703L\n"; 148 i.code += "const auto [" + dfn_n + ", " + dfn_ddp + "] = "; 149 i.code += call + ";\n"; 150 if (!requiresAnalyticalJacobian) { 151 i.code += "static_cast<void>(" + dfn_ddp + ");\n"; 152 } 153 i.code += "#else /* __cplusplus >= 201703L */\n"; 154 i.code += "auto " + dfn_n + " = real{};\n"; 155 if (requiresAnalyticalJacobian) { 156 i.code += "auto " + dfn_ddp + " = real{};\n"; 157 i.code += "std::tie(" + dfn_n + ", " + dfn_ddp + ") = "; 158 } else { 159 i.code += "std::tie(" + dfn_n + ", std::ignore) = "; 160 } 161 i.code += call + ";\n"; 162 i.code += "#endif /* __cplusplus >= 201703L */\n"; 163 i.code += "this->dfn" + id + " = " + dfn_n + ";\n"; 164 if (this->requiresLimitOnNucleationPorosity()) { 165 const auto fmax_n = 166 PorosityNucleationModel::getVariableId("fmax", id); 167 i.code += "if(this->fn" + id + " + this->dfn" + id + " >= this->" + 168 fmax_n + "){\n"; 169 i.code += "this->dfn" + id + " = this->" + fmax_n + " - this->fn;\n"; 170 i.code += "f" + f.name + " -= this->dfn" + id + ";\n"; 171 i.code += "} else {\n"; 172 } 173 i.code += "f" + f.name + " -= this->dfn" + id + ";\n"; 174 if (requiresAnalyticalJacobian) { 175 for (const auto& flow : iflows) { 176 i.code += "if(this->dp" + flow.first + " > real(0)){\n"; 177 i.code += "df" + f.name + "_ddp" + flow.first + // 178 " -= " + dfn_ddp + ";\n"; 179 i.code += "}\n"; 180 } 181 } 182 if (this->requiresLimitOnNucleationPorosity()) { 183 i.code += "}\n"; 184 } 185 if (this->porosity_evolution_algorithm == 186 PorosityEvolutionAlgorithm::STAGGERED_SCHEME) { 187 i.code += "}\n"; 188 } 189 bd.setCode(uh, BehaviourData::Integrator, i, 190 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 191 } // end of 192 // StrainBasedPorosityNucleationModelBase::endTreatment 193 194 std::string StrainBasedPorosityNucleationModelBase:: updateNextEstimateOfThePorosityIncrement(const BehaviourDescription & bd,const std::map<std::string,std::shared_ptr<bbrick::InelasticFlow>> & iflows,const std::string & id) const195 updateNextEstimateOfThePorosityIncrement( 196 const BehaviourDescription& bd, 197 const std::map<std::string, std::shared_ptr<bbrick::InelasticFlow>>& 198 iflows, 199 const std::string& id) const { 200 const auto mn = this->getModelName(); 201 const auto parameters = 202 PorosityNucleationModel::getVariableId("parameters", id); 203 // sum of the equivalent plastic strains 204 const auto p = [&iflows, &bd] { 205 auto first = true; 206 auto v = std::string{}; 207 for (const auto& flow : iflows) { 208 if (!first) { 209 v += " + "; 210 } 211 v += "this->p" + flow.first; 212 first = false; 213 } 214 return v; 215 }(); 216 // sum of the increments of the equivalent plastic strains 217 const auto dp = [&iflows, &bd] { 218 auto first = true; 219 auto v = std::string{}; 220 for (const auto& flow : iflows) { 221 if (!first) { 222 v += " + "; 223 } 224 v += "std::max(this->dp" + flow.first + ", real(0))"; 225 first = false; 226 } 227 return v; 228 }(); 229 // 230 auto c = std::string{}; 231 if (this->requiresLimitOnNucleationPorosity()) { 232 const auto dfn_n = PorosityNucleationModel::getVariableId("dfn", id); 233 const auto fmax_n = PorosityNucleationModel::getVariableId("fmax", id); 234 c += "const auto " + dfn_n + " = "; 235 c += "compute" + mn + "PorosityIncrement" + // 236 "(this->" + parameters + "," + p + "," + dp + ", " + 237 bd.getClassName() + "::theta);\n"; 238 c += "this->dfn" + id + " = "; 239 c += "std::min(" + dfn_n + ", this->" + fmax_n + " - this->fn);\n"; 240 } else { 241 c += "this->dfn" + id + " = compute" + mn + "PorosityIncrement" + // 242 "(this->" + parameters + "," + p + "," + dp + ", " + 243 bd.getClassName() + "::theta);\n"; 244 } 245 c += mfront::StandardElastoViscoPlasticityBrick:: 246 nextEstimateOfThePorosityIncrement; 247 c += " += this->dfn" + id + ";\n"; 248 return c; 249 } // end of 250 // StrainBasedPorosityNucleationModelBase::updateNextEstimateOfThePorosityIncrement 251 252 StrainBasedPorosityNucleationModelBase:: 253 ~StrainBasedPorosityNucleationModelBase() = default; 254 255 } // end of namespace bbrick 256 257 } // end of namespace mfront 258 259