1 /*!
2 * \file stensor_isotropic_function2.cxx
3 * \brief
4 * \author Thomas Helfer
5 * \date 25 juin 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 /* NDEBUG */
17
18 #include <cmath>
19 #include <fstream>
20 #include <cstdlib>
21 #include <iostream>
22
23 #include "TFEL/Tests/TestCase.hxx"
24 #include "TFEL/Tests/TestProxy.hxx"
25 #include "TFEL/Tests/TestManager.hxx"
26
27 #include "TFEL/Math/stensor.hxx"
28 #include "TFEL/Math/st2tost2.hxx"
29
30 struct StensorIsotropicFunctionDerivativeTest2 final
31 : public tfel::tests::TestCase {
StensorIsotropicFunctionDerivativeTest2StensorIsotropicFunctionDerivativeTest232 StensorIsotropicFunctionDerivativeTest2()
33 : tfel::tests::TestCase("TFEL/Math",
34 "StensorIsotropicFunctionDerivativeTest2") {
35 } // end of StensorIsotropicFunctionDerivativeTest
36
executeStensorIsotropicFunctionDerivativeTest237 tfel::tests::TestResult execute() override {
38 this->test<3>({-0.000733594, 0.000273806, 0.000273806, 0., 0., 0.});
39 this->test<3>({-0.002, 0.000746479, 0.000746479, 0., 0., 3.74031e-28});
40 return this->result;
41 } // end of execute
42 private:
43 template <unsigned short N>
testStensorIsotropicFunctionDerivativeTest244 void test(const tfel::math::stensor<N, double>& v) {
45 using size_type = typename tfel::math::stensor<N, double>::size_type;
46 constexpr const auto jes = tfel::math::stensor<N, double>::FSESJACOBIEIGENSOLVER;
47 constexpr const double eps = 1.e-8;
48 #if (defined __INTEL_COMPILER)
49 const auto K = tfel::math::st2tost2<N, double>::K();
50 #else
51 TFEL_CONSTEXPR const auto K = tfel::math::st2tost2<N, double>::K();
52 #endif
53 constexpr const auto E = 150e9;
54 constexpr const auto n = 0.3;
55 constexpr const auto m = E / (2 * (1 + n));
56 constexpr const auto rm = 0.6;
57 constexpr const double prec = 4.e-9*E;
58 auto f = [](const double x) -> double { return x > 0 ? x : 0; };
59 auto df = [](const double x) -> double {
60 return std::abs(x) < 1.e-12 ? 0.5 : ((x < 0) ? 0 : 1);
61 };
62 const auto s = deviator(v);
63 const auto fdf =
64 s.template computeIsotropicFunctionAndDerivative<jes>(f, df, 1.e-13);
65 const auto d = 2 * m * (K + (rm - 1) * fdf.second * K);
66 decltype(fdf.second) nd;
67 for (size_type i = 0; i != v.size(); ++i) {
68 auto v2 = v;
69 v2[i] += eps;
70 const auto s2_p = deviator(v2);
71 const auto nf_f =
72 2 * m *
73 (s2_p + (rm - 1) * s2_p.template computeIsotropicFunction<jes>(f));
74 v2[i] -= 2 * eps;
75 const auto s2_m = deviator(v2);
76 const auto nf_b =
77 2 * m *
78 (s2_m + (rm - 1) * s2_m.template computeIsotropicFunction<jes>(f));
79 for (size_type j = 0; j != v.size(); ++j) {
80 nd(j, i) = (nf_f[j] - nf_b[j]) / (2 * eps);
81 }
82 }
83 for (size_type i = 0; i != v.size(); ++i) {
84 for (size_type j = 0; j != v.size(); ++j) {
85 if (std::abs(d(i, j) - nd(i, j)) > prec) {
86 std::cout << i << " " << j << " " << d(i, j) << " " << nd(i, j) << " "
87 << std::abs(d(i, j) - nd(i, j)) << std::endl;
88 }
89 TFEL_TESTS_ASSERT(std::abs(d(i, j) - nd(i, j)) < prec);
90 }
91 }
92 }
93 };
94
95 TFEL_TESTS_GENERATE_PROXY(StensorIsotropicFunctionDerivativeTest2,
96 "StensorIsotropicFunctionDerivativeTest2");
97
98 /* coverity [UNCAUGHT_EXCEPT]*/
main()99 int main() {
100 auto& m = tfel::tests::TestManager::getTestManager();
101 m.addTestOutput(std::cout);
102 m.addXMLTestOutput("StensorIsotropicFunctionDerivative2.xml");
103 const auto r = m.execute();
104 return r.success() ? EXIT_SUCCESS : EXIT_FAILURE;
105 } // end of main
106