1 /*!
2  * \file  mtest/src/AbaqusStandardBehaviour.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 <sstream>
16 #include <algorithm>
17 
18 #include "TFEL/Raise.hxx"
19 #include "TFEL/Math/tmatrix.hxx"
20 #include "TFEL/Math/stensor.hxx"
21 #include "TFEL/Math/st2tost2.hxx"
22 #include "TFEL/System/ExternalLibraryManager.hxx"
23 #include "MFront/Abaqus/Abaqus.hxx"
24 #include "MFront/Abaqus/AbaqusComputeStiffnessTensor.hxx"
25 
26 #include "MTest/CurrentState.hxx"
27 #include "MTest/BehaviourWorkSpace.hxx"
28 #include "MTest/AbaqusStandardBehaviour.hxx"
29 
30 namespace mtest {
31 
getHypothesisSuffix(const Hypothesis h)32   std::string AbaqusStandardBehaviour::getHypothesisSuffix(const Hypothesis h) {
33     if (h == ModellingHypothesis::AXISYMMETRICAL) {
34       return "_AXIS";
35     } else if (h == ModellingHypothesis::PLANESTRAIN) {
36       return "_PSTRAIN";
37     } else if (h == ModellingHypothesis::PLANESTRESS) {
38       return "_PSTRESS";
39     } else if (h == ModellingHypothesis::TRIDIMENSIONAL) {
40       return "_3D";
41     }
42     tfel::raise(
43         "AbaqusStandardBehaviour::getHypothesisSuffix: "
44         "invalid hypothesis.");
45   } // end of AbaqusStandardBehaviour::getHypothesisSuffix
46 
extractBehaviourName(const std::string & b,const Hypothesis h)47   std::string AbaqusStandardBehaviour::extractBehaviourName(const std::string& b, const Hypothesis h) {
48     auto ends = [&b](const std::string& s) {
49       if (b.length() >= s.length()) {
50         return b.compare(b.length() - s.length(), s.length(), s) == 0;
51       }
52       return false;
53     };
54     const auto s = AbaqusStandardBehaviour::getHypothesisSuffix(h);
55     tfel::raise_if(!ends(s),
56                    "AbaqusStandardBehaviour::extractBehaviourName: "
57                    "invalid function name.");
58     return {b.begin(), b.begin() + b.length() - s.length()};
59   }
60 
AbaqusStandardBehaviour(const Hypothesis h,const std::string & l,const std::string & b)61   AbaqusStandardBehaviour::AbaqusStandardBehaviour(const Hypothesis h,
62                                                    const std::string& l,
63                                                    const std::string& b)
64       : StandardBehaviourBase(h, l, AbaqusStandardBehaviour::extractBehaviourName(b, h)) {
65     auto throw_if = [](const bool c, const std::string& m) {
66       tfel::raise_if(c,
67                      "AbaqusStandardBehaviour::"
68                      "AbaqusStandardBehaviour: " +
69                          m);
70     };
71     auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager();
72     const auto bn = AbaqusStandardBehaviour::extractBehaviourName(b, h);
73     throw_if(elm.getInterface(l, bn) != "Abaqus",
74              "invalid interface '" + elm.getInterface(l, bn) + "'");
75     this->fct = elm.getAbaqusExternalBehaviourFunction(l, b);
76     if (this->stype == 1u) {
77       this->omp = elm.getAbaqusOrthotropyManagementPolicy(l, bn);
78       if (this->omp == 2u) {
79         auto aivs = std::vector<std::string>{};
80         if ((h == ModellingHypothesis::PLANESTRESS) || (h == ModellingHypothesis::PLANESTRAIN) ||
81             (h == ModellingHypothesis::AXISYMMETRICAL) ||
82             (h == ModellingHypothesis::GENERALISEDPLANESTRAIN)) {
83           aivs = {"FirstOrthotropicDirection_1", "FirstOrthotropicDirection_2"};
84         } else if (h == ModellingHypothesis::TRIDIMENSIONAL) {
85           aivs = {"FirstOrthotropicDirection_1",  "FirstOrthotropicDirection_2",
86                   "FirstOrthotropicDirection_3",  "SecondOrthotropicDirection_1",
87                   "SecondOrthotropicDirection_2", "SecondOrthotropicDirection_3"};
88         } else {
89           throw_if(true, "unsupported modelling hypothesis");
90         }
91         for (const auto& iv : aivs) {
92           throw_if(std::find(this->ivnames.begin(), this->ivnames.end(), iv) != this->ivnames.end(),
93                    iv + " is a reserved name");
94           this->ivtypes.insert(this->ivtypes.begin(), 0);
95         }
96         this->ivnames.insert(this->ivnames.begin(), aivs.begin(), aivs.end());
97       }
98     }
99     auto tmp = std::vector<std::string>{};
100     if (this->etype == 0u) {
101       if (this->requiresStiffnessTensor) {
102         tmp.insert(tmp.end(), {"YoungModulus", "PoissonRatio"});
103       }
104       if (this->requiresThermalExpansionCoefficientTensor) {
105         tmp.push_back("ThermalExpansion");
106       }
107     } else if (this->etype == 1u) {
108       if ((h == ModellingHypothesis::PLANESTRESS) || (h == ModellingHypothesis::PLANESTRAIN) ||
109           (h == ModellingHypothesis::AXISYMMETRICAL)) {
110         if (this->requiresStiffnessTensor) {
111           tmp.insert(tmp.end(),
112                      {"YoungModulus1", "YoungModulus2", "YoungModulus3", "PoissonRatio12",
113                       "PoissonRatio23", "PoissonRatio13", "ShearModulus12"});
114         }
115         if (this->requiresThermalExpansionCoefficientTensor) {
116           tmp.insert(tmp.end(), {"ThermalExpansion1", "ThermalExpansion2", "ThermalExpansion3"});
117         }
118       } else if (h == ModellingHypothesis::TRIDIMENSIONAL) {
119         if (this->requiresStiffnessTensor) {
120           tmp.insert(tmp.end(), {"YoungModulus1", "YoungModulus2", "YoungModulus3",
121                                  "PoissonRatio12", "PoissonRatio23", "PoissonRatio13",
122                                  "ShearModulus12", "ShearModulus23", "ShearModulus13"});
123         }
124         if (this->requiresThermalExpansionCoefficientTensor) {
125           tmp.insert(tmp.end(), {"ThermalExpansion1", "ThermalExpansion2", "ThermalExpansion3"});
126         }
127       } else {
128         throw_if(true, "unsupported modelling hypothesis");
129       }
130     } else {
131       throw_if(true,
132                "unsupported behaviour type "
133                "(neither isotropic nor orthotropic)");
134     }
135     this->mpnames.insert(this->mpnames.begin(), tmp.begin(), tmp.end());
136   }
137 
getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> & r) const138   tfel::math::tmatrix<3u, 3u, real> AbaqusStandardBehaviour::getRotationMatrix(
139       const tfel::math::vector<real>&, const tfel::math::tmatrix<3u, 3u, real>& r) const {
140     return r;
141   }  // end of AbaqusStandardBehaviour::getRotationMatrix
142 
allocate(BehaviourWorkSpace & wk) const143   void AbaqusStandardBehaviour::allocate(BehaviourWorkSpace& wk) const {
144     const auto ndv = this->getGradientsSize();
145     const auto nth = this->getThermodynamicForcesSize();
146     const auto nstatev = this->getInternalStateVariablesSize();
147     wk.D.resize(nth, nth);
148     wk.kt.resize(nth, ndv);
149     wk.k.resize(nth, ndv);
150     wk.mps.resize(this->mpnames.size() == 0 ? 1u : this->mpnames.size(), real(0));
151     wk.ivs.resize(nstatev);
152     wk.nk.resize(nth, ndv);
153     wk.ne.resize(ndv);
154     wk.ns.resize(nth);
155     wk.nivs.resize(nstatev);
156     mtest::allocate(wk.cs, this->shared_from_this());
157   }  // end of AbaqusStandardBehaviour::allocate
158 
getDefaultStiffnessMatrixType() const159   StiffnessMatrixType AbaqusStandardBehaviour::getDefaultStiffnessMatrixType() const {
160     return StiffnessMatrixType::CONSISTENTTANGENTOPERATOR;
161   }
162 
computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const163   std::pair<bool, real> AbaqusStandardBehaviour::computePredictionOperator(
164       BehaviourWorkSpace& wk,
165       const CurrentState& s,
166       const StiffnessMatrixType ktype) const {
167     if (ktype == StiffnessMatrixType::ELASTICSTIFNESSFROMMATERIALPROPERTIES) {
168       return {false, real(-1)};
169     }
170     wk.cs = s;
171     return this->call_behaviour(wk.kt, wk.cs, wk, real(1), ktype, false);
172   }
173 
integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const174   std::pair<bool, real> AbaqusStandardBehaviour::integrate(
175       CurrentState& s,
176       BehaviourWorkSpace& wk,
177       const real dt,
178       const StiffnessMatrixType ktype) const {
179     return this->call_behaviour(wk.k, s, wk, dt, ktype, true);
180   }  // end of AbaqusStandardBehaviour::integrate
181 
182   AbaqusStandardBehaviour::~AbaqusStandardBehaviour() = default;
183 
184 }  // end of namespace mtest
185