1 /*! 2 * \file mtest/src/AsterCohesiveZoneModel.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 07 avril 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<cmath> 15 #include<ostream> 16 #include<algorithm> 17 18 #include"TFEL/Raise.hxx" 19 #include"TFEL/Math/tmatrix.hxx" 20 #include"TFEL/System/ExternalLibraryManager.hxx" 21 #include"MFront/MFrontLogStream.hxx" 22 #include"MFront/Aster/Aster.hxx" 23 #include"MTest/CurrentState.hxx" 24 #include"MTest/BehaviourWorkSpace.hxx" 25 #include"MTest/AsterCohesiveZoneModel.hxx" 26 27 namespace mtest 28 { 29 AsterCohesiveZoneModel(const Hypothesis h,const std::string & l,const std::string & b)30 AsterCohesiveZoneModel::AsterCohesiveZoneModel(const Hypothesis h, 31 const std::string& l, 32 const std::string& b) 33 : StandardBehaviourBase(h,l,b) 34 { 35 auto throw_if = [](const bool c, const std::string& m){ 36 tfel::raise_if(c,"AsterCohesiveZoneModel::AsterCohesiveZoneModel: "+m); 37 }; 38 auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager(); 39 throw_if(elm.getInterface(l,b)!="Aster", 40 "invalid interface '"+elm.getInterface(l,b)+"'"); 41 this->fct = elm.getAsterFunction(l,b); 42 this->emsg = elm.getAsterIntegrationErrorMessageFunction(l,b); 43 const auto& nh = ModellingHypothesis::toString(h); 44 this->mpnames = elm.getUMATMaterialPropertiesNames(l,b,nh); 45 throw_if(this->btype!=3u,"unsupported hypothesis"); 46 throw_if(this->stype!=0,"unsupported symmetry type"); 47 } 48 49 tfel::math::tmatrix<3u,3u,real> getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> &) const50 AsterCohesiveZoneModel::getRotationMatrix(const tfel::math::vector<real>&, 51 const tfel::math::tmatrix<3u,3u,real>&) const 52 { 53 tfel::raise("AsterCohesiveZoneModel::getRotationMatrix: invalid call"); 54 } // end of AsterCohesiveZoneModel::getRotationMatrix 55 56 void allocate(BehaviourWorkSpace & wk) const57 AsterCohesiveZoneModel::allocate(BehaviourWorkSpace& wk) const 58 { 59 const auto ndv = this->getGradientsSize(); 60 const auto nth = this->getThermodynamicForcesSize(); 61 const auto nstatev = this->getInternalStateVariablesSize(); 62 wk.kt.resize(nth,ndv); 63 wk.k.resize(nth,ndv); 64 wk.mps.resize(this->mpnames.size()==0 ? 1u : this->mpnames.size(),real(0)); 65 wk.ivs.resize(nstatev==0 ? 1u : nstatev,real(0)); 66 wk.nk.resize(nth,ndv); 67 wk.ne.resize(ndv); 68 wk.ns.resize(nth); 69 wk.nivs.resize(nstatev); 70 mtest::allocate(wk.cs,this->shared_from_this()); 71 } 72 73 void getGradientsDefaultInitialValues(tfel::math::vector<real> & v) const74 AsterCohesiveZoneModel::getGradientsDefaultInitialValues(tfel::math::vector<real>& v) const 75 { 76 std::fill(v.begin(),v.end(),real(0)); 77 } // end of AsterCohesiveZoneModel::setGradientsDefaultInitialValue 78 79 StiffnessMatrixType getDefaultStiffnessMatrixType() const80 AsterCohesiveZoneModel::getDefaultStiffnessMatrixType() const 81 { 82 return StiffnessMatrixType::CONSISTENTTANGENTOPERATOR; 83 } 84 85 std::pair<bool,real> computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const86 AsterCohesiveZoneModel::computePredictionOperator(BehaviourWorkSpace& wk, 87 const CurrentState& s, 88 const StiffnessMatrixType ktype) const 89 { 90 if(ktype==StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES){ 91 return {false,real(-1)}; 92 } 93 wk.cs = s; 94 return this->call_behaviour(wk.kt,wk.cs,wk,real(1),ktype,false); 95 } // end of AsterCohesiveZoneModel::computePredictionOperator 96 97 std::pair<bool,real> integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const98 AsterCohesiveZoneModel::integrate(CurrentState& s, 99 BehaviourWorkSpace& wk, 100 const real dt, 101 const StiffnessMatrixType ktype) const 102 { 103 return this->call_behaviour(wk.k,s,wk,dt,ktype,true); 104 } // end of AsterCohesiveZoneModel::integrate 105 106 std::pair<bool,real> call_behaviour(tfel::math::matrix<real> & Kt,CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype,const bool b) const107 AsterCohesiveZoneModel::call_behaviour(tfel::math::matrix<real>& Kt, 108 CurrentState& s, 109 BehaviourWorkSpace& wk, 110 const real dt, 111 const StiffnessMatrixType ktype, 112 const bool b) const 113 { 114 using namespace std; 115 using namespace tfel::math; 116 using namespace aster; 117 using tfel::math::vector; 118 AsterInt ntens; 119 AsterInt nprops = s.mprops1.size() == 0 ? 1 : static_cast<AsterInt>(s.mprops1.size()); 120 AsterInt nstatv; 121 AsterInt nummod; 122 const auto h = this->getHypothesis(); 123 if (h==ModellingHypothesis::AXISYMMETRICAL){ 124 ntens = 4; 125 nummod = 4u; 126 } else if (h==ModellingHypothesis::PLANESTRESS){ 127 ntens = 4; 128 nummod = 5u; 129 } else if (h==ModellingHypothesis::PLANESTRAIN){ 130 ntens = 4; 131 nummod = 6u; 132 } else if (h==ModellingHypothesis::TRIDIMENSIONAL){ 133 ntens = 6; 134 nummod = 3u; 135 } else { 136 tfel::raise("AsterCohesiveZoneModel::call_behaviour: " 137 "unsupported hypothesis"); 138 } 139 fill(Kt.begin(),Kt.end(),0.); 140 // choosing the type of stiffness matrix 141 if(b){ 142 if(ktype==StiffnessMatrixType::NOSTIFFNESS){ 143 // do nothing 144 } else if(ktype==StiffnessMatrixType::ELASTIC){ 145 Kt(0,0) = real(1); 146 } else if(ktype==StiffnessMatrixType::SECANTOPERATOR){ 147 Kt(0,0) = real(2); 148 } else if(ktype==StiffnessMatrixType::TANGENTOPERATOR){ 149 Kt(0,0) = real(3); 150 } else if(ktype==StiffnessMatrixType::CONSISTENTTANGENTOPERATOR){ 151 Kt(0,0) = real(4); 152 } else { 153 tfel::raise("AsterCohesiveZoneModel::call_behaviour: " 154 "invalid or unspecified stiffness matrix type"); 155 } 156 } else { 157 if(ktype==StiffnessMatrixType::ELASTIC){ 158 Kt(0,0) = real(-1); 159 } else if(ktype==StiffnessMatrixType::SECANTOPERATOR){ 160 Kt(0,0) = real(-2); 161 } else if(ktype==StiffnessMatrixType::TANGENTOPERATOR){ 162 Kt(0,0) = real(-3); 163 } else { 164 tfel::raise("AsterCohesiveZoneModel::call_behaviour : " 165 "invalid or unspecified stiffness matrix type"); 166 } 167 } 168 // using a local copy of material properties to handle the 169 // case where s.mprops1 is empty 170 copy(s.mprops1.begin(),s.mprops1.end(),wk.mps.begin()); 171 if(s.mprops1.empty()){ 172 wk.mps[0] = real(0); 173 } 174 // using a local copy of internal state variables to handle the 175 // case where s.iv0 is empty 176 copy(s.iv0.begin(),s.iv0.end(),wk.ivs.begin()); 177 if(s.iv0.empty()){ 178 wk.ivs[0] = real(0); 179 } 180 nstatv = static_cast<AsterInt>(wk.ivs.size()); 181 // rotation matrix 182 tmatrix<3u,3u,real> drot = transpose(s.r); 183 tvector<3u,real> ue0(real(0)); 184 tvector<3u,real> ude(real(0)); 185 copy(s.e0.begin(),s.e0.end(),ue0.begin()); 186 tmatrix<3u,3u,real>::size_type i; 187 for(i=0;i!=s.e1.size();++i){ 188 ude(i) = s.e1(i)-s.e0(i); 189 } 190 copy(s.s0.begin(),s.s0.end(),s.s1.begin()); 191 AsterReal ndt(1.); 192 (this->fct)(&(s.s1(0)),&(wk.ivs(0)),&Kt(0,0), 193 &ue0(0),&ude(0),&dt, 194 &(s.esv0(0)) ,&(s.desv(0)), 195 &(s.esv0(0))+1,&(s.desv(0))+1, 196 &ntens,&nstatv,&(wk.mps(0)), 197 &nprops,&drot(0,0),&ndt,&nummod); 198 if(ndt<1.){ 199 if(mfront::getVerboseMode()>=mfront::VERBOSE_LEVEL1){ 200 if(this->emsg!=nullptr){ 201 mfront::getLogStream() << this->emsg() << std::endl; 202 } 203 } 204 return {false,ndt}; 205 } 206 if(b){ 207 if(!s.iv0.empty()){ 208 copy_n(wk.ivs.begin(),s.iv1.size(),s.iv1.begin()); 209 } 210 } 211 return {true,ndt}; 212 } 213 214 void setOptionalMaterialPropertiesDefaultValues(EvolutionManager &,const EvolutionManager &) const215 AsterCohesiveZoneModel::setOptionalMaterialPropertiesDefaultValues(EvolutionManager&, 216 const EvolutionManager& ) const 217 {} // end of AsterCohesiveZoneModel::setOptionalMaterialPropertiesDefaultValues 218 219 AsterCohesiveZoneModel::~AsterCohesiveZoneModel() = default; 220 221 } // end of namespace mtest 222 223