1 /*! 2 * \file mtest/src/CastemCohesiveZoneModel.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 <limits> 16 #include <algorithm> 17 18 #include "TFEL/Raise.hxx" 19 #include "TFEL/Math/tmatrix.hxx" 20 #include "TFEL/System/ExternalLibraryManager.hxx" 21 #include "MFront/Castem/Castem.hxx" 22 #include "MTest/CurrentState.hxx" 23 #include "MTest/BehaviourWorkSpace.hxx" 24 #include "MTest/CastemCohesiveZoneModel.hxx" 25 26 namespace mtest { 27 CastemCohesiveZoneModel(const Hypothesis h,const std::string & l,const std::string & b)28 CastemCohesiveZoneModel::CastemCohesiveZoneModel(const Hypothesis h, 29 const std::string& l, 30 const std::string& b) 31 : StandardBehaviourBase(h, l, b) { 32 auto throw_if = [](const bool c, const std::string& m) { 33 tfel::raise_if(c, 34 "CastemCohesiveZoneModel::CastemCohesiveZoneModel: " + m); 35 }; 36 auto& elm = 37 tfel::system::ExternalLibraryManager::getExternalLibraryManager(); 38 throw_if(elm.getInterface(l, b) != "Castem", 39 "invalid interface '" + elm.getInterface(l, b) + "'"); 40 const auto& nh = ModellingHypothesis::toString(h); 41 this->fct = elm.getCastemExternalBehaviourFunction(l, b); 42 this->mpnames = elm.getUMATMaterialPropertiesNames(l, b, nh); 43 throw_if(this->btype != 3u, "invalid behaviour type"); 44 if (this->stype == 0) { 45 this->mpnames.insert(this->mpnames.begin(), "NormalThermalExpansion"); 46 this->mpnames.insert(this->mpnames.begin(), "MassDensity"); 47 // Those are the conventions used by Cast3M. The UMATInterface 48 // exchanges the 'NormalStiffness' and the 'TangentialStiffness' 49 // material properties to match MFront conventions 50 this->mpnames.insert(this->mpnames.begin(), "NormalStiffness"); 51 this->mpnames.insert(this->mpnames.begin(), "TangentialStiffness"); 52 } else { 53 throw_if(true, "unsupported symmetry type"); 54 } 55 } 56 getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> &) const57 tfel::math::tmatrix<3u, 3u, real> CastemCohesiveZoneModel::getRotationMatrix( 58 const tfel::math::vector<real>&, 59 const tfel::math::tmatrix<3u, 3u, real>&) const { 60 tfel::raise("CastemCohesiveZoneModel::getRotationMatrix: invalid call"); 61 } // end of CastemCohesiveZoneModel::getRotationMatrix 62 allocate(BehaviourWorkSpace & wk) const63 void CastemCohesiveZoneModel::allocate(BehaviourWorkSpace& wk) const { 64 const auto ndv = this->getGradientsSize(); 65 const auto nth = this->getThermodynamicForcesSize(); 66 const auto nstatev = this->getInternalStateVariablesSize(); 67 wk.D.resize(nth, ndv); 68 wk.kt.resize(nth, ndv); 69 wk.k.resize(nth, ndv); 70 wk.ivs.resize(nstatev == 0 ? 1u : nstatev, real(0)); 71 wk.nk.resize(nth, ndv); 72 wk.ne.resize(ndv); 73 wk.ns.resize(nth); 74 wk.nivs.resize(nstatev); 75 mtest::allocate(wk.cs, this->shared_from_this()); 76 } 77 getGradientsDefaultInitialValues(tfel::math::vector<real> & v) const78 void CastemCohesiveZoneModel::getGradientsDefaultInitialValues( 79 tfel::math::vector<real>& v) const { 80 std::fill(v.begin(), v.end(), real(0)); 81 } // end of CastemCohesiveZoneModel::setGradientsDefaultInitialValue 82 getDefaultStiffnessMatrixType() const83 StiffnessMatrixType CastemCohesiveZoneModel::getDefaultStiffnessMatrixType() 84 const { 85 return StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES; 86 } 87 computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const88 std::pair<bool, real> CastemCohesiveZoneModel::computePredictionOperator( 89 BehaviourWorkSpace& wk, 90 const CurrentState& s, 91 const StiffnessMatrixType ktype) const { 92 if (ktype == StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES) { 93 // rotation matrix 94 const auto drot = transpose(s.r); 95 this->computeElasticStiffness(wk.k, s.mprops1, drot); 96 return {true, 1}; 97 } 98 tfel::raise( 99 "CastemCohesiveZoneModel::computePredictionOperator: " 100 "computation of the tangent operator " 101 "is not supported"); 102 } // end of CastemCohesiveZoneModel::computePredictionOperator 103 integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const104 std::pair<bool, real> CastemCohesiveZoneModel::integrate( 105 CurrentState& s, 106 BehaviourWorkSpace& wk, 107 const real dt, 108 const StiffnessMatrixType ktype) const { 109 using namespace std; 110 using namespace tfel::math; 111 using namespace castem; 112 using tfel::math::vector; 113 CastemInt ntens; 114 CastemInt ndi; 115 CastemInt nprops = static_cast<CastemInt>(s.mprops1.size()); 116 CastemInt nstatv; 117 const auto h = this->getHypothesis(); 118 tfel::raise_if( 119 (h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN) || 120 (h == ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRESS) || 121 (h == ModellingHypothesis::AXISYMMETRICAL), 122 "CastemCohesiveZoneModel::integrate: " 123 "unsupported modelling hypothesis"); 124 if (h == ModellingHypothesis::PLANESTRESS) { 125 ndi = -2; 126 ntens = 2; 127 } else if (h == ModellingHypothesis::PLANESTRAIN) { 128 ndi = -1; 129 ntens = 2; 130 } else if (h == ModellingHypothesis::GENERALISEDPLANESTRAIN) { 131 ndi = -3; 132 ntens = 2; 133 } else if (h == ModellingHypothesis::TRIDIMENSIONAL) { 134 ndi = 2; 135 ntens = 3; 136 } else { 137 tfel::raise("CastemCohesiveZoneModel::integrate: unsupported hypothesis"); 138 } 139 tfel::raise_if((wk.D.getNbRows() != wk.k.getNbRows()) || 140 (wk.D.getNbCols() != wk.k.getNbCols()), 141 "CastemCohesiveZoneModel::integrate: " 142 "the memory has not been allocated correctly"); 143 tfel::raise_if(((s.iv0.size() == 0) && (wk.ivs.size() != 1u)) || 144 ((s.iv0.size() != 0) && (s.iv0.size() != wk.ivs.size())), 145 "CastemCohesiveZoneModel::integrate: " 146 "the memory has not been allocated correctly"); 147 std::fill(wk.D.begin(), wk.D.end(), 0.); 148 if (s.iv0.size() != 0) { 149 std::copy(s.iv0.begin(), s.iv0.end(), wk.ivs.begin()); 150 } 151 nstatv = static_cast<CastemInt>(wk.ivs.size()); 152 // rotation matrix 153 tmatrix<3u, 3u, real> drot = transpose(s.r); 154 CastemInt kinc(1); 155 tvector<3u, real> ue0(real(0)); 156 tvector<3u, real> ude(real(0)); 157 if (ntens == 2) { 158 ue0[0] = s.e0[1]; 159 ue0[1] = s.e0[0]; 160 ude[0] = s.e1[1] - s.e0[1]; 161 ude[1] = s.e1[0] - s.e0[0]; 162 s.s1[0] = s.s0[1]; 163 s.s1[1] = s.s0[0]; 164 } 165 if (ntens == 3) { 166 ue0[0] = s.e0[1]; 167 ue0[1] = s.e0[2]; 168 ue0[2] = s.e0[0]; 169 ude[0] = s.e1[1] - s.e0[1]; 170 ude[1] = s.e1[2] - s.e0[2]; 171 ude[2] = s.e1[0] - s.e0[0]; 172 s.s1[0] = s.s0[1]; 173 s.s1[1] = s.s0[2]; 174 s.s1[2] = s.s0[0]; 175 } 176 CastemReal ndt = std::numeric_limits<CastemReal>::max(); 177 (this->fct)(&s.s1(0), &wk.ivs(0), &(wk.D(0, 0)), nullptr, nullptr, nullptr, 178 nullptr, nullptr, nullptr, nullptr, &ue0(0), &ude(0), nullptr, 179 &dt, &(s.esv0(0)), &(s.desv(0)), &(s.esv0(0)) + 1, 180 &(s.desv(0)) + 1, nullptr, &ndi, nullptr, &ntens, &nstatv, 181 &(s.mprops1(0)), &nprops, nullptr, &drot(0, 0), &ndt, nullptr, 182 nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, 183 &kinc, 0); 184 if (kinc != 1) { 185 return {false, ndt}; 186 } 187 // tangent operator (...) 188 if (ktype != StiffnessMatrixType::NOSTIFFNESS) { 189 if (ktype == StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES) { 190 this->computeElasticStiffness(wk.k, s.mprops1, drot); 191 } else { 192 tfel::raise( 193 "CastemCohesiveZoneModel::integrate : " 194 "computation of the tangent operator " 195 "is not supported"); 196 } 197 } 198 if (!s.iv1.empty()) { 199 copy_n(wk.ivs.begin(), s.iv1.size(), s.iv1.begin()); 200 } 201 // turning things in standard conventions 202 if (ntens == 2) { 203 swap(s.s1[0], s.s1[1]); 204 } 205 if (ntens == 3) { 206 const real tmp = s.s1[0]; 207 s.s1[0] = s.s1[2]; 208 s.s1[2] = s.s1[1]; 209 s.s1[1] = tmp; 210 } 211 return {true, ndt}; 212 } // end of CastemCohesiveZoneModel::integrate 213 computeElasticStiffness(tfel::math::matrix<real> & Kt,const tfel::math::vector<real> & mp,const tfel::math::tmatrix<3u,3u,real> &) const214 void CastemCohesiveZoneModel::computeElasticStiffness( 215 tfel::math::matrix<real>& Kt, 216 const tfel::math::vector<real>& mp, 217 const tfel::math::tmatrix<3u, 3u, real>&) const { 218 using namespace tfel::math; 219 if (this->stype == 0u) { 220 const auto h = this->getHypothesis(); 221 const auto zero = real(0); 222 if ((h == ModellingHypothesis::PLANESTRESS) || 223 (h == ModellingHypothesis::PLANESTRAIN) || 224 (h == ModellingHypothesis::GENERALISEDPLANESTRAIN)) { 225 Kt(0, 0) = mp(1); 226 Kt(0, 1) = zero; 227 Kt(1, 0) = zero; 228 Kt(1, 1) = mp(0); 229 } else if (h == ModellingHypothesis::TRIDIMENSIONAL) { 230 Kt(0, 0) = mp(1); 231 Kt(1, 1) = mp(0); 232 Kt(2, 2) = mp(0); 233 Kt(0, 1) = Kt(0, 2) = zero; 234 Kt(1, 0) = Kt(1, 2) = zero; 235 Kt(2, 0) = Kt(2, 1) = zero; 236 } else { 237 tfel::raise( 238 "CastemCohesiveZoneModel::computeElasticStiffness: " 239 "unsupported hypothesis"); 240 } 241 } else if (this->stype == 1u) { 242 tfel::raise( 243 "CastemCohesiveZoneModel::computeElasticStiffness: " 244 "invalid behaviour type (orthotropic type " 245 "is not supported yet)"); 246 } else { 247 tfel::raise( 248 "CastemCohesiveZoneModel::computeElasticStiffness: " 249 "invalid behaviour type (neither " 250 "isotropic or orthotropic)"); 251 } 252 } 253 254 std::vector<std::string> getOptionalMaterialProperties() const255 CastemCohesiveZoneModel::getOptionalMaterialProperties() const { 256 return {"MassDensity","NormalThermalExpansion"}; 257 } // end of CastemCohesiveZoneModel::getOptionalMaterialProperties 258 setOptionalMaterialPropertiesDefaultValues(EvolutionManager & mp,const EvolutionManager & evm) const259 void CastemCohesiveZoneModel::setOptionalMaterialPropertiesDefaultValues( 260 EvolutionManager& mp, const EvolutionManager& evm) const { 261 Behaviour::setOptionalMaterialPropertyDefaultValue(mp, evm, "MassDensity", 262 0.); 263 if (this->stype == 0) { 264 Behaviour::setOptionalMaterialPropertyDefaultValue( 265 mp, evm, "NormalThermalExpansion", 0.); 266 } else { 267 tfel::raise( 268 "CastemCohesiveZoneModel::CastemCohesiveZoneModel : " 269 "unsupported symmetry type"); 270 } 271 } // end of 272 // CastemCohesiveZoneModel::setOptionalMaterialPropertiesDefaultValues 273 274 CastemCohesiveZoneModel::~CastemCohesiveZoneModel() = default; 275 276 } // end of namespace mtest 277