v0.11.0
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
constexpr double eps
Definition: forces_and_sources_getting_mult_H1_H1_atom_test.cpp:30
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312

◆ 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
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 };
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
matrix_function.cpp.

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 
279  // Test matrix againsst mathematica results
280  {
281  FTensor::Tensor2<double, 3, 3> t_A{
282 
283  1., 0.1, -0.5,
284 
285  0.1, 2., 0.,
286 
287  -0.5, 0., 3.};
288 
289  FTensor::Tensor2<double, 3, 3> t_N{
290 
291  0.969837, -0.0860972, 0.228042,
292 
293  0.0790574, 0.996073, 0.0398449,
294 
295  -0.230577, -0.0206147, 0.972836};
296 
297  FTensor::Tensor1<double, 3> t_L{0.873555, 2.00794, 3.11851};
298 
299  FTensor::Tensor2<double, 3, 3> t_S{
300 
301  1., 0., 0.,
302 
303  0., 1., 0.,
304 
305  0., 0., 1.};
306 
307  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff A";
308  print_mat(t_A);
309 
310  // Testing against values in mathematica for 0,2 directive
311  {
312  auto f = [](double v) { return exp(v); };
313  auto d_f = [](double v) { return exp(v); };
314  auto dd_f = [](double v) { return exp(v); };
315 
316  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
317  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Reconstruct mat";
318  print_mat(t_b);
319 
320  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
321 
322  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff";
323  print_ddg_direction(t_d, 0, 2);
324 
325  auto t_dd = EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S, 3);
326 
327  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff Diff";
328  print_ddg_direction(t_dd, 0, 2);
329 
330  auto norm2_t_b = t_b(i, j) * t_b(i, j);
331  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_b " << norm2_t_b;
332 
333  auto norm2_t_d = get_norm_t4(t_d, FTensor::Number<3>());
334  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_d " << norm2_t_d;
335 
336  auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
337  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
338 
339  constexpr double regression_t_b = 572.543;
340  constexpr double regression_t_d = 859.939;
341  constexpr double regression_t_dd = 859.939;
342 
343  constexpr double eps = 1e-2;
344  if (std::abs(norm2_t_b - regression_t_b) > eps)
345  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_b");
346  if (std::abs(norm2_t_d - regression_t_d) > eps)
347  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_d");
348  if (std::abs(norm2_t_dd - regression_t_dd) > eps)
349  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_dd");
350  }
351 
352  // Comparing with lapack calculated eigen values. Note resulst should be
353  // invarinat to the direction of eiegn vector. Eigen vector can be
354  // multiplied by -1 and result should be unchanged
355  {
356 
357  std::array<double, 9> a{1., 0.1, -0.5,
358 
359  0.1, 2., 0.,
360 
361  -0.5, 0., 3.};
362 
363  auto tuple = run_lapack(a);
364  auto &t_a = std::get<0>(tuple);
365  auto &t_eig_vec = std::get<1>(tuple);
366  auto &t_eig_vals = std::get<2>(tuple);
367 
368  auto t_eig_val_diff =
369  (t_eig_vals(i) - t_L(i)) * (t_eig_vals(i) - t_L(i));
370  MOFEM_LOG("ATOM_TEST", Sev::inform)
371  << "t_eig_val_diff " << t_eig_val_diff;
372 
373  auto f = [](double v) { return exp(v); };
374  auto d_f = [](double v) { return exp(v); };
375  auto dd_f = [](double v) { return exp(v); };
376 
377  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
378  auto t_c = EigenMatrix::getMat(t_eig_vals, t_eig_vec, f);
379  t_c(i, j) -= t_b(i, j);
380  print_mat(t_c);
381 
382  auto norm2_t_c = t_c(i, j) * t_c(i, j);
383  MOFEM_LOG("ATOM_TEST", Sev::inform)
384  << "Reconstruct mat difference with lapack eigen valyes and "
385  "vectors "
386  << norm2_t_c;
387 
388  constexpr double eps = 1e-8;
389  if (fabs(norm2_t_c) > eps)
390  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
391  "Matrix not reeconstructed");
392  }
393 
394  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
395  FTensor::Ddg<double, 3, 3> t_one;
396  t_one(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
397 
398  // Testsing linear function second second direvarive zero
399  {
400  auto f = [](double v) { return v; };
401  auto d_f = [](double v) { return 1; };
402  auto dd_f = [](double v) { return 0; };
403 
404  constexpr double eps = 1e-10;
405  {
406  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
407  t_b(i, j) -= (t_A(i, j) || t_A(j, i)) / 2;
408  auto norm2_t_b = t_b(i, j) * t_b(i, j);
409  MOFEM_LOG("ATOM_TEST", Sev::inform)
410  << "Result should be matrix itself " << norm2_t_b;
411  if (norm2_t_b > eps)
412  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
413  "This norm should be zero");
414  }
415 
416  {
417 
418  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
419  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
420 
421  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
422  print_ddg(t_d_a, "hand ");
423  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
424  print_ddg(t_d, "code ");
425 
426  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
427  MOFEM_LOG("ATOM_TEST", Sev::inform)
428  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
429  if (nrm2_t_d_a > eps)
430  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
431  "This norm should be zero");
432  }
433 
434  {
435  FTensor::Tensor2<double, 3, 3> t_S{
436 
437  1., 0., 0.,
438 
439  0., 1., 0.,
440 
441  0., 0., 1.};
442 
443  auto t_dd =
444  EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S, 3);
445 
446  auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
447  MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
448  if (norm2_t_dd > eps)
449  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
450  "This norm should be zero");
451  }
452  }
453 
454  // Testsing quadratic function second second direvarive zero
455  {
456  auto f = [](double v) { return v * v; };
457  auto d_f = [](double v) { return 2 * v; };
458  auto dd_f = [](double v) { return 2; };
459 
460  constexpr double eps = 1e-9;
461 
462  // check if multiplication gives right value
463  {
464  auto t_b = EigenMatrix::getMat(t_L, t_N, f);
465  FTensor::Tensor2<double, 3, 3> t_a;
466  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
467  print_mat(t_a);
468  auto norm2_t_a = t_a(i, j) * t_a(i, j);
469  MOFEM_LOG("ATOM_TEST", Sev::inform)
470  << "Result should be matrix times matrix " << norm2_t_a;
471  if (norm2_t_a > eps)
472  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
473  "This norm should be zero");
474  }
475 
476  // check first directive
477  {
478  auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
479  print_ddg_direction(t_d, 0, 2);
480  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<3>());
481  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
482  print_ddg(t_d_a, "hand ");
483  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
484  print_ddg(t_d, "code ");
485  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
486  MOFEM_LOG("ATOM_TEST", Sev::inform)
487  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
488  if (nrm2_t_d_a > eps)
489  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
490  "This norm should be zero");
491  }
492 
493  // check second directive
494  {
495  FTensor::Tensor2<double, 3, 3> t_S{
496 
497  1., 1. / 2., 1. / 3.,
498 
499  2. / 2., 1., 2. / 3.,
500 
501  3. / 2., 1., 3. / 3.};
502 
503  auto t_dd =
504  EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S, 3);
505  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<3>());
506 
507  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
508  print_ddg(t_dd_a, "hand ");
509  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
510  print_ddg(t_dd, "code ");
511 
512  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
513  MOFEM_LOG("ATOM_TEST", Sev::inform)
514  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
515  if (nrm2_t_dd_a > eps)
516  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
517  "This norm should be zero");
518  }
519  }
520  }
521 
522  // Testing two same eigen values
523  {
524 
525  std::array<double, 9> a{5., 4., 0,
526 
527  4., 5, 0.,
528 
529  0.0, 0., 9};
530 
531  auto tuple = run_lapack(a, swap01);
532  auto &t_a = std::get<0>(tuple);
533  auto &t_eig_vecs = std::get<1>(tuple);
534  auto &t_eig_vals = std::get<2>(tuple);
535 
536  auto f = [](double v) { return v; };
537  auto d_f = [](double v) { return 1; };
538  auto dd_f = [](double v) { return 0; };
539 
540  constexpr double eps = 1e-10;
541  {
542  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
543  t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
544  auto norm2_t_b = t_b(i, j) * t_b(i, j);
545  MOFEM_LOG("ATOM_TEST", Sev::inform)
546  << "Result should be matrix itself " << norm2_t_b;
547  if (norm2_t_b > eps)
548  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
549  "This norm should be zero");
550  }
551 
552  {
553  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
554  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
555  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
556  print_ddg(t_d_a, "hand ");
557  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
558  print_ddg(t_d, "code ");
559  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
560  MOFEM_LOG("ATOM_TEST", Sev::inform)
561  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
562  if (nrm2_t_d_a > eps)
563  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
564  "This norm should be zero");
565  }
566 
567  {
568  auto f = [](double v) { return v * v; };
569  auto d_f = [](double v) { return 2 * v; };
570  auto dd_f = [](double v) { return 2; };
571  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
572  auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
573  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
574  print_ddg(t_d_a, "hand ");
575  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
576  print_ddg(t_d, "code ");
577  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
578  MOFEM_LOG("ATOM_TEST", Sev::inform)
579  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
580  if (nrm2_t_d_a > eps)
581  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
582  "This norm should be zero");
583  }
584  }
585 
586  // Testing three same eigen values
587  {
588 
589  std::array<double, 9> a{4., 0., 0,
590 
591  0., 4., 0.,
592 
593  0.0, 0., 4.};
594 
595  auto f = [](double v) { return v; };
596  auto d_f = [](double v) { return 1; };
597  auto dd_f = [](double v) { return 0; };
598 
599  auto tuple = run_lapack(a);
600  auto &t_a = std::get<0>(tuple);
601  auto &t_eig_vecs = std::get<1>(tuple);
602  auto &t_eig_vals = std::get<2>(tuple);
603 
604  constexpr double eps = 1e-10;
605  {
606  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
607  t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
608  auto norm2_t_b = t_b(i, j) * t_b(i, j);
609  MOFEM_LOG("ATOM_TEST", Sev::inform)
610  << "Result should be matrix itself " << norm2_t_b;
611  if (norm2_t_b > eps)
612  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
613  "This norm should be zero");
614  }
615 
616  {
617  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
618  auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
619  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
620  print_ddg(t_d_a, "hand ");
621  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
622  print_ddg(t_d, "code ");
623  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
624  MOFEM_LOG("ATOM_TEST", Sev::inform)
625  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
626  if (nrm2_t_d_a > eps)
627  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
628  "This norm should be zero");
629  }
630 
631  {
632  auto f = [](double v) { return v * v; };
633  auto d_f = [](double v) { return 2 * v; };
634  auto dd_f = [](double v) { return 2; };
635  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
636  auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
637  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
638  print_ddg(t_d_a, "hand ");
639  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
640  print_ddg(t_d, "code ");
641  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
642  MOFEM_LOG("ATOM_TEST", Sev::inform)
643  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
644  if (nrm2_t_d_a > eps)
645  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
646  "This norm should be zero");
647  }
648  }
649 
650  // check second directive
651  {
652 
653  std::array<double, 9> a{0.1, 0., 0.,
654 
655  0., 0.1, 0.,
656 
657  0., 0., 0.1};
658 
659  auto tuple = run_lapack(a);
660  auto &t_a = std::get<0>(tuple);
661  auto &t_eig_vecs = std::get<1>(tuple);
662  auto &t_eig_vals = std::get<2>(tuple);
663 
664  t_eig_vals(0) -= 1e-4;
665  t_eig_vals(2) += 1e-4;
666 
667  constexpr double eps = 1e-10;
668 
669  auto f = [](double v) { return v; };
670  auto d_f = [](double v) { return 1; };
671  auto dd_f = [](double v) { return 0; };
672 
673  FTensor::Tensor2<double, 3, 3> t_S{
674 
675  1., 1. / 2., 1. / 3.,
676 
677  2. / 2., 1., 2. / 3.,
678 
679  3. / 2., 1., 3. / 3.};
680 
681  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
682  dd_f, t_S, 1);
683 
684  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
685  print_ddg(t_dd, "test ");
686 
687  double nrm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
688  MOFEM_LOG("ATOM_TEST", Sev::inform)
689  << "Direvarive hand calculation minus code " << nrm2_t_dd;
690  if (nrm2_t_dd > eps)
691  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
692  "This norm should be zero");
693  }
694 
695  // check second directive
696  {
697 
698  std::array<double, 9> a{2, 0., 0.,
699 
700  0., 2, 0.,
701 
702  0., 0., 2};
703 
704  auto tuple = run_lapack(a);
705  auto &t_a = std::get<0>(tuple);
706  auto &t_eig_vecs = std::get<1>(tuple);
707  auto &t_eig_vals = std::get<2>(tuple);
708 
709  constexpr double eps = 1e-10;
710 
711  auto f = [](double v) { return v * v; };
712  auto d_f = [](double v) { return 2 * v; };
713  auto dd_f = [](double v) { return 2; };
714  FTensor::Tensor2<double, 3, 3> t_S{
715 
716  1., 1. / 2., 1. / 3.,
717 
718  2. / 1., 1., 2. / 3.,
719 
720  3. / 1., 3. / 1., 1.};
721 
722  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
723  dd_f, t_S, 1);
724  // print_ddg(t_dd, "test ");
725 
726  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<3>());
727 
728  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
729  MOFEM_LOG("ATOM_TEST", Sev::inform)
730  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
731  if (nrm2_t_dd_a > eps)
732  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
733  "This norm should be zero");
734  }
735 
736  // check second directive two reapeating eiegn values
737  {
738 
739  std::array<double, 9> a{5., 4., 0.,
740 
741  4., 5., 0.,
742 
743  0., 0., 9};
744 
745  auto tuple = run_lapack(a, swap01);
746  auto &t_a = std::get<0>(tuple);
747  auto &t_eig_vecs = std::get<1>(tuple);
748  auto &t_eig_vals = std::get<2>(tuple);
749 
750  constexpr double eps = 1e-10;
751 
752  auto f = [](double v) { return v * v; };
753  auto d_f = [](double v) { return 2 * v; };
754  auto dd_f = [](double v) { return 2; };
755 
756  FTensor::Tensor2<double, 3, 3> t_S{
757 
758  1., 1. / 2., 1. / 3.,
759 
760  2. / 1., 1., 2. / 3.,
761 
762  3. / 1., 3. / 1., 1.};
763 
764  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
765  dd_f, t_S, 2);
766  print_ddg(t_dd, "test ");
767 
768  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<3>());
769 
770  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
771  MOFEM_LOG("ATOM_TEST", Sev::inform)
772  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
773  if (nrm2_t_dd_a > eps)
774  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
775  "This norm should be zero");
776  }
777 
778  // check second directive exponent
779  {
780 
781  std::array<double, 9> a{2, 0., 0.,
782 
783  0., 2, 0.,
784 
785  0., 0., 2};
786 
787  auto tuple = run_lapack(a);
788  auto &t_a = std::get<0>(tuple);
789  auto &t_eig_vecs = std::get<1>(tuple);
790  auto &t_eig_vals = std::get<2>(tuple);
791 
792  t_eig_vals(0) -= 1e-5;
793  t_eig_vals(2) += 1e-5;
794 
795  constexpr double eps = 1e-7;
796 
797  auto f = [](double v) { return exp(v); };
798  auto d_f = [](double v) { return exp(v); };
799  auto dd_f = [](double v) { return exp(v); };
800  FTensor::Tensor2<double, 3, 3> t_S{
801 
802  1., 1. / 2., 1. / 3.,
803 
804  2. / 1., 1., 2. / 3.,
805 
806  3. / 1., 3. / 1., 1.};
807 
808  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
809  dd_f, t_S, 3);
810  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
811  dd_f, t_S, 1);
812 
813  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
814  MOFEM_LOG("ATOM_TEST", Sev::verbose)
815  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
816 
817  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
818  MOFEM_LOG("ATOM_TEST", Sev::verbose)
819  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
820 
821  print_ddg(t_dd_1, "t_dd_1 ");
822  print_ddg(t_dd_2, "t_dd_2 ");
823 
824  FTensor::Ddg<double, 3, 3> t_dd_3;
825  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
826 
827  for (int ii = 0; ii != 3; ++ii)
828  for (int jj = 0; jj != 3; ++jj)
829  for (int kk = 0; kk != 3; ++kk)
830  for (int ll = 0; ll != 3; ++ll) {
831  constexpr double eps = 1e-4;
832  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
833  MOFEM_LOG("ATOM_TEST", Sev::error)
834  << "Error " << ii << " " << jj << " " << kk << " " << ll
835  << " " << t_dd_1(ii, jj, kk, ll) << " "
836  << t_dd_2(ii, jj, kk, ll);
837  }
838 
839  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
840  MOFEM_LOG("ATOM_TEST", Sev::inform)
841  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
842  if (nrm2_t_dd_3 > eps)
843  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
844  "This norm should be zero");
845  }
846 
847  // check second directive exponent agains perturned
848  {
849 
850  std::array<double, 9> a{5., 4., 0.,
851 
852  4., 5., 0.,
853 
854  0., 0., 9};
855 
856  auto tuple = run_lapack(a, swap01);
857  auto &t_a = std::get<0>(tuple);
858  auto &t_eig_vecs = std::get<1>(tuple);
859  auto &t_eig_vals = std::get<2>(tuple);
860 
861  t_eig_vals(0) -= 1e-4;
862  t_eig_vals(2) += 1e-4;
863 
864  constexpr double eps = 1e-4;
865 
866  auto f = [](double v) { return v * v; };
867  auto d_f = [](double v) { return 2 * v; };
868  auto dd_f = [](double v) { return 2; };
869 
870  FTensor::Tensor2<double, 3, 3> t_S{
871 
872  1., 1. / 2., 1. / 3.,
873 
874  2. / 1., 1., 2. / 3.,
875 
876  3. / 1., 3. / 1., 1.};
877 
878  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
879  dd_f, t_S, 3);
880  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
881  dd_f, t_S, 2);
882 
883  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
884  MOFEM_LOG("ATOM_TEST", Sev::verbose)
885  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
886 
887  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
888  MOFEM_LOG("ATOM_TEST", Sev::verbose)
889  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
890 
891  print_ddg(t_dd_1, "t_dd_1 ");
892  print_ddg(t_dd_2, "t_dd_2 ");
893 
894  FTensor::Ddg<double, 3, 3> t_dd_3;
895  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
896 
897  for (int ii = 0; ii != 3; ++ii)
898  for (int jj = 0; jj != 3; ++jj)
899  for (int kk = 0; kk != 3; ++kk)
900  for (int ll = 0; ll != 3; ++ll) {
901  constexpr double eps = 1e-3;
902  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
903  MOFEM_LOG("ATOM_TEST", Sev::error)
904  << "Error " << ii << " " << jj << " " << kk << " " << ll
905  << " " << t_dd_1(ii, jj, kk, ll) << " "
906  << t_dd_2(ii, jj, kk, ll);
907  }
908 
909  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
910  MOFEM_LOG("ATOM_TEST", Sev::inform)
911  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
912  if (nrm2_t_dd_3 > eps)
913  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
914  "This norm should be zero");
915  }
916 
917  // check second directive exponent
918  {
919 
920  std::array<double, 9> a{5., 4., 0.,
921 
922  4., 5., 0.,
923 
924  0., 0., 9};
925 
926  auto tuple = run_lapack(a, swap01);
927  auto &t_a = std::get<0>(tuple);
928  auto &t_eig_vecs = std::get<1>(tuple);
929  auto &t_eig_vals = std::get<2>(tuple);
930 
931  constexpr double eps = 1e-4;
932  constexpr int p = 3;
933 
934  auto f = [](double v) { return pow(v, p); };
935  auto d_f = [](double v) { return p * pow(v, p - 1); };
936  auto dd_f = [](double v) {
937  return p * (p - 1) * pow(v, std::max(0, p - 2));
938  };
939 
940  FTensor::Tensor2<double, 3, 3> t_S{
941 
942  1., 1. / 2., 1. / 3.,
943 
944  2. / 1., 1., 2. / 3.,
945 
946  3. / 1., 3. / 1., 1.};
947 
948  t_eig_vals(0) += 2e-5;
949  t_eig_vals(2) -= 2e-5;
950  auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
951  dd_f, t_S, 3);
952  auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
953  dd_f, t_S, 2);
954 
955  double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
956  MOFEM_LOG("ATOM_TEST", Sev::verbose)
957  << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
958 
959  double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
960  MOFEM_LOG("ATOM_TEST", Sev::verbose)
961  << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
962 
963  print_ddg(t_dd_1, "t_dd_1 ");
964  print_ddg(t_dd_2, "t_dd_2 ");
965 
966  FTensor::Ddg<double, 3, 3> t_dd_3;
967  t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
968 
969  for (int ii = 0; ii != 3; ++ii)
970  for (int jj = 0; jj != 3; ++jj)
971  for (int kk = 0; kk != 3; ++kk)
972  for (int ll = 0; ll != 3; ++ll) {
973  constexpr double eps = 1e-3;
974  if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
975  MOFEM_LOG("ATOM_TEST", Sev::error)
976  << "Error " << ii << " " << jj << " " << kk << " " << ll
977  << " " << t_dd_1(ii, jj, kk, ll) << " "
978  << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
979  }
980 
981  double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
982  MOFEM_LOG("ATOM_TEST", Sev::inform)
983  << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
984  if (nrm2_t_dd_3 > eps)
985  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
986  "This norm should be zero");
987  }
988 
989  // Speed
990  {
991 
992  std::array<double, 9> a{1., 0.1, -0.5,
993 
994  0.1, 2., 0.,
995 
996  -0.5, 0., 3.};
997 
998  auto tuple = run_lapack(a);
999  auto &t_a = std::get<0>(tuple);
1000  auto &t_eig_vecs = std::get<1>(tuple);
1001  auto &t_eig_vals = std::get<2>(tuple);
1002 
1003  auto f = [](double v) { return exp(v); };
1004  auto d_f = [](double v) { return exp(v); };
1005  auto dd_f = [](double v) { return exp(v); };
1006 
1007  FTensor::Tensor2<double, 3, 3> t_S{
1008 
1009  1., 1. / 2., 1. / 3.,
1010 
1011  2. / 1., 1., 2. / 3.,
1012 
1013  3. / 1., 3. / 1., 1.};
1014 
1015  MOFEM_LOG("ATOM_TEST", Sev::inform) << "Start";
1016  for (int ii = 0; ii != 1000; ++ii) {
1017  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 3);
1018  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1019  dd_f, t_S, 3);
1020  }
1021  MOFEM_LOG("ATOM_TEST", Sev::inform) << "End";
1022  }
1023 
1024  // 2d case
1025 
1026  auto run_lapack_2d = [](auto &a) {
1027  int info;
1028  double wkopt;
1029  double w[2];
1030 
1031  FTensor::Tensor2<double, 2, 2> t_a{
1032 
1033  a[0], a[1],
1034 
1035  a[2], a[3]};
1036 
1037  /* Query and allocate the optimal workspace */
1038  int lwork = -1;
1039  info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &wkopt, lwork);
1040  if (info > 0)
1041  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1042  lwork = (int)wkopt;
1043  std::vector<double> work(lwork);
1044  /* Solve eigenproblem */
1045  info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &*work.begin(), lwork);
1046  if (info > 0)
1047  THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1048 
1049  FTensor::Tensor2<double, 2, 2> t_eig_vecs{
1050 
1051  a[0 * 2 + 0], a[0 * 2 + 1],
1052 
1053  a[1 * 2 + 0], a[1 * 2 + 1]};
1054 
1055  FTensor::Tensor1<double, 2> t_eig_vals{w[0], w[1]};
1056 
1057  return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1058  };
1059 
1060  // Testsing quadratic function for 2d
1061  {
1062 
1063  std::array<double, 9> a{1., 0.1,
1064 
1065  0.1, 2.};
1066 
1067  auto tuple = run_lapack_2d(a);
1068  auto &t_A = std::get<0>(tuple);
1069  auto &t_eig_vecs = std::get<1>(tuple);
1070  auto &t_eig_vals = std::get<2>(tuple);
1071 
1072  auto f = [](double v) { return v * v; };
1073  auto d_f = [](double v) { return 2 * v; };
1074  auto dd_f = [](double v) { return 2; };
1075 
1076  constexpr double eps = 1e-6;
1077 
1078  FTensor::Index<'i', 2> i;
1079  FTensor::Index<'j', 2> j;
1080  FTensor::Index<'k', 2> k;
1081  FTensor::Index<'l', 2> l;
1082 
1083  // check if multiplication gives right value
1084  {
1085  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1086  FTensor::Tensor2<double, 2, 2> t_a;
1087  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1088  print_mat(t_a);
1089  auto norm2_t_a = t_a(i, j) * t_a(i, j);
1090  MOFEM_LOG("ATOM_TEST", Sev::inform)
1091  << "Result should be matrix times matrix " << norm2_t_a;
1092  if (norm2_t_a > eps)
1093  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1094  "This norm should be zero");
1095  }
1096 
1097  // check first directive
1098  {
1099  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
1100  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1101  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1102  MOFEM_LOG("ATOM_TEST", Sev::inform)
1103  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1104  if (nrm2_t_d_a > eps)
1105  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1106  "This norm should be zero");
1107  }
1108 
1109  // check second directive
1110  {
1111  FTensor::Tensor2<double, 2, 2> t_S{
1112 
1113  1., 1. / 2,
1114 
1115  2. / 2., 1.};
1116 
1117  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1118  dd_f, t_S, 2);
1119  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<2>());
1120 
1121  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1122  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1123 
1124  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1125  MOFEM_LOG("ATOM_TEST", Sev::inform)
1126  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1127  if (nrm2_t_dd_a > eps)
1128  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1129  "This norm should be zero");
1130  }
1131  }
1132 
1133  // Testsing quadratic function for repeating eigen valsues
1134  {
1135 
1136  std::array<double, 9> a{2., 0,
1137 
1138  0, 2.};
1139 
1140  auto tuple = run_lapack_2d(a);
1141  auto &t_A = std::get<0>(tuple);
1142  auto &t_eig_vecs = std::get<1>(tuple);
1143  auto &t_eig_vals = std::get<2>(tuple);
1144 
1145  auto f = [](double v) { return v * v; };
1146  auto d_f = [](double v) { return 2 * v; };
1147  auto dd_f = [](double v) { return 2; };
1148 
1149  constexpr double eps = 1e-6;
1150 
1151  FTensor::Index<'i', 2> i;
1152  FTensor::Index<'j', 2> j;
1153  FTensor::Index<'k', 2> k;
1154  FTensor::Index<'l', 2> l;
1155 
1156  // check if multiplication gives right value
1157  {
1158  auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1159  FTensor::Tensor2<double, 2, 2> t_a;
1160  t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1161  print_mat(t_a);
1162  auto norm2_t_a = t_a(i, j) * t_a(i, j);
1163  MOFEM_LOG("ATOM_TEST", Sev::inform)
1164  << "Result should be matrix times matrix " << norm2_t_a;
1165  if (norm2_t_a > eps)
1166  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1167  "This norm should be zero");
1168  }
1169 
1170  // check first directive
1171  {
1172  auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
1173  auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1174  double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1175  MOFEM_LOG("ATOM_TEST", Sev::inform)
1176  << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1177  if (nrm2_t_d_a > eps)
1178  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1179  "This norm should be zero");
1180  }
1181 
1182  // check second directive
1183  {
1184  FTensor::Tensor2<double, 2, 2> t_S{
1185 
1186  1., 1. / 2,
1187 
1188  2. / 2., 1.};
1189 
1190  auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1191  dd_f, t_S, 1);
1192  auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<2>());
1193 
1194  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1195  MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1196 
1197  double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1198  MOFEM_LOG("ATOM_TEST", Sev::inform)
1199  << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1200  if (nrm2_t_dd_a > eps)
1201  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1202  "This norm should be zero");
1203  }
1204  }
1205 
1206  }
1207  CATCH_ERRORS;
1208 
1209  CHKERR MoFEM::Core::Finalize();
1210 }
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:86
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
THROW_MESSAGE
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
print_mat
void print_mat(double *M, int m, int n)
print matric M
MOFEM_LOG_ATTRIBUTES
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:300
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:273
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)
Get the Diff Diff Mat object.
Definition: MatrixFunction.cpp:78
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:21
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:20
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:144
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:60
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100

