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/FiniteStrainSingleCrystalBrick.hxx"
32 
33 namespace mfront{
34 
35   const char* const FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption =
36     "shifted_elastic_deformation_gradient";
37 
38   const char* const FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute =
39     "FiniteStrainSingleCrystalBrick::shifted_elastic_deformation_gradient";
40 
FiniteStrainSingleCrystalBrick(AbstractBehaviourDSL & dsl_,BehaviourDescription & mb_,const Parameters & p,const DataMap & d)41   FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick(AbstractBehaviourDSL& dsl_,
42 								 BehaviourDescription& mb_,
43 								 const Parameters& p,
44 								 const DataMap& d)
45     : BehaviourBrickBase(dsl_,mb_)
46   {
47     auto throw_if = [](const bool b,const std::string& m){
48       tfel::raise_if(b,"FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick: "+m);
49     };
50     const auto h = ModellingHypothesis::TRIDIMENSIONAL;
51     // basic checks
52     if(this->bd.areModellingHypothesesDefined()){
53       const auto bmh = this->bd.getModellingHypotheses();
54       throw_if(bmh.size()!=1u,"the only supported hypothesis is 'Tridimensional'");
55       throw_if(*(bmh.begin())!=ModellingHypothesis::TRIDIMENSIONAL,
56 	       "the only supported hypothesis is 'Tridimensional'");
57     } else {
58       this->bd.setModellingHypotheses({ModellingHypothesis::TRIDIMENSIONAL});
59     }
60     throw_if(this->bd.getBehaviourType()!=BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR,
61 	     "this BehaviourBrick is only usable for small strain behaviours");
62     throw_if(this->bd.getIntegrationScheme()!=BehaviourDescription::IMPLICITSCHEME,
63 	     "this BehaviourBrick is only usable in implicit schemes");
64     throw_if(!p.empty(),"no parameters allowed");
65     // material symmetry
66     if(this->bd.isSymmetryTypeDefined()){
67       throw_if(this->bd.getSymmetryType()!=mfront::ORTHOTROPIC,
68 	       "material must be declared orthotropic");
69     } else {
70       this->bd.setSymmetryType(mfront::ORTHOTROPIC);
71     }
72     throw_if(this->bd.getElasticSymmetryType()!=mfront::ORTHOTROPIC,
73 	     "elastic symmetry type must be orthortropic");
74     // declaring the elastic strain as the first integration variable
75     throw_if(!this->bd.getBehaviourData(h).getIntegrationVariables().empty(),
76 	     "no integration variable shall be declared before declaring the "
77 	     "'FiniteStrainSingleCrystal' brick");
78     // options
79     this->checkOptionsNames(d,{FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption},
80 			    this->getName());
81     if(d.count(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption)){
82       const auto on = std::string(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientOption);
83       const auto an = std::string(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute);
84       const auto& b = d.at(on);
85       throw_if(!b.is<bool>(),"invalid type for option '"+on+"'");
86       throw_if(this->bd.hasAttribute(an),
87 	       "attribute '"+an+"' already used");
88       this->bd.setAttribute(an,b.get<bool>(),false);
89     }
90     // elastic strain
91     VariableDescription eel("StrainStensor","eel",1u,0u);
92     eel.description = "elastic strain";
93     this->bd.addIntegrationVariable(h,eel,BehaviourData::UNREGISTRED);
94     this->bd.setGlossaryName(h,"eel",tfel::glossary::Glossary::ElasticStrain);
95     // declaring the elastic part of the deformation gradient
96     VariableDescription Fe("DeformationGradientTensor","Fe",1u,0u);
97     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
98       Fe.description = "shifted elastic part of the deformation gradient";
99     } else {
100       Fe.description = "elastic part of the deformation gradient";
101     }
102     this->bd.addAuxiliaryStateVariable(h,Fe,BehaviourData::UNREGISTRED);
103     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
104       this->bd.setEntryName(h,"Fe","ShiftedElasticPartOfTheDeformationGradient");
105     } else {
106       this->bd.setEntryName(h,"Fe","ElasticPartOfTheDeformationGradient");
107     }
108     // additional includes
109     this->bd.appendToIncludes("#include\"TFEL/Math/General/CubicRoots.hxx\"");
110     // set that the tangent operator is computed
111     this->bd.setAttribute(h,BehaviourData::hasConsistentTangentOperator,true);
112     // reserve some specific variables
113     this->bd.reserveName(h,"ss");
114     this->bd.registerMemberName(h,"g");
115     this->bd.reserveName(h,"fg");
116     this->bd.reserveName(h,"fsscb_data");
117     this->bd.reserveName(h,"fsscb_tprd");
118     this->bd.reserveName(h,"fsscb_dfeel_dinv_dFp");
119     this->bd.reserveName(h,"fsscb_dC_dFe");
120     this->bd.reserveName(h,"fsscb_dS_dFe");
121     this->bd.reserveName(h,"fsscb_dtau_dFe");
122     this->bd.reserveName(h,"fsscb_dFe_dDF_tot");
123     this->bd.reserveName(h,"fsscb_dfeel_dDF");
124     this->bd.reserveName(h,"fsscb_Je");
125     this->bd.reserveName(h,"fsscb_Jg");
126     this->bd.reserveName(h,"fsscb_dinv_Fp_dDF");
127     this->bd.reserveName(h,"fsscb_dFe_dDF");
128   } // end of FiniteStrainSingleCrystalBrick::FiniteStrainSingleCrystalBrick
129 
completeVariableDeclaration() const130   void FiniteStrainSingleCrystalBrick::completeVariableDeclaration() const
131   {
132     using tfel::glossary::Glossary;
133     auto throw_if = [](const bool b,const std::string& m){
134       tfel::raise_if(b,"FiniteStrainSingleCrystalBrick:"
135 		     ":completeVariableDeclaration: "+m);
136     };
137     const auto h = ModellingHypothesis::TRIDIMENSIONAL;
138     if(getVerboseMode()>=VERBOSE_DEBUG){
139       getLogStream() << "FiniteStrainSingleCrystalBrick::completeVariableDeclaration: begin\n";
140     }
141     // defining the stiffness tensor, if necessary
142     if((!this->bd.getAttribute(BehaviourDescription::requiresStiffnessTensor,false))&&
143        (!this->bd.getAttribute(BehaviourDescription::computesStiffnessTensor,false))){
144       this->bd.setAttribute(BehaviourDescription::requiresStiffnessTensor,true,false);
145     }
146     LocalDataStructure d;
147     d.name = "fsscb_data";
148     // local data
149     d.addVariable(h,{"DeformationGradientTensor","dF"});
150     d.addVariable(h,{"DeformationGradientTensor","Fe_tr"});
151     d.addVariable(h,{"DeformationGradientTensor","Fe0"});
152     d.addVariable(h,{"StressStensor","S"});
153     d.addVariable(h,{"Tensor","inv_dFp"});
154     d.addVariable(h,{"real","J_inv_dFp"});
155     d.addVariable(h,{"StrainStensor","tmp"});
156     this->bd.addLocalDataStructure(d,BehaviourData::ALREADYREGISTRED);
157     // various checks
158     throw_if(!this->bd.hasCrystalStructure(),"no crystal structure defined");
159     throw_if(!this->bd.areSlipSystemsDefined(),"no slip systems defined");
160     const auto& ss = this->bd.getSlipSystems();
161     // declaring the plastic slip
162     VariableDescription g("strain","g",ss.getNumberOfSlipSystems(),0u);
163     g.description = "plastic slip";
164     this->bd.addStateVariable(h,g,BehaviourData::ALREADYREGISTRED);
165     this->bd.setEntryName(h,"g","PlasticSlip");
166     if(getVerboseMode()>=VERBOSE_DEBUG){
167       getLogStream() << "FiniteStrainSingleCrystalBrick::completeVariableDeclaration: end\n";
168     }
169   }
170 
endTreatment() const171   void FiniteStrainSingleCrystalBrick::endTreatment() const
172   {
173     const auto h  = ModellingHypothesis::TRIDIMENSIONAL;
174     const auto cn = this->bd.getClassName()+"SlipSystems<real>";
175     // local data values initialisation
176     CodeBlock init;
177     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
178       init.code = "this->Fe += DeformationGradientTensor::Id();\n";
179     }
180     init.code +=
181       "this->fsscb_data.dF    = (this->F1)*invert(this->F0);\n"
182       "this->fsscb_data.Fe0   = this->Fe;\n"
183       "this->fsscb_data.Fe_tr = (this->fsscb_data.dF)*(this->fsscb_data.Fe0);\n"
184       "this->eel = computeGreenLagrangeTensor(this->Fe);\n";
185     init.members  = {"F0","F1","Fe","eel"};
186     this->bd.setCode(h,BehaviourData::BeforeInitializeLocalVariables,
187     		     init,BehaviourData::CREATEORAPPEND,
188     		     BehaviourData::AT_END);
189     CodeBlock integrator;
190     integrator.code =
191       "const auto& ss = "+cn+"::getSlidingSystems();\n"
192       "this->fsscb_data.S = (this->D)*(this->eel+this->deel);\n"
193       "this->fsscb_data.tmp = StrainStensor::Id() + 2*(this->eel+this->deel);\n"
194       "// Mandel stress tensor\n"
195       "const auto M = eval(this->fsscb_data.tmp*(this->fsscb_data.S));\n";
196     const auto& idsl = dynamic_cast<const ImplicitDSLBase&>(this->dsl);
197     if((idsl.getSolver().usesJacobian())&&(!idsl.getSolver().requiresNumericalJacobian())){
198       integrator.code +=
199 	"// Mandel stress tensor derivative\n"
200 	"const auto dM_ddeel = eval(2*st2tot2<N,real>::tpld(this->fsscb_data.S)+\n"
201 	"			    st2tot2<N,real>::tprd(this->fsscb_data.tmp,this->D));\n";
202     }
203     integrator.code +=
204       "this->fsscb_data.inv_dFp = Tensor::Id();\n"
205       "for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
206       "  this->fsscb_data.inv_dFp -= (this->dg[i])*ss.mu[i];\n"
207       "}\n"
208       "this->fsscb_data.J_inv_dFp = det(this->fsscb_data.inv_dFp);\n"
209       "this->fsscb_data.inv_dFp /= CubicRoots::cbrt(this->fsscb_data.J_inv_dFp);\n"
210       "this->Fe = (this->fsscb_data.Fe_tr)*(this->fsscb_data.inv_dFp);\n"
211       "feel = this->eel+this->deel-computeGreenLagrangeTensor(this->Fe);\n";
212     if((idsl.getSolver().usesJacobian())&&(!idsl.getSolver().requiresNumericalJacobian())){
213       integrator.code +=
214 	"const auto fsscb_tprd = t2tot2<N,real>::tprd(this->fsscb_data.Fe_tr);\n"
215 	"const auto fsscb_dfeel_dinv_dFp = t2tost2<N,real>::dCdF(this->Fe)*fsscb_tprd;\n"
216 	"for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
217 	"  dfeel_ddg(i) = (fsscb_dfeel_dinv_dFp)*ss.mu[i]/2;\n"
218 	"}\n";
219     }
220     integrator.members  = {"eel","Fe","D"};
221     this->bd.setCode(h,BehaviourData::Integrator,
222     		     integrator,BehaviourData::CREATEORAPPEND,
223     		     BehaviourData::AT_BEGINNING);
224     CodeBlock fs;
225     fs.code =
226       "const auto& ss = "+cn+"::getSlidingSystems();\n"
227       "this->fsscb_data.inv_dFp = Tensor::Id();\n"
228       "for(unsigned short i=0;i!="+cn+"::Nss;++i){\n"
229       "  this->fsscb_data.inv_dFp -= dg[i]*ss.mu[i];\n"
230       "}\n"
231       "this->fsscb_data.J_inv_dFp = det(this->fsscb_data.inv_dFp);\n"
232       "this->fsscb_data.inv_dFp /= CubicRoots::cbrt(this->fsscb_data.J_inv_dFp);\n"
233       "this->Fe = (this->fsscb_data.Fe_tr)*(this->fsscb_data.inv_dFp);\n"
234       "this->fsscb_data.S = (this->D)*(this->eel);\n"
235       "this->sig = convertSecondPiolaKirchhoffStressToCauchyStress(this->fsscb_data.S,this->Fe);\n";
236     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
237       fs.code +=
238 	"this->Fe -= DeformationGradientTensor::Id();\n";
239     }
240     fs.members  = {"sig","Fe","D"};
241     this->bd.setCode(h,BehaviourData::ComputeFinalStress,
242     		     fs,BehaviourData::CREATE,
243     		     BehaviourData::AT_BEGINNING);
244     // tangent operator
245     CodeBlock to;
246     to.code =
247       "static_cast<void>(smt);\n"
248       "const auto& ss = "+cn+"::getSlidingSystems();\n";
249     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
250       to.code += "const auto fsscb_dC_dFe = t2tost2<N,real>::dCdF(this->Fe+DeformationGradientTensor::Id());\n";
251     } else {
252       to.code += "const auto fsscb_dC_dFe = t2tost2<N,real>::dCdF(this->Fe);\n";
253     }
254     to.code +=
255       "const auto fsscb_dS_dFe = eval((this->D)*fsscb_dC_dFe/2);\n";
256     if(this->bd.getAttribute<bool>(FiniteStrainSingleCrystalBrick::shiftedDeformationGradientAttribute,false)){
257       to.code +=
258 	"const auto fsscb_dtau_dFe = "
259 	"computePushForwardDerivative(fsscb_dS_dFe,this->fsscb_data.S,this->Fe+DeformationGradientTensor::Id());\n";
260     } else {
261       to.code +=
262 	"const auto fsscb_dtau_dFe = "
263 	"computePushForwardDerivative(fsscb_dS_dFe,this->fsscb_data.S,this->Fe);\n";
264     }
265     to.code +=
266       "const auto fsscb_dFe_dDF_tot = "
267       "t2tot2<N,real>::tpld(this->fsscb_data.inv_dFp,"
268       "                     t2tot2<N,real>::tpld(this->fsscb_data.Fe0));\n"
269       "const auto fsscb_dfeel_dDF = eval(-(fsscb_dC_dFe)*(fsscb_dFe_dDF_tot)/2);\n"
270       "st2tost2<N,real> fsscb_Je;\n"
271       "tvector<"+cn+"::Nss,Stensor> fsscb_Jg;\n"
272       "getPartialJacobianInvert(fsscb_Je,fsscb_Jg);\n"
273       "t2tot2<N,real> fsscb_dinv_Fp_dDF = (ss.mu[0])^(fsscb_Jg[0]|fsscb_dfeel_dDF);\n"
274       "for(unsigned short i=1;i!="+cn+"::Nss;++i){\n"
275       "  fsscb_dinv_Fp_dDF += (ss.mu[i])^(fsscb_Jg[i]|fsscb_dfeel_dDF);\n"
276       "}\n"
277       "const auto fsscb_dFe_dDF=\n"
278       "  fsscb_dFe_dDF_tot+t2tot2<N,real>::tprd(this->fsscb_data.Fe_tr,fsscb_dinv_Fp_dDF);\n"
279       "Dt = fsscb_dtau_dFe*fsscb_dFe_dDF;\n";
280     to.members  = {"Fe","D"};
281     const auto ton = std::string(BehaviourData::ComputeTangentOperator)+"-DTAU_DDF";
282     this->bd.setCode(h,ton,to,BehaviourData::CREATE,
283     		     BehaviourData::AT_BEGINNING);
284   } // end of FiniteStrainSingleCrystalBrick::endTreatment
285 
getName() const286   std::string FiniteStrainSingleCrystalBrick::getName() const{
287     return "FiniteStrainSingleCrystal";
288   }
289 
290   std::vector<tfel::material::ModellingHypothesis::Hypothesis>
getSupportedModellingHypotheses() const291   FiniteStrainSingleCrystalBrick::getSupportedModellingHypotheses() const
292   {
293     return {ModellingHypothesis::TRIDIMENSIONAL};
294   } // end of FiniteStrainSingleCrystalBrick::getSupportedModellingHypothesis
295 
296   FiniteStrainSingleCrystalBrick::~FiniteStrainSingleCrystalBrick() = default;
297 
298 } // end of namespace mfront
299