v0.13.2
Loading...
Searching...
No Matches
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>
14using namespace MoFEM;
15
16#include <MatrixFunction.hpp>
17
18FTensor::Index<'i', 3> i;
19FTensor::Index<'j', 3> j;
20FTensor::Index<'k', 3> k;
21FTensor::Index<'l', 3> l;
22
23template <typename T1, typename T2, int DIM>
24void 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
44template <typename T1, typename T2, int DIM>
45auto 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
81template <typename T1, typename T2, int DIM>
82auto 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
139template <typename T1, int DIM>
140auto 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
169template <typename T, int DIM>
170auto 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
180static char help[] = "...\n\n";
181
182int 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 auto d_f = [](double v) { return exp(v); };
378 auto dd_f = [](double v) { return exp(v); };
379
380 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
381 auto t_c = EigenMatrix::getMat(t_eig_vals, t_eig_vec, f);
382 t_c(i, j) -= t_b(i, j);
383 print_mat(t_c);
384
385 auto norm2_t_c = t_c(i, j) * t_c(i, j);
386 MOFEM_LOG("ATOM_TEST", Sev::inform)
387 << "Reconstruct mat difference with lapack eigen valyes and "
388 "vectors "
389 << norm2_t_c;
390
391 constexpr double eps = 1e-8;
392 if (fabs(norm2_t_c) > eps)
393 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
394 "Matrix not reeconstructed");
395 }
396
397 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
398 FTensor::Ddg<double, 3, 3> t_one;
399 t_one(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
400
401 // Testsing linear function second second direvarive zero
402 {
403 auto f = [](double v) { return v; };
404 auto d_f = [](double v) { return 1; };
405 auto dd_f = [](double v) { return 0; };
406
407 constexpr double eps = 1e-10;
408 {
409 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
410 t_b(i, j) -= (t_A(i, j) || t_A(j, i)) / 2;
411 auto norm2_t_b = t_b(i, j) * t_b(i, j);
412 MOFEM_LOG("ATOM_TEST", Sev::inform)
413 << "Result should be matrix itself " << norm2_t_b;
414 if (norm2_t_b > eps)
415 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
416 "This norm should be zero");
417 }
418
419 {
420
421 auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
422 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
423
424 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
425 print_ddg(t_d_a, "hand ");
426 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
427 print_ddg(t_d, "code ");
428
429 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
430 MOFEM_LOG("ATOM_TEST", Sev::inform)
431 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
432 if (nrm2_t_d_a > eps)
433 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
434 "This norm should be zero");
435 }
436
437 {
438 FTensor::Tensor2<double, 3, 3> t_S{
439
440 1., 0., 0.,
441
442 0., 1., 0.,
443
444 0., 0., 1.};
445
446 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
447 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
448
449 auto t_dd =
450 EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
451
452 auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
453 MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
454 if (norm2_t_dd > eps)
455 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
456 "This norm should be zero");
457 }
458 }
459
460 // Testsing quadratic function second second direvarive zero
461 {
462 auto f = [](double v) { return v * v; };
463 auto d_f = [](double v) { return 2 * v; };
464 auto dd_f = [](double v) { return 2; };
465
466 constexpr double eps = 1e-9;
467
468 // check if multiplication gives right value
469 {
470 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
471 FTensor::Tensor2<double, 3, 3> t_a;
472 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
473 print_mat(t_a);
474 auto norm2_t_a = t_a(i, j) * t_a(i, j);
475 MOFEM_LOG("ATOM_TEST", Sev::inform)
476 << "Result should be matrix times matrix " << norm2_t_a;
477 if (norm2_t_a > eps)
478 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
479 "This norm should be zero");
480 }
481
482 // check first directive
483 {
484 auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
485 print_ddg_direction(t_d, 0, 2);
486 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<3>());
487 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
488 print_ddg(t_d_a, "hand ");
489 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
490 print_ddg(t_d, "code ");
491 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
492 MOFEM_LOG("ATOM_TEST", Sev::inform)
493 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
494 if (nrm2_t_d_a > eps)
495 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
496 "This norm should be zero");
497 }
498
499 // check second directive
500 {
501 FTensor::Tensor2<double, 3, 3> t_S{
502
503 1., 1. / 2., 1. / 3.,
504
505 2. / 2., 1., 2. / 3.,
506
507 3. / 2., 1., 3. / 3.};
508
509 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
510 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
511
512 auto t_dd =
513 EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
514 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
515
516 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
517 print_ddg(t_dd_a, "hand ");
518 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
519 print_ddg(t_dd, "code ");
520
521 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
522 MOFEM_LOG("ATOM_TEST", Sev::inform)
523 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
524 if (nrm2_t_dd_a > eps)
525 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
526 "This norm should be zero");
527 }
528 }
529 }
530
531 // Testing two same eigen values
532 {
533
534 std::array<double, 9> a{5., 4., 0,
535
536 4., 5, 0.,
537
538 0.0, 0., 9};
539
540 auto tuple = run_lapack(a, swap01);
541 auto &t_a = std::get<0>(tuple);
542 auto &t_eig_vecs = std::get<1>(tuple);
543 auto &t_eig_vals = std::get<2>(tuple);
544
545 auto f = [](double v) { return v; };
546 auto d_f = [](double v) { return 1; };
547 auto dd_f = [](double v) { return 0; };
548
549 constexpr double eps = 1e-10;
550 {
551 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
552 t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
553 auto norm2_t_b = t_b(i, j) * t_b(i, j);
554 MOFEM_LOG("ATOM_TEST", Sev::inform)
555 << "Result should be matrix itself " << norm2_t_b;
556 if (norm2_t_b > eps)
557 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
558 "This norm should be zero");
559 }
560
561 {
562 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
563 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
564 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
565 print_ddg(t_d_a, "hand ");
566 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
567 print_ddg(t_d, "code ");
568 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
569 MOFEM_LOG("ATOM_TEST", Sev::inform)
570 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
571 if (nrm2_t_d_a > eps)
572 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
573 "This norm should be zero");
574 }
575
576 {
577 auto f = [](double v) { return v * v; };
578 auto d_f = [](double v) { return 2 * v; };
579 auto dd_f = [](double v) { return 2; };
580 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
581 auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
582 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
583 print_ddg(t_d_a, "hand ");
584 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
585 print_ddg(t_d, "code ");
586 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
587 MOFEM_LOG("ATOM_TEST", Sev::inform)
588 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
589 if (nrm2_t_d_a > eps)
590 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
591 "This norm should be zero");
592 }
593 }
594
595 // Testing three same eigen values
596 {
597
598 std::array<double, 9> a{4., 0., 0,
599
600 0., 4., 0.,
601
602 0.0, 0., 4.};
603
604 auto f = [](double v) { return v; };
605 auto d_f = [](double v) { return 1; };
606 auto dd_f = [](double v) { return 0; };
607
608 auto tuple = run_lapack(a);
609 auto &t_a = std::get<0>(tuple);
610 auto &t_eig_vecs = std::get<1>(tuple);
611 auto &t_eig_vals = std::get<2>(tuple);
612
613 constexpr double eps = 1e-10;
614 {
615 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
616 t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
617 auto norm2_t_b = t_b(i, j) * t_b(i, j);
618 MOFEM_LOG("ATOM_TEST", Sev::inform)
619 << "Result should be matrix itself " << norm2_t_b;
620 if (norm2_t_b > eps)
621 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
622 "This norm should be zero");
623 }
624
625 {
626 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
627 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
628 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
629 print_ddg(t_d_a, "hand ");
630 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
631 print_ddg(t_d, "code ");
632 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
633 MOFEM_LOG("ATOM_TEST", Sev::inform)
634 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
635 if (nrm2_t_d_a > eps)
636 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
637 "This norm should be zero");
638 }
639
640 {
641 auto f = [](double v) { return v * v; };
642 auto d_f = [](double v) { return 2 * v; };
643 auto dd_f = [](double v) { return 2; };
644 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
645 auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
646 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
647 print_ddg(t_d_a, "hand ");
648 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
649 print_ddg(t_d, "code ");
650 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
651 MOFEM_LOG("ATOM_TEST", Sev::inform)
652 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
653 if (nrm2_t_d_a > eps)
654 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
655 "This norm should be zero");
656 }
657 }
658
659 // check second directive
660 {
661
662 std::array<double, 9> a{0.1, 0., 0.,
663
664 0., 0.1, 0.,
665
666 0., 0., 0.1};
667
668 auto tuple = run_lapack(a);
669 auto &t_a = std::get<0>(tuple);
670 auto &t_eig_vecs = std::get<1>(tuple);
671 auto &t_eig_vals = std::get<2>(tuple);
672
673 t_eig_vals(0) -= 1e-4;
674 t_eig_vals(2) += 1e-4;
675
676 constexpr double eps = 1e-10;
677
678 auto f = [](double v) { return v; };
679 auto d_f = [](double v) { return 1; };
680 auto dd_f = [](double v) { return 0; };
681
682 FTensor::Tensor2<double, 3, 3> t_S{
683
684 1., 1. / 2., 1. / 3.,
685
686 2. / 2., 1., 2. / 3.,
687
688 3. / 2., 1., 3. / 3.};
689
690 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
691 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
692
693 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
694 dd_f, t_S_sym, 1);
695
696 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
697 print_ddg(t_dd, "test ");
698
699 double nrm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
700 MOFEM_LOG("ATOM_TEST", Sev::inform)
701 << "Direvarive hand calculation minus code " << nrm2_t_dd;
702 if (nrm2_t_dd > eps)
703 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
704 "This norm should be zero");
705 }
706
707 // check second directive
708 {
709
710 std::array<double, 9> a{2, 0., 0.,
711
712 0., 2, 0.,
713
714 0., 0., 2};
715
716 auto tuple = run_lapack(a);
717 auto &t_a = std::get<0>(tuple);
718 auto &t_eig_vecs = std::get<1>(tuple);
719 auto &t_eig_vals = std::get<2>(tuple);
720
721 constexpr double eps = 1e-10;
722
723 auto f = [](double v) { return v * v; };
724 auto d_f = [](double v) { return 2 * v; };
725 auto dd_f = [](double v) { return 2; };
726 FTensor::Tensor2<double, 3, 3> t_S{
727
728 1., 1. / 2., 1. / 3.,
729
730 2. / 1., 1., 2. / 3.,
731
732 3. / 1., 3. / 1., 1.};
733
734 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
735 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
736
737 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
738 dd_f, t_S_sym, 1);
739 // print_ddg(t_dd, "test ");
740
741 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
742
743 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
744 MOFEM_LOG("ATOM_TEST", Sev::inform)
745 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
746 if (nrm2_t_dd_a > eps)
747 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
748 "This norm should be zero");
749 }
750
751 // check second directive two reapeating eiegn values
752 {
753
754 std::array<double, 9> a{5., 4., 0.,
755
756 4., 5., 0.,
757
758 0., 0., 9};
759
760 auto tuple = run_lapack(a, swap01);
761 auto &t_a = std::get<0>(tuple);
762 auto &t_eig_vecs = std::get<1>(tuple);
763 auto &t_eig_vals = std::get<2>(tuple);
764
765 constexpr double eps = 1e-10;
766
767 auto f = [](double v) { return v * v; };
768 auto d_f = [](double v) { return 2 * v; };
769 auto dd_f = [](double v) { return 2; };
770
771 FTensor::Tensor2<double, 3, 3> t_S{
772
773 1., 1. / 2., 1. / 3.,
774
775 2. / 1., 1., 2. / 3.,
776
777 3. / 1., 3. / 1., 1.};
778
779 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
780 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
781
782 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
783 dd_f, t_S_sym, 2);
784 print_ddg(t_dd, "test ");
785
786 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
787
788 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
789 MOFEM_LOG("ATOM_TEST", Sev::inform)
790 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
791 if (nrm2_t_dd_a > eps)
792 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
793 "This norm should be zero");
794 }
795
796 // check second directive exponent
797 {
798
799 std::array<double, 9> a{2, 0., 0.,
800
801 0., 2, 0.,
802
803 0., 0., 2};
804
805 auto tuple = run_lapack(a);
806 auto &t_a = std::get<0>(tuple);
807 auto &t_eig_vecs = std::get<1>(tuple);
808 auto &t_eig_vals = std::get<2>(tuple);
809
810 t_eig_vals(0) -= 1e-5;
811 t_eig_vals(2) += 1e-5;
812
813 constexpr double eps = 1e-7;
814
815 auto f = [](double v) { return exp(v); };
816 auto d_f = [](double v) { return exp(v); };
817 auto dd_f = [](double v) { return exp(v); };
818 FTensor::Tensor2<double, 3, 3> t_S{
819
820 1., 1. / 2., 1. / 3.,
821
822 2. / 1., 1., 2. / 3.,
823
824 3. / 1., 3. / 1., 1.};
825
826 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
827 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
828
829 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
830 dd_f, t_S_sym, 3);
831 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
832 dd_f, t_S_sym, 1);
833
834 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
835 MOFEM_LOG("ATOM_TEST", Sev::verbose)
836 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
837
838 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
839 MOFEM_LOG("ATOM_TEST", Sev::verbose)
840 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
841
842 print_ddg(t_dd_1, "t_dd_1 ");
843 print_ddg(t_dd_2, "t_dd_2 ");
844
845 FTensor::Ddg<double, 3, 3> t_dd_3;
846 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
847
848 for (int ii = 0; ii != 3; ++ii)
849 for (int jj = 0; jj != 3; ++jj)
850 for (int kk = 0; kk != 3; ++kk)
851 for (int ll = 0; ll != 3; ++ll) {
852 constexpr double eps = 1e-4;
853 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
854 MOFEM_LOG("ATOM_TEST", Sev::error)
855 << "Error " << ii << " " << jj << " " << kk << " " << ll
856 << " " << t_dd_1(ii, jj, kk, ll) << " "
857 << t_dd_2(ii, jj, kk, ll);
858 }
859
860 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
861 MOFEM_LOG("ATOM_TEST", Sev::inform)
862 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
863 if (nrm2_t_dd_3 > eps)
864 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
865 "This norm should be zero");
866 }
867
868 // check second directive exponent agains perturned
869 {
870
871 std::array<double, 9> a{5., 4., 0.,
872
873 4., 5., 0.,
874
875 0., 0., 9};
876
877 auto tuple = run_lapack(a, swap01);
878 auto &t_a = std::get<0>(tuple);
879 auto &t_eig_vecs = std::get<1>(tuple);
880 auto &t_eig_vals = std::get<2>(tuple);
881
882 t_eig_vals(0) -= 1e-4;
883 t_eig_vals(2) += 1e-4;
884
885 constexpr double eps = 1e-4;
886
887 auto f = [](double v) { return v * v; };
888 auto d_f = [](double v) { return 2 * v; };
889 auto dd_f = [](double v) { return 2; };
890
891 FTensor::Tensor2<double, 3, 3> t_S{
892
893 1., 1. / 2., 1. / 3.,
894
895 2. / 1., 1., 2. / 3.,
896
897 3. / 1., 3. / 1., 1.};
898
899 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
900 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
901
902 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
903 dd_f, t_S_sym, 3);
904 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
905 dd_f, t_S_sym, 2);
906
907 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
908 MOFEM_LOG("ATOM_TEST", Sev::verbose)
909 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
910
911 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
912 MOFEM_LOG("ATOM_TEST", Sev::verbose)
913 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
914
915 print_ddg(t_dd_1, "t_dd_1 ");
916 print_ddg(t_dd_2, "t_dd_2 ");
917
918 FTensor::Ddg<double, 3, 3> t_dd_3;
919 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
920
921 for (int ii = 0; ii != 3; ++ii)
922 for (int jj = 0; jj != 3; ++jj)
923 for (int kk = 0; kk != 3; ++kk)
924 for (int ll = 0; ll != 3; ++ll) {
925 constexpr double eps = 1e-3;
926 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
927 MOFEM_LOG("ATOM_TEST", Sev::error)
928 << "Error " << ii << " " << jj << " " << kk << " " << ll
929 << " " << t_dd_1(ii, jj, kk, ll) << " "
930 << t_dd_2(ii, jj, kk, ll);
931 }
932
933 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
934 MOFEM_LOG("ATOM_TEST", Sev::inform)
935 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
936 if (nrm2_t_dd_3 > eps)
937 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
938 "This norm should be zero");
939 }
940
941 // check second directive exponent
942 {
943
944 std::array<double, 9> a{5., 4., 0.,
945
946 4., 5., 0.,
947
948 0., 0., 9};
949
950 auto tuple = run_lapack(a, swap01);
951 auto &t_a = std::get<0>(tuple);
952 auto &t_eig_vecs = std::get<1>(tuple);
953 auto &t_eig_vals = std::get<2>(tuple);
954
955 constexpr double eps = 1e-4;
956 constexpr int p = 3;
957
958 auto f = [](double v) { return pow(v, p); };
959 auto d_f = [](double v) { return p * pow(v, p - 1); };
960 auto dd_f = [](double v) {
961 return p * (p - 1) * pow(v, std::max(0, p - 2));
962 };
963
964 FTensor::Tensor2<double, 3, 3> t_S{
965
966 1., 1. / 2., 1. / 3.,
967
968 2. / 1., 1., 2. / 3.,
969
970 3. / 1., 3. / 1., 1.};
971
972 // FTensor::Tensor2_symmetric<double, 3> t_S_sym;
973 // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
974
975 t_eig_vals(0) += 2e-5;
976 t_eig_vals(2) -= 2e-5;
977 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
978 dd_f, t_S, 3);
979 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
980 dd_f, t_S, 2);
981
982 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
983 MOFEM_LOG("ATOM_TEST", Sev::verbose)
984 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
985
986 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
987 MOFEM_LOG("ATOM_TEST", Sev::verbose)
988 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
989
990 print_ddg(t_dd_1, "t_dd_1 ");
991 print_ddg(t_dd_2, "t_dd_2 ");
992
993 FTensor::Ddg<double, 3, 3> t_dd_3;
994 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
995
996 for (int ii = 0; ii != 3; ++ii)
997 for (int jj = 0; jj != 3; ++jj)
998 for (int kk = 0; kk != 3; ++kk)
999 for (int ll = 0; ll != 3; ++ll) {
1000 constexpr double eps = 1e-3;
1001 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
1002 MOFEM_LOG("ATOM_TEST", Sev::error)
1003 << "Error " << ii << " " << jj << " " << kk << " " << ll
1004 << " " << t_dd_1(ii, jj, kk, ll) << " "
1005 << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
1006 }
1007
1008 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
1009 MOFEM_LOG("ATOM_TEST", Sev::inform)
1010 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
1011 if (nrm2_t_dd_3 > eps)
1012 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1013 "This norm should be zero");
1014 }
1015
1016 // Speed
1017 {
1018
1019 std::array<double, 9> a{1., 0.1, -0.5,
1020
1021 0.1, 2., 0.,
1022
1023 -0.5, 0., 3.};
1024
1025 auto tuple = run_lapack(a);
1026 auto &t_a = std::get<0>(tuple);
1027 auto &t_eig_vecs = std::get<1>(tuple);
1028 auto &t_eig_vals = std::get<2>(tuple);
1029
1030 auto f = [](double v) { return exp(v); };
1031 auto d_f = [](double v) { return exp(v); };
1032 auto dd_f = [](double v) { return exp(v); };
1033
1034 FTensor::Tensor2<double, 3, 3> t_S{
1035
1036 1., 1. / 2., 1. / 3.,
1037
1038 2. / 1., 1., 2. / 3.,
1039
1040 3. / 1., 3. / 1., 1.};
1041
1042 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
1043 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1044
1045 MOFEM_LOG("ATOM_TEST", Sev::inform) << "Start";
1046 for (int ii = 0; ii != 1000; ++ii) {
1047 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 3);
1048 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1049 dd_f, t_S_sym, 3);
1050 }
1051 MOFEM_LOG("ATOM_TEST", Sev::inform) << "End";
1052 }
1053
1054 // 2d case
1055
1056 auto run_lapack_2d = [](auto &a) {
1057 int info;
1058 double wkopt;
1059 double w[2];
1060
1061 FTensor::Tensor2<double, 2, 2> t_a{
1062
1063 a[0], a[1],
1064
1065 a[2], a[3]};
1066
1067 /* Query and allocate the optimal workspace */
1068 int lwork = -1;
1069 info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &wkopt, lwork);
1070 if (info > 0)
1071 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1072 lwork = (int)wkopt;
1073 std::vector<double> work(lwork);
1074 /* Solve eigenproblem */
1075 info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &*work.begin(), lwork);
1076 if (info > 0)
1077 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1078
1079 FTensor::Tensor2<double, 2, 2> t_eig_vecs{
1080
1081 a[0 * 2 + 0], a[0 * 2 + 1],
1082
1083 a[1 * 2 + 0], a[1 * 2 + 1]};
1084
1085 FTensor::Tensor1<double, 2> t_eig_vals{w[0], w[1]};
1086
1087 return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1088 };
1089
1090 // Testsing quadratic function for 2d
1091 {
1092
1093 std::array<double, 9> a{1., 0.1,
1094
1095 0.1, 2.};
1096
1097 auto tuple = run_lapack_2d(a);
1098 auto &t_A = std::get<0>(tuple);
1099 auto &t_eig_vecs = std::get<1>(tuple);
1100 auto &t_eig_vals = std::get<2>(tuple);
1101
1102 auto f = [](double v) { return v * v; };
1103 auto d_f = [](double v) { return 2 * v; };
1104 auto dd_f = [](double v) { return 2; };
1105
1106 constexpr double eps = 1e-6;
1107
1108 FTensor::Index<'i', 2> i;
1109 FTensor::Index<'j', 2> j;
1110 FTensor::Index<'k', 2> k;
1111 FTensor::Index<'l', 2> l;
1112
1113 // check if multiplication gives right value
1114 {
1115 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1116 FTensor::Tensor2<double, 2, 2> t_a;
1117 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1118 print_mat(t_a);
1119 auto norm2_t_a = t_a(i, j) * t_a(i, j);
1120 MOFEM_LOG("ATOM_TEST", Sev::inform)
1121 << "Result should be matrix times matrix " << norm2_t_a;
1122 if (norm2_t_a > eps)
1123 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1124 "This norm should be zero");
1125 }
1126
1127 // check first directive
1128 {
1129 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
1130 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1131 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1132 MOFEM_LOG("ATOM_TEST", Sev::inform)
1133 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1134 if (nrm2_t_d_a > eps)
1135 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1136 "This norm should be zero");
1137 }
1138
1139 // check second directive
1140 {
1141 FTensor::Tensor2<double, 2, 2> t_S{
1142
1143 1., 1. / 2,
1144
1145 2. / 2., 1.};
1146
1147 // FTensor::Index<'i', 2> i;
1148 // FTensor::Index<'j', 2> j;
1149 // FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1150 // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1151
1152 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1153 dd_f, t_S, 2);
1154 auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<2>());
1155
1156 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1157 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1158
1159 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1160 MOFEM_LOG("ATOM_TEST", Sev::inform)
1161 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1162 if (nrm2_t_dd_a > eps)
1163 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1164 "This norm should be zero");
1165 }
1166 }
1167
1168 // Testsing quadratic function for repeating eigen valsues
1169 {
1170
1171 std::array<double, 9> a{2., 0,
1172
1173 0, 2.};
1174
1175 auto tuple = run_lapack_2d(a);
1176 auto &t_A = std::get<0>(tuple);
1177 auto &t_eig_vecs = std::get<1>(tuple);
1178 auto &t_eig_vals = std::get<2>(tuple);
1179
1180 auto f = [](double v) { return v * v; };
1181 auto d_f = [](double v) { return 2 * v; };
1182 auto dd_f = [](double v) { return 2; };
1183
1184 constexpr double eps = 1e-6;
1185
1186 FTensor::Index<'i', 2> i;
1187 FTensor::Index<'j', 2> j;
1188 FTensor::Index<'k', 2> k;
1189 FTensor::Index<'l', 2> l;
1190
1191 // check if multiplication gives right value
1192 {
1193 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1194 FTensor::Tensor2<double, 2, 2> t_a;
1195 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1196 print_mat(t_a);
1197 auto norm2_t_a = t_a(i, j) * t_a(i, j);
1198 MOFEM_LOG("ATOM_TEST", Sev::inform)
1199 << "Result should be matrix times matrix " << norm2_t_a;
1200 if (norm2_t_a > eps)
1201 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1202 "This norm should be zero");
1203 }
1204
1205 // check first directive
1206 {
1207 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
1208 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1209 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1210 MOFEM_LOG("ATOM_TEST", Sev::inform)
1211 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1212 if (nrm2_t_d_a > eps)
1213 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1214 "This norm should be zero");
1215 }
1216
1217 // check second directive
1218 {
1219 FTensor::Tensor2<double, 2, 2> t_S{
1220
1221 1., 1. / 2,
1222
1223 2. / 2., 1.};
1224
1225 FTensor::Index<'i', 2> i;
1226 FTensor::Index<'j', 2> j;
1227 FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1228 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1229
1230 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1231 dd_f, t_S_sym, 1);
1232 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<2>());
1233
1234 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1235 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1236
1237 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1238 MOFEM_LOG("ATOM_TEST", Sev::inform)
1239 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1240 if (nrm2_t_dd_a > eps)
1241 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1242 "This norm should be zero");
1243 }
1244 }
1245 }
1246 CATCH_ERRORS;
1247
1248 CHKERR MoFEM::Core::Finalize();
1249}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
FTensor.hpp
Tensors class implemented by Walter Landry.
main
int main()
Definition: adol-c_atom.cpp:46
a
constexpr double a
Definition: approx_sphere.cpp:30
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:14
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
print_mat
void print_mat(double *M, int m, int n)
print matric M
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:73
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:391
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MOFEM_LOG_ATTRIBUTES
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:296
MoFEM::LogManager::BitLineID
@ BitLineID
Definition: LogManager.hpp:48
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
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
get_diff_matrix
auto get_diff_matrix(T1 &t_d, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:140
diff_ddg
void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:24
help
static char help[]
Definition: matrix_function.cpp:180
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
get_diff2_matrix2
auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:82
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
get_norm_t4
auto get_norm_t4(T &t, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:170
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
get_diff_matrix2
auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:45
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
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
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
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
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::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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:300
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:342
Generated by Doxygen 1.9.5 and hosted at University of Glasgow