1 /*! 2 * \file mfront/src/Cazacu2001StressCriterion.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 "MFront/BehaviourBrick/BrickUtilities.hxx" 16 #include "MFront/BehaviourBrick/StressPotential.hxx" 17 #include "MFront/BehaviourBrick/OptionDescription.hxx" 18 #include "MFront/BehaviourBrick/Cazacu2001StressCriterion.hxx" 19 20 namespace mfront { 21 22 namespace bbrick { 23 getOptions() const24 std::vector<OptionDescription> Cazacu2001StressCriterion::getOptions() 25 const { 26 auto opts = StressCriterionBase::getOptions(); 27 opts.emplace_back("a", "coefficients of J20", 28 OptionDescription::ARRAYOFMATERIALPROPERTIES); 29 opts.emplace_back("b", "coefficients of J30", 30 OptionDescription::ARRAYOFMATERIALPROPERTIES); 31 opts.emplace_back("c", "", OptionDescription::MATERIALPROPERTY); 32 return opts; 33 } // end of Cazacu2001StressCriterion::getOptions() 34 35 std::vector<mfront::BehaviourSymmetryType> getSupportedBehaviourSymmetries() const36 Cazacu2001StressCriterion::getSupportedBehaviourSymmetries() const { 37 return {mfront::ISOTROPIC, mfront::ORTHOTROPIC}; 38 } // end of Cazacu2001StressCriterion::getSupportedBehaviourSymmetries() 39 initialize(BehaviourDescription & bd,AbstractBehaviourDSL & dsl,const std::string & id,const DataMap & d,const Role r)40 void Cazacu2001StressCriterion::initialize(BehaviourDescription& bd, 41 AbstractBehaviourDSL& dsl, 42 const std::string& id, 43 const DataMap& d, 44 const Role r) { 45 using ConstantMaterialProperty = 46 mfront::BehaviourDescription::ConstantMaterialProperty; 47 tfel::raise_if(bd.getSymmetryType() != mfront::ORTHOTROPIC, 48 "Cazacu2004OrthotropicStressCriterion::initialize: " 49 "the behaviour must be orthotropic"); 50 StressCriterionBase::initialize(bd, dsl, id, d, r); 51 const auto an = StressCriterion::getVariableId("a", id, r); 52 const auto bn = StressCriterion::getVariableId("b", id, r); 53 const auto cn = StressCriterion::getVariableId("c", id, r); 54 bd.appendToIncludes( 55 "#include\"TFEL/Material/Cazacu2001YieldCriterion.hxx\""); 56 tfel::raise_if(d.count("c") == 0, 57 "Cazacu2001StressCriterion::initialize: " 58 "material property 'c' is not defined"); 59 // J2O coefficients 60 if (d.count("a") == 0) { 61 tfel::raise("Cazacu2001StressCriterion::initialize: entry '" + 62 std::string("a") + "' is not defined"); 63 } 64 this->a = getArrayOfBehaviourDescriptionMaterialProperties<6u>(dsl, "a", 65 d.at("a")); 66 if (this->a[0].is<ConstantMaterialProperty>()) { 67 std::vector<double> values(6u); 68 for (unsigned short i = 0; i != 6; ++i) { 69 tfel::raise_if(!this->a[i].is<ConstantMaterialProperty>(), 70 "Cazacu2001StressCriterion::initialize: " 71 "if one component of '" + 72 std::string("a") + 73 "' is a constant value, all components must be " 74 "constant values"); 75 values[i] = this->a[i].get<ConstantMaterialProperty>().value; 76 } 77 const auto ean = [&r] { 78 if (r == StressCriterion::STRESSANDFLOWCRITERION) { 79 return "J2OCoefficients"; 80 } else if (r == StressCriterion::STRESSCRITERION) { 81 return "StressCriterionJ2OCoefficients"; 82 } 83 return "FlowCriterionJ2OCoefficients"; 84 }(); 85 addParameter(bd, an, ean, 6, values); 86 } else { 87 addLocalVariable(bd, "real", an, 6); 88 } 89 // J3O coefficients 90 if (d.count("b") == 0) { 91 tfel::raise("Cazacu2001StressCriterion::initialize: entry '" + 92 std::string("b") + "' is not defined"); 93 } 94 this->b = getArrayOfBehaviourDescriptionMaterialProperties<11u>( 95 dsl, "b", d.at("b")); 96 if (this->b[0].is<ConstantMaterialProperty>()) { 97 std::vector<double> values(11u); 98 for (unsigned short i = 0; i != 11; ++i) { 99 tfel::raise_if(!this->b[i].is<ConstantMaterialProperty>(), 100 "Cazacu2001StressCriterion::initialize: " 101 "if one component of '" + 102 std::string("b") + 103 "' is a constant value, all components must be " 104 "constant values"); 105 values[i] = this->b[i].get<ConstantMaterialProperty>().value; 106 } 107 const auto ebn = [&r] { 108 if (r == StressCriterion::STRESSANDFLOWCRITERION) { 109 return "J3OCoefficients"; 110 } else if (r == StressCriterion::STRESSCRITERION) { 111 return "StressCriterionJ3OCoefficients"; 112 } 113 return "FlowCriterionJ3OCoefficients"; 114 }(); 115 addParameter(bd, bn, ebn, 11, values); 116 } else { 117 addLocalVariable(bd, "real", bn, 11); 118 } 119 this->cp = getBehaviourDescriptionMaterialProperty(dsl, "c", d.at("c")); 120 declareParameterOrLocalVariable(bd, this->cp, "real", cn); 121 } // end of Cazacu2001StressCriterion::initialize 122 endTreatment(BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const std::string & id,const Role r)123 void Cazacu2001StressCriterion::endTreatment( 124 BehaviourDescription& bd, 125 const AbstractBehaviourDSL& dsl, 126 const std::string& id, 127 const Role r) { 128 constexpr const auto uh = 129 tfel::material::ModellingHypothesis::UNDEFINEDHYPOTHESIS; 130 const auto an = StressCriterion::getVariableId("a", id, r); 131 const auto bn = StressCriterion::getVariableId("b", id, r); 132 const auto cn = StressCriterion::getVariableId("c", id, r); 133 if (bd.getOrthotropicAxesConvention() == 134 tfel::material::OrthotropicAxesConvention::PIPE) { 135 for (const auto mh : bd.getModellingHypotheses()) { 136 tfel::raise_if( 137 ((mh != ModellingHypothesis::TRIDIMENSIONAL) && 138 (mh != ModellingHypothesis::AXISYMMETRICAL) && 139 (mh != 140 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN) && 141 (mh != 142 ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS)), 143 "Cazacu2001Criterion::endTreatment: " 144 "with the Pipe orthotropic axis convention, " 145 "only the following modelling hypotheses are supported:\n" 146 "- Tridimensional\n" 147 "- Axisymmetrical\n" 148 "- AxisymmetricalGeneralisedPlaneStrain\n" 149 "- AxisymmetricalGeneralisedPlaneStress."); 150 } 151 } 152 if (bd.getOrthotropicAxesConvention() == 153 tfel::material::OrthotropicAxesConvention::PLATE) { 154 for (const auto mh : bd.getModellingHypotheses()) { 155 tfel::raise_if( 156 ((mh != ModellingHypothesis::TRIDIMENSIONAL) && 157 (mh != ModellingHypothesis::PLANESTRESS)), 158 "Cazacu2001Criterion::endTreatment: " 159 "with the Plate orthotropic axis convention, " 160 "only the following modelling hypotheses are supported:\n" 161 "- Tridimensional\n" 162 "- Plane stress"); 163 } 164 } else { 165 tfel::raise_if( 166 bd.getOrthotropicAxesConvention() != 167 tfel::material::OrthotropicAxesConvention::DEFAULT, 168 "Cazacu2001Criterion::endTreatment: " 169 "internal error, unsupported orthotropic axes convention"); 170 for (const auto mh : bd.getModellingHypotheses()) { 171 tfel::raise_if(mh != ModellingHypothesis::TRIDIMENSIONAL, 172 "Cazacu2001Criterion::endTreatment: " 173 "an orthotropic axes convention must be choosen when " 174 "defining a stress free expansion in behaviours " 175 "which shall be valid in other modelling hypothesis " 176 "than 'Tridimensional'.\n" 177 "Either restrict the validity of the behaviour to " 178 "'Tridimensional' (see @ModellingHypothesis) or " 179 "choose and orthotropic axes convention as on option " 180 "to the @OrthotropicBehaviour keyword"); 181 } 182 } 183 auto c = std::string{}; 184 c += generateMaterialPropertiesInitializationCode(dsl, bd, an, this->a); 185 c += generateMaterialPropertiesInitializationCode(dsl, bd, bn, this->b); 186 c += generateMaterialPropertyInitializationCode(dsl, bd, cn, this->cp); 187 if (!c.empty()) { 188 CodeBlock i; 189 i.code = c; 190 bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, i, 191 BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING); 192 } 193 } // end of Cazacu2001StressCriterion::endTreatment 194 computeElasticPrediction(const std::string & id,const BehaviourDescription &,const StressPotential &) const195 std::string Cazacu2001StressCriterion::computeElasticPrediction( 196 const std::string& id, 197 const BehaviourDescription&, 198 const StressPotential&) const { 199 const auto an = StressCriterion::getVariableId( 200 "a", id, StressCriterion::STRESSCRITERION); 201 const auto bn = StressCriterion::getVariableId( 202 "b", id, StressCriterion::STRESSCRITERION); 203 const auto cn = StressCriterion::getVariableId( 204 "c", id, StressCriterion::STRESSCRITERION); 205 return "const auto seqel" + id + 206 " = computeCazacu2001StressCriterion(sel" + id + ",this->" + an + 207 ",this->" + bn + ",this->" + cn + ");\n"; 208 } // end of Cazacu2001StressCriterion::computeElasticPrediction 209 computeCriterion(const std::string & id,const BehaviourDescription &,const StressPotential &) const210 std::string Cazacu2001StressCriterion::computeCriterion( 211 const std::string& id, 212 const BehaviourDescription&, 213 const StressPotential&) const { 214 const auto an = StressCriterion::getVariableId( 215 "a", id, StressCriterion::STRESSCRITERION); 216 const auto bn = StressCriterion::getVariableId( 217 "b", id, StressCriterion::STRESSCRITERION); 218 const auto cn = StressCriterion::getVariableId( 219 "c", id, StressCriterion::STRESSCRITERION); 220 return "const auto seq" + id + " = computeCazacu2001StressCriterion(s" + 221 id + ",this->" + an + ",this->" + bn + ",this->" + cn + ");\n"; 222 } // end of Cazacu2001StressCriterion::computeCriterion 223 computeNormal(const std::string & id,const BehaviourDescription & bd,const StressPotential & sp,const Role r) const224 std::string Cazacu2001StressCriterion::computeNormal( 225 const std::string& id, 226 const BehaviourDescription& bd, 227 const StressPotential& sp, 228 const Role r) const { 229 const auto an = StressCriterion::getVariableId( 230 "a", id, StressCriterion::STRESSCRITERION); 231 const auto bn = StressCriterion::getVariableId( 232 "b", id, StressCriterion::STRESSCRITERION); 233 const auto cn = StressCriterion::getVariableId("c", id, r); 234 auto c = std::string{}; 235 if ((r == STRESSCRITERION) || (r == STRESSANDFLOWCRITERION)) { 236 #if __cplusplus >= 201703L 237 c += "const auto [seq" + id + ",dseq" + id + "_ds" + id + "] = "; 238 c += "computeCazacu2001StressCriterionNormal(s" + id + ", this->" + an + 239 ", this->" + bn + ", this->" + cn + "," + 240 sp.getEquivalentStressLowerBound(bd) + ");\n"; 241 #else /* __cplusplus >= 201703L */ 242 c += "stress seq" + id + ";\n"; 243 c += "Stensor dseq" + id + "_ds" + id + ";\n"; 244 c += "std::tie(seq" + id + ",dseq" + id + "_ds" + id + ") = "; 245 c += "computeCazacu2001StressCriterionNormal(s" + id + ", this->" + an + 246 ", this->" + bn + ", this->" + cn + "," + 247 sp.getEquivalentStressLowerBound(bd) + ");\n"; 248 #endif /* __cplusplus >= 201703L */ 249 } 250 if (r == STRESSANDFLOWCRITERION) { 251 c += "const auto& n" + id + " = dseq" + id + "_ds" + id + ";\n"; 252 } 253 if (r == FLOWCRITERION) { 254 #if __cplusplus >= 201703L 255 c += "const auto [seqf" + id + ", n" + id + "] = "; 256 c += "computeCazacu2001StressCriterionNormal(s" + id + ", this->" + an + 257 ", this->" + bn + ", this->" + cn + "," + 258 sp.getEquivalentStressLowerBound(bd) + ");\n"; 259 #else /* __cplusplus >= 201703L */ 260 c += "stress seqf" + id + ";\n"; 261 c += "Stensor n" + id + ";\n"; 262 c += "std::tie(seqf" + id + ",n" + id + ") = "; 263 c += "computeCazacu2001StressCriterionNormal(s" + id + ", this->" + an + 264 ", this->" + bn + ", this->" + cn + "," + 265 sp.getEquivalentStressLowerBound(bd) + ");\n"; 266 #endif /* __cplusplus >= 201703L */ 267 } 268 return c; 269 } // end of Cazacu2001StressCriterion::computeNormal 270 computeNormalDerivative(const std::string & id,const BehaviourDescription & bd,const StressPotential & sp,const Role r) const271 std::string Cazacu2001StressCriterion::computeNormalDerivative( 272 const std::string& id, 273 const BehaviourDescription& bd, 274 const StressPotential& sp, 275 const Role r) const { 276 const auto an = StressCriterion::getVariableId( 277 "a", id, StressCriterion::STRESSCRITERION); 278 const auto bn = StressCriterion::getVariableId( 279 "b", id, StressCriterion::STRESSCRITERION); 280 const auto cn = StressCriterion::getVariableId("c", id, r); 281 auto c = std::string{}; 282 if ((r == STRESSCRITERION) || (r == STRESSANDFLOWCRITERION)) { 283 #if __cplusplus >= 201703L 284 c += "const auto [seq" + id + ",dseq" + id + "_ds" + id + ",d2seq" + 285 id + "_ds" + id + "ds" + id + "] = "; 286 c += "computeCazacu2001StressCriterionSecondDerivative(s" + id + 287 ", this->" + an + ", this->" + bn + ", this->" + cn + "," + 288 sp.getEquivalentStressLowerBound(bd) + ");\n"; 289 #else /* __cplusplus >= 201703L */ 290 c += "stress seq" + id + ";\n"; 291 c += "Stensor dseq" + id + "_ds" + id + ";\n"; 292 c += "Stensor4 d2seq" + id + "_ds" + id + "ds" + id + ";\n"; 293 c += "std::tie(seq" + id + ",dseq" + id + "_ds" + id + ",d2seq" + id + 294 "_ds" + id + "ds" + id + ") = "; 295 c += "computeCazacu2001StressCriterionSecondDerivative(s" + id + 296 ", this->" + an + ", this->" + bn + ", this->" + cn + "," + 297 sp.getEquivalentStressLowerBound(bd) + ");\n"; 298 #endif /* __cplusplus >= 201703L */ 299 } 300 if (r == STRESSANDFLOWCRITERION) { 301 c += "const auto& n" + id + " = dseq" + id + "_ds" + id + ";\n"; 302 c += "const auto& dn" + id + "_ds" + id + " = "; 303 c += "d2seq" + id + "_ds" + id + "ds" + id + ";\n"; 304 } 305 if (r == FLOWCRITERION) { 306 #if __cplusplus >= 201703L 307 c += "const auto [seqf" + id + ", n" + id + ", dn" + id + "_ds" + id + 308 "] = "; 309 c += "computeCazacu2001StressCriterionSecondDerivative(s" + id + 310 ", this->" + an + ", this->" + bn + ",this->" + cn + "," + 311 sp.getEquivalentStressLowerBound(bd) + ");\n"; 312 #else /* __cplusplus >= 201703L */ 313 c += "stress seqf" + id + ";\n"; 314 c += "Stensor n" + id + ";\n"; 315 c += "Stensor4 dn" + id + "_ds" + id + ";\n"; 316 c += 317 "std::tie(seqf" + id + ",n" + id + ",dn" + id + "_ds" + id + ") = "; 318 c += "computeCazacu2001StressCriterionSecondDerivative(s" + id + 319 +", this->" + an + ", this->" + bn + ",this->" + cn + "," + 320 sp.getEquivalentStressLowerBound(bd) + ");\n"; 321 #endif /* __cplusplus >= 201703L */ 322 } 323 return c; 324 } // end of Cazacu2001StressCriterion::computeNormalDerivative 325 isCoupledWithPorosityEvolution() const326 bool Cazacu2001StressCriterion::isCoupledWithPorosityEvolution() const { 327 return false; 328 } // end of Cazacu2001StressCriterion::isCoupledWithPorosityEvolution 329 isNormalDeviatoric() const330 bool Cazacu2001StressCriterion::isNormalDeviatoric() const { 331 return true; 332 } // end of Cazacu2001StressCriterion::isNormalDeviatoric 333 334 StressCriterion::PorosityEffectOnFlowRule getPorosityEffectOnEquivalentPlasticStrain() const335 Cazacu2001StressCriterion::getPorosityEffectOnEquivalentPlasticStrain() const { 336 return StressCriterion::NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN; 337 } // end of Cazacu2001StressCriterion::getPorosityEffectOnEquivalentPlasticStrain() 338 339 Cazacu2001StressCriterion::~Cazacu2001StressCriterion() = default; 340 341 } // end of namespace bbrick 342 343 } // end of namespace mfront 344