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