1 /*! 2 * \file mfront/src/IsotropicStrainHardeningMisesCreepDSL.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 02 jui 2007 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<stdexcept> 16 #include"TFEL/Raise.hxx" 17 #include"MFront/DSLUtilities.hxx" 18 #include"MFront/MFrontDebugMode.hxx" 19 #include"MFront/DSLFactory.hxx" 20 #include"MFront/IsotropicStrainHardeningMisesCreepDSL.hxx" 21 22 namespace mfront{ 23 IsotropicStrainHardeningMisesCreepDSL()24 IsotropicStrainHardeningMisesCreepDSL::IsotropicStrainHardeningMisesCreepDSL() 25 { 26 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 27 this->mb.setDSLName("IsotropicStrainHardeningMisesCreep"); 28 // Default state vars 29 this->mb.addStateVariable(h,VariableDescription("StrainStensor","eel",1u,0u)); 30 this->mb.addStateVariable(h,VariableDescription("strain","p",1u,0u)); 31 this->mb.setGlossaryName(h,"eel","ElasticStrain"); 32 this->mb.setGlossaryName(h,"p","EquivalentViscoplasticStrain"); 33 // default local vars 34 this->reserveName("mu_3"); 35 this->mb.addLocalVariable(h,VariableDescription("DstrainDt","f",1u,0u)); 36 this->mb.addLocalVariable(h,VariableDescription("DF_DSEQ_TYPE","df_dseq",1u,0u)); 37 this->mb.addLocalVariable(h,VariableDescription("DstrainDt","df_dp",1u,0u)); 38 this->mb.addLocalVariable(h,VariableDescription("StressStensor","se",1u,0u)); 39 this->mb.addLocalVariable(h,VariableDescription("stress","seq",1u,0u)); 40 this->mb.addLocalVariable(h,VariableDescription("stress","seq_e",1u,0u)); 41 this->mb.addLocalVariable(h,VariableDescription("StrainStensor","n",1u,0u)); 42 this->mb.addLocalVariable(h,VariableDescription("strain","p_",1u,0u)); 43 this->mb.setAttribute(h,BehaviourData::hasConsistentTangentOperator,true); 44 this->mb.setAttribute(h,BehaviourData::isConsistentTangentOperatorSymmetric,true); 45 } 46 getName()47 std::string IsotropicStrainHardeningMisesCreepDSL::getName() 48 { 49 return "IsotropicStrainHardeningMisesCreep"; 50 } 51 getDescription()52 std::string IsotropicStrainHardeningMisesCreepDSL::getDescription() 53 { 54 return "this parser is used for standard strain hardening creep behaviours " 55 "of the form dp/dt=f(s,p) where p is the equivalent creep strain " 56 "and s the equivalent mises stress"; 57 } // end of IsotropicStrainHardeningMisesCreepDSL::getDescription 58 endsInputFileProcessing()59 void IsotropicStrainHardeningMisesCreepDSL::endsInputFileProcessing() 60 { 61 IsotropicBehaviourDSLBase::endsInputFileProcessing(); 62 const auto h = ModellingHypothesis::UNDEFINEDHYPOTHESIS; 63 tfel::raise_if(!this->mb.hasCode(h,BehaviourData::FlowRule), 64 "IsotropicMisesCreepDSL::endsInputFileProcessing: " 65 "no flow rule defined"); 66 } // end of IsotropicStrainHardeningMisesCreepDSL::endsInputFileProcessing 67 68 void writeBehaviourParserSpecificMembers(std::ostream & os,const Hypothesis h) const69 IsotropicStrainHardeningMisesCreepDSL::writeBehaviourParserSpecificMembers(std::ostream& os, 70 const Hypothesis h) const 71 { 72 this->checkBehaviourFile(os); 73 if(!this->mb.hasCode(h,BehaviourData::FlowRule)){ 74 this->throwRuntimeError("IsotropicStrainHardeningMisesCreepDSL::" 75 "writeBehaviourParserSpecificMembers", 76 "no flow rule declared " 77 "(use the @FlowRule directive)"); 78 } 79 os << "void computeFlow(){\n" 80 << "using namespace std;\n" 81 << "using namespace tfel::math;\n" 82 << "using namespace tfel::material;\n" 83 << "using std::vector;\n"; 84 writeMaterialLaws(os,this->mb.getMaterialLaws()); 85 os << this->mb.getCode(h,BehaviourData::FlowRule) 86 << "\n}\n\n" 87 << "bool NewtonIntegration(){\n" 88 << "using namespace std;\n" 89 << "using namespace tfel::math;\n" 90 << "bool converge=false;\n" 91 << "bool inversible=true;\n" 92 << "strain newton_f;\n" 93 << "strain newton_df;\n" 94 << "real newton_epsilon = 100*std::numeric_limits<real>::epsilon();\n" 95 << "stress mu_3 = 3*(this->mu);\n" 96 << "" 97 << "unsigned int iter = 0u;\n" 98 << "this->p_=this->p+this->dp;\n" 99 << "while((converge==false)&&\n" 100 << "(iter<(this->iterMax))&&\n" 101 << "(inversible==true)){\n" 102 << "this->seq=std::max(this->seq_e-mu_3*(this->theta)*(this->dp),real(0.f));\n" 103 << "this->computeFlow();\n" 104 << "newton_f = this->dp - (this->f)*(this->dt);\n" 105 << "newton_df = 1-(this->theta)*(this->dt)*((this->df_dp)-mu_3*(this->df_dseq));\n" 106 << "if(std::abs(base_cast(newton_df))" 107 << ">newton_epsilon){\n" 108 << "this->dp -= newton_f/newton_df;\n" 109 << "this->p_=this->p + (this->theta)*(this->dp);\n" 110 << "iter+=1;\n"; 111 if(getDebugMode()){ 112 os << "cout << \"" << this->mb.getClassName() 113 << "::NewtonIntegration() : iteration \" " 114 << "<< iter << \" : \" << std::abs(tfel::math::base_cast(newton_f)) << endl;\n"; 115 } 116 os << "converge = (std::abs(tfel::math::base_cast(newton_f))<" 117 << "(this->epsilon));\n" 118 << "} else {\n" 119 << "inversible=false;\n" 120 << "}\n" 121 << "}\n\n" 122 << "if(inversible==false){\n" 123 << "return false;\n" 124 << "}\n\n" 125 << "if(iter==this->iterMax){\n"; 126 if(getDebugMode()){ 127 os << "cout << \"" << this->mb.getClassName() 128 << "::NewtonIntegration() : no convergence after \" " 129 << "<< iter << \" iterations\"<< endl << endl;\n"; 130 os << "cout << *this << endl;\n"; 131 } 132 os << "return false;\n" 133 << "}\n\n"; 134 if(getDebugMode()){ 135 os << "cout << \"" << this->mb.getClassName() 136 << "::NewtonIntegration() : convergence after \" " 137 << "<< iter << \" iterations\"<< endl << endl;\n"; 138 } 139 os << "return true;\n" 140 << "}\n\n"; 141 } // end of writeBehaviourParserSpecificMembers 142 writeBehaviourIntegrator(std::ostream & os,const Hypothesis h) const143 void IsotropicStrainHardeningMisesCreepDSL::writeBehaviourIntegrator(std::ostream& os, 144 const Hypothesis h) const 145 { 146 const auto btype = this->mb.getBehaviourTypeFlag(); 147 const auto& d = this->mb.getBehaviourData(h); 148 this->checkBehaviourFile(os); 149 os << "/*!\n" 150 << "* \\brief Integrate behaviour law over the time step\n" 151 << "*/\n" 152 << "IntegrationResult\n" 153 << "integrate(const SMFlag smflag,const SMType smt) override{\n" 154 << "using namespace std;\n"; 155 if(this->mb.useQt()){ 156 os << "if(smflag!=MechanicalBehaviour<" << btype 157 << ",hypothesis,Type,use_qt>::STANDARDTANGENTOPERATOR){\n" 158 << "throw(runtime_error(\"invalid tangent operator flag\"));\n" 159 << "}\n"; 160 } else { 161 os << "if(smflag!=MechanicalBehaviour<" << btype 162 << ",hypothesis,Type,false>::STANDARDTANGENTOPERATOR){\n" 163 << "throw(runtime_error(\"invalid tangent operator flag\"));\n" 164 << "}\n"; 165 } 166 os << "if(!this->NewtonIntegration()){\n"; 167 if(this->mb.useQt()){ 168 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::FAILURE;\n"; 169 } else { 170 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::FAILURE;\n"; 171 } 172 os << "}\n" 173 << "if(smt!=NOSTIFFNESSREQUESTED){\n" 174 << "if(!this->computeConsistentTangentOperator(smt)){\n"; 175 if(this->mb.useQt()){ 176 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::FAILURE;\n"; 177 } else { 178 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::FAILURE;\n"; 179 } 180 os << "}\n" 181 << "}\n" 182 << "this->deel = this->deto-(this->dp)*(this->n);\n" 183 << "this->updateStateVariables();\n" 184 << "this->sig = (this->lambda_tdt)*trace(this->eel)*StrainStensor::Id()+2*(this->mu_tdt)*(this->eel);\n" 185 << "this->updateAuxiliaryStateVariables();\n"; 186 for(const auto& v : d.getPersistentVariables()){ 187 this->writePhysicalBoundsChecks(os,v,false); 188 } 189 for(const auto& v : d.getPersistentVariables()){ 190 this->writeBoundsChecks(os,v,false); 191 } 192 if(this->mb.useQt()){ 193 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,use_qt>::SUCCESS;\n"; 194 } else { 195 os << "return MechanicalBehaviour<" << btype << ",hypothesis,Type,false>::SUCCESS;\n"; 196 } 197 os << "}\n\n"; 198 } 199 200 void writeBehaviourComputeTangentOperator(std::ostream & os,const Hypothesis) const201 IsotropicStrainHardeningMisesCreepDSL::writeBehaviourComputeTangentOperator(std::ostream& os, 202 const Hypothesis) const 203 { 204 os << "bool computeConsistentTangentOperator(const SMType smt){\n" 205 << "using tfel::material::computeElasticStiffness;\n" 206 << "using tfel::math::st2tost2;\n" 207 << "if(smt==CONSISTENTTANGENTOPERATOR){\n" 208 << "computeElasticStiffness<N,Type>::exe(this->Dt,this->lambda_tdt,this->mu_tdt);\n" 209 << "if(this->seq_e>(0.01*(this->young))*std::numeric_limits<stress>::epsilon()){\n" 210 << "const real ccto_tmp_1 = this->dp/this->seq_e;\n" 211 << "const auto& M = st2tost2<N,Type>::M();\n" 212 << "this->Dt += -4*(this->mu_tdt)*(this->mu)*(this->theta)*(ccto_tmp_1*M-(ccto_tmp_1-this->df_dseq*(this->dt)/(Type(1)+(this->theta)*(this->dt)*(Type(3)*(this->mu)*this->df_dseq-(this->df_dp))))*((this->n)^(this->n)));\n" 213 << "}\n" 214 << "} else if((smt==ELASTIC)||(smt==SECANTOPERATOR)){\n" 215 << "computeElasticStiffness<N,Type>::exe(this->Dt,this->lambda_tdt,this->mu_tdt);\n" 216 << "} else {\n" 217 << "return false;" 218 << "}\n" 219 << "return true;\n" 220 << "}\n\n"; 221 } 222 223 void writeBehaviourParserSpecificInitializeMethodPart(std::ostream & os,const Hypothesis) const224 IsotropicStrainHardeningMisesCreepDSL::writeBehaviourParserSpecificInitializeMethodPart(std::ostream& os, 225 const Hypothesis) const 226 { 227 this->checkBehaviourFile(os); 228 os << "this->se=2*(this->mu)*(tfel::math::deviator(this->eel+(" 229 << this->mb.getClassName() 230 << "::theta)*(this->deto)));\n" 231 << "this->seq_e = sigmaeq(this->se);\n" 232 << "if(this->seq_e>(0.01*(this->young))*std::numeric_limits<stress>::epsilon()){\n" 233 << "this->n = 1.5f*(this->se)/(this->seq_e);\n" 234 << "} else {\n" 235 << "this->n = StrainStensor(strain(0));\n" 236 << "}\n"; 237 } // end of IsotropicStrainHardeningMisesCreepDSL::writeBehaviourParserSpecificInitializeMethodPart 238 239 IsotropicStrainHardeningMisesCreepDSL::~IsotropicStrainHardeningMisesCreepDSL() = default; 240 241 } // end of namespace mfront 242 243 244