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