1 /*!
2  * \file   Chaboche2012KinematicHardeningRule.cxx
3  * \brief
4  * \author Thomas Helfer
5  * \date   04/04/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 <sstream>
15 #include "TFEL/Raise.hxx"
16 #include "MFront/BehaviourBrick/BrickUtilities.hxx"
17 #include "MFront/BehaviourBrick/StressPotential.hxx"
18 #include "MFront/BehaviourBrick/StressCriterion.hxx"
19 #include "MFront/BehaviourBrick/OptionDescription.hxx"
20 #include "MFront/BehaviourBrick/Chaboche2012KinematicHardeningRule.hxx"
21 
22 namespace mfront {
23 
24   namespace bbrick {
25 
26     std::vector<OptionDescription>
getOptions() const27     Chaboche2012KinematicHardeningRule::getOptions() const {
28       auto opts = KinematicHardeningRuleBase::getOptions();
29       opts.emplace_back("D", "back-strain callback coefficient",
30                         OptionDescription::MATERIALPROPERTY);
31       opts.emplace_back("m", "", OptionDescription::MATERIALPROPERTY);
32       opts.emplace_back("Phi_inf", "", OptionDescription::MATERIALPROPERTY,
33                         std::vector<std::string>(1u, "b"),
34                         std::vector<std::string>{});
35       opts.emplace_back("b", "", OptionDescription::MATERIALPROPERTY,
36                         std::vector<std::string>(1u, "Phi_inf"),
37                         std::vector<std::string>{});
38       opts.emplace_back("w", "", OptionDescription::MATERIALPROPERTY);
39       return opts;
40     }  // end of Chaboche2012KinematicHardeningRule::getOptions()
41 
initialize(BehaviourDescription & bd,AbstractBehaviourDSL & dsl,const std::string & fid,const std::string & kid,const DataMap & d)42     void Chaboche2012KinematicHardeningRule::initialize(
43         BehaviourDescription& bd,
44         AbstractBehaviourDSL& dsl,
45         const std::string& fid,
46         const std::string& kid,
47         const DataMap& d) {
48       constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS;
49       // `this` must be captured in gcc 4.7.2
50       auto get_mp = [&bd, &dsl, &d, &fid, &kid, this](
51           BehaviourDescription::MaterialProperty& mp, const std::string& n) {
52         const auto nid = KinematicHardeningRule::getVariableId(n, fid, kid);
53         if (d.count(n) == 0) {
54           tfel::raise(
55               "Chaboche2012KinematicHardeningRule::initialize: "
56               "material property '" +
57               n + "' is not defined");
58         }
59         mp = getBehaviourDescriptionMaterialProperty(dsl, n, d.at(n));
60         declareParameterOrLocalVariable(bd, mp, "real", nid);
61       };
62       KinematicHardeningRuleBase::initialize(bd, dsl, fid, kid, d);
63       get_mp(this->D, "D");
64       get_mp(this->m, "m");
65       if ((d.count("b") != 0) || (d.count("Phi_inf") != 0)) {
66         get_mp(this->bp, "b");
67         get_mp(this->Phi_inf, "Phi_inf");
68       }
69       get_mp(this->w, "w");
70       // reserved named
71       const auto DJan = KinematicHardeningRule::getVariableId("DJa", fid, kid);
72       const auto iDJan =
73           KinematicHardeningRule::getVariableId("iDJa", fid, kid);
74       const auto dDJan =
75           KinematicHardeningRule::getVariableId("dDJa", fid, kid);
76       const auto Psin = KinematicHardeningRule::getVariableId("Psi", fid, kid);
77       const auto dPsin =
78           KinematicHardeningRule::getVariableId("dPsi", fid, kid) + "_d" +
79           KinematicHardeningRule::getVariableId("DJa", fid, kid);
80       bd.reserveName(uh, Psin);
81       bd.reserveName(uh, dPsin);
82       bd.reserveName(uh, DJan);
83       bd.reserveName(uh, iDJan);
84       bd.reserveName(uh, dDJan);
85       bd.reserveName(uh, "r" + DJan + "_m_1");
86       bd.reserveName(uh, "r" + DJan + "_m");
87       bd.reserveName(uh, DJan + "_m_1");
88       bd.reserveName(uh, DJan + "_m");
89     }  // end of Chaboche2012KinematicHardeningRule::initialize
90 
endTreatment(BehaviourDescription & bd,const AbstractBehaviourDSL & dsl,const std::string & fid,const std::string & kid) const91     void Chaboche2012KinematicHardeningRule::endTreatment(
92         BehaviourDescription& bd,
93         const AbstractBehaviourDSL& dsl,
94         const std::string& fid,
95         const std::string& kid) const {
96       KinematicHardeningRuleBase::endTreatment(bd, dsl, fid, kid);
97       constexpr const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS;
98       const auto Dn = KinematicHardeningRule::getVariableId("D", fid, kid);
99       const auto mn = KinematicHardeningRule::getVariableId("m", fid, kid);
100       const auto wn = KinematicHardeningRule::getVariableId("w", fid, kid);
101       const auto bn = KinematicHardeningRule::getVariableId("b", fid, kid);
102       const auto Phi_infn =
103           KinematicHardeningRule::getVariableId("Phi_inf", fid, kid);
104       auto c = generateMaterialPropertyInitializationCode(dsl, bd, Dn, this->D);
105       c += generateMaterialPropertyInitializationCode(dsl, bd, mn, this->m);
106       c += generateMaterialPropertyInitializationCode(dsl, bd, wn, this->w);
107       if (!this->bp.empty()) {
108         c += generateMaterialPropertyInitializationCode(dsl, bd, bn, this->bp);
109         c += generateMaterialPropertyInitializationCode(dsl, bd, Phi_infn,
110                                                         this->Phi_inf);
111       }
112       if (!c.empty()) {
113         CodeBlock i;
114         i.code = c;
115         bd.setCode(uh, BehaviourData::BeforeInitializeLocalVariables, i,
116                    BehaviourData::CREATEORAPPEND, BehaviourData::AT_BEGINNING);
117       }
118     }  // end of Chaboche2012KinematicHardeningRule::endTreatment
119 
120     std::string
buildBackStrainImplicitEquations(const BehaviourDescription & bd,const StressPotential & sp,const StressCriterion & fc,const std::vector<std::shared_ptr<KinematicHardeningRule>> & khrs,const std::string & fid,const std::string & kid,const bool b) const121     Chaboche2012KinematicHardeningRule::buildBackStrainImplicitEquations(
122         const BehaviourDescription& bd,
123         const StressPotential& sp,
124         const StressCriterion& fc,
125         const std::vector<std::shared_ptr<KinematicHardeningRule>>& khrs,
126         const std::string& fid,
127         const std::string& kid,
128         const bool b) const {
129       const auto an = KinematicHardeningRule::getVariableId("a", fid, kid);
130       const auto Cn = KinematicHardeningRule::getVariableId("C", fid, kid);
131       const auto Dn = KinematicHardeningRule::getVariableId("D", fid, kid);
132       const auto mn = KinematicHardeningRule::getVariableId("m", fid, kid);
133       const auto wn = KinematicHardeningRule::getVariableId("w", fid, kid);
134       const auto bn = KinematicHardeningRule::getVariableId("b", fid, kid);
135       const auto Phi_infn =
136           KinematicHardeningRule::getVariableId("Phi_inf", fid, kid);
137       const auto DJan = KinematicHardeningRule::getVariableId("DJa", fid, kid);
138       const auto dDJan =
139           KinematicHardeningRule::getVariableId("dDJa", fid, kid);
140       const auto Phin = KinematicHardeningRule::getVariableId("Phi", fid, kid);
141       const auto dPhin =
142           KinematicHardeningRule::getVariableId("dPhi", fid, kid) + "_ddp" +
143           fid;
144       const auto Psin = KinematicHardeningRule::getVariableId("Psi", fid, kid);
145       const auto dPsin =
146           KinematicHardeningRule::getVariableId("dPsi", fid, kid) + "_d" +
147           KinematicHardeningRule::getVariableId("Ja", fid, kid);
148       const auto n = "n" + fid;
149       auto c = std::string{};
150       c += "const auto " + DJan + " = ";
151       c += "(this->" + Dn + ")*sigmaeq(" + an + "_);\n";
152       c += "if((2 * " + DJan + " > 3 * (this->" + wn + "))&&";
153       c += "(this->dp" + fid + ">0)){\n";
154       if (b) {
155 #if __cplusplus >= 201703L
156         c += "const auto [" + Psin + "," + dPsin + "] = ";
157 #else  /* __cplusplus >= 201703L */
158         c += "auto " + Psin + " = real{};\n";
159         c += "auto " + dPsin + " = real{};\n";
160         c += "std::tie(" + Psin + "," + dPsin + ") = ";
161 #endif /* __cplusplus >= 201703L */
162         c += "[this," + DJan + "]() -> std::pair<real,real>{\n";
163         c += "const auto r" + DJan + " = (" + DJan + "-3*(this->" + wn +
164              ")/2)/((1 - (this->" + wn + "))*" + DJan + ");\n";
165         c += "const auto r" + DJan + "_m = ";
166         c += "pow(r" + DJan + ",this->" + mn + ");\n";
167         c += "return {r" + DJan + "_m,";
168         c += "r" + DJan + "_m * 3 * (this->" + wn + ") * (this->" + Dn +
169              ") *(this->" + mn + ")/((2 * " + DJan + "-3 * (this->" + wn +
170              ")) * " + DJan + ")};\n";
171         c += "}();\n";
172       } else {
173         c += "const auto " + Psin + " = ";
174         c += "pow((" + DJan + "-3*(this->" + wn + ")/2)/((1-this->" + wn +
175              ")*" + DJan + "),this->" + mn + ");\n";
176       }
177       if (!this->bp.empty()) {
178         const auto exp_mbpn =
179             KinematicHardeningRule::getVariableId("exp_mbp", fid, kid);
180         c += "const auto " + exp_mbpn + " = ";
181         c += "exp(-(this->" + bn + ")*(this->p" + fid +
182              " + ((this->theta) * (this->dp" + fid + "))));\n";
183         c += "const auto " + Phin + " = ";
184         c += Phi_infn + " + ";
185         c += "(1 - " + Phi_infn + ") * " + exp_mbpn + ";\n";
186         c += "const auto " + dPhin + " = ";
187         c += "(this->" + bn + ") * (this->theta) *" + Phi_infn + " * " +
188              exp_mbpn + ";\n";
189         c += "f" + an + " -= ";
190         c += "(this->dp" + fid + ") * (n" + fid + " - (this->" + Dn + ") * " +
191              Phin + " * " + Psin + " * " + an + "_);\n";
192       } else {
193         c += "f" + an + " -= ";
194         c += "(this->dp" + fid + ") * (n" + fid + " - (this->" + Dn + ") * " +
195              Psin + " * " + an + "_);\n";
196       }
197       if (b) {
198         c += "const auto dJ" + an + "_d" + an + " = ";
199         c += "3 * (this->" + Dn + ") * deviator(" + an + "_) / ";
200         c += "(2 * " + DJan + ");\n";
201         if (!this->bp.empty()) {
202           c += "df" + an + "_ddp" + fid + " = - (";
203           c += "n" + fid + "- (this->" + Dn + ") * " + Phin + " * " + Psin +
204                " * " + an + "_);\n";
205           c += "df" + an + "_ddp" + fid + " += ";
206           c += "(this->dp" + fid + ") *";
207           c += "(this->" + Dn + ") * " + dPhin + " * " + Psin + " * " + an +
208                "_);\n";
209         } else {
210           c += "df" + an + "_ddp" + fid + " = -(";
211           c += "n" + fid + "- (this->" + Dn + ") * " + Psin + " * " + an +
212                "_);\n";
213         }
214         // opposite of dfa_ds
215         const auto mdf_ds = "(this->dp" + fid + ") * dn" + fid + "_ds" + fid;
216         c += sp.generateImplicitEquationDerivatives(
217             bd, "StrainStensor", an, "-" + mdf_ds, fc.isNormalDeviatoric());
218         auto kid2 = decltype(khrs.size()){};
219         for (const auto& khr : khrs) {
220           c += khr->generateImplicitEquationDerivatives(an, mdf_ds, fid,
221                                                         std::to_string(kid2));
222           ++kid2;
223         }
224         c += "df" + an + "_dd" + an + " += ";
225         if (!this->bp.empty()) {
226           c += Phin + " * ";
227         }
228         c += "(this->" + Dn + ") * (this->theta) * (this->dp" + fid + ") * ";
229         c += "(" + Psin + " * Stensor4::Id() + " + dPsin + " * (" + an +
230              "_^dJ" + an + "_d" + an + "));\n";
231       }
232       c += "} else {\n";
233       // simple Prager hardening rule
234       c += "f" + an + " -= ";
235       c += "(this->dp" + fid + ") * n" + fid + ";\n";
236       if (b) {
237         // opposite of dfa_ds
238         const auto mdf_ds = "(this->dp" + fid + ")*dn" + fid + "_ds" + fid;
239         c += "df" + an + "_ddp" + fid + " = - n" + fid + ";\n";
240         c += sp.generateImplicitEquationDerivatives(
241             bd, "StrainStensor", an, "-" + mdf_ds, fc.isNormalDeviatoric());
242         auto kid2 = decltype(khrs.size()){};
243         for (const auto& khr : khrs) {
244           c += khr->generateImplicitEquationDerivatives(an, mdf_ds, fid,
245                                                         std::to_string(kid2));
246           ++kid2;
247         }
248       }
249       c += "}\n";
250       return c;
251     }  // end of
252     // Chaboche2012KinematicHardeningRule::buildBackStrainImplicitEquations
253 
254     Chaboche2012KinematicHardeningRule::~Chaboche2012KinematicHardeningRule() =
255         default;
256 
257   }  // end of namespace bbrick
258 
259 }  // end of namespace mfront
260