v0.14.0
matrix_function.cpp
Go to the documentation of this file.
1 /**
2  * @file matrix_function.cpp
3  * @example matrix_function.cpp
4  * @brief Test and example for matrix function
5  *
6  * For reference see \cite miehe2001algorithms
7  *
8  */
9 
10 //#define FTENSOR_DEBUG
11 #include <FTensor.hpp>
12 
13 #include <MoFEM.hpp>
14 using namespace MoFEM;
15 
16 #include <MatrixFunction.hpp>
17 
18 FTensor::Index<'i', 3> i;
19 FTensor::Index<'j', 3> j;
20 FTensor::Index<'k', 3> k;
21 FTensor::Index<'l', 3> l;
22 
23 template <typename T1, typename T2, int DIM>
24 void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number<DIM> &) {
25  constexpr double eps = 1e-4;
26  for (int ii = 0; ii != DIM; ++ii)
27  for (int jj = 0; jj != DIM; ++jj)
28  for (int kk = 0; kk != DIM; ++kk)
29  for (int ll = 0; ll != DIM; ++ll) {
30 
31  if (std::abs(t_1(ii, jj, kk, ll) - t_2(ii, jj, kk, ll)) > eps)
32  MOFEM_LOG("ATOM_TEST", Sev::error)
33  << "Error " << ii << " " << jj << " " << kk << " " << ll << " "
34  << t_1(ii, jj, kk, ll) << " " << t_2(ii, jj, kk, ll);
35  }
36 
37  for (int ii = 0; ii != DIM; ++ii)
38  for (int jj = 0; jj != DIM; ++jj)
39  for (int kk = 0; kk != DIM; ++kk)
40  for (int ll = 0; ll != DIM; ++ll)
41  t_1(ii, jj, kk, ll) -= t_2(ii, jj, kk, ll);
42 };
43 
44 template <typename T1, typename T2, int DIM>
45 auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number<DIM> &) {
46  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
47  FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_d_a;
48 
49  FTensor::Index<'i', DIM> i;
50  FTensor::Index<'j', DIM> j;
51  FTensor::Index<'k', DIM> k;
52  FTensor::Index<'l', DIM> l;
53 
54  t_d_a(i, j, k, l) = 0;
55 
56  for (int ii = 0; ii != DIM; ++ii)
57  for (int jj = 0; jj != DIM; ++jj)
58  for (int kk = 0; kk != DIM; ++kk)
59  for (int ll = 0; ll != DIM; ++ll)
60  for (int zz = 0; zz != DIM; ++zz) {
61 
62  auto diff = [&](auto ii, auto jj, auto kk, auto ll, int zz) {
63  return
64 
65  t_a(ii, zz) * t_kd(zz, kk) * t_kd(jj, ll)
66 
67  +
68 
69  t_kd(ii, kk) * t_kd(zz, ll) * t_a(zz, jj);
70  };
71 
72  t_d_a(ii, jj, kk, ll) +=
73  (diff(ii, jj, kk, ll, zz) + diff(ii, jj, ll, kk, zz)) / 2.;
74  }
75 
76  diff_ddg(t_d_a, t_d, FTensor::Number<DIM>());
77 
78  return t_d_a;
79 };
80 
81 template <typename T1, typename T2, int DIM>
82 auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number<DIM> &) {
83  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
84  FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_dd_a;
85 
86  FTensor::Index<'i', DIM> i;
87  FTensor::Index<'j', DIM> j;
88  FTensor::Index<'k', DIM> k;
89  FTensor::Index<'l', DIM> l;
90 
91  t_dd_a(i, j, k, l) = 0;
92 
93  for (int ii = 0; ii != DIM; ++ii)
94  for (int jj = 0; jj != DIM; ++jj)
95  for (int kk = 0; kk != DIM; ++kk)
96  for (int ll = 0; ll != DIM; ++ll)
97  for (int mm = 0; mm != DIM; ++mm)
98  for (int nn = 0; nn != DIM; ++nn)
99  for (int zz = 0; zz != DIM; ++zz) {
100 
101  auto diff = [&](auto ii, auto jj, auto kk, auto ll, int mm,
102  int nn, int zz) {
103  return
104 
105  t_s(ii, jj) *
106  (t_kd(ii, mm) * t_kd(zz, nn) * t_kd(zz, kk) * t_kd(jj, ll)
107 
108  +
109 
110  t_kd(ii, kk) * t_kd(zz, ll) * t_kd(zz, mm) *
111  t_kd(jj, nn));
112  };
113 
114  t_dd_a(kk, ll, mm, nn) += (
115 
116  diff(ii, jj, kk, ll, mm, nn, zz)
117 
118  +
119 
120  diff(ii, jj, ll, kk, mm, nn, zz)
121 
122  +
123 
124  diff(ii, jj, kk, ll, nn, mm, zz)
125 
126  +
127 
128  diff(ii, jj, ll, kk, nn, mm, zz)
129 
130  ) /
131  4.;
132  }
133 
134  diff_ddg(t_dd_a, t_dd, FTensor::Number<DIM>());
135 
136  return t_dd_a;
137 };
138 
139 template <typename T1, int DIM>
140 auto get_diff_matrix(T1 &t_d, const FTensor::Number<DIM> &) {
141  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
142  FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_d_a;
143 
144  FTensor::Index<'i', DIM> i;
145  FTensor::Index<'j', DIM> j;
146  FTensor::Index<'k', DIM> k;
147  FTensor::Index<'l', DIM> l;
148 
149  t_d_a(i, j, k, l) = 0;
150 
151  for (int ii = 0; ii != DIM; ++ii)
152  for (int jj = 0; jj != DIM; ++jj)
153  for (int kk = 0; kk != DIM; ++kk)
154  for (int ll = 0; ll != DIM; ++ll) {
155 
156  auto diff = [&](auto ii, auto jj, auto kk, auto ll) {
157  return t_kd(ii, kk) * t_kd(jj, ll);
158  };
159 
160  t_d_a(ii, jj, kk, ll) =
161  (diff(ii, jj, kk, ll) + diff(ii, jj, ll, kk)) / 2.;
162  }
163 
164  diff_ddg(t_d_a, t_d, FTensor::Number<3>());
165 
166  return t_d_a;
167 };
168 
169 template <typename T, int DIM>
170 auto get_norm_t4(T &t, const FTensor::Number<DIM> &) {
171  double r = 0;
172  for (int ii = 0; ii != DIM; ++ii)
173  for (int jj = 0; jj != DIM; ++jj)
174  for (int kk = 0; kk != DIM; ++kk)
175  for (int ll = 0; ll != DIM; ++ll)
176  r += t(ii, jj, kk, ll) * t(ii, jj, kk, ll);
177  return r;
178 };
179 
180 static char help[] = "...\n\n";
181 
182 int main(int argc, char *argv[]) {
183 
184  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
185 
186  auto core_log = logging::core::get();
187  core_log->add_sink(
188  LogManager::createSink(LogManager::getStrmSelf(), "ATOM_TEST"));
189  LogManager::setLog("ATOM_TEST");
190  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
191  MOFEM_LOG_ATTRIBUTES("ATOM_TEST", LogManager::BitLineID);
192 
193  try {
194 
195  auto print_ddg = [](auto &t, auto str = "") {
196  constexpr double eps = 1e-6;
197  for (int ii = 0; ii != 3; ++ii)
198  for (int jj = 0; jj != 3; ++jj)
199  for (int kk = 0; kk != 3; ++kk)
200  for (int ll = 0; ll != 3; ++ll) {
201  double v = t(ii, jj, kk, ll);
202  double w = std::abs(v) < eps ? 0 : v;
203  MOFEM_LOG("ATOM_TEST", Sev::noisy)
204  << str << std::fixed << std::setprecision(3) << std::showpos
205  << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
206  << " : " << w;
207  }
208  };
209 
210  auto print_ddg_direction = [](auto &t, auto kk, int ll) {
211  for (int ii = 0; ii != 3; ++ii)
212  for (int jj = 0; jj <= ii; ++jj)
213  MOFEM_LOG("ATOM_TEST", Sev::noisy)
214  << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
215  << " : " << t(ii, jj, kk, ll);
216  };
217 
218  auto print_mat = [](auto &t) {
219  for (int ii = 0; ii != 3; ++ii)
220  for (int jj = 0; jj != 3; ++jj)
221  MOFEM_LOG("ATOM_TEST", Sev::noisy)
222  << ii + 1 << " " << jj + 1 << " : " << t(ii, jj);
223  };
224 
225  enum swap { swap12, swap01 };
226  auto run_lapack = [](auto &a, swap s = swap12) {
227  int info;
228  double wkopt;
229  double w[3];
230 
231  FTensor::Tensor2<double, 3, 3> t_a{
232 
233  a[0], a[1], a[2],
234 
235  a[3], a[4], a[5],
236 
237  a[6], a[7], a[8]};
238 
239  /* Query and allocate the optimal workspace */
240  int lwork = -1;
241  info = lapack_dsyev('V', 'U', 3, a.data(), 3, w, &wkopt, lwork);
242  if (info > 0)
243  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
244  lwork = (int)wkopt;
245  std::vector<double> work(lwork);
246  /* Solve eigenproblem */
247  info = lapack_dsyev('V', 'U', 3, a.data(), 3, w, &*work.begin(), lwork);
248  if (info > 0)
249  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
250 
251  if (s == swap12) {
252  FTensor::Tensor2<double, 3, 3> t_eig_vec{
253 
254  a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
255 
256  a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
257 
258  a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2]};
259 
260  FTensor::Tensor1<double, 3> t_eig_vals{w[0], w[2], w[1]};
261  return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
262  }
263 
264  FTensor::Tensor2<double, 3, 3> t_eig_vec{
265 
266  a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2],
267 
268  a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
269 
270  a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
271 
272  };
273 
274  FTensor::Tensor1<double, 3> t_eig_vals{w[1], w[0], w[2]};
275  return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
276  };
277 
278  // Test matrix againsst mathematica results
279  {
280  FTensor::Tensor2<double, 3, 3> t_A{
281 
282  1., 0.1, -0.5,
283 
284  0.1, 2., 0.,
285 
286  -0.5, 0., 3.};
287 
288  FTensor::Tensor2<double, 3, 3> t_N{
289 
290  0.969837, -0.0860972, 0.228042,
291 
292  0.0790574, 0.996073, 0.0398449,
293 
294  -0.230577, -0.0206147, 0.972836};
295 
296  FTensor::Tensor1<double, 3> t_L{0.873555, 2.00794, 3.11851};
297 
298  FTensor::Tensor2<double, 3, 3> t_S{
299 
300  1., 0., 0.,
301 
302  0., 1., 0.,
303 
304  0., 0., 1.};
305 
306  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
307  t_S_sym(i, j) = (t_S(i, j) || t_S(j, i)) / 2.;
308 
309  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff A";
310  print_mat(t_A);
311 
312  // Testing against values in mathematica for 0,2 directive
313  {
314  auto f = [](double v) { return exp(v); };
315  auto d_f = [](double v) { return exp(v); };
316  auto dd_f = [](double v) { return exp(v); };
317 
318  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
319  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Reconstruct mat";
320  print_mat(t_b);
321 
322  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
323 
324  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff";
325  print_ddg_direction(t_d, 0, 2);
326 
327  auto t_dd =
328  EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
329 
330  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff Diff";
331  print_ddg_direction(t_dd, 0, 2);
332 
333  auto norm2_t_b = t_b(i, j) * t_b(i, j);
334  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_b " << norm2_t_b;
335 
336  auto norm2_t_d = get_norm_t4(t_d, FTensor::Number<3>());
337  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_d " << norm2_t_d;
338 
339  auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
340  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
341 
342  constexpr double regression_t_b = 572.543;
343  constexpr double regression_t_d = 859.939;
344  constexpr double regression_t_dd = 859.939;
345 
346  constexpr double eps = 1e-2;
347  if (std::abs(norm2_t_b - regression_t_b) > eps)
348  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_b");
349  if (std::abs(norm2_t_d - regression_t_d) > eps)
350  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_d");
351  if (std::abs(norm2_t_dd - regression_t_dd) > eps)
352  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_dd");
353  }
354 
355  // Comparing with lapack calculated eigen values. Note resulst should be
356  // invarinat to the direction of eiegn vector. Eigen vector can be
357  // multiplied by -1 and result should be unchanged
358  {
359 
360  std::array<double, 9> a{1., 0.1, -0.5,
361 
362  0.1, 2., 0.,
363 
364  -0.5, 0., 3.};
365 
366  auto tuple = run_lapack(a);
367  // auto &t_a = std::get<0>(tuple);
368  auto &t_eig_vec = std::get<1>(tuple);
369  auto &t_eig_vals = std::get<2>(tuple);
370 
371  auto t_eig_val_diff =
372  (t_eig_vals(i) - t_L(i)) * (t_eig_vals(i) - t_L(i));
373  MOFEM_LOG("ATOM_TEST", Sev::inform)
374  << "t_eig_val_diff " << t_eig_val_diff;
375 
376  auto f = [](double v) { return exp(v); };
377 
378  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
379  auto t_c = EigenMatrix::getMat(t_eig_vals, t_eig_vec, f);
380  t_c(i, j) -= t_b(i, j);
381  print_mat(t_c);
382 
383  auto norm2_t_c = t_c(i, j) * t_c(i, j);
384  MOFEM_LOG("ATOM_TEST", Sev::inform)
385  << "Reconstruct mat difference with lapack eigen valyes and "
386  "vectors "
387  << norm2_t_c;
388 
389  constexpr double eps = 1e-8;
390  if (fabs(norm2_t_c) > eps)
391  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
392  "Matrix not reeconstructed");
393  }
394 
395  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
396  FTensor::Ddg<double, 3, 3> t_one;
397  t_one(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
398 
399  // Testsing linear function second second direvarive zero
400  {
401  auto f = [](double v) { return v; };
402  auto d_f = [](double v) { return 1; };
403  auto dd_f = [](double v) { return 0; };
404 
405  constexpr double eps = 1e-10;
406  {
407  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
408  t_b(i, j) -= (t_A(i, j) || t_A(j, i)) / 2;
409  auto norm2_t_b = t_b(i, j) * t_b(i, j);
410  MOFEM_LOG("ATOM_TEST", Sev::inform)
411  << "Result should be matrix itself " << norm2_t_b;
412  if (norm2_t_b > eps)
413  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
414  "This norm should be zero");
415  }
416 
417  {
418 
419  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
420  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
421 
422  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
423  print_ddg(t_d_a, "hand ");
424  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
425  print_ddg(t_d, "code ");
426 
427  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
428  MOFEM_LOG("ATOM_TEST", Sev::inform)
429  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
430  if (nrm2_t_d_a > eps)
431  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
432  "This norm should be zero");
433  }
434 
435  {
436  FTensor::Tensor2<double, 3, 3> t_S{
437 
438  1., 0., 0.,
439 
440  0., 1., 0.,
441 
442  0., 0., 1.};
443 
444  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
445  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
446 
447  auto t_dd =
448  EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
449 
450  auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
451  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
452  if (norm2_t_dd > eps)
453  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
454  "This norm should be zero");
455  }
456  }
457 
458  // Testsing quadratic function second second direvarive zero
459  {
460  auto f = [](double v) { return v * v; };
461  auto d_f = [](double v) { return 2 * v; };
462  auto dd_f = [](double v) { return 2; };
463 
464  constexpr double eps = 1e-9;
465 
466  // check if multiplication gives right value
467  {
468  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
469  FTensor::Tensor2<double, 3, 3> t_a;
470  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
471  print_mat(t_a);
472  auto norm2_t_a = t_a(i, j) * t_a(i, j);
473  MOFEM_LOG("ATOM_TEST", Sev::inform)
474  << "Result should be matrix times matrix " << norm2_t_a;
475  if (norm2_t_a > eps)
476  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
477  "This norm should be zero");
478  }
479 
480  // check first directive
481  {
482  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
483  print_ddg_direction(t_d, 0, 2);
484  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<3>());
485  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
486  print_ddg(t_d_a, "hand ");
487  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
488  print_ddg(t_d, "code ");
489  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
490  MOFEM_LOG("ATOM_TEST", Sev::inform)
491  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
492  if (nrm2_t_d_a > eps)
493  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
494  "This norm should be zero");
495  }
496 
497  // check second directive
498  {
499  FTensor::Tensor2<double, 3, 3> t_S{
500 
501  1., 1. / 2., 1. / 3.,
502 
503  2. / 2., 1., 2. / 3.,
504 
505  3. / 2., 1., 3. / 3.};
506 
507  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
508  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
509 
510  auto t_dd =
511  EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
512  auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
513 
514  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
515  print_ddg(t_dd_a, "hand ");
516  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
517  print_ddg(t_dd, "code ");
518 
519  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
520  MOFEM_LOG("ATOM_TEST", Sev::inform)
521  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
522  if (nrm2_t_dd_a > eps)
523  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
524  "This norm should be zero");
525  }
526  }
527  }
528 
529  // Testing two same eigen values
530  {
531 
532  std::array<double, 9> a{5., 4., 0,
533 
534  4., 5, 0.,
535 
536  0.0, 0., 9};
537 
538  auto tuple = run_lapack(a, swap01);
539  auto &t_a = std::get<0>(tuple);
540  auto &t_eig_vecs = std::get<1>(tuple);
541  auto &t_eig_vals = std::get<2>(tuple);
542 
543  auto f = [](double v) { return v; };
544  auto d_f = [](double v) { return 1; };
545 
546  constexpr double eps = 1e-10;
547  {
548  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
549  t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
550  auto norm2_t_b = t_b(i, j) * t_b(i, j);
551  MOFEM_LOG("ATOM_TEST", Sev::inform)
552  << "Result should be matrix itself " << norm2_t_b;
553  if (norm2_t_b > eps)
554  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
555  "This norm should be zero");
556  }
557 
558  {
559  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
560  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
561  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
562  print_ddg(t_d_a, "hand ");
563  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
564  print_ddg(t_d, "code ");
565  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
566  MOFEM_LOG("ATOM_TEST", Sev::inform)
567  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
568  if (nrm2_t_d_a > eps)
569  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
570  "This norm should be zero");
571  }
572 
573  {
574  auto f = [](double v) { return v * v; };
575  auto d_f = [](double v) { return 2 * v; };
576  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
577  auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
578  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
579  print_ddg(t_d_a, "hand ");
580  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
581  print_ddg(t_d, "code ");
582  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
583  MOFEM_LOG("ATOM_TEST", Sev::inform)
584  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
585  if (nrm2_t_d_a > eps)
586  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
587  "This norm should be zero");
588  }
589  }
590 
591  // Testing three same eigen values
592  {
593 
594  std::array<double, 9> a{4., 0., 0,
595 
596  0., 4., 0.,
597 
598  0.0, 0., 4.};
599 
600  auto f = [](double v) { return v; };
601  auto d_f = [](double v) { return 1; };
602 
603  auto tuple = run_lapack(a);
604  auto &t_a = std::get<0>(tuple);
605  auto &t_eig_vecs = std::get<1>(tuple);
606  auto &t_eig_vals = std::get<2>(tuple);
607 
608  constexpr double eps = 1e-10;
609  {
610  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
611  t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
612  auto norm2_t_b = t_b(i, j) * t_b(i, j);
613  MOFEM_LOG("ATOM_TEST", Sev::inform)
614  << "Result should be matrix itself " << norm2_t_b;
615  if (norm2_t_b > eps)
616  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
617  "This norm should be zero");
618  }
619 
620  {
621  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
622  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
623  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
624  print_ddg(t_d_a, "hand ");
625  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
626  print_ddg(t_d, "code ");
627  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
628  MOFEM_LOG("ATOM_TEST", Sev::inform)
629  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
630  if (nrm2_t_d_a > eps)
631  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
632  "This norm should be zero");
633  }
634 
635  {
636  auto f = [](double v) { return v * v; };
637  auto d_f = [](double v) { return 2 * v; };
638  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
639  auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
640  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
641  print_ddg(t_d_a, "hand ");
642  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
643  print_ddg(t_d, "code ");
644  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
645  MOFEM_LOG("ATOM_TEST", Sev::inform)
646  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
647  if (nrm2_t_d_a > eps)
648  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
649  "This norm should be zero");
650  }
651  }
652 
653  // check second directive
654  {
655 
656  std::array<double, 9> a{0.1, 0., 0.,
657 
658  0., 0.1, 0.,
659 
660  0., 0., 0.1};
661 
662  auto tuple = run_lapack(a);
663  // auto &t_a = std::get<0>(tuple);
664  auto &t_eig_vecs = std::get<1>(tuple);
665  auto &t_eig_vals = std::get<2>(tuple);
666 
667  t_eig_vals(0) -= 1e-4;
668  t_eig_vals(2) += 1e-4;
669 
670  constexpr double eps = 1e-10;
671 
672  auto f = [](double v) { return v; };
673  auto d_f = [](double v) { return 1; };
674  auto dd_f = [](double v) { return 0; };
675 
676  FTensor::Tensor2<double, 3, 3> t_S{
677 
678  1., 1. / 2., 1. / 3.,
679 
680  2. / 2., 1., 2. / 3.,
681 
682  3. / 2., 1., 3. / 3.};
683 
684  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
685  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
686 
687  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
688  dd_f, t_S_sym, 1);
689 
690  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
691  print_ddg(t_dd, "test ");
692 
693  double nrm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
694  MOFEM_LOG("ATOM_TEST", Sev::inform)
695  << "Direvarive hand calculation minus code " << nrm2_t_dd;
696  if (nrm2_t_dd > eps)
697  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
698  "This norm should be zero");
699  }
700 
701  // check second directive
702  {
703 
704  std::array<double, 9> a{2, 0., 0.,
705 
706  0., 2, 0.,
707 
708  0., 0., 2};
709 
710  auto tuple = run_lapack(a);
711  // auto &t_a = std::get<0>(tuple);
712  auto &t_eig_vecs = std::get<1>(tuple);
713  auto &t_eig_vals = std::get<2>(tuple);
714 
715  constexpr double eps = 1e-10;
716 
717  auto f = [](double v) { return v * v; };
718  auto d_f = [](double v) { return 2 * v; };
719  auto dd_f = [](double v) { return 2; };
720  FTensor::Tensor2<double, 3, 3> t_S{
721 
722  1., 1. / 2., 1. / 3.,
723 
724  2. / 1., 1., 2. / 3.,
725 
726  3. / 1., 3. / 1., 1.};
727 
728  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
729  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
730 
731  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
732  dd_f, t_S_sym, 1);
733  // print_ddg(t_dd, "test ");
734 
735  auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
736 
737  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
738  MOFEM_LOG("ATOM_TEST", Sev::inform)
739  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
740  if (nrm2_t_dd_a > eps)
741  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
742  "This norm should be zero");
743  }
744 
745  // check second directive two reapeating eiegn values
746  {
747 
748  std::array<double, 9> a{5., 4., 0.,
749 
750  4., 5., 0.,
751 
752  0., 0., 9};
753 
754  auto tuple = run_lapack(a, swap01);
755  // auto &t_a = std::get<0>(tuple);
756  auto &t_eig_vecs = std::get<1>(tuple);
757  auto &t_eig_vals = std::get<2>(tuple);
758 
759  constexpr double eps = 1e-10;
760 
761  auto f = [](double v) { return v * v; };
762  auto d_f = [](double v) { return 2 * v; };
763  auto dd_f = [](double v) { return 2; };
764 
765  FTensor::Tensor2<double, 3, 3> t_S{
766 
767  1., 1. / 2., 1. / 3.,
768 
769  2. / 1., 1., 2. / 3.,
770 
771  3. / 1., 3. / 1., 1.};
772 
773  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
774  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
775 
776  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
777  dd_f, t_S_sym, 2);
778  print_ddg(t_dd, "test ");
779 
780  auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
781 
782  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
783  MOFEM_LOG("ATOM_TEST", Sev::inform)
784  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
785  if (nrm2_t_dd_a > eps)
786  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
787  "This norm should be zero");
788  }
789 
790  // check second directive exponent
791  {
792 
793  std::array<double, 9> a{2, 0., 0.,
794 
795  0., 2, 0.,
796 
797  0., 0., 2};
798 
799  auto tuple = run_lapack(a);
800  // auto &t_a = std::get<0>(tuple);
801  auto &t_eig_vecs = std::get<1>(tuple);
802  auto &t_eig_vals = std::get<2>(tuple);
803 
804  t_eig_vals(0) -= 1e-5;
805  t_eig_vals(2) += 1e-5;
806 
807  constexpr double eps = 1e-7;
808 
809  auto f = [](double v) { return exp(v); };
810  auto d_f = [](double v) { return exp(v); };
811  auto dd_f = [](double v) { return exp(v); };
812  FTensor::Tensor2<double, 3, 3> t_S{
813 
814  1., 1. / 2., 1. / 3.,
815 
816  2. / 1., 1., 2. / 3.,
817 
818  3. / 1., 3. / 1., 1.};
819 
820  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
821  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
822 
823  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
824  dd_f, t_S_sym, 3);
825  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
826  dd_f, t_S_sym, 1);
827 
828  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
829  MOFEM_LOG("ATOM_TEST", Sev::verbose)
830  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
831 
832  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
833  MOFEM_LOG("ATOM_TEST", Sev::verbose)
834  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
835 
836  print_ddg(t_dd_1, "t_dd_1 ");
837  print_ddg(t_dd_2, "t_dd_2 ");
838 
839  FTensor::Ddg<double, 3, 3> t_dd_3;
840  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
841 
842  for (int ii = 0; ii != 3; ++ii)
843  for (int jj = 0; jj != 3; ++jj)
844  for (int kk = 0; kk != 3; ++kk)
845  for (int ll = 0; ll != 3; ++ll) {
846  constexpr double eps = 1e-4;
847  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
848  MOFEM_LOG("ATOM_TEST", Sev::error)
849  << "Error " << ii << " " << jj << " " << kk << " " << ll
850  << " " << t_dd_1(ii, jj, kk, ll) << " "
851  << t_dd_2(ii, jj, kk, ll);
852  }
853 
854  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
855  MOFEM_LOG("ATOM_TEST", Sev::inform)
856  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
857  if (nrm2_t_dd_3 > eps)
858  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
859  "This norm should be zero");
860  }
861 
862  // check second directive exponent agains perturned
863  {
864 
865  std::array<double, 9> a{5., 4., 0.,
866 
867  4., 5., 0.,
868 
869  0., 0., 9};
870 
871  auto tuple = run_lapack(a, swap01);
872  // auto &t_a = std::get<0>(tuple);
873  auto &t_eig_vecs = std::get<1>(tuple);
874  auto &t_eig_vals = std::get<2>(tuple);
875 
876  t_eig_vals(0) -= 1e-4;
877  t_eig_vals(2) += 1e-4;
878 
879  constexpr double eps = 1e-4;
880 
881  auto f = [](double v) { return v * v; };
882  auto d_f = [](double v) { return 2 * v; };
883  auto dd_f = [](double v) { return 2; };
884 
885  FTensor::Tensor2<double, 3, 3> t_S{
886 
887  1., 1. / 2., 1. / 3.,
888 
889  2. / 1., 1., 2. / 3.,
890 
891  3. / 1., 3. / 1., 1.};
892 
893  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
894  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
895 
896  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
897  dd_f, t_S_sym, 3);
898  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
899  dd_f, t_S_sym, 2);
900 
901  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
902  MOFEM_LOG("ATOM_TEST", Sev::verbose)
903  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
904 
905  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
906  MOFEM_LOG("ATOM_TEST", Sev::verbose)
907  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
908 
909  print_ddg(t_dd_1, "t_dd_1 ");
910  print_ddg(t_dd_2, "t_dd_2 ");
911 
912  FTensor::Ddg<double, 3, 3> t_dd_3;
913  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
914 
915  for (int ii = 0; ii != 3; ++ii)
916  for (int jj = 0; jj != 3; ++jj)
917  for (int kk = 0; kk != 3; ++kk)
918  for (int ll = 0; ll != 3; ++ll) {
919  constexpr double eps = 1e-3;
920  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
921  MOFEM_LOG("ATOM_TEST", Sev::error)
922  << "Error " << ii << " " << jj << " " << kk << " " << ll
923  << " " << t_dd_1(ii, jj, kk, ll) << " "
924  << t_dd_2(ii, jj, kk, ll);
925  }
926 
927  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
928  MOFEM_LOG("ATOM_TEST", Sev::inform)
929  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
930  if (nrm2_t_dd_3 > eps)
931  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
932  "This norm should be zero");
933  }
934 
935  // check second directive exponent
936  {
937 
938  std::array<double, 9> a{5., 4., 0.,
939 
940  4., 5., 0.,
941 
942  0., 0., 9};
943 
944  auto tuple = run_lapack(a, swap01);
945  // auto &t_a = std::get<0>(tuple);
946  auto &t_eig_vecs = std::get<1>(tuple);
947  auto &t_eig_vals = std::get<2>(tuple);
948 
949  constexpr double eps = 1e-4;
950  constexpr int p = 3;
951 
952  auto f = [](double v) { return pow(v, p); };
953  auto d_f = [](double v) { return p * pow(v, p - 1); };
954  auto dd_f = [](double v) {
955  return p * (p - 1) * pow(v, std::max(0, p - 2));
956  };
957 
958  FTensor::Tensor2<double, 3, 3> t_S{
959 
960  1., 1. / 2., 1. / 3.,
961 
962  2. / 1., 1., 2. / 3.,
963 
964  3. / 1., 3. / 1., 1.};
965 
966  // FTensor::Tensor2_symmetric<double, 3> t_S_sym;
967  // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
968 
969  t_eig_vals(0) += 2e-5;
970  t_eig_vals(2) -= 2e-5;
971  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
972  dd_f, t_S, 3);
973  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
974  dd_f, t_S, 2);
975 
976  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
977  MOFEM_LOG("ATOM_TEST", Sev::verbose)
978  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
979 
980  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
981  MOFEM_LOG("ATOM_TEST", Sev::verbose)
982  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
983 
984  print_ddg(t_dd_1, "t_dd_1 ");
985  print_ddg(t_dd_2, "t_dd_2 ");
986 
987  FTensor::Ddg<double, 3, 3> t_dd_3;
988  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
989 
990  for (int ii = 0; ii != 3; ++ii)
991  for (int jj = 0; jj != 3; ++jj)
992  for (int kk = 0; kk != 3; ++kk)
993  for (int ll = 0; ll != 3; ++ll) {
994  constexpr double eps = 1e-3;
995  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
996  MOFEM_LOG("ATOM_TEST", Sev::error)
997  << "Error " << ii << " " << jj << " " << kk << " " << ll
998  << " " << t_dd_1(ii, jj, kk, ll) << " "
999  << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
1000  }
1001 
1002  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
1003  MOFEM_LOG("ATOM_TEST", Sev::inform)
1004  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
1005  if (nrm2_t_dd_3 > eps)
1006  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1007  "This norm should be zero");
1008  }
1009 
1010  // Speed
1011  {
1012 
1013  std::array<double, 9> a{1., 0.1, -0.5,
1014 
1015  0.1, 2., 0.,
1016 
1017  -0.5, 0., 3.};
1018 
1019  auto tuple = run_lapack(a);
1020  // auto &t_a = std::get<0>(tuple);
1021  auto &t_eig_vecs = std::get<1>(tuple);
1022  auto &t_eig_vals = std::get<2>(tuple);
1023 
1024  auto f = [](double v) { return exp(v); };
1025  auto d_f = [](double v) { return exp(v); };
1026  auto dd_f = [](double v) { return exp(v); };
1027 
1028  FTensor::Tensor2<double, 3, 3> t_S{
1029 
1030  1., 1. / 2., 1. / 3.,
1031 
1032  2. / 1., 1., 2. / 3.,
1033 
1034  3. / 1., 3. / 1., 1.};
1035 
1036  FTensor::Tensor2_symmetric<double, 3> t_S_sym;
1037  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1038 
1039  MOFEM_LOG("ATOM_TEST", Sev::inform) << "Start";
1040  for (int ii = 0; ii != 1000; ++ii) {
1041  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 3);
1042  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1043  dd_f, t_S_sym, 3);
1044  std::ignore = t_d;
1045  std::ignore = t_dd;
1046  }
1047  MOFEM_LOG("ATOM_TEST", Sev::inform) << "End";
1048  }
1049 
1050  // 2d case
1051 
1052  auto run_lapack_2d = [](auto &a) {
1053  int info;
1054  double wkopt;
1055  double w[2];
1056 
1057  FTensor::Tensor2<double, 2, 2> t_a{
1058 
1059  a[0], a[1],
1060 
1061  a[2], a[3]};
1062 
1063  /* Query and allocate the optimal workspace */
1064  int lwork = -1;
1065  info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &wkopt, lwork);
1066  if (info > 0)
1067  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1068  lwork = (int)wkopt;
1069  std::vector<double> work(lwork);
1070  /* Solve eigenproblem */
1071  info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &*work.begin(), lwork);
1072  if (info > 0)
1073  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1074 
1075  FTensor::Tensor2<double, 2, 2> t_eig_vecs{
1076 
1077  a[0 * 2 + 0], a[0 * 2 + 1],
1078 
1079  a[1 * 2 + 0], a[1 * 2 + 1]};
1080 
1081  FTensor::Tensor1<double, 2> t_eig_vals{w[0], w[1]};
1082 
1083  return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1084  };
1085 
1086  // Testing quadratic function for 2d
1087  {
1088 
1089  std::array<double, 9> a{1., 0.1,
1090 
1091  0.1, 2.};
1092 
1093  auto tuple = run_lapack_2d(a);
1094  auto &t_A = std::get<0>(tuple);
1095  auto &t_eig_vecs = std::get<1>(tuple);
1096  auto &t_eig_vals = std::get<2>(tuple);
1097 
1098  auto f = [](double v) { return v * v; };
1099  auto d_f = [](double v) { return 2 * v; };
1100  auto dd_f = [](double v) { return 2; };
1101 
1102  constexpr double eps = 1e-6;
1103 
1104  FTensor::Index<'i', 2> i;
1105  FTensor::Index<'j', 2> j;
1106  FTensor::Index<'k', 2> k;
1107  FTensor::Index<'l', 2> l;
1108 
1109  // check if multiplication gives right value
1110  {
1111  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1112  FTensor::Tensor2<double, 2, 2> t_a;
1113  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1114  print_mat(t_a);
1115  auto norm2_t_a = t_a(i, j) * t_a(i, j);
1116  MOFEM_LOG("ATOM_TEST", Sev::inform)
1117  << "Result should be matrix times matrix " << norm2_t_a;
1118  if (norm2_t_a > eps)
1119  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1120  "This norm should be zero");
1121  }
1122 
1123  // check first directive
1124  {
1125  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
1126  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1127  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1128  MOFEM_LOG("ATOM_TEST", Sev::inform)
1129  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1130  if (nrm2_t_d_a > eps)
1131  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1132  "This norm should be zero");
1133  }
1134 
1135  // check second directive
1136  {
1137  FTensor::Tensor2<double, 2, 2> t_S{
1138 
1139  1., 1. / 2,
1140 
1141  2. / 2., 1.};
1142 
1143  // FTensor::Index<'i', 2> i;
1144  // FTensor::Index<'j', 2> j;
1145  // FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1146  // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1147 
1148  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1149  dd_f, t_S, 2);
1150  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<2>());
1151 
1152  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1153  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1154 
1155  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1156  MOFEM_LOG("ATOM_TEST", Sev::inform)
1157  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1158  if (nrm2_t_dd_a > eps)
1159  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1160  "This norm should be zero");
1161  }
1162  }
1163 
1164  // Testing quadratic function for repeating eigen values
1165  {
1166 
1167  std::array<double, 9> a{2., 0,
1168 
1169  0, 2.};
1170 
1171  auto tuple = run_lapack_2d(a);
1172  auto &t_A = std::get<0>(tuple);
1173  auto &t_eig_vecs = std::get<1>(tuple);
1174  auto &t_eig_vals = std::get<2>(tuple);
1175 
1176  auto f = [](double v) { return v * v; };
1177  auto d_f = [](double v) { return 2 * v; };
1178  auto dd_f = [](double v) { return 2; };
1179 
1180  constexpr double eps = 1e-6;
1181 
1182  FTensor::Index<'i', 2> i;
1183  FTensor::Index<'j', 2> j;
1184  FTensor::Index<'k', 2> k;
1185  FTensor::Index<'l', 2> l;
1186 
1187  // check if multiplication gives right value
1188  {
1189  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1190  FTensor::Tensor2<double, 2, 2> t_a;
1191  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1192  print_mat(t_a);
1193  auto norm2_t_a = t_a(i, j) * t_a(i, j);
1194  MOFEM_LOG("ATOM_TEST", Sev::inform)
1195  << "Result should be matrix times matrix " << norm2_t_a;
1196  if (norm2_t_a > eps)
1197  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1198  "This norm should be zero");
1199  }
1200 
1201  // check first directive
1202  {
1203  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
1204  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1205  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1206  MOFEM_LOG("ATOM_TEST", Sev::inform)
1207  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1208  if (nrm2_t_d_a > eps)
1209  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1210  "This norm should be zero");
1211  }
1212 
1213  // check second directive
1214  {
1215  FTensor::Tensor2<double, 2, 2> t_S{
1216 
1217  1., 1. / 2,
1218 
1219  2. / 2., 1.};
1220 
1221  FTensor::Index<'i', 2> i;
1222  FTensor::Index<'j', 2> j;
1223  FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1224  t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1225 
1226  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1227  dd_f, t_S_sym, 1);
1228  auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<2>());
1229 
1230  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1231  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1232 
1233  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1234  MOFEM_LOG("ATOM_TEST", Sev::inform)
1235  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1236  if (nrm2_t_dd_a > eps)
1237  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1238  "This norm should be zero");
1239  }
1240  }
1241  }
1242  CATCH_ERRORS;
1243 
1244  CHKERR MoFEM::Core::Finalize();
1245 }
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
FTensor::Tensor1< double, 3 >
print_mat
void print_mat(double *M, int m, int n)
print matric M
main
int main(int argc, char *argv[])
Definition: matrix_function.cpp:182
help
static char help[]
Definition: matrix_function.cpp:180
get_diff2_matrix2
auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:82
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MOFEM_LOG_ATTRIBUTES
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:296
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
MoFEM::LogManager::BitLineID
@ BitLineID
Definition: LogManager.hpp:48
FTensor::Number
Definition: Number.hpp:11
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MatrixFunction.hpp
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
get_norm_t4
auto get_norm_t4(T &t, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:170
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
FTensor::Index< 'i', 3 >
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
diff_ddg
void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:24
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
FTensor::Ddg< double, 3, 3 >
lapack_dsyev
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:261
EigenMatrix::getDiffDiffMat
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
Definition: MatrixFunction.cpp:78
FTensor.hpp
Tensors class implemented by Walter Landry.
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
get_diff_matrix
auto get_diff_matrix(T1 &t_d, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:140
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
get_diff_matrix2
auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:45
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
convert.int
int
Definition: convert.py:64
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
Generated by Doxygen 1.8.17 and hosted at University of Glasgow