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