1 /*!
2 * \file OrthotropicPlasticityTest.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/OrthotropicPlasticity.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 tfel::math::tvector<11,real> & c,const real e)32 getNumericalApproximation(const tfel::math::stensor<N,real>& sig,
33 const tfel::math::tvector<11,real>& c,
34 const real e){
35 tfel::math::stensor<N,real> r;
36 for(unsigned short j=0;j!=sig.size();++j){
37 auto sigp = sig;
38 auto sigm = sig;
39 sigp(j)+=e;
40 sigm(j)-=e;
41 const auto Jp = tfel::material::computeJ3O(sigp,c);
42 const auto Jm = tfel::material::computeJ3O(sigm,c);
43 const auto dJ = (Jp-Jm)/(2*e);
44 r(j)=dJ;
45 }
46 return r;
47 } // end of getNumericalApproximation
48
49 struct OrthotropicPlasticityTest final
50 : public tfel::tests::TestCase
51 {
OrthotropicPlasticityTestOrthotropicPlasticityTest52 OrthotropicPlasticityTest()
53 : tfel::tests::TestCase("TFEL/Material",
54 "OrthotropicPlasticityTest")
55 {} // end of OrthotropicPlasticityTest
executeOrthotropicPlasticityTest56 tfel::tests::TestResult execute() override
57 {
58 const double sqrt2 = std::sqrt(2.);
59 const double v1[6] = {8.2e-2,-4.5e-2,7.2e-2,2.3e-2*sqrt2,0,0.};
60 const double v2[6] = {-8.2e-2,4.5e-2,7.2e-2,0e-2,2.3e-2*sqrt2,0.};
61 const double v3[6] = {8.2e-2,4.5e-2,-7.2e-2,0e-2,0.e-2,2.3e-2*sqrt2};
62 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};
63 for(const auto v : {v1,v2,v3,v4}){
64 this->check(tfel::math::stensor<1u,double>(v));
65 this->check(tfel::math::stensor<2u,double>(v));
66 this->check(tfel::math::stensor<3u,double>(v));
67 }
68 return this->result;
69 } // end of execute
70 private:
71 template<unsigned short N>
checkOrthotropicPlasticityTest72 void check(const tfel::math::stensor<N,double>& s){
73 this->test0(s);
74 this->test1(s);
75 this->test2(s);
76 // this->test3(s);
77 }
78 template<unsigned short N>
test0OrthotropicPlasticityTest79 void test0(const tfel::math::stensor<N,double>& s){
80 constexpr const auto prec = 1e-14;
81 constexpr const auto one = double(1);
82 const auto J3 = det(deviator(s));
83 const auto J3O = tfel::material::computeJ3O(s,one,one,one,
84 one,one,one,
85 one,one,one,
86 one,one);
87 TFEL_TESTS_ASSERT(std::abs(J3-J3O)<prec);
88 }
89 template<unsigned short N>
test1OrthotropicPlasticityTest90 void test1(const tfel::math::stensor<N,double>& s){
91 constexpr const auto eps = 1e-6;
92 constexpr const auto prec = 2e-13;
93 constexpr const auto one = double(1);
94 const auto c = tfel::math::tvector<11,double>{one,one,one,
95 one,one,one,
96 one,one,one,
97 one,one};
98 const auto ndJ = getNumericalApproximation(s,c,eps);
99 const auto dJ = tfel::material::computeJ3ODerivative(s,c);
100 for(unsigned short i=0;i!=s.size();++i){
101 if(std::abs(dJ(i)-ndJ(i))>prec){
102 std::cout << "dJ : " << i << " " << dJ(i) << " " << ndJ(i)
103 << " " << std::abs(dJ(i)-ndJ(i)) << std::endl;
104 }
105 TFEL_TESTS_ASSERT(std::abs(dJ(i)-ndJ(i))<prec);
106 }
107 }
108 template<unsigned short N>
test2OrthotropicPlasticityTest109 void test2(const tfel::math::stensor<N,double>& s){
110 using namespace tfel::math;
111 using namespace tfel::material;
112 using size_type = typename stensor<N,double>::size_type;
113 constexpr const auto one = double(1);
114 const auto c = tfel::math::tvector<11,double>{one,one,one,
115 one,one,one,
116 one,one,one,
117 one,one};
118 const double eps = 1.e-5;
119 const double prec = 1.e-10;
120 const auto d2J3O = computeJ3OSecondDerivative(s,c);
121 st2tost2<N,double> nd2J3O;
122 for(size_type i=0;i!=s.size();++i){
123 stensor<N> s_e(s);
124 s_e[i] += eps;
125 const stensor<N> dJ3Op = computeJ3ODerivative(s_e,c);
126 s_e[i] -= 2*eps;
127 const stensor<N> dJ3Om = computeJ3ODerivative(s_e,c);
128 for(size_type j=0;j!=s.size();++j){
129 nd2J3O(j,i) = (dJ3Op(j)-dJ3Om(j))/(2*eps);
130 }
131 }
132 for(size_type i=0;i!=s.size();++i){
133 for(size_type j=0;j!=s.size();++j){
134 TFEL_TESTS_ASSERT(std::abs(nd2J3O(i,j)-d2J3O(i,j))<prec);
135 if(std::abs(nd2J3O(i,j)-d2J3O(i,j))>prec){
136 std::cout << "Error " << N << " (" << i << ", " << j << ") "
137 << "[" << i*StensorDimeToSize<N>::value+j << "]: "
138 << nd2J3O(i,j) << " vs " << d2J3O(i,j) << " "
139 << std::abs(nd2J3O(i,j)-d2J3O(i,j)) << std::endl;
140 }
141 }
142 }
143 }
144 // template<unsigned short N>
145 // void test3(const tfel::math::stensor<N,double> s){
146 // using namespace tfel::math;
147 // using Stensor = stensor<N,double>;
148 // using Stensor4 = st2tost2<N,double>;
149 // using size_type = typename stensor<N,double>::size_type;
150 // const auto id = Stensor::Id();
151 // const double prec = 1.e-10;
152 // const auto i1 = trace(s);
153 // const auto di2 = i1*id-s;
154 // const auto d2i2 = (id^id)-Stensor4::Id();
155 // const auto d2i3 = computeDeterminantSecondDerivative(s);
156 // const auto dJ3 = tfel::material::computeJ3SecondDerivative(s);
157 // const st2tost2<N,double> dJ3_2 =
158 // (4*i1/9)*(id^id)-((id^di2)+(di2^id)+i1*d2i2)/3+d2i3;
159 // for(size_type i=0;i!=s.size();++i){
160 // for(size_type j=0;j!=s.size();++j){
161 // TFEL_TESTS_ASSERT(std::abs(dJ3_2(i,j)-dJ3(i,j))<prec);
162 // if(std::abs(dJ3_2(i,j)-dJ3(i,j))>prec){
163 // std::cout << "Error " << N << " (" << i << ", " << j << ") "
164 // << "[" << i*StensorDimeToSize<N>::value+j << "]: "
165 // << dJ3_2(i,j) << " vs " << dJ3(i,j) << " "
166 // << std::abs(dJ3_2(i,j)-dJ3(i,j)) << std::endl;
167 // }
168 // }
169 // }
170 // }
171 };
172
173 TFEL_TESTS_GENERATE_PROXY(OrthotropicPlasticityTest,
174 "OrthotropicPlasticityTest");
175
176 /* coverity [UNCAUGHT_EXCEPT]*/
main()177 int main()
178 {
179 auto& m = tfel::tests::TestManager::getTestManager();
180 m.addTestOutput(std::cout);
181 m.addXMLTestOutput("OrthotropicPlasticityTest.xml");
182 return m.execute().success() ? EXIT_SUCCESS : EXIT_FAILURE;
183 }
184