1 /*! 2 * \file Chaboche2012KinematicHardeningRule.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 04/04/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 <sstream> 15 #include "TFEL/Raise.hxx" 16 #include "MFront/BehaviourBrick/BrickUtilities.hxx" 17 #include "MFront/BehaviourBrick/StressPotential.hxx" 18 #include "MFront/BehaviourBrick/StressCriterion.hxx" 19 #include "MFront/BehaviourBrick/OptionDescription.hxx" 20 #include "MFront/BehaviourBrick/Chaboche2012KinematicHardeningRule.hxx" 21 22 namespace mfront { 23 24 namespace bbrick { 25 26 std::vector<OptionDescription> getOptions() const27 Chaboche2012KinematicHardeningRule::getOptions() const { 28 auto opts = KinematicHardeningRuleBase::getOptions(); 29 opts.emplace_back("D", "back-strain callback coefficient", 30 OptionDescription::MATERIALPROPERTY); 31 opts.emplace_back("m", "", OptionDescription::MATERIALPROPERTY); 32 opts.emplace_back("Phi_inf", "", OptionDescription::MATERIALPROPERTY, 33 std::vector<std::string>(1u, "b"), 34 std::vector<std::string>{}); 35 opts.emplace_back("b", "", OptionDescription::MATERIALPROPERTY, 36 std::vector<std::string>(1u, "Phi_inf"), 37 std::vector<std::string>{}); 38 opts.emplace_back("w", "", OptionDescription::MATERIALPROPERTY); 39 return opts; 40 } // end of Chaboche2012KinematicHardeningRule::getOptions() 41 initialize(BehaviourDescription & bd,AbstractBehaviourDSL & dsl,const std::string & fid,const std::string & kid,const DataMap & d)42 void Chaboche2012KinematicHardeningRule::initialize( 43 BehaviourDescription& bd, 44 AbstractBehaviourDSL& dsl, 45 const std::string& fid, 46 const std::string& kid, 47 const DataMap& d) { 48 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 49 // `this` must be captured in gcc 4.7.2 50 auto get_mp = [&bd, &dsl, &d, &fid, &kid, this]( 51 BehaviourDescription::MaterialProperty& mp, const std::string& n) { 52 const auto nid = KinematicHardeningRule::getVariableId(n, fid, kid); 53 if (d.count(n) == 0) { 54 tfel::raise( 55 "Chaboche2012KinematicHardeningRule::initialize: " 56 "material property '" + 57 n + "' is not defined"); 58 } 59 mp = getBehaviourDescriptionMaterialProperty(dsl, n, d.at(n)); 60 declareParameterOrLocalVariable(bd, mp, "real", nid); 61 }; 62 KinematicHardeningRuleBase::initialize(bd, dsl, fid, kid, d); 63 get_mp(this->D, "D"); 64 get_mp(this->m, "m"); 65 if ((d.count("b") != 0) || (d.count("Phi_inf") != 0)) { 66 get_mp(this->bp, "b"); 67 get_mp(this->Phi_inf, "Phi_inf"); 68 } 69 get_mp(this->w, "w"); 70 // reserved named 71 const auto DJan = KinematicHardeningRule::getVariableId("DJa", fid, kid); 72 const auto iDJan = 73 KinematicHardeningRule::getVariableId("iDJa", fid, kid); 74 const auto dDJan = 75 KinematicHardeningRule::getVariableId("dDJa", fid, kid); 76 const auto Psin = KinematicHardeningRule::getVariableId("Psi", fid, kid); 77 const auto dPsin = 78 KinematicHardeningRule::getVariableId("dPsi", fid, kid) + "_d" + 79 KinematicHardeningRule::getVariableId("DJa", fid, kid); 80 bd.reserveName(uh, Psin); 81 bd.reserveName(uh, dPsin); 82 bd.reserveName(uh, DJan); 83 bd.reserveName(uh, iDJan); 84 bd.reserveName(uh, dDJan); 85 bd.reserveName(uh, "r" + DJan + "_m_1"); 86 bd.reserveName(uh, "r" + DJan + "_m"); 87 bd.reserveName(uh, DJan + "_m_1"); 88 bd.reserveName(uh, DJan + "_m"); 89 } // end of Chaboche2012KinematicHardeningRule::initialize 90 endTreatment(BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const std::string & fid,const std::string & kid) const91 void Chaboche2012KinematicHardeningRule::endTreatment( 92 BehaviourDescription& bd, 93 const AbstractBehaviourDSL& dsl, 94 const std::string& fid, 95 const std::string& kid) const { 96 KinematicHardeningRuleBase::endTreatment(bd, dsl, fid, kid); 97 constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 98 const auto Dn = KinematicHardeningRule::getVariableId("D", fid, kid); 99 const auto mn = KinematicHardeningRule::getVariableId("m", fid, kid); 100 const auto wn = KinematicHardeningRule::getVariableId("w", fid, kid); 101 const auto bn = KinematicHardeningRule::getVariableId("b", fid, kid); 102 const auto Phi_infn = 103 KinematicHardeningRule::getVariableId("Phi_inf", fid, kid); 104 auto c = generateMaterialPropertyInitializationCode(dsl, bd, Dn, this->D); 105 c += generateMaterialPropertyInitializationCode(dsl, bd, mn, this->m); 106 c += generateMaterialPropertyInitializationCode(dsl, bd, wn, this->w); 107 if (!this->bp.empty()) { 108 c += generateMaterialPropertyInitializationCode(dsl, bd, bn, this->bp); 109 c += generateMaterialPropertyInitializationCode(dsl, bd, Phi_infn, 110 this->Phi_inf); 111 } 112 if (!c.empty()) { 113 CodeBlock i; 114 i.code = c; 115 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, i, 116 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 117 } 118 } // end of Chaboche2012KinematicHardeningRule::endTreatment 119 120 std::string buildBackStrainImplicitEquations(const BehaviourDescription & bd,const StressPotential & sp,const StressCriterion & fc,const std::vector<std::shared_ptr<KinematicHardeningRule>> & khrs,const std::string & fid,const std::string & kid,const bool b) const121 Chaboche2012KinematicHardeningRule::buildBackStrainImplicitEquations( 122 const BehaviourDescription& bd, 123 const StressPotential& sp, 124 const StressCriterion& fc, 125 const std::vector<std::shared_ptr<KinematicHardeningRule>>& khrs, 126 const std::string& fid, 127 const std::string& kid, 128 const bool b) const { 129 const auto an = KinematicHardeningRule::getVariableId("a", fid, kid); 130 const auto Cn = KinematicHardeningRule::getVariableId("C", fid, kid); 131 const auto Dn = KinematicHardeningRule::getVariableId("D", fid, kid); 132 const auto mn = KinematicHardeningRule::getVariableId("m", fid, kid); 133 const auto wn = KinematicHardeningRule::getVariableId("w", fid, kid); 134 const auto bn = KinematicHardeningRule::getVariableId("b", fid, kid); 135 const auto Phi_infn = 136 KinematicHardeningRule::getVariableId("Phi_inf", fid, kid); 137 const auto DJan = KinematicHardeningRule::getVariableId("DJa", fid, kid); 138 const auto dDJan = 139 KinematicHardeningRule::getVariableId("dDJa", fid, kid); 140 const auto Phin = KinematicHardeningRule::getVariableId("Phi", fid, kid); 141 const auto dPhin = 142 KinematicHardeningRule::getVariableId("dPhi", fid, kid) + "_ddp" + 143 fid; 144 const auto Psin = KinematicHardeningRule::getVariableId("Psi", fid, kid); 145 const auto dPsin = 146 KinematicHardeningRule::getVariableId("dPsi", fid, kid) + "_d" + 147 KinematicHardeningRule::getVariableId("Ja", fid, kid); 148 const auto n = "n" + fid; 149 auto c = std::string{}; 150 c += "const auto " + DJan + " = "; 151 c += "(this->" + Dn + ")*sigmaeq(" + an + "_);\n"; 152 c += "if((2 * " + DJan + " > 3 * (this->" + wn + "))&&"; 153 c += "(this->dp" + fid + ">0)){\n"; 154 if (b) { 155 #if __cplusplus >= 201703L 156 c += "const auto [" + Psin + "," + dPsin + "] = "; 157 #else /* __cplusplus >= 201703L */ 158 c += "auto " + Psin + " = real{};\n"; 159 c += "auto " + dPsin + " = real{};\n"; 160 c += "std::tie(" + Psin + "," + dPsin + ") = "; 161 #endif /* __cplusplus >= 201703L */ 162 c += "[this," + DJan + "]() -> std::pair<real,real>{\n"; 163 c += "const auto r" + DJan + " = (" + DJan + "-3*(this->" + wn + 164 ")/2)/((1 - (this->" + wn + "))*" + DJan + ");\n"; 165 c += "const auto r" + DJan + "_m = "; 166 c += "pow(r" + DJan + ",this->" + mn + ");\n"; 167 c += "return {r" + DJan + "_m,"; 168 c += "r" + DJan + "_m * 3 * (this->" + wn + ") * (this->" + Dn + 169 ") *(this->" + mn + ")/((2 * " + DJan + "-3 * (this->" + wn + 170 ")) * " + DJan + ")};\n"; 171 c += "}();\n"; 172 } else { 173 c += "const auto " + Psin + " = "; 174 c += "pow((" + DJan + "-3*(this->" + wn + ")/2)/((1-this->" + wn + 175 ")*" + DJan + "),this->" + mn + ");\n"; 176 } 177 if (!this->bp.empty()) { 178 const auto exp_mbpn = 179 KinematicHardeningRule::getVariableId("exp_mbp", fid, kid); 180 c += "const auto " + exp_mbpn + " = "; 181 c += "exp(-(this->" + bn + ")*(this->p" + fid + 182 " + ((this->theta) * (this->dp" + fid + "))));\n"; 183 c += "const auto " + Phin + " = "; 184 c += Phi_infn + " + "; 185 c += "(1 - " + Phi_infn + ") * " + exp_mbpn + ";\n"; 186 c += "const auto " + dPhin + " = "; 187 c += "(this->" + bn + ") * (this->theta) *" + Phi_infn + " * " + 188 exp_mbpn + ";\n"; 189 c += "f" + an + " -= "; 190 c += "(this->dp" + fid + ") * (n" + fid + " - (this->" + Dn + ") * " + 191 Phin + " * " + Psin + " * " + an + "_);\n"; 192 } else { 193 c += "f" + an + " -= "; 194 c += "(this->dp" + fid + ") * (n" + fid + " - (this->" + Dn + ") * " + 195 Psin + " * " + an + "_);\n"; 196 } 197 if (b) { 198 c += "const auto dJ" + an + "_d" + an + " = "; 199 c += "3 * (this->" + Dn + ") * deviator(" + an + "_) / "; 200 c += "(2 * " + DJan + ");\n"; 201 if (!this->bp.empty()) { 202 c += "df" + an + "_ddp" + fid + " = - ("; 203 c += "n" + fid + "- (this->" + Dn + ") * " + Phin + " * " + Psin + 204 " * " + an + "_);\n"; 205 c += "df" + an + "_ddp" + fid + " += "; 206 c += "(this->dp" + fid + ") *"; 207 c += "(this->" + Dn + ") * " + dPhin + " * " + Psin + " * " + an + 208 "_);\n"; 209 } else { 210 c += "df" + an + "_ddp" + fid + " = -("; 211 c += "n" + fid + "- (this->" + Dn + ") * " + Psin + " * " + an + 212 "_);\n"; 213 } 214 // opposite of dfa_ds 215 const auto mdf_ds = "(this->dp" + fid + ") * dn" + fid + "_ds" + fid; 216 c += sp.generateImplicitEquationDerivatives( 217 bd, "StrainStensor", an, "-" + mdf_ds, fc.isNormalDeviatoric()); 218 auto kid2 = decltype(khrs.size()){}; 219 for (const auto& khr : khrs) { 220 c += khr->generateImplicitEquationDerivatives(an, mdf_ds, fid, 221 std::to_string(kid2)); 222 ++kid2; 223 } 224 c += "df" + an + "_dd" + an + " += "; 225 if (!this->bp.empty()) { 226 c += Phin + " * "; 227 } 228 c += "(this->" + Dn + ") * (this->theta) * (this->dp" + fid + ") * "; 229 c += "(" + Psin + " * Stensor4::Id() + " + dPsin + " * (" + an + 230 "_^dJ" + an + "_d" + an + "));\n"; 231 } 232 c += "} else {\n"; 233 // simple Prager hardening rule 234 c += "f" + an + " -= "; 235 c += "(this->dp" + fid + ") * n" + fid + ";\n"; 236 if (b) { 237 // opposite of dfa_ds 238 const auto mdf_ds = "(this->dp" + fid + ")*dn" + fid + "_ds" + fid; 239 c += "df" + an + "_ddp" + fid + " = - n" + fid + ";\n"; 240 c += sp.generateImplicitEquationDerivatives( 241 bd, "StrainStensor", an, "-" + mdf_ds, fc.isNormalDeviatoric()); 242 auto kid2 = decltype(khrs.size()){}; 243 for (const auto& khr : khrs) { 244 c += khr->generateImplicitEquationDerivatives(an, mdf_ds, fid, 245 std::to_string(kid2)); 246 ++kid2; 247 } 248 } 249 c += "}\n"; 250 return c; 251 } // end of 252 // Chaboche2012KinematicHardeningRule::buildBackStrainImplicitEquations 253 254 Chaboche2012KinematicHardeningRule::~Chaboche2012KinematicHardeningRule() = 255 default; 256 257 } // end of namespace bbrick 258 259 } // end of namespace mfront 260