1 /*!
2  * \file   FiniteStrainSingleCrystalBrick.cxx
3  * \brief
4  * \author Thomas Helfer
5  * \date   04/10/2016
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<fstream>
15 #include<sstream>
16 #include<stdexcept>
17 #include "TFEL/Raise.hxx"
18 #include "TFEL/Utilities/Data.hxx"
19 #include "TFEL/Utilities/StringAlgorithms.hxx"
20 #include "TFEL/Glossary/Glossary.hxx"
21 #include "TFEL/Glossary/GlossaryEntry.hxx"
22 #include "TFEL/System/System.hxx"
23 #include "MFront/DSLUtilities.hxx"
24 #include "MFront/MFrontHeader.hxx"
25 #include "MFront/MFrontLogStream.hxx"
26 #include "MFront/AbstractBehaviourDSL.hxx"
27 #include "MFront/LocalDataStructure.hxx"
28 #include "MFront/BehaviourDescription.hxx"
29 #include "MFront/ImplicitDSLBase.hxx"
30 #include "MFront/NonLinearSystemSolver.hxx"
31 #include "MFront/BehaviourBrick/BrickUtilities.hxx"
32 #include "MFront/BehaviourBrick/OptionDescription.hxx"
33 #include "MFront/BehaviourBrick/HookeStressPotentialBase.hxx"
34 #include "MFront/FiniteStrainSingleCrystalBrick.hxx"
35 
36 namespace mfront{
37 
38   const char* const FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption =
39     "shifted_elastic_deformation_gradient";
40 
41   const char* const FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute =
42     "FiniteStrainSingleCrystalBrick::shifted_elastic_deformation_gradient";
43 
FiniteStrainSingleCrystalBrick(AbstractBehaviourDSL & dsl_,BehaviourDescription & mb_)44   FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick(
45       AbstractBehaviourDSL& dsl_, BehaviourDescription& mb_)
46       : BehaviourBrickBase(dsl_, mb_) {
47   }  // end of FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick
48 
getName() const49   std::string FiniteStrainSingleCrystalBrick::getName() const{
50     return "FiniteStrainSingleCrystal";
51   }
52 
getDescription() const53   BehaviourBrickDescription FiniteStrainSingleCrystalBrick::getDescription() const {
54     auto d = BehaviourBrickDescription{};
55     d.behaviourType = tfel::material::MechanicalBehaviourBase::STANDARDFINITESTRAINBEHAVIOUR;
56     d.integrationScheme = IntegrationScheme::IMPLICITSCHEME;
57     d.supportedModellingHypotheses = {ModellingHypothesis::TRIDIMENSIONAL};
58     d.supportedBehaviourSymmetries = {mfront::ORTHOTROPIC};
59     d.requireCrystalStructureDefinition = true;
60     d.requireStiffnessTensorDefinition = true;
61     d.managedIntegrationVariables.push_back("eel");
62     d.managedCodeBlocks = {BehaviourData::ComputePredictionOperator,
63                            BehaviourData::ComputeThermodynamicForces,
64                            BehaviourData::ComputeTangentOperator};
65     return d;
66   } // end of FiniteStrainSingleCrystalBrick::getDescription
67 
68   std::vector<bbrick::OptionDescription>
getOptions(const bool) const69   FiniteStrainSingleCrystalBrick::getOptions(const bool) const {
70     auto opts = mfront::bbrick::HookeStressPotentialBase::
71         getOrthotropicBehaviourElasticMaterialPropertiesOptions();
72     opts.emplace_back(
73         FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption, "",
74         bbrick::OptionDescription::BOOLEAN);
75     return opts;
76   }  // end of FiniteStrainSingleCrystalBrick::getOptions
77 
initialize(const Parameters & p,const DataMap & d)78   void FiniteStrainSingleCrystalBrick::initialize(const Parameters& p,
79                                                   const DataMap& d) {
80     auto throw_if = [](const bool b,const std::string& m){
81       tfel::raise_if(b,"FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick: "+m);
82     };
83     const auto uh = ModellingHypothesis::UNDEFINEDHYPOTHESIS;
84     const auto h = ModellingHypothesis::TRIDIMENSIONAL;
85     // basic checks
86     if(this->bd.areModellingHypothesesDefined()){
87       const auto bmh = this->bd.getModellingHypotheses();
88       throw_if(bmh.size()!=1u,"the only supported hypothesis is 'Tridimensional'");
89       throw_if(*(bmh.begin())!=ModellingHypothesis::TRIDIMENSIONAL,
90 	       "the only supported hypothesis is 'Tridimensional'");
91     } else {
92       this->bd.setModellingHypotheses({ModellingHypothesis::TRIDIMENSIONAL});
93     }
94     throw_if(this->bd.getBehaviourType()!=BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR,
95 	     "this BehaviourBrick is only usable for small strain behaviours");
96     throw_if(this->bd.getIntegrationScheme()!=BehaviourDescription::IMPLICITSCHEME,
97 	     "this BehaviourBrick is only usable in implicit schemes");
98     throw_if(!p.empty(),"no parameters allowed");
99     // material symmetry
100     if(this->bd.isSymmetryTypeDefined()){
101       throw_if(this->bd.getSymmetryType()!=mfront::ORTHOTROPIC,
102 	       "material must be declared orthotropic");
103     } else {
104       this->bd.setSymmetryType(mfront::ORTHOTROPIC);
105     }
106     throw_if(this->bd.getElasticSymmetryType()!=mfront::ORTHOTROPIC,
107 	     "elastic symmetry type must be orthortropic");
108     // declaring the elastic strain as the first integration variable
109     throw_if(!this->bd.getBehaviourData(h).getIntegrationVariables().empty(),
110 	     "no integration variable shall be declared before declaring the "
111 	     "'FiniteStrainSingleCrystal' brick");
112     // options
113     this->checkOptionsNames(d);
114     // elastic properties
115     auto get_mp = [this, &d](const char* const n) {
116       if (d.count(n) == 0) {
117         tfel::raise(
118             "FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick: "
119             "option '" +
120             std::string(n) + "' is not defined");
121       }
122       return mfront::bbrick::getBehaviourDescriptionMaterialProperty(
123           this->dsl, n, d.at(n));
124     };
125     if ((d.count("young_modulus1") != 0) || (d.count("young_modulus2") != 0) ||
126         (d.count("young_modulus3") != 0) || (d.count("poisson_ratio12") != 0) ||
127         (d.count("poisson_ratio23") != 0) ||
128         (d.count("poisson_ratio13") != 0) ||
129         (d.count("shear_modulus12") != 0) ||
130         (d.count("shear_modulus23") != 0) ||
131         (d.count("shear_modulus13") != 0)) {
132       std::vector<BehaviourDescription::MaterialProperty> emps = {
133           get_mp("young_modulus1"),  get_mp("young_modulus2"),
134           get_mp("young_modulus3"),  get_mp("poisson_ratio12"),
135           get_mp("poisson_ratio23"), get_mp("poisson_ratio13"),
136           get_mp("shear_modulus12"), get_mp("shear_modulus23"),
137           get_mp("shear_modulus13")};
138       bd.setElasticMaterialProperties(emps);
139       bd.setAttribute(BehaviourDescription::computesStiffnessTensor, true,
140                       false);
141       bd.setAttribute(BehaviourDescription::requiresUnAlteredStiffnessTensor,
142                       true, false);
143     }
144     // shifted deformation gradient
145     if (d.count(
146             FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption)) {
147       const auto on = std::string(
148           FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption);
149       const auto an = std::string(
150           FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute);
151       const auto& b = d.at(on);
152       throw_if(!b.is<bool>(), "invalid type for option '" + on + "'");
153       throw_if(this->bd.hasAttribute(an),
154                "attribute '" + an + "' already used");
155       this->bd.setAttribute(an, b.get<bool>(), false);
156     }
157     // elastic strain
158     VariableDescription eel("StrainStensor", "εᵉˡ", "eel", 1u, 0u);
159     eel.description = "elastic strain";
160     this->bd.addIntegrationVariable(uh, eel, BehaviourData::UNREGISTRED);
161     this->bd.setGlossaryName(h, "eel", tfel::glossary::Glossary::ElasticStrain);
162     // declaring the plastic slip
163     auto add_plastic_slips = [this, throw_if, uh, h] {
164       throw_if(
165           this->bd.getBehaviourData(h).getIntegrationVariables().size() != 1u,
166           "no integration variable shall be declared before declaring the "
167           "'FiniteStrainSingleCrystal' brick");
168       const auto& ss = this->bd.getSlipSystems();
169       VariableDescription g("strain", "g", ss.getNumberOfSlipSystems(), 0u);
170       g.description = "plastic slip";
171       this->bd.addStateVariable(uh, g, BehaviourData::UNREGISTRED);
172       this->bd.setEntryName(h, "g", "PlasticSlip");
173     };
174     if ((!this->bd.hasCrystalStructure()) ||
175         (!this->bd.areSlipSystemsDefined())) {
176       // no slip systems defined, yet
177       auto& idsl = dynamic_cast<ImplicitDSLBase&>(this->dsl);
178       idsl.addHook("@SlipSystem", add_plastic_slips);
179       idsl.addHook("@GlidingSystem", add_plastic_slips);
180       idsl.addHook("@SlidingSystem", add_plastic_slips);
181       idsl.addHook("@SlipSystems", add_plastic_slips);
182       idsl.addHook("@GlidingSystems", add_plastic_slips);
183       idsl.addHook("@SlidingSystems", add_plastic_slips);
184       } else {
185         add_plastic_slips();
186       }
187       // declaring the elastic part of the deformation gradient
188       VariableDescription Fe("DeformationGradientTensor", "Fe", 1u, 0u);
189       if (this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::
190                                           shiftedDeformationGradientAttribute,
191                                       false)) {
192         Fe.description = "shifted elastic part of the deformation gradient";
193       } else {
194         Fe.description = "elastic part of the deformation gradient";
195       }
196       this->bd.addAuxiliaryStateVariable(h, Fe, BehaviourData::UNREGISTRED);
197       if (this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::
198                                           shiftedDeformationGradientAttribute,
199                                       false)) {
200         this->bd.setEntryName(h, "Fe",
201                               "ShiftedElasticPartOfTheDeformationGradient");
202       } else {
203         this->bd.setEntryName(h, "Fe", "ElasticPartOfTheDeformationGradient");
204       }
205       // additional includes
206       this->bd.appendToIncludes("#include\"TFEL/Math/General/CubicRoots.hxx\"");
207       // set that the tangent operator is computed
208       this->bd.setAttribute(h, BehaviourData::hasConsistentTangentOperator,
209                             true);
210       // reserve some specific variables
211       this->bd.reserveName(h, "ss");
212       this->bd.reserveName(h, "fsscb_data");
213       this->bd.reserveName(h, "fsscb_tprd");
214       this->bd.reserveName(h, "fsscb_dfeel_dinv_dFp");
215       this->bd.reserveName(h, "fsscb_dC_dFe");
216       this->bd.reserveName(h, "fsscb_dS_dFe");
217       this->bd.reserveName(h, "fsscb_dtau_dFe");
218       this->bd.reserveName(h, "fsscb_dFe_dDF_tot");
219       this->bd.reserveName(h, "fsscb_dfeel_dDF");
220       this->bd.reserveName(h, "fsscb_Je");
221       this->bd.reserveName(h, "fsscb_Jg");
222       this->bd.reserveName(h, "fsscb_dinv_Fp_dDF");
223       this->bd.reserveName(h, "fsscb_dFe_dDF");
224   } // end of FiniteStrainSingleCrystalBrick::initialize
225 
completeVariableDeclaration() const226   void FiniteStrainSingleCrystalBrick::completeVariableDeclaration() const {
227     using tfel::glossary::Glossary;
228     const auto h = ModellingHypothesis::TRIDIMENSIONAL;
229     auto throw_if = [](const bool b, const std::string& m) {
230       tfel::raise_if(b,
231                      "FiniteStrainSingleCrystalBrick:"
232                      ":completeVariableDeclaration: " +
233                          m);
234     };
235     if (getVerboseMode() >= VERBOSE_DEBUG) {
236       getLogStream() << "FiniteStrainSingleCrystalBrick::"
237                         "completeVariableDeclaration: begin\n";
238     }
239     // various checks
240     throw_if(!this->bd.hasCrystalStructure(), "no crystal structure defined");
241     throw_if(!this->bd.areSlipSystemsDefined(),"no slip systems defined");
242    // defining the stiffness tensor, if necessary
243     if((!this->bd.getAttribute(BehaviourDescription::requiresStiffnessTensor,false))&&
244        (!this->bd.getAttribute(BehaviourDescription::computesStiffnessTensor,false))){
245       this->bd.setAttribute(BehaviourDescription::requiresStiffnessTensor,true,false);
246     }
247     LocalDataStructure d;
248     d.name = "fsscb_data";
249     // local data
250     d.addVariable(h,{"DeformationGradientTensor","dF"});
251     d.addVariable(h,{"DeformationGradientTensor","Fe_tr"});
252     d.addVariable(h,{"DeformationGradientTensor","Fe0"});
253     d.addVariable(h,{"StressStensor","S"});
254     d.addVariable(h,{"Tensor","inv_dFp"});
255     d.addVariable(h,{"real","J_inv_dFp"});
256     d.addVariable(h,{"StrainStensor","tmp"});
257     this->bd.addLocalDataStructure(d,BehaviourData::ALREADYREGISTRED);
258     if (getVerboseMode() >= VERBOSE_DEBUG) {
259       getLogStream() << "FiniteStrainSingleCrystalBrick::"
260                         "completeVariableDeclaration: end\n";
261     }
262   }
263 
endTreatment() const264   void FiniteStrainSingleCrystalBrick::endTreatment() const {
265     const auto h  = ModellingHypothesis::TRIDIMENSIONAL;
266     const auto cn = this->bd.getClassName()+"SlipSystems<real>";
267     // local data values initialisation
268     CodeBlock init;
269     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
270       init.code = "this->Fe += DeformationGradientTensor::Id();\n";
271     }
272     init.code +=
273       "this->fsscb_data.dF    = (this->F1)*invert(this->F0);\n"
274       "this->fsscb_data.Fe0   = this->Fe;\n"
275       "this->fsscb_data.Fe_tr = (this->fsscb_data.dF)*(this->fsscb_data.Fe0);\n"
276       "this->eel = computeGreenLagrangeTensor(this->Fe);\n";
277     init.members  = {"F0","F1","Fe","eel"};
278     this->bd.setCode(h,BehaviourData::BeforeInitializeLocalVariables,
279     		     init,BehaviourData::CREATEORAPPEND,
280     		     BehaviourData::AT_END);
281     CodeBlock integrator;
282     integrator.code =
283       "const auto& ss = "+cn+"::getSlidingSystems();\n"
284       "this->fsscb_data.S = (this->D)*(this->eel+this->deel);\n"
285       "this->fsscb_data.tmp = StrainStensor::Id() + 2*(this->eel+this->deel);\n"
286       "// Mandel stress tensor\n"
287       "const auto M = eval(this->fsscb_data.tmp*(this->fsscb_data.S));\n";
288     const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(this->dsl);
289     if((idsl.getSolver().usesJacobian())&&(!idsl.getSolver().requiresNumericalJacobian())){
290       integrator.code +=
291 	"// Mandel stress tensor derivative\n"
292 	"const auto dM_ddeel = eval(2*st2tot2<N,real>::tpld(this->fsscb_data.S)+\n"
293 	"			    st2tot2<N,real>::tprd(this->fsscb_data.tmp,this->D));\n";
294     }
295     integrator.code +=
296       "this->fsscb_data.inv_dFp = Tensor::Id();\n"
297       "for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
298       "  this->fsscb_data.inv_dFp -= (this->dg[i])*ss.mu[i];\n"
299       "}\n"
300       "this->fsscb_data.J_inv_dFp = det(this->fsscb_data.inv_dFp);\n"
301       "this->fsscb_data.inv_dFp /= CubicRoots::cbrt(this->fsscb_data.J_inv_dFp);\n"
302       "this->Fe = (this->fsscb_data.Fe_tr)*(this->fsscb_data.inv_dFp);\n"
303       "feel = this->eel+this->deel-computeGreenLagrangeTensor(this->Fe);\n";
304     if((idsl.getSolver().usesJacobian())&&(!idsl.getSolver().requiresNumericalJacobian())){
305       integrator.code +=
306 	"const auto fsscb_tprd = t2tot2<N,real>::tprd(this->fsscb_data.Fe_tr);\n"
307 	"const auto fsscb_dfeel_dinv_dFp = t2tost2<N,real>::dCdF(this->Fe)*fsscb_tprd;\n"
308 	"for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
309 	"  dfeel_ddg(i) = (fsscb_dfeel_dinv_dFp)*ss.mu[i]/2;\n"
310 	"}\n";
311     }
312     integrator.members  = {"eel","Fe","D"};
313     this->bd.setCode(h,BehaviourData::Integrator,
314     		     integrator,BehaviourData::CREATEORAPPEND,
315     		     BehaviourData::AT_BEGINNING);
316     CodeBlock fs;
317     fs.code =
318       "const auto& ss = "+cn+"::getSlidingSystems();\n"
319       "this->fsscb_data.inv_dFp = Tensor::Id();\n"
320       "for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
321       "  this->fsscb_data.inv_dFp -= dg[i]*ss.mu[i];\n"
322       "}\n"
323       "this->fsscb_data.J_inv_dFp = det(this->fsscb_data.inv_dFp);\n"
324       "this->fsscb_data.inv_dFp /= CubicRoots::cbrt(this->fsscb_data.J_inv_dFp);\n"
325       "this->Fe = (this->fsscb_data.Fe_tr)*(this->fsscb_data.inv_dFp);\n"
326       "this->fsscb_data.S = (this->D)*(this->eel);\n"
327       "this->sig = convertSecondPiolaKirchhoffStressToCauchyStress(this->fsscb_data.S,this->Fe);\n";
328     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
329       fs.code +=
330 	"this->Fe -= DeformationGradientTensor::Id();\n";
331     }
332     fs.members  = {"sig","Fe","D"};
333     this->bd.setCode(h,BehaviourData::ComputeFinalThermodynamicForces,
334     		     fs,BehaviourData::CREATE,
335     		     BehaviourData::AT_BEGINNING);
336     // tangent operator
337     CodeBlock to;
338     if (idsl.getSolver().usesJacobian()) {
339       to.attributes["requires_jacobian_decomposition"] = true;
340       to.attributes["uses_get_partial_jacobian_invert"] = true;
341     }
342     to.code =
343       "static_cast<void>(smt);\n"
344       "const auto& ss = "+cn+"::getSlidingSystems();\n";
345     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
346       to.code += "const auto fsscb_dC_dFe = t2tost2<N,real>::dCdF(this->Fe+DeformationGradientTensor::Id());\n";
347     } else {
348       to.code += "const auto fsscb_dC_dFe = t2tost2<N,real>::dCdF(this->Fe);\n";
349     }
350     to.code +=
351       "const auto fsscb_dS_dFe = eval((this->D)*fsscb_dC_dFe/2);\n";
352     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
353       to.code +=
354 	"const auto fsscb_dtau_dFe = "
355 	"computePushForwardDerivative(fsscb_dS_dFe,this->fsscb_data.S,this->Fe+DeformationGradientTensor::Id());\n";
356     } else {
357       to.code +=
358 	"const auto fsscb_dtau_dFe = "
359 	"computePushForwardDerivative(fsscb_dS_dFe,this->fsscb_data.S,this->Fe);\n";
360     }
361     to.code +=
362       "const auto fsscb_dFe_dDF_tot = "
363       "t2tot2<N,real>::tpld(this->fsscb_data.inv_dFp,"
364       "                     t2tot2<N,real>::tpld(this->fsscb_data.Fe0));\n"
365       "const auto fsscb_dfeel_dDF = eval(-(fsscb_dC_dFe)*(fsscb_dFe_dDF_tot)/2);\n"
366       "st2tost2<N,real> fsscb_Je;\n"
367       "tvector<"+cn+"::Nss,Stensor> fsscb_Jg;\n"
368       "getPartialJacobianInvert(fsscb_Je,fsscb_Jg);\n"
369       "t2tot2<N,real> fsscb_dinv_Fp_dDF = (ss.mu[0])^(fsscb_Jg[0]|fsscb_dfeel_dDF);\n"
370       "for(unsigned short i=1;i!="+cn+"::Nss;++i){\n"
371       "  fsscb_dinv_Fp_dDF += (ss.mu[i])^(fsscb_Jg[i]|fsscb_dfeel_dDF);\n"
372       "}\n"
373       "const auto fsscb_dFe_dDF=\n"
374       "  fsscb_dFe_dDF_tot+t2tot2<N,real>::tprd(this->fsscb_data.Fe_tr,fsscb_dinv_Fp_dDF);\n"
375       "Dt = fsscb_dtau_dFe*fsscb_dFe_dDF;\n";
376     to.members  = {"Fe","D"};
377     const auto ton = std::string(BehaviourData::ComputeTangentOperator)+"-DTAU_DDF";
378     this->bd.setCode(h,ton,to,BehaviourData::CREATE,
379     		     BehaviourData::AT_BEGINNING);
380   } // end of FiniteStrainSingleCrystalBrick::endTreatment
381 
382   std::vector<tfel::material::ModellingHypothesis::Hypothesis>
getSupportedModellingHypotheses() const383   FiniteStrainSingleCrystalBrick::getSupportedModellingHypotheses() const
384   {
385     return {ModellingHypothesis::TRIDIMENSIONAL};
386   } // end of FiniteStrainSingleCrystalBrick::getSupportedModellingHypothesis
387 
388   FiniteStrainSingleCrystalBrick::~FiniteStrainSingleCrystalBrick() = default;
389 
390 } // end of namespace mfront
391