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