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