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