Variable Documentation

◆ help

char help[] = "...\n\n"
static
Examples
matrix_function.cpp.

Definition at line 180 of file matrix_function.cpp.

◆ i

FTensor::Index<'i', 3> i
Examples
ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, EshelbianPlasticity.cpp, HookeElement.cpp, HookeElement.hpp, MagneticElement.hpp, NavierStokesElement.cpp, NeoHookean.hpp, Remodeling.cpp, analytical_poisson_field_split.cpp, bernstein_bezier_generate_base.cpp, cell_forces.cpp, contact.cpp, continuity_check_on_skeleton_3d.cpp, continuity_check_on_skeleton_with_simple_2d_for_hcurl.cpp, continuity_check_on_skeleton_with_simple_2d_for_hdiv.cpp, dynamic_elastic.cpp, edge_and_bubble_shape_functions_on_quad.cpp, eigen_elastic.cpp, elasticity.cpp, field_blas_set_vertex_dofs.cpp, gauss_points_on_quad.cpp, hcurl_check_approx_in_2d.cpp, hcurl_curl_operator.cpp, hcurl_divergence_operator_2d.cpp, heat_method.cpp, lorentz_force.cpp, matrix_function.cpp, mixed_poisson.cpp, navier_stokes.cpp, plastic.cpp, prism_elements_from_surface.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, reaction_diffusion.cpp, scalar_check_approximation_2d.cpp, simple_contact.cpp, simple_elasticity.cpp, simple_interface.cpp, and testing_jacobian_of_hook_scaled_with_density_element.cpp.

