1/*! 2 * \file Barlat.ixx 3 * \brief 4 * \author Thomas Helfer 5 * \date 17 nov. 2017 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#ifndef LIB_TFEL_MATERIAL_BARLAT_IXX 15#define LIB_TFEL_MATERIAL_BARLAT_IXX 16 17#include"TFEL/Material/OrthotropicStressLinearTransformation.hxx" 18 19namespace tfel{ 20 21 namespace material{ 22 23 namespace internals{ 24 25 /*! 26 * \brief add the terms relative to the eigenvectors derivatives 27 * \param[out] d2Phi_ds2: second derivative of the Barlat equivalent stress 28 * \param[in] dPhi_dvp: first derivative of the Barlat 29 * equivalent stress with respect to the eigenvalues 30 * \param[in] d2Phi_dvp2: second derivative of the Barlat 31 * equivalent stress with respect to the eigenvalues 32 * \param[in] vp: eigen values 33 * \param[in] m: matrix for the eigen vectors 34 * \param[in] e: criterion used to check if two eigenvalues are equal 35 */ 36 template<typename StressStensor> 37 typename std::enable_if<tfel::math::StensorTraits<StressStensor>::dime==1u, 38 void>::type 39 completeBaralatStressSecondDerivative(tfel::material::BarlatStressSecondDerivativeType<StressStensor>&, 40 const tfel::math::tvector<3u,tfel::material::BarlatBaseType<StressStensor>>&, 41 const tfel::math::tvector<6u,tfel::material::BarlatInvertStressType<StressStensor>>&, 42 const tfel::math::tvector<3u,BarlatStressType<StressStensor>>&, 43 const tfel::math::tmatrix<3u,3u,BarlatBaseType<StressStensor>>&, 44 const tfel::material::BarlatStressType<StressStensor>) 45 {} // end of completeBaralatStressSecondDerivative 46 /*! 47 * \brief add the terms relative to the eigenvectors derivatives 48 * \param[out] d2Phi_ds2: second derivative of the Barlat equivalent stress 49 * \param[in] dPhi_dvp: first derivative of the Barlat 50 * equivalent stress with respect to the eigenvalues 51 * \param[in] d2Phi_dvp2: second derivative of the Barlat 52 * equivalent stress with respect to the eigenvalues 53 * \param[in] vp: eigen values 54 * \param[in] m: matrix for the eigen vectors 55 * \param[in] e: criterion used to check if two eigenvalues are equal 56 */ 57 template<typename StressStensor> 58 typename std::enable_if<tfel::math::StensorTraits<StressStensor>::dime==2u, 59 void>::type 60 completeBaralatStressSecondDerivative(tfel::material::BarlatStressSecondDerivativeType<StressStensor>& d2Phi_ds2, 61 const tfel::math::tvector<3u,tfel::material::BarlatBaseType<StressStensor>>& dPhi_dvp, 62 const tfel::math::tvector<6u,tfel::material::BarlatInvertStressType<StressStensor>>& d2Phi_dvp2, 63 const tfel::math::tvector<3u,BarlatStressType<StressStensor>>& vp, 64 const tfel::math::tmatrix<3u,3u,BarlatBaseType<StressStensor>>& m, 65 const tfel::material::BarlatStressType<StressStensor> e) 66 { 67 using namespace tfel::math; 68 using base = tfel::material::BarlatBaseType<StressStensor>; 69 constexpr const auto icste = Cste<base>::isqrt2; 70 const tvector<3u,base> v0 = m.template column_view<0u>(); 71 const tvector<3u,base> v1 = m.template column_view<1u>(); 72 const stensor<2u,base> n01 = stensor<2u,base>::buildFromVectorsSymmetricDiadicProduct(v0,v1)*icste; 73 if(std::abs(vp[0]-vp[1])<e){ 74 d2Phi_ds2 += ((d2Phi_dvp2[0]+d2Phi_dvp2[1]-2*d2Phi_dvp2[3])/2)*(n01^n01); 75 } else { 76 d2Phi_ds2 += (dPhi_dvp[0]-dPhi_dvp[1])/(vp[0]-vp[1])*(n01^n01); 77 } 78 } // end of completeBaralatStressSecondDerivative 79 /*! 80 * \brief add the terms relative to the eigenvectors derivatives 81 * \param[out] d2Phi_ds2: second derivative of the Barlat equivalent stress 82 * \param[in] dPhi_dvp: first derivative of the Barlat 83 * equivalent stress with respect to the eigenvalues 84 * \param[in] d2Phi_dvp2: second derivative of the Barlat 85 * equivalent stress with respect to the eigenvalues 86 * \param[in] vp: eigen values 87 * \param[in] m: matrix for the eigen vectors 88 * \param[in] e: criterion used to check if two eigenvalues are equal 89 */ 90 template<typename StressStensor> 91 typename std::enable_if<tfel::math::StensorTraits<StressStensor>::dime==3u, 92 void>::type 93 completeBaralatStressSecondDerivative(tfel::material::BarlatStressSecondDerivativeType<StressStensor>& d2Phi_ds2, 94 const tfel::math::tvector<3u,tfel::material::BarlatBaseType<StressStensor>>& dPhi_dvp, 95 const tfel::math::tvector<6u,tfel::material::BarlatInvertStressType<StressStensor>>& d2Phi_dvp2, 96 const tfel::math::tvector<3u,BarlatStressType<StressStensor>>& vp, 97 const tfel::math::tmatrix<3u,3u,BarlatBaseType<StressStensor>>& m, 98 const tfel::material::BarlatStressType<StressStensor> e) 99 { 100 using namespace tfel::math; 101 using base = tfel::material::BarlatBaseType<StressStensor>; 102 constexpr const auto cste = Cste<base>::isqrt2; 103 const tvector<3u,base> v0 = m.template column_view<0u>(); 104 const tvector<3u,base> v1 = m.template column_view<1u>(); 105 const tvector<3u,base> v2 = m.template column_view<2u>(); 106 const stensor<3u,base> n01 = stensor<3u,base>::buildFromVectorsSymmetricDiadicProduct(v0,v1)*cste; 107 const stensor<3u,base> n02 = stensor<3u,base>::buildFromVectorsSymmetricDiadicProduct(v0,v2)*cste; 108 const stensor<3u,base> n12 = stensor<3u,base>::buildFromVectorsSymmetricDiadicProduct(v1,v2)*cste; 109 if(std::abs(vp[0]-vp[1])<e){ 110 d2Phi_ds2 += ((d2Phi_dvp2[0]+d2Phi_dvp2[1]-2*d2Phi_dvp2[3])/2)*(n01^n01); 111 } else { 112 d2Phi_ds2 += (dPhi_dvp[0]-dPhi_dvp[1])/(vp[0]-vp[1])*(n01^n01); 113 } 114 if(std::abs(vp[0]-vp[2])<e){ 115 d2Phi_ds2 += ((d2Phi_dvp2[0]+d2Phi_dvp2[2]-2*d2Phi_dvp2[4])/2)*(n02^n02); 116 } else { 117 d2Phi_ds2 += (dPhi_dvp[0]-dPhi_dvp[2])/(vp[0]-vp[2])*(n02^n02); 118 } 119 if(std::abs(vp[1]-vp[2])<e){ 120 d2Phi_ds2 += ((d2Phi_dvp2[1]+d2Phi_dvp2[2]-2*d2Phi_dvp2[5])/2)*(n12^n12); 121 } else { 122 d2Phi_ds2 += (dPhi_dvp[1]-dPhi_dvp[2])/(vp[1]-vp[2])*(n12^n12); 123 } 124 } // end of completeBaralatStressSecondDerivative 125 126 } // end namespace internals 127 128 template<unsigned short N, typename real> 129 tfel::math::st2tost2<N,real> 130 makeBarlatLinearTransformation(const real c12,const real c21, const real c13, 131 const real c31,const real c23, const real c32, 132 const real c44,const real c55, const real c66) 133 { 134 return makeOrthotropicStressLinearTransformation<N,real>(c12,c21,c13,c31,c23, 135 c32,c44,c55,c66); 136 } // end of makeBarlatLinearTransformation 137 138 template<unsigned short N, typename real> 139 tfel::math::st2tost2<N,real> 140 makeBarlatLinearTransformation(const tfel::math::tvector<9u,real>& c) 141 { 142 return makeOrthotropicStressLinearTransformation<N,real>(c); 143 } // end of makeBarlatLinearTransformation 144 145 template<ModellingHypothesis::Hypothesis H, 146 OrthotropicAxesConvention c,typename real> 147 tfel::math::st2tost2<ModellingHypothesisToSpaceDimension<H>::value,real> 148 makeBarlatLinearTransformation(const real c12,const real c21, const real c13, 149 const real c31,const real c23, const real c32, 150 const real c44,const real c55, const real c66) 151 { 152 return makeOrthotropicStressLinearTransformation<H,c,real>(c12,c21,c13,c31,c23, 153 c32,c44,c55,c66); 154 } // end of makeBarlatLinearTransformation 155 156 template<ModellingHypothesis::Hypothesis H, 157 OrthotropicAxesConvention oac,typename real> 158 tfel::math::st2tost2<ModellingHypothesisToSpaceDimension<H>::value,real> 159 makeBarlatLinearTransformation(const tfel::math::tvector<9u,real>& c) 160 { 161 return makeOrthotropicStressLinearTransformation<H,oac,real>(c); 162 } // end of makeBarlatLinearTransformation 163 164 template<typename StressStensor, 165 typename BarlatExponentType, 166 tfel::math::stensor_common::EigenSolver es> 167 BarlatStressType<StressStensor> 168 computeBarlatStress(const StressStensor& s, 169 const BarlatLinearTransformationType<StressStensor>& l1, 170 const BarlatLinearTransformationType<StressStensor>& l2, 171 const BarlatExponentType a, 172 const BarlatStressType<StressStensor> e) 173 { 174 using real = BarlatBaseType<StressStensor>; 175 const auto seq = sigmaeq(s); 176 if(seq<e){ 177 return seq*0; 178 } 179 const auto iseq = 1/seq; 180 const auto s1 = eval(l1*s); 181 const auto s2 = eval(l2*s); 182 const auto vp1 = s1.template computeEigenValues<es>()*iseq; 183 const auto vp2 = s2.template computeEigenValues<es>()*iseq; 184 return seq*std::pow((std::pow(std::abs(vp1[0]-vp2[0]),a)+ 185 std::pow(std::abs(vp1[0]-vp2[1]),a)+ 186 std::pow(std::abs(vp1[0]-vp2[2]),a)+ 187 std::pow(std::abs(vp1[1]-vp2[0]),a)+ 188 std::pow(std::abs(vp1[1]-vp2[1]),a)+ 189 std::pow(std::abs(vp1[1]-vp2[2]),a)+ 190 std::pow(std::abs(vp1[2]-vp2[0]),a)+ 191 std::pow(std::abs(vp1[2]-vp2[1]),a)+ 192 std::pow(std::abs(vp1[2]-vp2[2]),a))/4,1/real(a)); 193 } // end of computeBarlatStress 194 195 template<typename StressStensor, 196 typename BarlatExponentType, 197 tfel::math::stensor_common::EigenSolver es> 198 std::tuple<BarlatStressType<StressStensor>, 199 BarlatStressNormalType<StressStensor>> 200 computeBarlatStressNormal(const StressStensor& s, 201 const BarlatLinearTransformationType<StressStensor>& l1, 202 const BarlatLinearTransformationType<StressStensor>& l2, 203 const BarlatExponentType a, 204 const BarlatStressType<StressStensor> e) 205 { 206 constexpr const auto N = tfel::math::StensorTraits<StressStensor>::dime; 207 using stress = BarlatStressType<StressStensor>; 208 using real = BarlatBaseType<StressStensor>; 209 using normal = BarlatStressNormalType<StressStensor>; 210 // Von Mises stress used to normalise the stress eigenvalues 211 const auto seq = sigmaeq(s); 212 if(seq<e){ 213 return std::make_tuple(stress{0},normal{real{0}}); 214 } 215 const auto iseq = 1/seq; 216 // linear transformations 217 const auto s1 = eval(l1*s); 218 const auto s2 = eval(l2*s); 219 // no structured bindings yet 220 tfel::math::tvector<3u,stress> vp1; 221 tfel::math::tmatrix<3u,3u,real> m1; 222 std::tie(vp1,m1) = s1.template computeEigenVectors<es>(); 223 const auto n1 = tfel::math::stensor<N,stress>::computeEigenTensors(m1); 224 tfel::math::tvector<3u,stress> vp2; 225 tfel::math::tmatrix<3u,3u,real> m2; 226 std::tie(vp2,m2) = s2.template computeEigenVectors<es>(); 227 const auto n2 = tfel::math::stensor<N,stress>::computeEigenTensors(m2); 228 // Barlat stress 229 const auto rvp1 = vp1*iseq; 230 const auto rvp2 = vp2*iseq; 231 const auto Phi_a = (std::pow(std::abs(rvp1[0]-rvp2[0]),a)+ 232 std::pow(std::abs(rvp1[0]-rvp2[1]),a)+ 233 std::pow(std::abs(rvp1[0]-rvp2[2]),a)+ 234 std::pow(std::abs(rvp1[1]-rvp2[0]),a)+ 235 std::pow(std::abs(rvp1[1]-rvp2[1]),a)+ 236 std::pow(std::abs(rvp1[1]-rvp2[2]),a)+ 237 std::pow(std::abs(rvp1[2]-rvp2[0]),a)+ 238 std::pow(std::abs(rvp1[2]-rvp2[1]),a)+ 239 std::pow(std::abs(rvp1[2]-rvp2[2]),a))/4; 240 // Barlat equivalent stress 241 const auto Phi = seq*std::pow(Phi_a,1/real(a)); 242 // For the derivatives, the stress eigenvalues are normalised by 243 // the Barlat equivalent stress 244 const auto rvp1b = tfel::math::eval(vp1/Phi); 245 const auto rvp2b = tfel::math::eval(vp2/Phi); 246 const tfel::math::tvector<9u,real> drvpb_am2 = 247 {real(std::pow(std::abs(rvp1b[0]-rvp2b[0]),a-2)), 248 real(std::pow(std::abs(rvp1b[0]-rvp2b[1]),a-2)), 249 real(std::pow(std::abs(rvp1b[0]-rvp2b[2]),a-2)), 250 real(std::pow(std::abs(rvp1b[1]-rvp2b[0]),a-2)), 251 real(std::pow(std::abs(rvp1b[1]-rvp2b[1]),a-2)), 252 real(std::pow(std::abs(rvp1b[1]-rvp2b[2]),a-2)), 253 real(std::pow(std::abs(rvp1b[2]-rvp2b[0]),a-2)), 254 real(std::pow(std::abs(rvp1b[2]-rvp2b[1]),a-2)), 255 real(std::pow(std::abs(rvp1b[2]-rvp2b[2]),a-2))}; 256 const tfel::math::tvector<3u,real> dPhi_dsvp1 = 257 {((rvp1b[0]-rvp2b[0])*drvpb_am2[0]+ 258 (rvp1b[0]-rvp2b[1])*drvpb_am2[1]+ 259 (rvp1b[0]-rvp2b[2])*drvpb_am2[2])/4, 260 ((rvp1b[1]-rvp2b[0])*drvpb_am2[3]+ 261 (rvp1b[1]-rvp2b[1])*drvpb_am2[4]+ 262 (rvp1b[1]-rvp2b[2])*drvpb_am2[5])/4, 263 ((rvp1b[2]-rvp2b[0])*drvpb_am2[6]+ 264 (rvp1b[2]-rvp2b[1])*drvpb_am2[7]+ 265 (rvp1b[2]-rvp2b[2])*drvpb_am2[8])/4}; 266 const tfel::math::tvector<3u,real> dPhi_dsvp2 = 267 {((rvp2b[0]-rvp1b[0])*drvpb_am2[0]+ 268 (rvp2b[0]-rvp1b[1])*drvpb_am2[3]+ 269 (rvp2b[0]-rvp1b[2])*drvpb_am2[6])/4, 270 ((rvp2b[1]-rvp1b[0])*drvpb_am2[1]+ 271 (rvp2b[1]-rvp1b[1])*drvpb_am2[4]+ 272 (rvp2b[1]-rvp1b[2])*drvpb_am2[7])/4, 273 ((rvp2b[2]-rvp1b[0])*drvpb_am2[2]+ 274 (rvp2b[2]-rvp1b[1])*drvpb_am2[5]+ 275 (rvp2b[2]-rvp1b[2])*drvpb_am2[8])/4}; 276 const auto dPhi_ds1 = dPhi_dsvp1[0]*std::get<0>(n1)+ 277 dPhi_dsvp1[1]*std::get<1>(n1)+ 278 dPhi_dsvp1[2]*std::get<2>(n1); 279 const auto dPhi_ds2 = dPhi_dsvp2[0]*std::get<0>(n2)+ 280 dPhi_dsvp2[1]*std::get<1>(n2)+ 281 dPhi_dsvp2[2]*std::get<2>(n2); 282 return std::make_tuple(Phi,eval(dPhi_ds1*l1+dPhi_ds2*l2)); 283 } // end of computeBarlatStressNormal 284 285 template<typename StressStensor,typename BarlatExponentType> 286 BarlatStressAndDerivativesWithRespectToEigenvalues<StressStensor> 287 computeBarlatStressSecondDerivative(const tfel::math::tvector<3u,BarlatStressType<StressStensor>>& vp1, 288 const tfel::math::tvector<3u,BarlatStressType<StressStensor>>& vp2, 289 const BarlatStressType<StressStensor> seq, 290 const BarlatExponentType a) 291 { 292 using real = BarlatBaseType<StressStensor>; 293 BarlatStressAndDerivativesWithRespectToEigenvalues<StressStensor> d; 294 // Barlat equivalent stress 295 const auto iseq = 1/seq; 296 const auto rvp1 = vp1*iseq; 297 const auto rvp2 = vp2*iseq; 298 const auto Phi_a = (std::pow(std::abs(rvp1[0]-rvp2[0]),a)+ 299 std::pow(std::abs(rvp1[0]-rvp2[1]),a)+ 300 std::pow(std::abs(rvp1[0]-rvp2[2]),a)+ 301 std::pow(std::abs(rvp1[1]-rvp2[0]),a)+ 302 std::pow(std::abs(rvp1[1]-rvp2[1]),a)+ 303 std::pow(std::abs(rvp1[1]-rvp2[2]),a)+ 304 std::pow(std::abs(rvp1[2]-rvp2[0]),a)+ 305 std::pow(std::abs(rvp1[2]-rvp2[1]),a)+ 306 std::pow(std::abs(rvp1[2]-rvp2[2]),a))/4; 307 d.Phi = seq*std::pow(Phi_a,1/real(a)); 308 // For the derivatives, the stress eigenvalues are normalised by 309 // the Barlat equivalent stress 310 const auto rvp1b = tfel::math::eval(vp1/d.Phi); 311 const auto rvp2b = tfel::math::eval(vp2/d.Phi); 312 const tfel::math::tvector<9u,real> drvpb_am2 = 313 {real(std::pow(std::abs(rvp1b[0]-rvp2b[0]),a-2)), // s'1-s''1 314 real(std::pow(std::abs(rvp1b[0]-rvp2b[1]),a-2)), // s'1-s''2 315 real(std::pow(std::abs(rvp1b[0]-rvp2b[2]),a-2)), // s'1-s''3 316 real(std::pow(std::abs(rvp1b[1]-rvp2b[0]),a-2)), // s'2-s''1 317 real(std::pow(std::abs(rvp1b[1]-rvp2b[1]),a-2)), // s'2-s''2 318 real(std::pow(std::abs(rvp1b[1]-rvp2b[2]),a-2)), // s'2-s''3 319 real(std::pow(std::abs(rvp1b[2]-rvp2b[0]),a-2)), // s'3-s''1 320 real(std::pow(std::abs(rvp1b[2]-rvp2b[1]),a-2)), // s'3-s''2 321 real(std::pow(std::abs(rvp1b[2]-rvp2b[2]),a-2))}; // s'3-s''3 322 d.dPhi_dsvp1 = {((rvp1b[0]-rvp2b[0])*drvpb_am2[0]+ 323 (rvp1b[0]-rvp2b[1])*drvpb_am2[1]+ 324 (rvp1b[0]-rvp2b[2])*drvpb_am2[2])/4, 325 ((rvp1b[1]-rvp2b[0])*drvpb_am2[3]+ 326 (rvp1b[1]-rvp2b[1])*drvpb_am2[4]+ 327 (rvp1b[1]-rvp2b[2])*drvpb_am2[5])/4, 328 ((rvp1b[2]-rvp2b[0])*drvpb_am2[6]+ 329 (rvp1b[2]-rvp2b[1])*drvpb_am2[7]+ 330 (rvp1b[2]-rvp2b[2])*drvpb_am2[8])/4}; 331 d.dPhi_dsvp2 = {((rvp2b[0]-rvp1b[0])*drvpb_am2[0]+ 332 (rvp2b[0]-rvp1b[1])*drvpb_am2[3]+ 333 (rvp2b[0]-rvp1b[2])*drvpb_am2[6])/4, 334 ((rvp2b[1]-rvp1b[0])*drvpb_am2[1]+ 335 (rvp2b[1]-rvp1b[1])*drvpb_am2[4]+ 336 (rvp2b[1]-rvp1b[2])*drvpb_am2[7])/4, 337 ((rvp2b[2]-rvp1b[0])*drvpb_am2[2]+ 338 (rvp2b[2]-rvp1b[1])*drvpb_am2[5]+ 339 (rvp2b[2]-rvp1b[2])*drvpb_am2[8])/4}; 340 // second derivative 341 const auto cste = (a-1)/d.Phi; 342 d.d2Phi_dsvp12 = {cste*((drvpb_am2[0]+drvpb_am2[1]+drvpb_am2[2])/4- 343 d.dPhi_dsvp1[0]*d.dPhi_dsvp1[0]), // s1s1 344 cste*((drvpb_am2[3]+drvpb_am2[4]+drvpb_am2[5])/4- 345 d.dPhi_dsvp1[1]*d.dPhi_dsvp1[1]), // s2s2 346 cste*((drvpb_am2[6]+drvpb_am2[7]+drvpb_am2[8])/4- 347 d.dPhi_dsvp1[2]*d.dPhi_dsvp1[2]), // s3s3 348 -cste*d.dPhi_dsvp1[0]*d.dPhi_dsvp1[1], // s1s2 349 -cste*d.dPhi_dsvp1[0]*d.dPhi_dsvp1[2], // s1s3 350 -cste*d.dPhi_dsvp1[1]*d.dPhi_dsvp1[2]}; // s2s3 351 d.d2Phi_dsvp22 = {cste*((drvpb_am2[0]+drvpb_am2[3]+drvpb_am2[6])/4- 352 d.dPhi_dsvp2[0]*d.dPhi_dsvp2[0]), // s1s1 353 cste*((drvpb_am2[1]+drvpb_am2[4]+drvpb_am2[7])/4- 354 d.dPhi_dsvp2[1]*d.dPhi_dsvp2[1]), // s2s2 355 cste*((drvpb_am2[2]+drvpb_am2[5]+drvpb_am2[8])/4- 356 d.dPhi_dsvp2[2]*d.dPhi_dsvp2[2]), // s3s3 357 -cste*d.dPhi_dsvp2[0]*d.dPhi_dsvp2[1], // s1s2 358 -cste*d.dPhi_dsvp2[0]*d.dPhi_dsvp2[2], // s1s3 359 -cste*d.dPhi_dsvp2[1]*d.dPhi_dsvp2[2]}; // s2s3 360 d.d2Phi_dsvp1dsvp2 = {-cste*(drvpb_am2[0]/4+d.dPhi_dsvp1[0]*d.dPhi_dsvp2[0]), //s'1 s''1 361 -cste*(drvpb_am2[1]/4+d.dPhi_dsvp1[0]*d.dPhi_dsvp2[1]), //s'1 s''2 362 -cste*(drvpb_am2[2]/4+d.dPhi_dsvp1[0]*d.dPhi_dsvp2[2]), //s'1 s''3 363 -cste*(drvpb_am2[3]/4+d.dPhi_dsvp1[1]*d.dPhi_dsvp2[0]), //s'2 s''1 364 -cste*(drvpb_am2[4]/4+d.dPhi_dsvp1[1]*d.dPhi_dsvp2[1]), //s'2 s''2 365 -cste*(drvpb_am2[5]/4+d.dPhi_dsvp1[1]*d.dPhi_dsvp2[2]), //s'2 s''3 366 -cste*(drvpb_am2[6]/4+d.dPhi_dsvp1[2]*d.dPhi_dsvp2[0]), //s'3 s''1 367 -cste*(drvpb_am2[7]/4+d.dPhi_dsvp1[2]*d.dPhi_dsvp2[1]), //s'3 s''2 368 -cste*(drvpb_am2[8]/4+d.dPhi_dsvp1[2]*d.dPhi_dsvp2[2])}; //s'3 s''3 369 return d; 370 } // end of BarlatStressSecondDerivative 371 372 template<typename StressStensor, 373 typename BarlatExponentType, 374 tfel::math::stensor_common::EigenSolver es> 375 std::tuple<BarlatStressType<StressStensor>, 376 BarlatStressNormalType<StressStensor>, 377 BarlatStressSecondDerivativeType<StressStensor>> 378 computeBarlatStressSecondDerivative(const StressStensor& s, 379 const BarlatLinearTransformationType<StressStensor>& l1, 380 const BarlatLinearTransformationType<StressStensor>& l2, 381 const BarlatExponentType a, 382 const BarlatStressType<StressStensor> e) 383 { 384 constexpr const auto N = tfel::math::StensorTraits<StressStensor>::dime; 385 using stress = BarlatStressType<StressStensor>; 386 using real = BarlatBaseType<StressStensor>; 387 using normal = BarlatStressNormalType<StressStensor>; 388 using istress = BarlatInvertStressType<StressStensor>; 389 using second_derivative = BarlatStressSecondDerivativeType<StressStensor>; 390 // Von Mises stress used to normalise the stress eigenvalues 391 const auto seq = sigmaeq(s); 392 if(seq<e){ 393 return std::make_tuple(stress{0},normal{real{0}}, 394 second_derivative{istress{0}}); 395 } 396 // linear transformations 397 const auto tl1 = transpose(l1); 398 const auto tl2 = transpose(l2); 399 const auto s1 = eval(l1*s); 400 const auto s2 = eval(l2*s); 401 // no structured bindings yet 402 tfel::math::tvector<3u,stress> vp1; 403 tfel::math::tmatrix<3u,3u,real> m1; 404 std::tie(vp1,m1) = s1.template computeEigenVectors<es>(); 405 const auto n1 = tfel::math::stensor<N,stress>::computeEigenTensors(m1); 406 tfel::math::tvector<3u,stress> vp2; 407 tfel::math::tmatrix<3u,3u,real> m2; 408 std::tie(vp2,m2) = s2.template computeEigenVectors<es>(); 409 const auto n2 = tfel::math::stensor<N,stress>::computeEigenTensors(m2); 410 // derivatives of the Barlat potentiel with respect of the eigen 411 // values of the s1 and s2 412 const auto d = computeBarlatStressSecondDerivative<StressStensor>(vp1,vp2,seq,a); 413 // Barlat normal 414 const auto dPhi_ds1 = (d.dPhi_dsvp1[0]*std::get<0>(n1)+ 415 d.dPhi_dsvp1[1]*std::get<1>(n1)+ 416 d.dPhi_dsvp1[2]*std::get<2>(n1)); 417 const auto dPhi_ds2 = (d.dPhi_dsvp2[0]*std::get<0>(n2)+ 418 d.dPhi_dsvp2[1]*std::get<1>(n2)+ 419 d.dPhi_dsvp2[2]*std::get<2>(n2)); 420 const auto dPhi_ds = eval(dPhi_ds1*l1+dPhi_ds2*l2); 421 // second derivatives with respect to s' 422 auto d2Phi_ds1ds1 = 423 eval(d.d2Phi_dsvp12[0]*(std::get<0>(n1)^std::get<0>(n1))+ // s'1s'1 424 d.d2Phi_dsvp12[3]*(std::get<0>(n1)^std::get<1>(n1))+ // s'1s'2 425 d.d2Phi_dsvp12[4]*(std::get<0>(n1)^std::get<2>(n1))+ // s'1s'3 426 d.d2Phi_dsvp12[3]*(std::get<1>(n1)^std::get<0>(n1))+ // s'2s'1 427 d.d2Phi_dsvp12[1]*(std::get<1>(n1)^std::get<1>(n1))+ // s'2s'2 428 d.d2Phi_dsvp12[5]*(std::get<1>(n1)^std::get<2>(n1))+ // s'2s'3 429 d.d2Phi_dsvp12[4]*(std::get<2>(n1)^std::get<0>(n1))+ // s'3s'1 430 d.d2Phi_dsvp12[5]*(std::get<2>(n1)^std::get<1>(n1))+ // s'3s'2 431 d.d2Phi_dsvp12[2]*(std::get<2>(n1)^std::get<2>(n1))); // s'3s'3 432 // a now, terms related to the eigen vectors derivatives 433 internals::completeBaralatStressSecondDerivative<StressStensor>(d2Phi_ds1ds1, 434 d.dPhi_dsvp1, 435 d.d2Phi_dsvp12, 436 vp1,m1,e); 437 // second derivatives with respect to s's'' 438 const auto d2Phi_ds1ds2 = 439 eval(d.d2Phi_dsvp1dsvp2[0]*(std::get<0>(n1)^std::get<0>(n2))+ // s'1s''1 440 d.d2Phi_dsvp1dsvp2[1]*(std::get<0>(n1)^std::get<1>(n2))+ // s'1s''2 441 d.d2Phi_dsvp1dsvp2[2]*(std::get<0>(n1)^std::get<2>(n2))+ // s'1s''3 442 d.d2Phi_dsvp1dsvp2[3]*(std::get<1>(n1)^std::get<0>(n2))+ // s'2s''1 443 d.d2Phi_dsvp1dsvp2[4]*(std::get<1>(n1)^std::get<1>(n2))+ // s'2s''2 444 d.d2Phi_dsvp1dsvp2[5]*(std::get<1>(n1)^std::get<2>(n2))+ // s'2s''3 445 d.d2Phi_dsvp1dsvp2[6]*(std::get<2>(n1)^std::get<0>(n2))+ // s'3s''1 446 d.d2Phi_dsvp1dsvp2[7]*(std::get<2>(n1)^std::get<1>(n2))+ // s'3s''2 447 d.d2Phi_dsvp1dsvp2[8]*(std::get<2>(n1)^std::get<2>(n2))); // s'3s''3 448 // second derivatives with respect to s'' 449 auto d2Phi_ds2ds2 = 450 eval(d.d2Phi_dsvp22[0]*(std::get<0>(n2)^std::get<0>(n2))+ // s''1s''1 451 d.d2Phi_dsvp22[3]*(std::get<0>(n2)^std::get<1>(n2))+ // s''1s''2 452 d.d2Phi_dsvp22[4]*(std::get<0>(n2)^std::get<2>(n2))+ // s''1s''3 453 d.d2Phi_dsvp22[3]*(std::get<1>(n2)^std::get<0>(n2))+ // s''2s''1 454 d.d2Phi_dsvp22[1]*(std::get<1>(n2)^std::get<1>(n2))+ // s''2s''2 455 d.d2Phi_dsvp22[5]*(std::get<1>(n2)^std::get<2>(n2))+ // s''2s''3 456 d.d2Phi_dsvp22[4]*(std::get<2>(n2)^std::get<0>(n2))+ // s''3s''1 457 d.d2Phi_dsvp22[5]*(std::get<2>(n2)^std::get<1>(n2))+ // s''3s''2 458 d.d2Phi_dsvp22[2]*(std::get<2>(n2)^std::get<2>(n2))); // s''3s''3 459 // a now, terms related to the eigen vectors derivatives 460 internals::completeBaralatStressSecondDerivative<StressStensor>(d2Phi_ds2ds2, 461 d.dPhi_dsvp2, 462 d.d2Phi_dsvp22, 463 vp2,m2,e); 464 // second derivatives with respect to s''s' 465 const auto d2Phi_ds2ds1 = 466 eval(d.d2Phi_dsvp1dsvp2[0]*(std::get<0>(n2)^std::get<0>(n1))+ // s''1s'1 467 d.d2Phi_dsvp1dsvp2[3]*(std::get<0>(n2)^std::get<1>(n1))+ // s''1s'2 468 d.d2Phi_dsvp1dsvp2[6]*(std::get<0>(n2)^std::get<2>(n1))+ // s''1s'3 469 d.d2Phi_dsvp1dsvp2[1]*(std::get<1>(n2)^std::get<0>(n1))+ // s''2s'1 470 d.d2Phi_dsvp1dsvp2[4]*(std::get<1>(n2)^std::get<1>(n1))+ // s''2s'2 471 d.d2Phi_dsvp1dsvp2[7]*(std::get<1>(n2)^std::get<2>(n1))+ // s''2s'3 472 d.d2Phi_dsvp1dsvp2[2]*(std::get<2>(n2)^std::get<0>(n1))+ // s''3s'1 473 d.d2Phi_dsvp1dsvp2[5]*(std::get<2>(n2)^std::get<1>(n1))+ // s''3s'2 474 d.d2Phi_dsvp1dsvp2[8]*(std::get<2>(n2)^std::get<2>(n1))); // s''3s'3 475 const auto d2Phi_ds2 = eval(tl1*(d2Phi_ds1ds1*l1+d2Phi_ds1ds2*l2)+ 476 tl2*(d2Phi_ds2ds1*l1+d2Phi_ds2ds2*l2)); 477 return std::make_tuple(d.Phi,dPhi_ds,d2Phi_ds2); 478 } // end of computeBarlatSecondDerivative 479 480 } // end of namespace material 481 482} // end of namespace tfel 483 484#endif /* LIB_TFEL_MATERIAL_BARLAT_IXX */ 485