1 /*! 2 * \file mfront/src/MultipleIsotropicMisesFlowsDSL.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 31 jan 2008 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<string> 15 #include<sstream> 16 #include<stdexcept> 17 #include"TFEL/Raise.hxx" 18 #include"MFront/DSLUtilities.hxx" 19 #include"MFront/MFrontDebugMode.hxx" 20 #include"MFront/DSLFactory.hxx" 21 #include"MFront/MultipleIsotropicMisesFlowsDSL.hxx" 22 23 namespace mfront{ 24 25 MultipleIsotropicMisesFlowsDSL::FlowHandler::~FlowHandler() noexcept = default; 26 MultipleIsotropicMisesFlowsDSL()27 MultipleIsotropicMisesFlowsDSL::MultipleIsotropicMisesFlowsDSL() 28 { 29 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 30 this->mb.setDSLName("MultipleIsotropicMisesFlows"); 31 // Default state vars 32 this->mb.addStateVariable(h,VariableDescription("StrainStensor","eel",1u,0u)); 33 this->mb.addStateVariable(h,VariableDescription("strain","p",1u,0u)); 34 this->mb.setGlossaryName(h,"eel","ElasticStrain"); 35 this->mb.setGlossaryName(h,"p","EquivalentStrain"); 36 // default local vars 37 this->reserveName("mu_3_theta"); 38 this->reserveName("surf"); 39 this->mb.addLocalVariable(h,VariableDescription("StressStensor","se",1u,0u)); 40 this->mb.addLocalVariable(h,VariableDescription("stress","seq",1u,0u)); 41 this->mb.addLocalVariable(h,VariableDescription("stress","seq_e",1u,0u)); 42 this->mb.addLocalVariable(h,VariableDescription("StrainStensor","n",1u,0u)); 43 this->mb.addLocalVariable(h,VariableDescription("strain","p_",1u,0u)); 44 } 45 getName()46 std::string MultipleIsotropicMisesFlowsDSL::getName() 47 { 48 return "MultipleIsotropicMisesFlows"; 49 } 50 getDescription()51 std::string MultipleIsotropicMisesFlowsDSL::getDescription() 52 { 53 return "this parser is used to define behaviours combining several " 54 "isotropic flows. Supported flow type are 'Creep' (dp/dt=f(s)) " 55 "'StrainHardeningCreep' (dp/dt=f(s,p)) and 'Plasticity' (f(p,s)=0) " 56 "where p is the equivalent plastic strain and s the equivalent " 57 "mises stress"; 58 } // end of MultipleIsotropicMisesFlowsDSL::getDescription 59 writeBehaviourParserSpecificIncludes(std::ostream & os) const60 void MultipleIsotropicMisesFlowsDSL::writeBehaviourParserSpecificIncludes(std::ostream& os) const 61 { 62 this->checkBehaviourFile(os); 63 os << "#include\"TFEL/Math/General/BaseCast.hxx\"\n" 64 << "#include\"TFEL/Math/TinyMatrixSolve.hxx\"\n\n"; 65 } 66 67 void writeBehaviourParserSpecificMembers(std::ostream & os,const Hypothesis) const68 MultipleIsotropicMisesFlowsDSL::writeBehaviourParserSpecificMembers(std::ostream& os, 69 const Hypothesis) const 70 { 71 using namespace std; 72 vector<FlowHandler>::const_iterator p; 73 vector<FlowHandler>::const_iterator p2; 74 unsigned short n; 75 bool genericTheta; 76 this->checkBehaviourFile(os); 77 tfel::raise_if(this->flows.empty(),"MultipleIsotropicMisesFlowsDSL::" 78 "writeBehaviourParserSpecificMembers : " 79 "no flow rule defined"); 80 for(p=this->flows.begin(),n=0;p!=this->flows.end();++p,++n){ 81 if(p->flow==FlowHandler::PlasticFlow){ 82 os << "void computeFlow" << n << "(stress& f,\n" 83 << "real& df_dseq,\n" 84 << "stress& df_dp){\n"; 85 } else if(p->flow==FlowHandler::CreepFlow){ 86 os << "void computeFlow" << n << "(DstrainDt& f,\n" 87 << "DF_DSEQ_TYPE& df_dseq){\n"; 88 } else if(p->flow==FlowHandler::StrainHardeningCreepFlow){ 89 os << "void computeFlow" << n << "(DstrainDt& f,\n" 90 << "DF_DSEQ_TYPE& df_dseq,\n" 91 << "DstrainDt& df_dp){\n"; 92 } 93 os << "using namespace std;\n"; 94 os << "using namespace tfel::math;\n"; 95 os << "using namespace tfel::material;\n"; 96 os << "using std::vector;\n"; 97 writeMaterialLaws(os,this->mb.getMaterialLaws()); 98 os << p->flowRule << endl; 99 os << "}\n\n"; 100 } 101 os << "bool NewtonIntegration(){\n" 102 << "using namespace std;\n" 103 << "using namespace tfel::math;\n" 104 << "tvector<" << this->flows.size() << ",strain> vdp(strain(real(0.)));\n" 105 << "tvector<" << this->flows.size() << ",strain> newton_f;\n" 106 << "tmatrix<" << this->flows.size() 107 << "," << this->flows.size() << ",strain> newton_df;\n"; 108 109 genericTheta=false; 110 for(p=this->flows.begin(),n=0;p!=this->flows.end();++p,++n){ 111 if(p->hasSpecificTheta){ 112 ostringstream otheta; 113 otheta << "mu_3_theta" << n; 114 os << "stress "+otheta.str()+" = 3*(real("; 115 os << p->theta << "))*(this->mu);\n"; 116 } else { 117 genericTheta=true; 118 } 119 } 120 121 if(genericTheta){ 122 os << "stress mu_3_theta = 3*("; 123 os << this->mb.getClassName() << "::theta)*(this->mu);\n"; 124 } 125 bool found=false; 126 for(p=this->flows.begin();(p!=this->flows.end())&&!(found);++p){ 127 if(p->flow==FlowHandler::PlasticFlow){ 128 os << "real surf;\n"; 129 os << "real newton_epsilon = 100*std::numeric_limits<real>::epsilon();\n"; 130 found = true; 131 } 132 } 133 os << "unsigned int iter=0u;\n" 134 << "bool converge=false;\n" 135 << "while((converge==false)&&\n" 136 << "(iter<(" << this->mb.getClassName() << "::iterMax))){\n"; 137 for(p=this->flows.begin(),n=0;p!=this->flows.end();++p,++n){ 138 os << "this->p_ = this->p" << n << " + ("; 139 if(p->hasSpecificTheta){ 140 os << "real(" << p->theta << ")"; 141 } else { 142 os << this->mb.getClassName() << "::theta"; 143 } 144 os << ")*(vdp(" << n << "));\n"; 145 if(p->hasSpecificTheta){ 146 ostringstream otheta; 147 os << "this->seq = std::max(this->seq_e" << n << "-"; 148 otheta << "mu_3_theta" << n << "*("; 149 os << otheta.str(); 150 } else { 151 os << "this->seq = std::max(this->seq_e-"; 152 os << "mu_3_theta*("; 153 } 154 p2=this->flows.begin(); 155 unsigned short n2 = 0u; 156 while(p2!=this->flows.end()){ 157 os << "vdp(" << n2 << ")"; 158 ++p2; 159 ++n2; 160 if(p2!=this->flows.end()){ 161 os << "+"; 162 } 163 } 164 os << "),real(0.f));\n"; 165 if(p->flow==FlowHandler::PlasticFlow){ 166 os << "this->computeFlow" << n << "(" 167 << "this->f" << n << "," 168 << "this->df_dseq" << n << "," 169 << "this->df_dp" << n << ");\n"; 170 os << "surf = (this->f" << n << ")/(this->young);\n"; 171 os << "if(((surf>newton_epsilon)&&((vdp(" << n << "))>=0))||" 172 << "((vdp(" << n << "))>newton_epsilon)){"; 173 os << "newton_f(" << n << ") = surf;\n"; 174 for(p2=this->flows.begin(),n2=0;p2!=this->flows.end();++p2,++n2){ 175 if(p2==p){ 176 os << "newton_df(" << n << "," << n << ")"; 177 if(p->hasSpecificTheta){ 178 os << " = ((real("<< p->theta << "))*(this->df_dp" << n << ")" 179 << "-mu_3_theta" << n << "*(this->df_dseq" << n <<"))/(this->young);\n"; 180 } else { 181 os << " = (("<< this->mb.getClassName() << "::theta)*(this->df_dp" << n << ")" 182 << "-mu_3_theta*(this->df_dseq" << n <<"))/(this->young);\n"; 183 } 184 } else { 185 os << "newton_df(" << n << "," << n2 << ")"; 186 if(p->hasSpecificTheta){ 187 os << " = -mu_3_theta" << n << "*(this->df_dseq" << n <<")/(this->young);\n"; 188 } else { 189 os << " = -mu_3_theta*(this->df_dseq" << n <<")/(this->young);\n"; 190 } 191 } 192 } 193 os << "} else {\n"; 194 os << "newton_f(" << n << ") =(vdp(" << n << "));\n"; 195 os << "newton_df(" << n << "," << n << ") = real(1.);\n"; 196 for(p2=this->flows.begin(),n2=0;p2!=this->flows.end();++p2,++n2){ 197 if(p2!=p){ 198 os << "newton_df(" << n << "," << n2 << ") = real(0.);\n"; 199 } 200 } 201 os << "}\n"; 202 } else if (p->flow==FlowHandler::CreepFlow){ 203 os << "this->computeFlow" << n << "(" 204 << "this->f" << n << "," 205 << "this->df_dseq" << n << "" 206 << ");\n"; 207 os << "newton_f(" << n << ") = vdp(" << n << ") - (this->f" << n << ")*(this->dt);\n"; 208 os << "newton_df(" << n << "," << n << ") = 1+"; 209 if(p->hasSpecificTheta){ 210 os << "mu_3_theta" << n; 211 } else { 212 os << "mu_3_theta"; 213 } 214 os << "*(this->df_dseq" << n << ")*(this->dt);\n"; 215 for(p2=this->flows.begin(),n2=0;p2!=this->flows.end();++p2,++n2){ 216 if(p2!=p){ 217 os << "newton_df(" << n << "," << n2 << ") = "; 218 if(p->hasSpecificTheta){ 219 os << "mu_3_theta" << n; 220 } else { 221 os << "mu_3_theta"; 222 } 223 os << "*(this->df_dseq" << n << ")*(this->dt);\n"; 224 } 225 } 226 } else { 227 os << "this->computeFlow" << n << "(" 228 << "this->f" << n << "," 229 << "this->df_dseq" << n << "," 230 << "this->df_dp" << n << ");\n"; 231 os << "newton_f(" << n << ") = vdp(" << n << ") - (this->f" << n << ")*(this->dt);\n"; 232 os << "newton_df(" << n << "," << n << ") = 1-(this->dt)*("; 233 if(p->hasSpecificTheta){ 234 os << "(real(" << p->theta << "))"; 235 } else { 236 os << "(" << this->mb.getClassName() << "::theta)"; 237 } 238 os << "*(this->df_dp" << n << ")-"; 239 if(p->hasSpecificTheta){ 240 os << "mu_3_theta" << n; 241 } else { 242 os << "mu_3_theta"; 243 } 244 os << "*(this->df_dseq" << n << "));\n"; 245 for(p2=this->flows.begin(),n2=0;p2!=this->flows.end();++p2,++n2){ 246 if(p2!=p){ 247 os << "newton_df(" << n << "," << n2 << ") = "; 248 if(p->hasSpecificTheta){ 249 os << "mu_3_theta" << n; 250 } else { 251 os << "mu_3_theta"; 252 } 253 os << "*(this->df_dseq" << n << ")*(this->dt);\n"; 254 } 255 } 256 } 257 } 258 os << "real error=static_cast<real>(0.);\n"; 259 for(p=this->flows.begin(),n=0;p!=this->flows.end();++p,++n){ 260 os << "error+=std::abs(tfel::math::base_cast(newton_f(" << n << ")));\n"; 261 } 262 os << "try{" << endl 263 << "TinyMatrixSolve<" << this->flows.size() << "," << "real>::exe(newton_df,newton_f);\n" 264 << "} catch(LUException&){" << endl 265 << "return false;" << endl 266 << "}" << endl 267 << "vdp -= newton_f;\n" 268 << "iter+=1;\n"; 269 if(getDebugMode()){ 270 os << "cout << \"" << this->mb.getClassName() 271 << "::NewtonIntegration() : iteration \" " 272 << "<< iter << \" : \" << (error/(real(" << this->flows.size() << "))) << endl;\n"; 273 } 274 os << "converge = ((error)/(real(" << this->flows.size() << "))<" 275 << "(" << this->mb.getClassName() << "::epsilon));\n" 276 << "}\n\n" 277 << "if(iter==" << this->mb.getClassName() << "::iterMax){\n"; 278 if(getDebugMode()){ 279 os << "cout << \"" << this->mb.getClassName() 280 << "::NewtonIntegration() : no convergence after \" " 281 << "<< iter << \" iterations\"<< endl << endl;\n"; 282 os << "cout << *this << endl;\n"; 283 } 284 os << "return false;" << endl 285 << "}\n\n"; 286 for(p=this->flows.begin(),n=0;p!=this->flows.end();++p,++n){ 287 os << "this->dp"<< n << " = " << "vdp(" << n<< ");\n"; 288 } 289 if(getDebugMode()){ 290 os << "cout << \"" << this->mb.getClassName() 291 << "::NewtonIntegration() : convergence after \" " 292 << "<< iter << \" iterations\"<< endl << endl;\n"; 293 } 294 os << "return true;" << endl 295 << "\n}\n\n"; 296 } // end of writeBehaviourParserSpecificMembers 297 writeBehaviourIntegrator(std::ostream & os,const Hypothesis h) const298 void MultipleIsotropicMisesFlowsDSL::writeBehaviourIntegrator(std::ostream& os, 299 const Hypothesis h) const 300 { 301 const auto btype = this->mb.getBehaviourTypeFlag(); 302 const auto& d = this->mb.getBehaviourData(h); 303 unsigned short n; 304 this->checkBehaviourFile(os); 305 os << "/*!\n" 306 << "* \\brief Integrate behaviour law over the time step\n" 307 << "*/\n" 308 << "IntegrationResult\n" 309 << "integrate(const SMFlag smflag,const SMType smt) override{\n" 310 << "using namespace std;\n"; 311 if(this->mb.useQt()){ 312 os << "if(smflag!=MechanicalBehaviour<" << btype 313 << ",hypothesis,Type,use_qt>::STANDARDTANGENTOPERATOR){\n" 314 << "throw(runtime_error(\"invalid tangent operator flag\"));\n" 315 << "}\n"; 316 } else { 317 os << "if(smflag!=MechanicalBehaviour<" << btype 318 << ",hypothesis,Type,false>::STANDARDTANGENTOPERATOR){\n" 319 << "throw(runtime_error(\"invalid tangent operator flag\"));\n" 320 << "}\n"; 321 } 322 os << "if(!this->NewtonIntegration()){\n"; 323 if(this->mb.useQt()){ 324 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::FAILURE;\n"; 325 } else { 326 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::FAILURE;\n"; 327 } 328 os << "}\n" 329 << "this->dp = "; 330 auto p2=this->flows.begin(); 331 n=0; 332 while(p2!=this->flows.end()){ 333 os << "this->dp" << n << ""; 334 ++n; 335 ++p2; 336 if(p2!=this->flows.end()){ 337 os << "+"; 338 } 339 } 340 os << ";\n" 341 << "if(smt!=NOSTIFFNESSREQUESTED){\n" 342 << "if(!this->computeConsistentTangentOperator(smt)){\n"; 343 if(this->mb.useQt()){ 344 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::FAILURE;\n"; 345 } else { 346 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::FAILURE;\n"; 347 } 348 os << "}\n" 349 << "}\n" 350 << "this->deel = this->deto-dp*(this->n);\n" 351 << "this->updateStateVariables();\n" 352 << "this->sig = (this->lambda)*trace(this->eel)*StrainStensor::Id()+2*(this->mu)*(this->eel);\n" 353 << "this->updateAuxiliaryStateVariables();\n"; 354 for(const auto& v : d.getPersistentVariables()){ 355 this->writePhysicalBoundsChecks(os,v,false); 356 } 357 for(const auto& v : d.getPersistentVariables()){ 358 this->writeBoundsChecks(os,v,false); 359 } 360 if(this->mb.useQt()){ 361 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::SUCCESS;\n"; 362 } else { 363 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::SUCCESS;\n"; 364 } 365 os << "}\n\n"; 366 } 367 treatFlowRule()368 void MultipleIsotropicMisesFlowsDSL::treatFlowRule() 369 { 370 using namespace std; 371 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 372 FlowHandler flow; 373 this->checkNotEndOfFile("MultipleIsotropicMisesFlowsDSL::treatFlowRule", 374 "Expected flow rule name."); 375 if(this->current->value=="Plasticity"){ 376 ostringstream p; 377 ostringstream f; 378 ostringstream df_dseq; 379 ostringstream df_dp; 380 p << "p" << this->flows.size(); 381 f << "f" << this->flows.size(); 382 df_dseq << "df_dseq" << this->flows.size(); 383 df_dp << "df_dp" << this->flows.size(); 384 this->mb.addStateVariable(h,VariableDescription("strain",p.str(),1u,0u), 385 BehaviourData::FORCEREGISTRATION); 386 this->mb.addLocalVariable(h,VariableDescription("stress",f.str(),1u,0u)); 387 this->mb.addLocalVariable(h,VariableDescription("real",df_dseq.str(),1u,0u)); 388 this->mb.addLocalVariable(h,VariableDescription("stress",df_dp.str(),1u,0u)); 389 flow.flow = FlowHandler::PlasticFlow; 390 } else if(this->current->value=="Creep"){ 391 ostringstream p; 392 ostringstream f; 393 ostringstream df_dseq; 394 p << "p" << this->flows.size(); 395 f << "f" << this->flows.size(); 396 df_dseq << "df_dseq" << this->flows.size(); 397 this->mb.addStateVariable(h,VariableDescription("strain",p.str(),1u,0u), 398 BehaviourData::FORCEREGISTRATION); 399 this->mb.addLocalVariable(h,VariableDescription("DstrainDt",f.str(),1u,0u)); 400 this->mb.addLocalVariable(h,VariableDescription("DF_DSEQ_TYPE",df_dseq.str(),1u,0u)); 401 flow.flow = FlowHandler::CreepFlow; 402 } else if(this->current->value=="StrainHardeningCreep"){ 403 ostringstream p; 404 ostringstream f; 405 ostringstream df_dseq; 406 ostringstream df_dp; 407 p << "p" << this->flows.size(); 408 f << "f" << this->flows.size(); 409 df_dseq << "df_dseq" << this->flows.size(); 410 df_dp << "df_dp" << this->flows.size(); 411 this->mb.addStateVariable(h,VariableDescription("strain",p.str(),1u,0u), 412 BehaviourData::FORCEREGISTRATION); 413 this->mb.addLocalVariable(h,VariableDescription("DstrainDt",f.str(),1u,0u)); 414 this->mb.addLocalVariable(h,VariableDescription("DF_DSEQ_TYPE",df_dseq.str(),1u,0u)); 415 this->mb.addLocalVariable(h,VariableDescription("DstrainDt",df_dp.str(),1u,0u)); 416 flow.flow = FlowHandler::StrainHardeningCreepFlow; 417 } else { 418 this->throwRuntimeError("MultipleIsotropicMisesFlowsDSL::treatFlowRule", 419 "Unknown flow rule (read '"+this->current->value+ 420 "').Valid flow rules are 'Plasticity', 'Creep' and 'StrainHardeningCreep'."); 421 } 422 ++(this->current); 423 this->checkNotEndOfFile("MultipleIsotropicMisesFlowsDSL::treatFlowRule", 424 "Expected the beginning of a block or a specific theta value."); 425 if(this->current->value!="{"){ 426 istringstream converter(this->current->value); 427 ostringstream otheta; 428 ostringstream ose; 429 ostringstream oseq_e; 430 flow.hasSpecificTheta = true; 431 converter >> flow.theta; 432 if(!converter||(!converter.eof())){ 433 this->throwRuntimeError("MultipleIsotropicMisesFlowsDSL::treatFlowRule", 434 "Could not read theta value (read '"+this->current->value+"')."); 435 } 436 otheta << "mu_3_theta" << this->flows.size(); 437 ose << "se" << this->flows.size(); 438 oseq_e << "seq_e" << this->flows.size(); 439 this->reserveName(ose.str()); 440 this->mb.addLocalVariable(h,VariableDescription("stress",oseq_e.str(),1u,0u)); 441 ++(this->current); 442 } else { 443 flow.hasSpecificTheta = false; 444 } 445 ostringstream cname; 446 cname << BehaviourData::FlowRule << flows.size() << '\n'; 447 this->readCodeBlock(*this,cname.str(), 448 &MultipleIsotropicMisesFlowsDSL::flowRuleVariableModifier,true,false); 449 flow.flowRule = this->mb.getCode(ModellingHypothesis::UNDEFINEDHYPOTHESIS,cname.str()); 450 this->flows.push_back(flow); 451 } // end of MultipleIsotropicMisesFlowsDSL::treatFlowRule 452 453 void writeBehaviourParserSpecificInitializeMethodPart(std::ostream & os,const Hypothesis) const454 MultipleIsotropicMisesFlowsDSL::writeBehaviourParserSpecificInitializeMethodPart(std::ostream& os, 455 const Hypothesis) const 456 { 457 this->checkBehaviourFile(os); 458 os << "this->se=2*(this->mu)*(tfel::math::deviator(this->eel+(" 459 << this->mb.getClassName() 460 << "::theta)*(this->deto)));\n" 461 << "this->seq_e = sigmaeq(this->se);\n"; 462 unsigned short n = 0; 463 for(const auto& f : this->flows){ 464 if(f.hasSpecificTheta){ 465 os << "StressStensor se" << n << "=2*(this->mu)*(tfel::math::deviator(this->eel+("; 466 os << f.theta << ")*(this->deto)));\n"; 467 os << "this->seq_e" << n << " = sigmaeq(se" << n << ");\n"; 468 } 469 ++n; 470 } 471 os << "if(this->seq_e>100*std::numeric_limits<stress>::epsilon()){\n" 472 << "this->n = 1.5f*(this->se)/(this->seq_e);\n" 473 << "} else {\n" 474 << "this->n = StrainStensor(strain(0));\n" 475 << "}\n"; 476 } // end of MultipleIsotropicMisesFlowsDSL::writeBehaviourParserSpecificInitializeMethodPart 477 writeBehaviourComputeTangentOperator(std::ostream & os,const Hypothesis) const478 void MultipleIsotropicMisesFlowsDSL::writeBehaviourComputeTangentOperator(std::ostream& os, 479 const Hypothesis) const 480 { 481 os << "bool computeConsistentTangentOperator(const SMType smt){\n" 482 << "using namespace std;\n" 483 << "using tfel::material::computeElasticStiffness;\n" 484 << "if((smt==ELASTIC)||(smt==SECANTOPERATOR)){\n" 485 << "computeElasticStiffness<N,Type>::exe(this->Dt,this->lambda,this->mu);\n" 486 << "return true;\n" 487 << "}\n" 488 << "return false;\n" 489 << "}\n\n"; 490 } 491 492 MultipleIsotropicMisesFlowsDSL::~MultipleIsotropicMisesFlowsDSL() = default; 493 494 } // end of namespace mfront 495 496