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