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