1 /*! 2 * \file mfront/src/LSDYNAFiniteStrain.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 21/08/2016 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<iostream> 17 #include<stdexcept> 18 19 #include"TFEL/Raise.hxx" 20 #include"TFEL/Math/tensor.hxx" 21 #include"TFEL/Math/stensor.hxx" 22 #include"TFEL/Math/tmatrix.hxx" 23 #include"TFEL/Math/stensor.hxx" 24 #include"TFEL/Math/st2tost2.hxx" 25 #include"TFEL/Math/Stensor/StensorView.hxx" 26 #include"TFEL/Math/Stensor/ConstStensorView.hxx" 27 #include"TFEL/Math/ST2toST2/ST2toST2View.hxx" 28 #include"TFEL/Math/ST2toST2/ConstST2toST2View.hxx" 29 #include"TFEL/Math/ST2toST2/UmatNormaliseTangentOperator.hxx" 30 #include"TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx" 31 #include"MFront/LSDYNA/LSDYNATangentOperator.hxx" 32 #include"MFront/LSDYNA/LSDYNAFiniteStrain.hxx" 33 34 namespace lsdyna 35 { 36 37 void computeGreenLagrangeStrain(LSDYNAReal * const e,const LSDYNAReal * const F,const LSDYNAInt NTENS,const bool ps)38 LSDYNAFiniteStrain::computeGreenLagrangeStrain(LSDYNAReal* const e, 39 const LSDYNAReal* const F, 40 const LSDYNAInt NTENS, 41 const bool ps) 42 { 43 // warning : F is given in the fortran convention 44 // (F0 F3 F6 45 // F1 F4 F7 46 // F2 F5 F8) 47 // maxima : 48 // (0.5*(F2^2+F1^2+F0^2-1) 0.5*(F2*F5+F1*F4+F0*F3) 0.5*(F2*F8+F1*F7+F0*F6) 49 // 0.5*(F2*F5+F1*F4+F0*F3) 0.5*(F5^2+F4^2+F3^2-1) 0.5*(F5*F8+F4*F7+F3*F6) 50 // 0.5*(F2*F8+F1*F7+F0*F6) 0.5*(F5*F8+F4*F7+F3*F6) 0.5*(F8^2+F7^2+F6^2-1) 51 e[0]=0.5*(F[0]*F[0]+F[1]*F[1]+F[2]*F[2]-1); 52 e[1]=0.5*(F[3]*F[3]+F[4]*F[4]+F[5]*F[5]-1); 53 if(!ps){ 54 e[2]=0.5*(F[6]*F[6]+F[7]*F[7]+F[8]*F[8]-1); 55 } else { 56 e[2]=0.; 57 } 58 // warning : e must be computed using Voigt 59 // notations. Extradiagonals terms are mulitplied by two 60 if(NTENS==4){ 61 e[3]= F[2]*F[5]+F[1]*F[4]+F[0]*F[3]; 62 } else if (NTENS==6){ 63 e[3]=F[2]*F[5]+F[1]*F[4]+F[0]*F[3]; 64 e[4]=F[2]*F[8]+F[1]*F[7]+F[0]*F[6]; 65 e[5]=F[5]*F[8]+F[4]*F[7]+F[3]*F[6]; 66 } 67 } // end of LSDYNAFiniteStrain::computeGreenLagrangeStrain 68 69 70 void computeSecondPiolaKirchhoffStressFromCauchyStress(LSDYNAReal * const STRESS,const LSDYNAReal * const F,const LSDYNAInt NTENS,const bool ps,const LSDYNAReal Fzz)71 LSDYNAFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(LSDYNAReal* const STRESS, 72 const LSDYNAReal* const F, 73 const LSDYNAInt NTENS, 74 const bool ps, 75 const LSDYNAReal Fzz) 76 { 77 // warning F is given in the fortran convention 78 // (F0 F3 F6 79 // F1 F4 F7 80 // F2 F5 F8) 81 const LSDYNAReal F0 = F[0]; 82 const LSDYNAReal F1 = F[1]; 83 const LSDYNAReal F2 = F[2]; 84 const LSDYNAReal F3 = F[3]; 85 const LSDYNAReal F4 = F[4]; 86 const LSDYNAReal F5 = F[5]; 87 const LSDYNAReal F6 = F[6]; 88 const LSDYNAReal F7 = F[7]; 89 const LSDYNAReal F8 = !ps ? F[8] : Fzz; 90 // determinant 91 const LSDYNAReal J = F0*(F4*F8-F5*F7)-F3*(F1*F8-F2*F7)+(F1*F5-F2*F4)*F6; 92 // inverse of the determinant 93 const LSDYNAReal iJ = 1/J; 94 // inverse of F 95 // maxima : F4*F8-F5*F7 96 const LSDYNAReal iF0 = iJ*(F4*F8-F5*F7); 97 // maxima : F2*F7-F1*F8 98 const LSDYNAReal iF1 = iJ*(F2*F7-F1*F8); 99 // maxima : F1*F5-F2*F4 100 const LSDYNAReal iF2 = iJ*(F1*F5-F2*F4); 101 // maxima : F5*F6-F3*F8 102 const LSDYNAReal iF3 = iJ*(F5*F6-F3*F8); 103 // maxima : F0*F8-F2*F6 104 const LSDYNAReal iF4 = iJ*(F0*F8-F2*F6); 105 // maxima : F2*F3-F0*F5 106 const LSDYNAReal iF5 = iJ*(F2*F3-F0*F5); 107 // maxima : F3*F7-F4*F6 108 const LSDYNAReal iF6 = iJ*(F3*F7-F4*F6); 109 // maxima : F1*F6-F0*F7 110 const LSDYNAReal iF7 = iJ*(F1*F6-F0*F7); 111 // maxima : F0*F4-F1*F3 112 const LSDYNAReal iF8 = iJ*(F0*F4-F1*F3); 113 // sk2 114 const LSDYNAReal p0 = STRESS[0]; 115 const LSDYNAReal p1 = STRESS[1]; 116 const LSDYNAReal p2 = STRESS[2]; 117 const LSDYNAReal p3 = (NTENS>=4) ? STRESS[3] : LSDYNAReal(0); 118 const LSDYNAReal p4 = (NTENS==6) ? STRESS[4] : LSDYNAReal(0); 119 const LSDYNAReal p5 = (NTENS==6) ? STRESS[5] : LSDYNAReal(0); 120 STRESS[0] = J *(iF3*(p5*iF6+p1*iF3+p3*iF0)+iF0*(p4*iF6+p3*iF3+p0*iF0)+iF6*(p2*iF6+p5*iF3+p4*iF0)); 121 STRESS[1] = J*(iF4*(p5*iF7+p1*iF4+p3*iF1)+iF1*(p4*iF7+p3*iF4+p0*iF1)+iF7*(p2*iF7+p5*iF4+p4*iF1)); 122 STRESS[2] = J*(iF5*(p5*iF8+p1*iF5+p3*iF2)+iF2*(p4*iF8+p3*iF5+p0*iF2)+iF8*(p2*iF8+p5*iF5+p4*iF2)); 123 if(NTENS==4){ 124 if(ps){ 125 STRESS[2] = 0; 126 } 127 STRESS[3] = J*(iF3*(p5*iF7+p1*iF4+p3*iF1)+iF0*(p4*iF7+p3*iF4+p0*iF1)+iF6*(p2*iF7+p5*iF4+p4*iF1)); 128 } else if(NTENS==6){ 129 STRESS[3] = J*(iF3*(p5*iF7+p1*iF4+p3*iF1)+iF0*(p4*iF7+p3*iF4+p0*iF1)+iF6*(p2*iF7+p5*iF4+p4*iF1)); 130 STRESS[4] = J*(iF3*(p5*iF8+p1*iF5+p3*iF2)+iF0*(p4*iF8+p3*iF5+p0*iF2)+iF6*(p2*iF8+p5*iF5+p4*iF2)); 131 STRESS[5] = J*(iF4*(p5*iF8+p1*iF5+p3*iF2)+iF1*(p4*iF8+p3*iF5+p0*iF2)+iF7*(p2*iF8+p5*iF5+p4*iF2)); 132 } 133 } // end of LSDYNAFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress 134 135 void computeSecondPiolaKirchhoffStressFromCauchyStress(LSDYNAReal * const sk2,const LSDYNAReal * const STRESS,const LSDYNAReal * const F,const LSDYNAInt NTENS,const bool ps,const LSDYNAReal Fzz)136 LSDYNAFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(LSDYNAReal* const sk2, 137 const LSDYNAReal* const STRESS, 138 const LSDYNAReal* const F, 139 const LSDYNAInt NTENS, 140 const bool ps, 141 const LSDYNAReal Fzz) 142 { 143 sk2[0] = STRESS[0]; 144 sk2[1] = STRESS[1]; 145 sk2[2] = STRESS[2]; 146 if(NTENS==4){ 147 sk2[3] = STRESS[3]; 148 } else if(NTENS==6){ 149 sk2[3] = STRESS[3]; 150 sk2[4] = STRESS[4]; 151 sk2[5] = STRESS[5]; 152 } 153 LSDYNAFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(sk2,F,NTENS,ps,Fzz); 154 } // end of LSDYNAFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress 155 156 void computeCauchyStressFromSecondPiolaKirchhoffStress(LSDYNAReal * const s,const LSDYNAReal * const F,const LSDYNAInt NTENS,const bool ps,const LSDYNAReal Fzz)157 LSDYNAFiniteStrain::computeCauchyStressFromSecondPiolaKirchhoffStress(LSDYNAReal* const s, 158 const LSDYNAReal* const F, 159 const LSDYNAInt NTENS, 160 const bool ps, 161 const LSDYNAReal Fzz) 162 { 163 // warning F is given in the fortran convention 164 // (F0 F3 F6 165 // F1 F4 F7 166 // F2 F5 F8) 167 const LSDYNAReal F0 = F[0]; 168 const LSDYNAReal F1 = F[1]; 169 const LSDYNAReal F2 = F[2]; 170 const LSDYNAReal F3 = F[3]; 171 const LSDYNAReal F4 = F[4]; 172 const LSDYNAReal F5 = F[5]; 173 const LSDYNAReal F6 = F[6]; 174 const LSDYNAReal F7 = F[7]; 175 const LSDYNAReal F8 = !ps ? F[8] : Fzz; 176 // determinant 177 const LSDYNAReal inv_J = 1/(F0*(F4*F8-F5*F7)-F3*(F1*F8-F2*F7)+(F1*F5-F2*F4)*F6); 178 const LSDYNAReal p0 = s[0]; 179 const LSDYNAReal p1 = s[1]; 180 const LSDYNAReal p2 = s[2]; 181 const LSDYNAReal p3 = (NTENS>=4) ? s[3] : LSDYNAReal(0); 182 const LSDYNAReal p4 = (NTENS==6) ? s[4] : LSDYNAReal(0); 183 const LSDYNAReal p5 = (NTENS==6) ? s[5] : LSDYNAReal(0); 184 // maxima : F3*(p5*F6+p1*F3+p3*F0)+F0*(p4*F6+p3*F3+p0*F0)+F6*(p2*F6+p5*F3+p4*F0) 185 s[0] = inv_J*(F3*(p5*F6+p1*F3+p3*F0)+F0*(p4*F6+p3*F3+p0*F0)+F6*(p2*F6+p5*F3+p4*F0)); 186 // maxima : F4*(p5*F7+p1*F4+p3*F1)+F1*(p4*F7+p3*F4+p0*F1)+F7*(p2*F7+p5*F4+p4*F1) 187 s[1] = inv_J*(F4*(p5*F7+p1*F4+p3*F1)+F1*(p4*F7+p3*F4+p0*F1)+F7*(p2*F7+p5*F4+p4*F1)); 188 // maxima : F5*(p5*F8+p1*F5+p3*F2)+F2*(p4*F8+p3*F5+p0*F2)+F8*(p2*F8+p5*F5+p4*F2) 189 s[2] = inv_J*(F5*(p5*F8+p1*F5+p3*F2)+F2*(p4*F8+p3*F5+p0*F2)+F8*(p2*F8+p5*F5+p4*F2)); 190 if(NTENS==4){ 191 if(ps){ 192 s[2] = 0; 193 } 194 // maxima : F3*(p5*F7+p1*F4+p3*F1)+F0*(p4*F7+p3*F4+p0*F1)+F6*(p2*F7+p5*F4+p4*F1) 195 s[3] = inv_J*(F3*(p5*F7+p1*F4+p3*F1)+F0*(p4*F7+p3*F4+p0*F1)+F6*(p2*F7+p5*F4+p4*F1)); 196 } else if(NTENS==6){ 197 // maxima : F3*(p5*F7+p1*F4+p3*F1)+F0*(p4*F7+p3*F4+p0*F1)+F6*(p2*F7+p5*F4+p4*F1) 198 s[3] = inv_J*(F3*(p5*F7+p1*F4+p3*F1)+F0*(p4*F7+p3*F4+p0*F1)+F6*(p2*F7+p5*F4+p4*F1)); 199 // maxima : F3*(p5*F8+p1*F5+p3*F2)+F0*(p4*F8+p3*F5+p0*F2)+F6*(p2*F8+p5*F5+p4*F2) 200 s[4] = inv_J*(F3*(p5*F8+p1*F5+p3*F2)+F0*(p4*F8+p3*F5+p0*F2)+F6*(p2*F8+p5*F5+p4*F2)); 201 // maxima : F4*(p5*F8+p1*F5+p3*F2)+F1*(p4*F8+p3*F5+p0*F2)+F7*(p2*F8+p5*F5+p4*F2) 202 s[5] = inv_J*(F4*(p5*F8+p1*F5+p3*F2)+F1*(p4*F8+p3*F5+p0*F2)+F7*(p2*F8+p5*F5+p4*F2)); 203 } 204 } // end of LSDYNAFiniteStrain::computeCauchyStressFromSecondPiolaKirchhoffStress 205 206 } // end of namespace lsdyna 207