1 /*!
2 * \file tests/Math/t2tost2/SaintVenantKirchhoffTangentOperator.cxx
3 * \brief
4 * \author Thomas Helfer
5 * \brief 08 juin 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 #ifdef NDEBUG
15 #undef NDEBUG
16 #endif /* NDEBUG */
17
18 #include<cmath>
19 #include<string>
20 #include<cstdlib>
21 #include<cassert>
22 #include<iostream>
23
24 #include"TFEL/Math/tensor.hxx"
25 #include"TFEL/Math/stensor.hxx"
26 #include"TFEL/Math/t2tost2.hxx"
27 #include"TFEL/Math/st2tost2.hxx"
28 #include"TFEL/Material/Lame.hxx"
29
30 #include"TFEL/Tests/TestCase.hxx"
31 #include"TFEL/Tests/TestProxy.hxx"
32 #include"TFEL/Tests/TestManager.hxx"
33
34 // fixing a bug on current glibc++ cygwin versions (19/08/2015)
35 #if defined __CYGWIN__ && (!defined _GLIBCXX_USE_C99)
36 #include<sstream>
37 namespace std{
38 template<typename T>
to_string(const T & v)39 std::string to_string(const T& v){
40 std::ostringstream s;
41 s << v;
42 return s.str();
43 }
44 }
45 #endif /* defined __CYGWIN__ && (!defined _GLIBCXX_USE_C99) */
46
47 template<unsigned short N>
48 struct SaintVenantKirchhoffTangentOperator final
49 : public tfel::tests::TestCase
50 {
SaintVenantKirchhoffTangentOperatorSaintVenantKirchhoffTangentOperator51 SaintVenantKirchhoffTangentOperator()
52 : tfel::tests::TestCase("TFEL/Math",
53 "SaintVenantKirchhoffTangentOperator"+
54 std::to_string(static_cast<unsigned int>(N)))
55 {} // end of SaintVenantKirchhoffTangentOperator
executeSaintVenantKirchhoffTangentOperator56 tfel::tests::TestResult execute() override
57 {
58 using namespace std;
59 using namespace tfel::math;
60 using namespace tfel::material;
61 const double eps = 1.e-7;
62 const double E = 150e9;
63 const double prec1 = E*1.e-7;
64 const double prec2 = 1e-8;
65 const double nu = 0.33;
66 const double lambda = computeLambda(E,nu);
67 const double mu = computeMu(E,nu);
68 const st2tost2<N,double> D = lambda*st2tost2<N,double>::IxI()+2*mu*st2tost2<N,double>::Id();
69 tensor<N> F;
70 for(typename tensor<N>::size_type i=0;i!=F.size();++i){
71 if(i<3){
72 F[i] = 1+0.1*i*i;
73 } else {
74 F[i] = 0.03*(i-2);
75 }
76 }
77 const stensor<N> C = computeRightCauchyGreenTensor(F);
78 const stensor<N> e_gl = 0.5*(C-stensor<N>::Id());
79 const stensor<N> S = D*e_gl;
80 stensor<N> sig = convertSecondPiolaKirchhoffStressToCauchyStress(S,F);
81 const t2tost2<N> dC = t2tost2<N>::dCdF(F);
82 const t2tost2<N> dS = 0.5*D*dC;
83 t2tost2<N> dtau;
84 computePushForwardDerivative(dtau,dS,S,F);
85 tensor<N> dJ;
86 computeDeterminantDerivative(dJ,F);
87 t2tost2<N> dsig;
88 computeCauchyStressDerivativeFromKirchhoffStressDerivative(dsig,dtau,sig,F);
89 t2tost2<N> ndtau;
90 t2tost2<N> ndsig;
91 t2tost2<N> ndS;
92 tensor<N> ndJ;
93 for(typename tensor<N>::size_type i=0;i!=F.size();++i){
94 tensor<N> Fe(F);
95 Fe[i] += eps;
96 const stensor<N> Cef = computeRightCauchyGreenTensor(Fe);
97 const stensor<N> e_glef = 0.5*(Cef-stensor<N>::Id());
98 const stensor<N> Sef = D*e_glef;
99 const double Jef = det(Fe);
100 const stensor<N> tauef = pushForward(Sef,Fe);
101 const stensor<N> sigef = convertSecondPiolaKirchhoffStressToCauchyStress(Sef,Fe);
102 Fe[i] -= 2*eps;
103 const stensor<N> Ceb = computeRightCauchyGreenTensor(Fe);
104 const stensor<N> e_gleb = 0.5*(Ceb-stensor<N>::Id());
105 const stensor<N> Seb = D*e_gleb;
106 const double Jeb = det(Fe);
107 const stensor<N> taueb = pushForward(Seb,Fe);
108 const stensor<N> sigeb = convertSecondPiolaKirchhoffStressToCauchyStress(Seb,Fe);
109 for(typename tensor<N>::size_type j=0;j!=sig.size();++j){
110 ndS(j,i) = (Sef(j)-Seb(j))/(2*eps);
111 ndtau(j,i) = (tauef(j)-taueb(j))/(2*eps);
112 ndsig(j,i) = (sigef(j)-sigeb(j))/(2*eps);
113 }
114 ndJ(i) = (Jef-Jeb)/(2*eps);
115 }
116 for(typename tensor<N>::size_type i=0;i!=sig.size();++i){
117 for(typename tensor<N>::size_type j=0;j!=F.size();++j){
118 TFEL_TESTS_ASSERT(abs(ndS(i,j)-dS(i,j))<prec1);
119 if(abs(ndS(i,j)-dS(i,j))>prec1){
120 cout << "dS error (" << N << "): " << i << " " << j << " " << ndS(i,j) << " " << dS(i,j) << " " << abs(ndS(i,j)-dS(i,j))/(max(abs(abs(ndS(i,j))),abs(dS(i,j)))) << endl;
121 }
122 TFEL_TESTS_ASSERT(abs(ndtau(i,j)-dtau(i,j))<prec1);
123 if(abs(ndtau(i,j)-dtau(i,j))>prec1){
124 cout << "dtau error (" << N << "): " << i << " " << j << " " << ndtau(i,j) << " " << dtau(i,j) << " " << abs(ndtau(i,j)-dtau(i,j))/(max(abs(abs(ndtau(i,j))),abs(dtau(i,j)))) << endl;
125 }
126 TFEL_TESTS_ASSERT(abs(ndsig(i,j)-dsig(i,j))<prec1);
127 if(abs(ndsig(i,j)-dsig(i,j))>prec1){
128 cout << "dsig error (" << N << "): " << i << " " << j << " " << ndsig(i,j) << " " << dsig(i,j) << " " << abs(ndsig(i,j)-dsig(i,j))/(max(abs(abs(ndsig(i,j))),abs(dsig(i,j)))) << endl;
129 }
130 }
131 }
132 for(typename tensor<N>::size_type i=0;i!=F.size();++i){
133 TFEL_TESTS_ASSERT(abs(ndJ(i)-dJ(i))<prec2);
134 if(abs(ndJ(i)-dJ(i))>prec2){
135 cout << "dJ error (" << N << "): " << i << " " << ndJ(i) << " " << dJ(i) << " " << abs(ndJ(i)-dJ(i)) << endl;
136 }
137 }
138 return this->result;
139 } // end of execute
140 };
141
142 using SaintVenantKirchhoffTangentOperator_1D = SaintVenantKirchhoffTangentOperator<1U>;
143 using SaintVenantKirchhoffTangentOperator_2D = SaintVenantKirchhoffTangentOperator<2U>;
144 using SaintVenantKirchhoffTangentOperator_3D = SaintVenantKirchhoffTangentOperator<3U>;
145 TFEL_TESTS_GENERATE_PROXY(SaintVenantKirchhoffTangentOperator_1D,"SaintVenantKirchhoffTangentOperator-1D");
146 TFEL_TESTS_GENERATE_PROXY(SaintVenantKirchhoffTangentOperator_2D,"SaintVenantKirchhoffTangentOperator-2D");
147 TFEL_TESTS_GENERATE_PROXY(SaintVenantKirchhoffTangentOperator_3D,"SaintVenantKirchhoffTangentOperator-3D");
148
149 /* coverity [UNCAUGHT_EXCEPT]*/
main()150 int main(){
151 using namespace tfel::tests;
152 auto& manager = TestManager::getTestManager();
153 manager.addTestOutput(std::cout);
154 manager.addXMLTestOutput("SaintVenantKirchhoffTangentOperator.xml");
155 TestResult r = manager.execute();
156 if(!r.success()){
157 return EXIT_FAILURE;
158 }
159 return EXIT_SUCCESS;
160 }
161
162