1 /*!
2 * \file IsotropicPlasticityTest.cxx
3 * \brief
4 * \author Thomas Helfer
5 * \date 14/01/2018
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
17
18 #include<iostream>
19 #include<cstdlib>
20 #include<cmath>
21
22 #include<cassert>
23
24 #include"TFEL/Tests/TestCase.hxx"
25 #include"TFEL/Tests/TestProxy.hxx"
26 #include"TFEL/Tests/TestManager.hxx"
27
28 #include"TFEL/Material/IsotropicPlasticity.hxx"
29
30 template<unsigned short N,typename real>
31 static tfel::math::stensor<N,real>
getNumericalApproximation(const tfel::math::stensor<N,real> & sig,const real e)32 getNumericalApproximation(const tfel::math::stensor<N,real>& sig,
33 const real e){
34 tfel::math::stensor<N,real> r;
35 for(unsigned short j=0;j!=sig.size();++j){
36 auto sigp = sig;
37 auto sigm = sig;
38 sigp(j)+=e;
39 sigm(j)-=e;
40 const auto Jp = det(deviator(sigp));
41 const auto Jm = det(deviator(sigm));
42 const auto dJ = (Jp-Jm)/(2*e);
43 r(j)=dJ;
44 }
45 return r;
46 } // end of getNumericalApproximation
47
48 struct IsotropicPlasticityTest final
49 : public tfel::tests::TestCase
50 {
IsotropicPlasticityTestIsotropicPlasticityTest51 IsotropicPlasticityTest()
52 : tfel::tests::TestCase("TFEL/Material",
53 "IsotropicPlasticityTest")
54 {} // end of IsotropicPlasticityTest
executeIsotropicPlasticityTest55 tfel::tests::TestResult execute() override
56 {
57 const double sqrt2 = std::sqrt(2.);
58 const double v1[6] = {8.2e-2,-4.5e-2,7.2e-2,2.3e-2*sqrt2,0,0.};
59 const double v2[6] = {-8.2e-2,4.5e-2,7.2e-2,0e-2,2.3e-2*sqrt2,0.};
60 const double v3[6] = {8.2e-2,4.5e-2,-7.2e-2,0e-2,0.e-2,2.3e-2*sqrt2};
61 const double v4[6] = {8.2e-2,4.5e-2,-7.2e-2,3.14e-2*sqrt2,-6.42e-2*sqrt2,2.3e-2*sqrt2};
62 for(const auto v : {v1,v2,v3,v4}){
63 this->check(tfel::math::stensor<1u,double>(v));
64 this->check(tfel::math::stensor<2u,double>(v));
65 this->check(tfel::math::stensor<3u,double>(v));
66 }
67 return this->result;
68 } // end of execute
69 private:
70 template<unsigned short N>
checkIsotropicPlasticityTest71 void check(const tfel::math::stensor<N,double>& s){
72 this->test1(s);
73 this->test2(s);
74 this->test3(s);
75 }
76 template<unsigned short N>
test1IsotropicPlasticityTest77 void test1(const tfel::math::stensor<N,double>& s){
78 TFEL_CONSTEXPR const auto eps = 1e-6;
79 TFEL_CONSTEXPR const auto prec = 2e-13;
80 const auto ndJ = getNumericalApproximation(s,eps);
81 const auto dJ = tfel::material::computeJ3Derivative(s);
82 for(unsigned short i=0;i!=s.size();++i){
83 if(std::abs(dJ(i)-ndJ(i))>prec){
84 std::cout << "dJ : " << i << " " << dJ(i) << " " << ndJ(i)
85 << " " << std::abs(dJ(i)-ndJ(i)) << std::endl;
86 }
87 TFEL_TESTS_ASSERT(std::abs(dJ(i)-ndJ(i))<prec);
88 }
89 }
90 template<unsigned short N>
test2IsotropicPlasticityTest91 void test2(const tfel::math::stensor<N,double>& s){
92 using namespace tfel::math;
93 using namespace tfel::material;
94 using size_type = typename stensor<N,double>::size_type;
95 const double eps = 1.e-5;
96 const double prec = 1.e-10;
97 const auto dJ3 = computeJ3SecondDerivative(s);
98 st2tost2<N,double> ndJ3;
99 for(size_type i=0;i!=s.size();++i){
100 stensor<N> s_e(s);
101 s_e[i] += eps;
102 const stensor<N> dJp = computeJ3Derivative(s_e);
103 s_e[i] -= 2*eps;
104 const stensor<N> dJm = computeJ3Derivative(s_e);
105 for(size_type j=0;j!=s.size();++j){
106 ndJ3(j,i) = (dJp(j)-dJm(j))/(2*eps);
107 }
108 }
109 for(size_type i=0;i!=s.size();++i){
110 for(size_type j=0;j!=s.size();++j){
111 TFEL_TESTS_ASSERT(std::abs(ndJ3(i,j)-dJ3(i,j))<prec);
112 if(std::abs(ndJ3(i,j)-dJ3(i,j))>prec){
113 std::cout << "Error " << N << " (" << i << ", " << j << ") "
114 << "[" << i*StensorDimeToSize<N>::value+j << "]: "
115 << ndJ3(i,j) << " vs " << dJ3(i,j) << " "
116 << std::abs(ndJ3(i,j)-dJ3(i,j)) << std::endl;
117 }
118 }
119 }
120 }
121 template<unsigned short N>
test3IsotropicPlasticityTest122 void test3(const tfel::math::stensor<N,double> s){
123 using namespace tfel::math;
124 using Stensor = stensor<N,double>;
125 using Stensor4 = st2tost2<N,double>;
126 using size_type = typename stensor<N,double>::size_type;
127 const auto id = Stensor::Id();
128 const double prec = 1.e-10;
129 const auto i1 = trace(s);
130 const auto di2 = i1*id-s;
131 const auto d2i2 = (id^id)-Stensor4::Id();
132 const auto d2i3 = computeDeterminantSecondDerivative(s);
133 const auto dJ3 = tfel::material::computeJ3SecondDerivative(s);
134 const st2tost2<N,double> dJ3_2 =
135 (4*i1/9)*(id^id)-((id^di2)+(di2^id)+i1*d2i2)/3+d2i3;
136 for(size_type i=0;i!=s.size();++i){
137 for(size_type j=0;j!=s.size();++j){
138 TFEL_TESTS_ASSERT(std::abs(dJ3_2(i,j)-dJ3(i,j))<prec);
139 if(std::abs(dJ3_2(i,j)-dJ3(i,j))>prec){
140 std::cout << "Error " << N << " (" << i << ", " << j << ") "
141 << "[" << i*StensorDimeToSize<N>::value+j << "]: "
142 << dJ3_2(i,j) << " vs " << dJ3(i,j) << " "
143 << std::abs(dJ3_2(i,j)-dJ3(i,j)) << std::endl;
144 }
145 }
146 }
147 }
148 };
149
150 TFEL_TESTS_GENERATE_PROXY(IsotropicPlasticityTest,
151 "IsotropicPlasticityTest");
152
153 /* coverity [UNCAUGHT_EXCEPT]*/
main()154 int main()
155 {
156 auto& m = tfel::tests::TestManager::getTestManager();
157 m.addTestOutput(std::cout);
158 m.addXMLTestOutput("IsotropicPlasticityTest.xml");
159 return m.execute().success() ? EXIT_SUCCESS : EXIT_FAILURE;
160 }
161