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