1 /*! 2 * \file mfront/src/AnsysFiniteStrain.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/Ansys/AnsysTangentOperator.hxx" 32 #include"MFront/Ansys/AnsysFiniteStrain.hxx" 33 34 namespace ansys 35 { 36 37 void computeGreenLagrangeStrain(AnsysReal * const e,const AnsysReal * const F,const AnsysInt NTENS,const bool ps)38 AnsysFiniteStrain::computeGreenLagrangeStrain(AnsysReal* const e, 39 const AnsysReal* const F, 40 const AnsysInt 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]=(F[0]*F[0]+F[1]*F[1]+F[2]*F[2]-1)/2; 52 e[1]=(F[3]*F[3]+F[4]*F[4]+F[5]*F[5]-1)/2; 53 if(!ps){ 54 e[2]=(F[6]*F[6]+F[7]*F[7]+F[8]*F[8]-1)/2; 55 } else { 56 e[2]=0.; 57 } 58 if(NTENS==4){ 59 e[3]= F[2]*F[5]+F[1]*F[4]+F[0]*F[3]; 60 } else if (NTENS==6){ 61 e[3]= F[2]*F[5]+F[1]*F[4]+F[0]*F[3]/2; 62 // Ansys ordering conventions differs from TFEL' ones 63 e[5]= F[2]*F[8]+F[1]*F[7]+F[0]*F[6]; 64 e[4]= F[5]*F[8]+F[4]*F[7]+F[3]*F[6]; 65 } 66 } // end of AnsysFiniteStrain::computeGreenLagrangeStrain 67 68 69 void computeSecondPiolaKirchhoffStressFromCauchyStress(AnsysReal * const STRESS,const AnsysReal * const F,const AnsysInt NTENS,const bool ps,const AnsysReal Fzz)70 AnsysFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(AnsysReal* const STRESS, 71 const AnsysReal* const F, 72 const AnsysInt NTENS, 73 const bool ps, 74 const AnsysReal Fzz) 75 { 76 // warning F is given in the fortran convention 77 // (F0 F3 F6 78 // F1 F4 F7 79 // F2 F5 F8) 80 const AnsysReal F0 = F[0]; 81 const AnsysReal F1 = F[1]; 82 const AnsysReal F2 = F[2]; 83 const AnsysReal F3 = F[3]; 84 const AnsysReal F4 = F[4]; 85 const AnsysReal F5 = F[5]; 86 const AnsysReal F6 = F[6]; 87 const AnsysReal F7 = F[7]; 88 const AnsysReal F8 = !ps ? F[8] : Fzz; 89 // determinant 90 const AnsysReal J = F0*(F4*F8-F5*F7)-F3*(F1*F8-F2*F7)+(F1*F5-F2*F4)*F6; 91 // inverse of the determinant 92 const AnsysReal iJ = 1/J; 93 // inverse of F 94 // maxima : F4*F8-F5*F7 95 const AnsysReal iF0 = iJ*(F4*F8-F5*F7); 96 // maxima : F2*F7-F1*F8 97 const AnsysReal iF1 = iJ*(F2*F7-F1*F8); 98 // maxima : F1*F5-F2*F4 99 const AnsysReal iF2 = iJ*(F1*F5-F2*F4); 100 // maxima : F5*F6-F3*F8 101 const AnsysReal iF3 = iJ*(F5*F6-F3*F8); 102 // maxima : F0*F8-F2*F6 103 const AnsysReal iF4 = iJ*(F0*F8-F2*F6); 104 // maxima : F2*F3-F0*F5 105 const AnsysReal iF5 = iJ*(F2*F3-F0*F5); 106 // maxima : F3*F7-F4*F6 107 const AnsysReal iF6 = iJ*(F3*F7-F4*F6); 108 // maxima : F1*F6-F0*F7 109 const AnsysReal iF7 = iJ*(F1*F6-F0*F7); 110 // maxima : F0*F4-F1*F3 111 const AnsysReal iF8 = iJ*(F0*F4-F1*F3); 112 // sk2 113 const AnsysReal p0 = STRESS[0]; 114 const AnsysReal p1 = STRESS[1]; 115 const AnsysReal p2 = STRESS[2]; 116 const AnsysReal p3 = (NTENS>=4) ? STRESS[3] : AnsysReal(0); 117 const AnsysReal p4 = (NTENS==6) ? STRESS[5] : AnsysReal(0); 118 const AnsysReal p5 = (NTENS==6) ? STRESS[4] : AnsysReal(0); 119 STRESS[0] = J *(iF3*(p5*iF6+p1*iF3+p3*iF0)+iF0*(p4*iF6+p3*iF3+p0*iF0)+iF6*(p2*iF6+p5*iF3+p4*iF0)); 120 STRESS[1] = J*(iF4*(p5*iF7+p1*iF4+p3*iF1)+iF1*(p4*iF7+p3*iF4+p0*iF1)+iF7*(p2*iF7+p5*iF4+p4*iF1)); 121 STRESS[2] = J*(iF5*(p5*iF8+p1*iF5+p3*iF2)+iF2*(p4*iF8+p3*iF5+p0*iF2)+iF8*(p2*iF8+p5*iF5+p4*iF2)); 122 if(NTENS==4){ 123 if(ps){ 124 STRESS[2] = 0; 125 } 126 STRESS[3] = J*(iF3*(p5*iF7+p1*iF4+p3*iF1)+iF0*(p4*iF7+p3*iF4+p0*iF1)+iF6*(p2*iF7+p5*iF4+p4*iF1)); 127 } else if(NTENS==6){ 128 STRESS[3] = J*(iF3*(p5*iF7+p1*iF4+p3*iF1)+iF0*(p4*iF7+p3*iF4+p0*iF1)+iF6*(p2*iF7+p5*iF4+p4*iF1)); 129 STRESS[5] = J*(iF3*(p5*iF8+p1*iF5+p3*iF2)+iF0*(p4*iF8+p3*iF5+p0*iF2)+iF6*(p2*iF8+p5*iF5+p4*iF2)); 130 STRESS[4] = J*(iF4*(p5*iF8+p1*iF5+p3*iF2)+iF1*(p4*iF8+p3*iF5+p0*iF2)+iF7*(p2*iF8+p5*iF5+p4*iF2)); 131 } 132 } // end of AnsysFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress 133 134 void computeSecondPiolaKirchhoffStressFromCauchyStress(AnsysReal * const sk2,const AnsysReal * const STRESS,const AnsysReal * const F,const AnsysInt NTENS,const bool ps,const AnsysReal Fzz)135 AnsysFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(AnsysReal* const sk2, 136 const AnsysReal* const STRESS, 137 const AnsysReal* const F, 138 const AnsysInt NTENS, 139 const bool ps, 140 const AnsysReal Fzz) 141 { 142 sk2[0] = STRESS[0]; 143 sk2[1] = STRESS[1]; 144 sk2[2] = STRESS[2]; 145 if(NTENS==4){ 146 sk2[3] = STRESS[3]; 147 } else if(NTENS==6){ 148 sk2[3] = STRESS[3]; 149 sk2[4] = STRESS[4]; 150 sk2[5] = STRESS[5]; 151 } 152 AnsysFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress(sk2,F,NTENS,ps,Fzz); 153 } // end of AnsysFiniteStrain::computeSecondPiolaKirchhoffStressFromCauchyStress 154 155 void computeCauchyStressFromSecondPiolaKirchhoffStress(AnsysReal * const s,const AnsysReal * const F,const AnsysInt NTENS,const bool ps,const AnsysReal Fzz)156 AnsysFiniteStrain::computeCauchyStressFromSecondPiolaKirchhoffStress(AnsysReal* const s, 157 const AnsysReal* const F, 158 const AnsysInt NTENS, 159 const bool ps, 160 const AnsysReal Fzz) 161 { 162 // warning F is given in the fortran convention 163 // (F0 F3 F6 164 // F1 F4 F7 165 // F2 F5 F8) 166 const AnsysReal F0 = F[0]; 167 const AnsysReal F1 = F[1]; 168 const AnsysReal F2 = F[2]; 169 const AnsysReal F3 = F[3]; 170 const AnsysReal F4 = F[4]; 171 const AnsysReal F5 = F[5]; 172 const AnsysReal F6 = F[6]; 173 const AnsysReal F7 = F[7]; 174 const AnsysReal F8 = !ps ? F[8] : Fzz; 175 // determinant 176 const AnsysReal inv_J = 1/(F0*(F4*F8-F5*F7)-F3*(F1*F8-F2*F7)+(F1*F5-F2*F4)*F6); 177 // pk2 values using standard TFEL ordering 178 const AnsysReal p0 = s[0]; 179 const AnsysReal p1 = s[1]; 180 const AnsysReal p2 = s[2]; 181 const AnsysReal p3 = (NTENS>=4) ? s[3] : AnsysReal(0); 182 const AnsysReal p4 = (NTENS==6) ? s[5] : AnsysReal(0); 183 const AnsysReal p5 = (NTENS==6) ? s[4] : AnsysReal(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 // Ansys ordering conventions differs from TFEL' ones 200 // maxima : F3*(p5*F8+p1*F5+p3*F2)+F0*(p4*F8+p3*F5+p0*F2)+F6*(p2*F8+p5*F5+p4*F2) 201 s[5] = inv_J*(F3*(p5*F8+p1*F5+p3*F2)+F0*(p4*F8+p3*F5+p0*F2)+F6*(p2*F8+p5*F5+p4*F2)); 202 // maxima : F4*(p5*F8+p1*F5+p3*F2)+F1*(p4*F8+p3*F5+p0*F2)+F7*(p2*F8+p5*F5+p4*F2) 203 s[4] = inv_J*(F4*(p5*F8+p1*F5+p3*F2)+F1*(p4*F8+p3*F5+p0*F2)+F7*(p2*F8+p5*F5+p4*F2)); 204 } 205 } // end of AnsysFiniteStrain::computeCauchyStressFromSecondPiolaKirchhoffStress 206 207 void computeAnsysTangentOperatorFromCSE(AnsysReal * const D,const AnsysReal * const CSE,const AnsysReal * const F,const AnsysReal * const STRESS,const AnsysInt NTENS,const bool ps)208 AnsysFiniteStrain::computeAnsysTangentOperatorFromCSE(AnsysReal* const D, 209 const AnsysReal* const CSE, 210 const AnsysReal* const F, 211 const AnsysReal* const STRESS, 212 const AnsysInt NTENS, 213 const bool ps){ 214 using namespace tfel::math; 215 using namespace tfel::material; 216 using TangentOperator = FiniteStrainBehaviourTangentOperatorBase; 217 if(ps){ 218 std::cout << "AnsysFiniteStrain::computeAnsysTangentOperatorFromCSE: " 219 << "planeStress support is not yet implemented" << std::endl; 220 std::exit(-1); 221 } 222 if(NTENS==4){ 223 using tensor = tensor<2u,AnsysReal>; 224 stensor<2u,AnsysReal> s; 225 s.importTab(STRESS); 226 auto F1 = tensor::buildFromFortranMatrix(F); 227 const auto C = AnsysTangentOperator<AnsysReal>::convert2D(CSE); 228 ST2toST2View<2u,AnsysReal> Ca(D); 229 Ca = convert<TangentOperator::ABAQUS, 230 TangentOperator::DS_DEGL>(C,tensor::Id(),F1,s); 231 // now convert to ansys 232 AnsysTangentOperator<AnsysReal>::normalize(Ca); 233 } else if(NTENS==6){ 234 using tensor = tensor<3u,AnsysReal>; 235 stensor<3u,AnsysReal> s; 236 s.importTab(STRESS); 237 auto F1 = tensor::buildFromFortranMatrix(F); 238 const auto C = AnsysTangentOperator<AnsysReal>::convert3D(CSE); 239 ST2toST2View<3u,AnsysReal> Ca(D); 240 Ca = convert<TangentOperator::ABAQUS, 241 TangentOperator::DS_DEGL>(C,tensor::Id(),F1,s); 242 // now convert to ansys 243 AnsysTangentOperator<AnsysReal>::normalize(Ca); 244 } 245 } // end of AnsysFiniteStrain::computeAnsysTangentOperatorFromCSE 246 247 } // end of namespace ansys 248