1 /*! 2 * \file mfront/src/CastemStressFreeExpansionHandler.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 04 mars 2014 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 <stdexcept> 17 #include "TFEL/Raise.hxx" 18 #include "TFEL/Math/General/MathConstants.hxx" 19 #include "MFront/Castem/CastemStressFreeExpansionHandler.hxx" 20 21 namespace castem { 22 CastemStandardSmallStrainStressFreeExpansionHandler(CastemReal * const e,CastemReal * const de,const CastemReal * const s0,const CastemReal * const s1,const CastemInt d)23 void CastemStandardSmallStrainStressFreeExpansionHandler( 24 CastemReal *const e, 25 CastemReal *const de, 26 const CastemReal *const s0, 27 const CastemReal *const s1, 28 const CastemInt d) { 29 constexpr const auto cste = tfel::math::Cste<CastemReal>::sqrt2; 30 e[0] -= s0[0]; 31 e[1] -= s0[1]; 32 e[2] -= s0[2]; 33 de[0] -= (s1[0] - s0[0]); 34 de[1] -= (s1[1] - s0[1]); 35 de[2] -= (s1[2] - s0[2]); 36 if (d == 2) { 37 e[3] -= s0[3] * cste; 38 de[3] -= (s1[3] - s0[3]) * cste; 39 } 40 if (d == 3) { 41 e[3] -= s0[3] * cste; 42 e[4] -= s0[3] * cste; 43 e[5] -= s0[5] * cste; 44 de[3] -= (s1[3] - s0[3]) * cste; 45 de[4] -= (s1[4] - s0[4]) * cste; 46 de[5] -= (s1[5] - s0[5]) * cste; 47 } 48 } // end of CastemStandardSmallStrainStressFreeExpansionHandler 49 CastemLogarithmicStrainStressFreeExpansionHandler(CastemReal * const e,CastemReal * const de,const CastemReal * const s0,const CastemReal * const s1,const CastemInt d)50 void CastemLogarithmicStrainStressFreeExpansionHandler( 51 CastemReal *const e, 52 CastemReal *const de, 53 const CastemReal *const s0, 54 const CastemReal *const s1, 55 const CastemInt d) { 56 constexpr const auto cste = tfel::math::Cste<CastemReal>::sqrt2; 57 CastemReal log_s0[6]; 58 CastemReal log_s1[6]; 59 if (d == 1) { 60 log_s0[0] = std::log1p(s0[0]); 61 log_s0[1] = std::log1p(s0[1]); 62 log_s0[2] = std::log1p(s0[2]); 63 log_s1[0] = std::log1p(s1[0]); 64 log_s1[1] = std::log1p(s1[1]); 65 log_s1[2] = std::log1p(s1[2]); 66 } else if (d == 2) { 67 tfel::raise_if( 68 (std::abs(s0[3]) > 10 * std::numeric_limits<CastemReal>::min()) || 69 (std::abs(s1[3]) > 10 * std::numeric_limits<CastemReal>::min()), 70 "CastemLogarithmicStrainStressFreeExpansionHandler: " 71 "stress free expansion is assumed diagonal"); 72 log_s0[0] = std::log1p(s0[0]); 73 log_s0[1] = std::log1p(s0[1]); 74 log_s0[2] = std::log1p(s0[2]); 75 log_s0[3] = CastemReal(0); 76 log_s1[0] = std::log1p(s1[0]); 77 log_s1[1] = std::log1p(s1[1]); 78 log_s1[2] = std::log1p(s1[2]); 79 log_s1[3] = CastemReal(0); 80 } else if (d == 3) { 81 tfel::raise_if( 82 (std::abs(s0[3]) > 10 * std::numeric_limits<CastemReal>::min()) || 83 (std::abs(s1[3]) > 10 * std::numeric_limits<CastemReal>::min()) || 84 (std::abs(s0[4]) > 10 * std::numeric_limits<CastemReal>::min()) || 85 (std::abs(s1[4]) > 10 * std::numeric_limits<CastemReal>::min()) || 86 (std::abs(s0[5]) > 10 * std::numeric_limits<CastemReal>::min()) || 87 (std::abs(s1[5]) > 10 * std::numeric_limits<CastemReal>::min()), 88 "CastemLogarithmicStrainStressFreeExpansionHandler: " 89 "stress free expansion is assumed diagonal"); 90 log_s0[0] = std::log1p(s0[0]); 91 log_s0[1] = std::log1p(s0[1]); 92 log_s0[2] = std::log1p(s0[2]); 93 log_s0[3] = CastemReal(0); 94 log_s0[4] = CastemReal(0); 95 log_s0[5] = CastemReal(0); 96 log_s1[0] = std::log1p(s1[0]); 97 log_s1[1] = std::log1p(s1[1]); 98 log_s1[2] = std::log1p(s1[2]); 99 log_s1[3] = CastemReal(0); 100 log_s1[4] = CastemReal(0); 101 log_s1[5] = CastemReal(0); 102 } else { 103 tfel::raise( 104 "CastemLogarithmicStrainStressFreeExpansionHandler: " 105 "invalid dimension"); 106 } 107 e[0] -= log_s0[0]; 108 e[1] -= log_s0[1]; 109 e[2] -= log_s0[2]; 110 de[0] -= (log_s1[0] - log_s0[0]); 111 de[1] -= (log_s1[1] - log_s0[1]); 112 de[2] -= (log_s1[2] - log_s0[2]); 113 if (d == 2) { 114 e[3] -= log_s0[3] * cste; 115 de[3] -= (log_s1[3] - log_s0[3]) * cste; 116 } 117 if (d == 3) { 118 e[3] -= log_s0[3] * cste; 119 e[4] -= log_s0[3] * cste; 120 e[5] -= log_s0[5] * cste; 121 de[3] -= (log_s1[3] - log_s0[3]) * cste; 122 de[4] -= (log_s1[4] - log_s0[4]) * cste; 123 de[5] -= (log_s1[5] - log_s0[5]) * cste; 124 } 125 } // end of CastemLogarithmicStrainStressFreeExpansionHandler 126 127 } // end of namespace castem 128