Definition at line 18 of file matrix_function.cpp.

◆ j

FTensor::Index<'j', 3> j
Examples
ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, EshelbianPlasticity.cpp, HookeElement.cpp, HookeElement.hpp, NavierStokesElement.cpp, NeoHookean.hpp, Remodeling.cpp, analytical_poisson_field_split.cpp, bases_on_reference_rectangle.cpp, bernstein_bezier_generate_base.cpp, cell_forces.cpp, contact.cpp, dynamic_elastic.cpp, eigen_elastic.cpp, elasticity.cpp, hcurl_check_approx_in_2d.cpp, hcurl_curl_operator.cpp, heat_method.cpp, lorentz_force.cpp, matrix_function.cpp, plastic.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, scalar_check_approximation_2d.cpp, and simple_elasticity.cpp.

Definition at line 19 of file matrix_function.cpp.

◆ k

FTensor::Index<'k', 3> k
Examples
ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, EshelbianPlasticity.cpp, HookeElement.cpp, HookeElement.hpp, NavierStokesElement.cpp, Remodeling.cpp, bernstein_bezier_generate_base.cpp, contact.cpp, dynamic_elastic.cpp, eigen_elastic.cpp, helmholtz.cpp, lorentz_force.cpp, matrix_function.cpp, plastic.cpp, prism_polynomial_approximation.cpp, reaction_diffusion.cpp, and simple_elasticity.cpp.

Definition at line 20 of file matrix_function.cpp.

◆ l

FTensor::Index<'l', 3> l
Examples
ElasticityMixedFormulation.hpp, EshelbianOperators.cpp, HookeElement.cpp, HookeElement.hpp, contact.cpp, dynamic_elastic.cpp, eigen_elastic.cpp, matrix_function.cpp, plastic.cpp, and simple_elasticity.cpp.

Definition at line 21 of file matrix_function.cpp.

Generated on Sun Apr 11 2021 23:15:43 for MoFEM by Doxygen 1.9.1 and hosted at University of Glasgow