1 /*! 2 * \file mtest/src/CastemStandardBehaviour.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 18 november 2013 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"TFEL/Raise.hxx" 15 #include"TFEL/System/ExternalLibraryManager.hxx" 16 #include"MFront/Castem/Castem.hxx" 17 #include"MFront/MFrontLogStream.hxx" 18 #include"MTest/Evolution.hxx" 19 #include"MTest/BehaviourWorkSpace.hxx" 20 #include"MTest/CastemStandardBehaviour.hxx" 21 22 namespace mtest 23 { 24 setMaterialProperties(StandardBehaviourDescription & umb,const tfel::material::ModellingHypothesis::Hypothesis h)25 static void setMaterialProperties(StandardBehaviourDescription& umb, 26 const tfel::material::ModellingHypothesis::Hypothesis h){ 27 using tfel::material::ModellingHypothesis; 28 auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager(); 29 umb.mpnames = elm.getUMATMaterialPropertiesNames(umb.library,umb.behaviour, 30 ModellingHypothesis::toString(h)); 31 if(umb.stype==0){ 32 if(h==ModellingHypothesis::PLANESTRESS){ 33 umb.mpnames.insert(umb.mpnames.begin(),{ 34 "YoungModulus","PoissonRatio","MassDensity", 35 "ThermalExpansion","PlateWidth"}); 36 } else { 37 umb.mpnames.insert(umb.mpnames.begin(),{ 38 "YoungModulus","PoissonRatio","MassDensity","ThermalExpansion"}); 39 } 40 } else { 41 if(h==ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN){ 42 umb.mpnames.insert(umb.mpnames.begin(),{ 43 "YoungModulus1","YoungModulus2","YoungModulus3", 44 "PoissonRatio12","PoissonRatio23","PoissonRatio13", 45 "MassDensity", 46 "ThermalExpansion1","ThermalExpansion2","ThermalExpansion3"}); 47 } else if(h==ModellingHypothesis::PLANESTRESS){ 48 umb.mpnames.insert(umb.mpnames.begin(),{ 49 "YoungModulus1","YoungModulus2","PoissonRatio12", 50 "ShearModulus12","V1X","V1Y","YoungModulus3", 51 "PoissonRatio23","PoissonRatio13","MassDensity", 52 "ThermalExpansion1","ThermalExpansion2","PlateWidth"}); 53 } else if((h==ModellingHypothesis::PLANESTRAIN)|| 54 (h==ModellingHypothesis::AXISYMMETRICAL)|| 55 (h==ModellingHypothesis::GENERALISEDPLANESTRAIN)){ 56 umb.mpnames.insert(umb.mpnames.begin(),{ 57 "YoungModulus1","YoungModulus2","YoungModulus3", 58 "PoissonRatio12","PoissonRatio23","PoissonRatio13", 59 "ShearModulus12","V1X","V1Y","MassDensity", 60 "ThermalExpansion1","ThermalExpansion2","ThermalExpansion3"}); 61 } else if(h==ModellingHypothesis::TRIDIMENSIONAL){ 62 umb.mpnames.insert(umb.mpnames.begin(),{ 63 "YoungModulus1","YoungModulus2","YoungModulus3", 64 "PoissonRatio12","PoissonRatio23","PoissonRatio13", 65 "ShearModulus12","ShearModulus23","ShearModulus13", 66 "V1X","V1Y","V1Z","V2X","V2Y","V2Z","MassDensity", 67 "ThermalExpansion1","ThermalExpansion2","ThermalExpansion3"}); 68 } else { 69 tfel::raise("setMaterialProperties: unsupported hypothesis"); 70 } 71 } 72 } // end of setMaterialProperties 73 CastemStandardBehaviour(const Hypothesis h,const std::string & l,const std::string & b)74 CastemStandardBehaviour::CastemStandardBehaviour(const Hypothesis h, 75 const std::string& l, 76 const std::string& b) 77 : StandardBehaviourBase(h,l,b) 78 { 79 auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager(); 80 tfel::raise_if(elm.getInterface(l,b)!="Castem", 81 "CastemStandardBehaviour::CastemStandardBehaviour: " 82 "invalid interface '"+elm.getInterface(l,b)+"'"); 83 this->fct = elm.getCastemExternalBehaviourFunction(l,b); 84 setMaterialProperties(*this,h); 85 } // end of CastemStandardBehaviour::CastemStandardBehaviour 86 CastemStandardBehaviour(const StandardBehaviourDescription & umb)87 CastemStandardBehaviour::CastemStandardBehaviour(const StandardBehaviourDescription& umb) 88 : StandardBehaviourBase(umb) 89 { 90 auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager(); 91 this->fct = elm.getCastemExternalBehaviourFunction(this->library,this->behaviour); 92 } // end of CastemStandardBehaviour::CastemStandardBehaviour 93 94 tfel::math::tmatrix<3u,3u,real> getRotationMatrix(const tfel::math::vector<real> & mp,const tfel::math::tmatrix<3u,3u,real> & r) const95 CastemStandardBehaviour::getRotationMatrix(const tfel::math::vector<real>& mp, 96 const tfel::math::tmatrix<3u,3u,real>& r) const 97 { 98 tfel::math::tmatrix<3u,3u,real> nr(0.); 99 const auto h = this->getHypothesis(); 100 if(h==ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN){ 101 nr(0,0) = nr(1,1) = nr(2,2) = 1.; // identité 102 } else if((h==ModellingHypothesis::AXISYMMETRICAL)|| 103 (h==ModellingHypothesis::GENERALISEDPLANESTRAIN)|| 104 (h==ModellingHypothesis::PLANESTRAIN)|| 105 (h==ModellingHypothesis::PLANESTRESS)){ 106 real V1X = 0; 107 real V1Y = 0; 108 if(h==ModellingHypothesis::PLANESTRESS){ 109 V1X = mp[4]; 110 V1Y = mp[5]; 111 } else { 112 V1X = mp[7]; 113 V1Y = mp[8]; 114 } 115 nr(0,0) = r(0,0)*V1X+r(0,1)*V1Y; 116 nr(1,0) = r(1,0)*V1X+r(1,1)*V1Y; 117 nr(0,1) =-nr(1,0); 118 nr(1,1) = nr(0,0); 119 nr(2,2) = 1; 120 } else if(h==ModellingHypothesis::TRIDIMENSIONAL){ 121 const real V1X = mp[9]; 122 const real V1Y = mp[10]; 123 const real V1Z = mp[11]; 124 const real V2X = mp[12]; 125 const real V2Y = mp[13]; 126 const real V2Z = mp[14]; 127 nr(0,0)=r(0,0)*V1X+r(0,1)*V1Y+r(0,2)*V1Z; 128 nr(1,0)=r(1,0)*V1X+r(1,1)*V1Y+r(1,2)*V1Z; 129 nr(2,0)=r(2,0)*V1X+r(2,1)*V1Y+r(2,2)*V1Z; 130 nr(0,1)=r(0,0)*V2X+r(0,1)*V2Y+r(0,2)*V2Z; 131 nr(1,1)=r(1,0)*V2X+r(1,1)*V2Y+r(1,2)*V2Z; 132 nr(2,1)=r(2,0)*V2X+r(2,1)*V2Y+r(2,2)*V2Z; 133 nr(0,2)=nr(1,0)*nr(2,1)-nr(1,1)*nr(2,0); 134 nr(1,2)=nr(2,0)*nr(0,1)-nr(2,1)*nr(0,0); 135 nr(2,2)=nr(0,0)*nr(1,1)-nr(0,1)*nr(1,0); 136 } 137 return nr; 138 } // end of CastemStandardBehaviour::getRotationMatrix 139 allocate(BehaviourWorkSpace & wk) const140 void CastemStandardBehaviour::allocate(BehaviourWorkSpace& wk) const{ 141 const auto ndv = this->getGradientsSize(); 142 const auto nth = this->getThermodynamicForcesSize(); 143 const auto nstatev = this->getInternalStateVariablesSize(); 144 wk.D.resize(nth,nth); 145 wk.k.resize(nth,ndv); 146 wk.kt.resize(nth,ndv); 147 wk.ivs.resize(nstatev==0 ? 1u : nstatev,real(0)); 148 wk.nk.resize(nth,ndv); 149 wk.ne.resize(ndv); 150 wk.ns.resize(nth); 151 wk.nivs.resize(nstatev); 152 if((this->usesGenericPlaneStressAlgorithm)&&(this->stype==0u)){ 153 wk.mps.resize(this->mpnames.size()+1); 154 } else { 155 wk.mps.resize(this->mpnames.size()); 156 } 157 mtest::allocate(wk.cs,this->shared_from_this()); 158 } // end of CastemStandardBehaviour::allocate 159 160 StiffnessMatrixType getDefaultStiffnessMatrixType() const161 CastemStandardBehaviour::getDefaultStiffnessMatrixType() const 162 { 163 return StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES; 164 } // end of CastemStandardBehaviour::getDefaultStiffnessMatrixType 165 166 void setOptionalMaterialPropertiesDefaultValues(EvolutionManager & mp,const EvolutionManager & evm) const167 CastemStandardBehaviour::setOptionalMaterialPropertiesDefaultValues(EvolutionManager& mp, 168 const EvolutionManager& evm) const 169 { 170 const auto h = this->getHypothesis(); 171 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"MassDensity",0.); 172 if(this->stype==1){ 173 if(h!=ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN){ 174 // orthotropic behaviour 175 const bool bv1x = evm.find("V1X")!=evm.end(); 176 const bool bv1y = evm.find("V1Y")!=evm.end(); 177 if((h==ModellingHypothesis::GENERALISEDPLANESTRAIN)|| 178 (h==ModellingHypothesis::AXISYMMETRICAL)|| 179 (h==ModellingHypothesis::PLANESTRESS)|| 180 (h==ModellingHypothesis::PLANESTRAIN)){ 181 const bool b = bv1x&&bv1y; 182 tfel::raise_if(((bv1x)&&(!b))||((bv1y)&&(!b)), 183 "Behaviour::setOptionalMaterialPropertiesDefaultValues : " 184 "if one component of the orthotropic basis is defined, all " 185 "the components must be defined."); 186 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V1X",1.); 187 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V1Y",0.); 188 } else if(h==ModellingHypothesis::TRIDIMENSIONAL){ 189 const bool bv1z = evm.find("V1Z")!=evm.end(); 190 const bool bv2x = evm.find("V2X")!=evm.end(); 191 const bool bv2y = evm.find("V2Y")!=evm.end(); 192 const bool bv2z = evm.find("V2Z")!=evm.end(); 193 const bool b = bv1x&&bv1y&&bv1z&&bv2x&&bv2y&&bv2z; 194 tfel::raise_if(((bv1x)&&(!b))||((bv1y)&&(!b))||((bv1z)&&(!b))|| 195 ((bv2x)&&(!b))||((bv2y)&&(!b))||((bv2z)&&(!b)), 196 "Behaviour::setOptionalMaterialPropertiesDefaultValues : " 197 "if one component of the orthotropic basis is defined, all " 198 "the components must be defined."); 199 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V1X",1.); 200 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V1Y",0.); 201 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V1Z",0.); 202 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V2X",0.); 203 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V2Y",1.); 204 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"V2Z",0.); 205 } else { 206 tfel::raise("Behaviour::setOptionalMaterialPropertiesDefaultValues : " 207 "unsupported hypothesis"); 208 } 209 } 210 } 211 Behaviour::setOptionalMaterialPropertyDefaultValue(mp,evm,"PlateWidth",1.); 212 } // end of CastemStandardBehaviour::setOptionalMaterialPropertiesDefaultValues 213 buildMaterialProperties(BehaviourWorkSpace & wk,const CurrentState & s) const214 void CastemStandardBehaviour::buildMaterialProperties(BehaviourWorkSpace& wk, 215 const CurrentState& s) const{ 216 auto throw_if = [](const bool c, const std::string& m){ 217 tfel::raise_if(c,"CastemSmallStrainBehaviour::buildMaterialProperties: "+m); 218 }; 219 if(this->usesGenericPlaneStressAlgorithm){ 220 if(this->stype==0u){ 221 throw_if(wk.mps.size()!=s.mprops1.size()+1, 222 "temporary material properties vector was not allocated properly"); 223 throw_if(s.mprops1.size()<3,"invalid number of material properties"); 224 wk.mps[0] = s.mprops1[0]; 225 wk.mps[1] = s.mprops1[1]; 226 wk.mps[2] = s.mprops1[2]; 227 wk.mps[3] = s.mprops1[3]; 228 // plate width 229 wk.mps[4] = castem::CastemReal(1); 230 std::copy(s.mprops1.begin()+4u,s.mprops1.end(),wk.mps.begin()+5u); 231 } else if(this->stype==1u){ 232 throw_if(wk.mps.size()!=s.mprops1.size(), 233 "temporary material properties vector was not allocated properly"); 234 throw_if(s.mprops1.size()<13,"invalid number of material properties"); 235 // YoungModulus1 236 wk.mps[0] = s.mprops1[0]; 237 // YoungModulus2 238 wk.mps[1] = s.mprops1[1]; 239 // PoissonRatio12 240 wk.mps[2] = s.mprops1[3]; 241 // ShearModulus12 242 wk.mps[3] = s.mprops1[6]; 243 // V1X 244 wk.mps[4] = s.mprops1[7]; 245 // V1Y 246 wk.mps[5] = s.mprops1[8]; 247 // YoungModulus3 248 wk.mps[6] = s.mprops1[2]; 249 // PoissonRatio23 250 wk.mps[7] = s.mprops1[4]; 251 // PoissonRatio13 252 wk.mps[8] = s.mprops1[5]; 253 // MassDensity 254 wk.mps[9] = s.mprops1[9]; 255 // ThermalExpansion1 256 wk.mps[10] = s.mprops1[10]; 257 // ThermalExpansion2 258 wk.mps[11] = s.mprops1[11]; 259 // ThermalExpansion3 (does not exists in mps) 260 // plate width 261 wk.mps[12] = castem::CastemReal(1); 262 std::copy(s.mprops1.begin()+13u,s.mprops1.end(),wk.mps.begin()+13u); 263 } else { 264 throw_if(true,"unsupported symmetry type"); 265 } 266 } else { 267 throw_if(wk.mps.size()!=s.mprops1.size(), 268 "temporary material properties vector was not allocated properly"); 269 std::copy(s.mprops1.begin(),s.mprops1.end(),wk.mps.begin()); 270 } 271 } // end of CastemStandardBehaviour::buildMaterialProperties 272 273 CastemStandardBehaviour::~CastemStandardBehaviour() = default; 274 275 } 276