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