1 /*!
2  * \file   mtest/src/DianaFEASmallStrainBehaviour.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 <cstring>
17 #include <algorithm>
18 
19 #include "TFEL/Raise.hxx"
20 #include "TFEL/Math/tmatrix.hxx"
21 #include "TFEL/Math/stensor.hxx"
22 #include "TFEL/Math/st2tost2.hxx"
23 #include "TFEL/Math/ST2toST2/ST2toST2View.hxx"
24 #include "TFEL/System/ExternalLibraryManager.hxx"
25 #include "MFront/DianaFEA/DianaFEA.hxx"
26 #include "MFront/MFrontLogStream.hxx"
27 #include "MFront/DianaFEA/DianaFEAComputeStiffnessTensor.hxx"
28 
29 #include "MTest/Evolution.hxx"
30 #include "MTest/CurrentState.hxx"
31 #include "MTest/BehaviourWorkSpace.hxx"
32 #include "MTest/DianaFEASmallStrainBehaviour.hxx"
33 #include "MTest/UmatNormaliseTangentOperator.hxx"
34 
35 namespace mtest {
36 
DianaFEASmallStrainBehaviour(const Hypothesis h,const std::string & l,const std::string & b)37   DianaFEASmallStrainBehaviour::DianaFEASmallStrainBehaviour(
38       const Hypothesis h, const std::string& l, const std::string& b)
39       : StandardBehaviourBase(h, l, b) {
40     auto& elm = tfel::system::ExternalLibraryManager::getExternalLibraryManager();
41     this->fct = elm.getDianaFEAExternalBehaviourFunction(l, b);
42     if (this->stype != 0) {
43       tfel::raise(
44           "DianaFEASmallStrainBehaviour::DianaFEASmallStrainBehaviour: "
45           "only isotropic behaviours are supported");
46     }
47     auto tmp = std::vector<std::string>{};
48     if (this->requiresStiffnessTensor) {
49       tmp.insert(tmp.end(), {"YoungModulus", "PoissonRatio"});
50     }
51     if (this->requiresThermalExpansionCoefficientTensor) {
52       tmp.push_back("ThermalExpansion");
53     }
54     this->mpnames.insert(this->mpnames.begin(),tmp.begin(),tmp.end());
55   }  // end of DianaFEASmallStrainBehaviour::DianaFEASmallStrainBehaviour
56 
allocate(BehaviourWorkSpace & wk) const57   void DianaFEASmallStrainBehaviour::allocate(BehaviourWorkSpace& wk) const{
58     const auto ndv     = this->getGradientsSize();
59     const auto nth     = this->getThermodynamicForcesSize();
60     const auto nstatev = this->getInternalStateVariablesSize();
61     wk.D.resize(nth,nth);
62     wk.k.resize(nth,ndv);
63     wk.kt.resize(nth,ndv);
64     wk.ivs.resize(nstatev);
65     wk.nk.resize(nth,ndv);
66     wk.ne.resize(ndv);
67     wk.ns.resize(nth);
68     wk.nivs.resize(nstatev);
69     wk.mps.resize(this->mpnames.size());
70     mtest::allocate(wk.cs,this->shared_from_this());
71   } // end of DianaFEASmallStrainBehaviour::allocate
72 
getRotationMatrix(const tfel::math::vector<real> &,const tfel::math::tmatrix<3u,3u,real> & r) const73   tfel::math::tmatrix<3u, 3u, real> DianaFEASmallStrainBehaviour::getRotationMatrix(
74       const tfel::math::vector<real>&, const tfel::math::tmatrix<3u, 3u, real>& r) const {
75     return r;
76   }  // end of DianaFEASmallStrainBehaviour::getRotationMatrix
77 
getDefaultStiffnessMatrixType() const78   StiffnessMatrixType DianaFEASmallStrainBehaviour::getDefaultStiffnessMatrixType() const {
79     return StiffnessMatrixType::CONSISTENTTANGENTOPERATOR;
80   } // end of DianaFEASmallStrainBehaviour::getDefaultStiffnessMatrixType
81 
getGradientsDefaultInitialValues(tfel::math::vector<real> & v) const82   void DianaFEASmallStrainBehaviour::getGradientsDefaultInitialValues(
83       tfel::math::vector<real>& v) const {
84     std::fill(v.begin(), v.end(), real(0));
85   }  // end of DianaFEASmallStrainBehaviour::setGradientsDefaultInitialValue
86 
computePredictionOperator(BehaviourWorkSpace & wk,const CurrentState & s,const StiffnessMatrixType ktype) const87   std::pair<bool, real> DianaFEASmallStrainBehaviour::computePredictionOperator(
88       BehaviourWorkSpace& wk,
89       const CurrentState& s,
90       const StiffnessMatrixType ktype) const {
91     wk.cs = s;
92     return this->call_behaviour(wk.kt, wk.cs, wk, real(1), ktype, false);
93   }  // end of DianaFEASmallStrainBehaviour::computePredictionOperator
94 
integrate(CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype) const95   std::pair<bool, real> DianaFEASmallStrainBehaviour::integrate(
96       CurrentState& s,
97       BehaviourWorkSpace& wk,
98       const real dt,
99       const StiffnessMatrixType ktype) const {
100     return this->call_behaviour(wk.k, s, wk, dt, ktype, true);
101   }  // end of DianaFEASmallStrainBehaviour::integrate
102 
call_behaviour(tfel::math::matrix<real> & Kt,CurrentState & s,BehaviourWorkSpace & wk,const real dt,const StiffnessMatrixType ktype,const bool b) const103   std::pair<bool, real> DianaFEASmallStrainBehaviour::call_behaviour(
104       tfel::math::matrix<real>& Kt,
105       CurrentState& s,
106       BehaviourWorkSpace& wk,
107       const real dt,
108       const StiffnessMatrixType ktype,
109       const bool b) const {
110     using namespace std;
111     using namespace tfel::math;
112     using namespace dianafea;
113     using tfel::math::vector;
114     using tfel::material::getSpaceDimension;
115     constexpr const auto sqrt2 = Cste<real>::sqrt2;
116     auto throw_if = [](const bool c, const std::string& m) {
117       tfel::raise_if(c, "DianaFEASmallStrainBehaviour::call_behaviour: " + m);
118     };
119     throw_if((wk.D.getNbRows() != Kt.getNbRows()) ||
120                  (wk.D.getNbCols() != Kt.getNbCols()),
121              "the memory has not been allocated correctly");
122     throw_if(s.iv0.size() != wk.ivs.size(),
123              "the memory has not been allocated correctly (2)");
124     std::fill(wk.D.begin(), wk.D.end(), 0.);
125     // choosing the type of stiffness matrix
126     StandardBehaviourBase::initializeTangentOperator(wk.D, ktype, b);
127     // state variable initial values
128     std::copy(s.iv0.begin(), s.iv0.end(), wk.ivs.begin());
129     // strain
130     stensor<3u, real> ue0(real(0));
131     stensor<3u, real> ude(real(0));
132     std::copy(s.e0.begin(), s.e0.end(), ue0.begin());
133     for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) {
134       ude(i) = s.e1(i) - s.e0(i);
135     }
136     std::copy(s.s0.begin(), s.s0.end(), s.s1.begin());
137     // thermal strain
138     for (decltype(s.e1.size()) i = 0; i != s.e1.size(); ++i) {
139       ue0(i) -= s.e_th0(i);
140       ude(i) -= s.e_th1(i) - s.e_th0(i);
141     }
142     // convert from umat conventions
143     for (decltype(s.e1.size()) i = 3; i != this->getThermodynamicForcesSize(); ++i) {
144       s.s1(i) /= sqrt2;
145       ue0(i) *= sqrt2;
146       ude(i) *= sqrt2;
147     }
148     const auto nprops = static_cast<DianaFEAInt>(s.mprops1.size());
149     const auto ntens =
150         static_cast<DianaFEAInt>(this->getThermodynamicForcesSize());
151     const auto nstatv = static_cast<DianaFEAInt>(wk.ivs.size());
152     (this->fct)(&(s.s1(0)), &(wk.D(0, 0)),
153                 wk.ivs.empty() ? nullptr : &wk.ivs(0), &ntens, &nprops, &nstatv,
154                 &ue0(0), &ude(0), &dt,
155                 s.mprops1.empty() ? nullptr : &s.mprops1(0), &(s.esv0(0)),
156                 &(s.desv(0)));
157     // tangent operator
158     if (ktype != StiffnessMatrixType::NOSTIFFNESS) {
159       UmatNormaliseTangentOperator::exe(
160           &Kt(0, 0), wk.D, getSpaceDimension(this->getHypothesis()));
161     }
162     std::copy(wk.ivs.begin(), wk.ivs.end(), s.iv1.begin());
163     // turning things in standard conventions
164     for (decltype(s.e1.size()) i = 3; i != this->getThermodynamicForcesSize();
165          ++i) {
166       s.s1(i) *= sqrt2;
167     }
168     return {true, std::numeric_limits<DianaFEAReal>::max()};
169   }  // end of DianaFEASmallStrainBehaviour::integrate
170 
171   std::vector<std::string>
getOptionalMaterialProperties() const172   DianaFEASmallStrainBehaviour::getOptionalMaterialProperties() const {
173     if (this->stype != 0) {
174       tfel::raise(
175           "DianaFEASmallStrainBehaviour::"
176           "setOptionalMaterialProperties: "
177           "unsupported symmetry type");
178     }
179     return std::vector<std::string>(1u, "ThermalExpansion");
180   }  // end of DianaFEASmallStrainBehaviour::getOptionalMaterialProperties
181 
setOptionalMaterialPropertiesDefaultValues(EvolutionManager & mp,const EvolutionManager & evm) const182   void DianaFEASmallStrainBehaviour::setOptionalMaterialPropertiesDefaultValues(
183       EvolutionManager& mp, const EvolutionManager& evm) const {
184     if (this->stype == 0) {
185       Behaviour::setOptionalMaterialPropertyDefaultValue(
186           mp, evm, "ThermalExpansion", 0.);
187     } else {
188       tfel::raise(
189           "DianaFEASmallStrainBehaviour::"
190           "setOptionalMaterialPropertiesDefaultValues: "
191           "unsupported symmetry type");
192     }
193   }  // end of
194   // DianaFEASmallStrainBehaviour::setOptionalMaterialPropertiesDefaultValues
195 
196   DianaFEASmallStrainBehaviour::~DianaFEASmallStrainBehaviour() = default;
197 
198 }  // end of namespace mtest
199