v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
matrix_function.cpp File Reference
#include <FTensor.hpp>
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>

Go to the source code of this file.

Functions

template<typename T1 , typename T2 , int DIM>
void diff_ddg (T1 &t_1, T2 &t_2, const FTensor::Number< DIM > &)
 
template<typename T1 , typename T2 , int DIM>
auto get_diff_matrix2 (T1 &t_a, T2 &t_d, const FTensor::Number< DIM > &)
 
template<typename T1 , typename T2 , int DIM>
auto get_diff2_matrix2 (T1 &t_s, T2 &t_dd, const FTensor::Number< DIM > &)
 
template<typename T1 , int DIM>
auto get_diff_matrix (T1 &t_d, const FTensor::Number< DIM > &)
 
template<typename T , int DIM>
auto get_norm_t4 (T &t, const FTensor::Number< DIM > &)
 
int main (int argc, char *argv[])
 

Variables

FTensor::Index< 'i', 3 > i
 
FTensor::Index< 'j', 3 > j
 
FTensor::Index< 'k', 3 > k
 
FTensor::Index< 'l', 3 > l
 
static char help [] = "...\n\n"
 

Function Documentation

◆ diff_ddg()

template<typename T1 , typename T2 , int DIM>
void diff_ddg ( T1 &  t_1,
T2 &  t_2,
const FTensor::Number< DIM > &   
)
Examples
matrix_function.cpp.

Definition at line 24 of file matrix_function.cpp.

24 {
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};
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308

◆ get_diff2_matrix2()

template<typename T1 , typename T2 , int DIM>
auto get_diff2_matrix2 ( T1 &  t_s,
T2 &  t_dd,
const FTensor::Number< DIM > &   
)
Examples
matrix_function.cpp.

Definition at line 82 of file matrix_function.cpp.

82 {
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};
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
diff_ddg
void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:24
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20

◆ get_diff_matrix()

template<typename T1 , int DIM>
auto get_diff_matrix ( T1 &  t_d,
const FTensor::Number< DIM > &   
)
Examples
matrix_function.cpp.

Definition at line 140 of file matrix_function.cpp.

140 {
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};

◆ get_diff_matrix2()

template<typename T1 , typename T2 , int DIM>
auto get_diff_matrix2 ( T1 &  t_a,
T2 &  t_d,
const FTensor::Number< DIM > &   
)
Examples
matrix_function.cpp.

Definition at line 45 of file matrix_function.cpp.

45 {
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};

◆ get_norm_t4()

template<typename T , int DIM>
auto get_norm_t4 ( T &  t,
const FTensor::Number< DIM > &   
)
Examples
matrix_function.cpp.

Definition at line 170 of file matrix_function.cpp.

170 {
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};
sdf.r
int r
Definition: sdf.py:8
t
constexpr double t
plate stiffness
Definition: plate.cpp:59

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 182 of file matrix_function.cpp.

182 {
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}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
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
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
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
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
get_norm_t4
auto get_norm_t4(T &t, const FTensor::Number< DIM > &)
Definition: matrix_function.cpp:170
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
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:109
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:298
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 180 of file matrix_function.cpp.

◆ i

FTensor::Index<'i', 3> i

Definition at line 18 of file matrix_function.cpp.

◆ j

FTensor::Index<'j', 3> j
Examples
ContactOps.hpp, ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, HookeElement.cpp, HookeElement.hpp, HookeInternalStressElement.hpp, MagneticElement.hpp, NavierStokesElement.cpp, PlasticOps.hpp, PlasticOpsGeneric.hpp, PlasticOpsLargeStrains.hpp, PlasticOpsSmallStrains.hpp, Remodeling.cpp, ThermoElasticOps.hpp, analytical_poisson_field_split.cpp, approx_sphere.cpp, bernstein_bezier_generate_base.cpp, cell_forces.cpp, dynamic_first_order_con_law.cpp, elasticity.cpp, free_surface.cpp, hcurl_curl_operator.cpp, higher_derivatives.cpp, level_set.cpp, lorentz_force.cpp, matrix_function.cpp, mortar_contact_thermal.cpp, plastic.cpp, plate.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, scalar_check_approximation.cpp, shallow_wave.cpp, simple_contact_thermal.cpp, simple_elasticity.cpp, and thermo_elastic.cpp.

Definition at line 19 of file matrix_function.cpp.

◆ k

FTensor::Index<'k', 3> k
Examples
ContactOps.hpp, ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, HookeElement.cpp, HookeElement.hpp, HookeInternalStressElement.hpp, MagneticElement.hpp, NavierStokesElement.cpp, PlasticOps.hpp, PlasticOpsGeneric.hpp, PlasticOpsLargeStrains.hpp, PlasticOpsSmallStrains.hpp, Remodeling.cpp, ThermoElasticOps.hpp, approx_sphere.cpp, bernstein_bezier_generate_base.cpp, dynamic_first_order_con_law.cpp, eigen_elastic.cpp, free_surface.cpp, hcurl_curl_operator.cpp, heat_equation.cpp, helmholtz.cpp, higher_derivatives.cpp, level_set.cpp, lorentz_force.cpp, matrix_function.cpp, phase.cpp, plastic.cpp, plate.cpp, prism_polynomial_approximation.cpp, scalar_check_approximation.cpp, simple_elasticity.cpp, and thermo_elastic.cpp.

Definition at line 20 of file matrix_function.cpp.

◆ l

FTensor::Index<'l', 3> l
Examples
ContactOps.hpp, ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, HookeElement.cpp, HookeElement.hpp, HookeInternalStressElement.hpp, PlasticOps.hpp, PlasticOpsGeneric.hpp, PlasticOpsLargeStrains.hpp, PlasticOpsSmallStrains.hpp, ThermoElasticOps.hpp, approx_sphere.cpp, dynamic_first_order_con_law.cpp, free_surface.cpp, hanging_node_approx.cpp, higher_derivatives.cpp, level_set.cpp, matrix_function.cpp, plastic.cpp, plate.cpp, shallow_wave.cpp, simple_elasticity.cpp, and thermo_elastic.cpp.

Definition at line 21 of file matrix_function.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow