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