v0.14.0
MatrixFunctionTemplate.hpp
Go to the documentation of this file.
1 /** \file FunctionMatrix.hpp
2  \brief Get function from matrix
3  \ingroup ftensor
4 
5  For reference see \cite miehe2001algorithms
6 
7  Usage:
8 
9  To calculate exponent of matrix, first and second derivatives
10  \code
11  auto f = [](double v) { return exp(v); };
12  auto d_f = [](double v) { return exp(v); };
13  auto dd_f = [](double v) { return exp(v); };
14  \endcode
15 
16  Calculate matrix here t_L are vector of eigen values, and t_N is matrix of
17  eigen vectors.
18  \code
19  auto t_A = EigenMatrixImp<double, double, 3, 3>::getMat(t_L, t_N, f);
20  \endcode
21  Return t_A is symmetric tensor rank two.
22 
23  Calculate derivarive
24  \code
25  auto t_P = EigenMatrixImp<double, double, 3, 3>::getDiffMat(t_L, t_N, f, d_f);
26  \endcode
27  where return t_SL is 4th order tensor (symmetry on first two and
28  second to indices, i.e. minor symmetrise)
29 
30  Calculate second derivative, L, such that S:L, for given S,
31  \code
32  FTensor::Tensor2<double, 3, 3> t_S{
33 
34  1., 0., 0.,
35 
36  0., 1., 0.,
37 
38  0., 0., 1.};
39 
40  auto t_SL = EigenMatrixImp<double, double, 3, 3>::getDiffDiffMat(
41  t_L, t_N, f, d_f, dd_f, t_S)
42  \endcode
43  where return t_SL is 4th order tensor (symmetry on first two and
44  second to indices, i.e. minor symmetrise)
45 
46  You can calculate eigen values using lapack.
47 
48  *
49  */
50 
51 
52 
53 namespace EigenMatrix {
54 
55 template <int N> using Number = FTensor::Number<N>;
56 
57 template <int N1, int N2, int Dim>
58 inline auto get_sym_index(const Number<N1> &, const Number<N2> &,
59  const Number<Dim> &) {
60  if constexpr (N1 > N2)
61  return N1 + (N2 * (2 * Dim - N2 - 1)) / 2;
62  else
63  return N2 + (N1 * (2 * Dim - N1 - 1)) / 2;
64 }
65 
66 inline auto get_sym_index(const int N1, const int N2, const int Dim) {
67  if (N1 > N2)
68  return N1 + (N2 * (2 * Dim - N2 - 1)) / 2;
69  else
70  return N2 + (N1 * (2 * Dim - N1 - 1)) / 2;
71 }
72 
73 template <int N1, int N2, int Dim>
74 inline auto get_nodiag_index(const Number<N1> &, const Number<N2> &,
75  const Number<Dim> &) {
76  static_assert(N1 != N2, "Bad index");
77  if constexpr (N2 > N1)
78  return (Dim - 1) * N1 + N2 - 1;
79  else
80  return (Dim - 1) * N1 + N2;
81 }
82 
83 inline auto get_nodiag_index(const int N1, const int N2, int Dim) {
84  if (N2 > N1)
85  return (Dim - 1) * N1 + N2 - 1;
86  else
87  return (Dim - 1) * N1 + N2;
88 }
89 
90 template <typename E, typename C> struct d2MCoefficients {
91  using Val = typename E::Val;
92  using Vec = typename E::Vec;
93  using Fun = typename E::Fun;
94 
95  using NumberNb = typename E::NumberNb;
96  using NumberDim = typename E::NumberDim;
97 
98  template <int N> using Number = FTensor::Number<N>;
99 
101  E &e;
102 
103  template <int a, int b>
104  inline auto get(const Number<a> &, const Number<b> &, const int i,
105  const int j, const int k, const int l,
106  const Number<3> &) const {
107  return e.aS[get_sym_index(Number<a>(), Number<b>(), NumberDim())](i, j, k,
108  l) *
109  e.fVal(a) * e.aF(a, b);
110  }
111 
112  template <int a, int b>
113  inline auto get(const Number<a> &, const Number<b> &, const int i,
114  const int j, const int k, const int l,
115  const Number<2> &) const {
116  if constexpr (a == 1 || b == 1)
117  return get(Number<a>(), Number<b>(), i, j, k, l, Number<3>());
118  else
119  return get(Number<a>(), Number<b>(), i, j, k, l, Number<1>());
120  }
121 
122  template <int a, int b>
123  inline auto get(const Number<a> &, const Number<b> &, const int i,
124  const int j, const int k, const int l,
125  const Number<1>) const {
126  return e.aS[get_sym_index(Number<a>(), Number<b>(), NumberDim())](i, j, k,
127  l) *
128  e.dfVal(a) / static_cast<C>(2);
129  }
130 };
131 
132 template <typename E, typename C, typename G> struct d2MImpl {
133  using Val = typename E::Val;
134  using Vec = typename E::Vec;
135  using Fun = typename E::Fun;
136 
137  template <int N> using Number = FTensor::Number<N>;
138  d2MImpl(E &e) : g(e), e(e) {}
139  G g;
140  E &e;
141 
142  template <int b, int a>
143  inline C term(const int i, const int j, const int k, const int l) const {
144  if constexpr (a != b) {
145  return g.get(Number<a>(), Number<b>(), i, j, k, l,
146  typename E::NumberNb());
147  }
148  return 0;
149  }
150 
151  template <int nb, int a>
152  inline C eval(const Number<nb> &, const Number<a> &, const int i, const int j,
153  const int k, const int l) const {
154  return term<nb - 1, a>(i, j, k, l) +
155  eval(Number<nb - 1>(), Number<a>(), i, j, k, l);
156  }
157 
158  template <int a>
159  inline C eval(const Number<1> &, const Number<a> &, const int i, const int j,
160  const int k, const int l) const {
161  return term<0, a>(i, j, k, l);
162  }
163 };
164 
165 template <typename E, typename C> struct Fdd4MImpl {
166  using Val = typename E::Val;
167  using Vec = typename E::Vec;
168  using Fun = typename E::Fun;
169 
170  Fdd4MImpl(E &e) : e(e) {}
171  E &e;
172 
173  using NumberNb = typename E::NumberNb;
174  using NumberDim = typename E::NumberDim;
175 
176 
177  template <int A, int a, int b, int I, int J, int K, int L>
178  inline auto fd2M() const {
179  return e.d2MType1[b][get_sym_index(Number<A>(), Number<a>(), NumberDim())](
181  }
182 
183  template <int a, int b, int i, int j, int k, int l, int m, int n>
184  inline auto term_fd2S(const Number<a> &, const Number<b> &, const Number<i> &,
185  const Number<j> &, const Number<k> &, const Number<l> &,
186  const Number<m> &, const Number<n> &) const {
187  if constexpr (i == j && k == l) {
188  return 4 *
189  (
190 
191  fd2M<a, a, b, i, k, m, n>() * e.aM[b](Number<j>(), Number<l>())
192 
193  +
194 
195  fd2M<b, a, b, i, k, m, n>() * e.aM[a](Number<j>(), Number<l>())
196 
197  );
198 
199  } else if constexpr (i == j)
200  return 2 *
201  (
202 
203  fd2M<a, a, b, i, k, m, n>() *
204  e.aM[b](Number<j>(), Number<l>()) +
205  fd2M<b, a, b, j, l, m, n>() * e.aM[a](Number<i>(), Number<k>())
206 
207  +
208 
209  fd2M<b, a, b, i, k, m, n>() *
210  e.aM[a](Number<j>(), Number<l>()) +
211  fd2M<a, a, b, j, l, m, n>() * e.aM[b](Number<i>(), Number<k>())
212 
213  );
214  else if constexpr (k == l)
215  return 2 *
216  (
217 
218  fd2M<a, a, b, i, k, m, n>() *
219  e.aM[b](Number<j>(), Number<l>()) +
220  fd2M<b, a, b, j, l, m, n>() *
221  e.aM[a](Number<i>(), Number<k>()) +
222 
223  +
224 
225  fd2M<b, a, b, i, k, m, n>() *
226  e.aM[a](Number<j>(), Number<l>()) +
227  fd2M<a, a, b, j, l, m, n>() * e.aM[b](Number<i>(), Number<k>())
228 
229  );
230  else
231  return fd2M<a, a, b, i, k, m, n>() * e.aM[b](Number<j>(), Number<l>()) +
232  fd2M<b, a, b, j, l, m, n>() * e.aM[a](Number<i>(), Number<k>()) +
233  fd2M<a, a, b, i, l, m, n>() * e.aM[b](Number<j>(), Number<k>()) +
234  fd2M<b, a, b, j, k, m, n>() * e.aM[a](Number<i>(), Number<l>())
235 
236  +
237 
238  fd2M<b, a, b, i, k, m, n>() * e.aM[a](Number<j>(), Number<l>()) +
239  fd2M<a, a, b, j, l, m, n>() * e.aM[b](Number<i>(), Number<k>()) +
240  fd2M<b, a, b, i, l, m, n>() * e.aM[a](Number<j>(), Number<k>()) +
241  fd2M<a, a, b, j, k, m, n>() * e.aM[b](Number<i>(), Number<l>());
242  }
243 
244  template <int NB, int a, int b, int i, int j, int k, int l, int m, int n>
245  inline C term_SM(const Number<a> &, const Number<b> &, const Number<NB> &,
246  const Number<i> &, const Number<j> &, const Number<k> &,
247  const Number<l> &, const Number<m> &,
248  const Number<n> &) const {
249 
250  if constexpr (NB == 1) {
251  return 0;
252 
253  } else if constexpr (NB == 2) {
254 
255  if constexpr (a == 1 || b == 1) {
256  return
257 
258  e.aF2(Number<a>(), Number<b>()) *
261  Number<i>(), Number<j>(), Number<k>(), Number<l>()) -
264  Number<i>(), Number<j>(), Number<k>(), Number<l>()));
265 
266  } else {
267  return 0;
268  }
269 
270  } else {
271 
272  return
273 
274  e.aF2(Number<a>(), Number<b>()) *
277  Number<i>(), Number<j>(), Number<k>(), Number<l>()) -
280  Number<i>(), Number<j>(), Number<k>(), Number<l>()));
281  }
282 
283  return 0;
284  }
285 
286  template <int nb, int a, int i, int j, int k, int l, int m, int n>
287  inline C eval_fdS2(const Number<nb> &, const Number<a> &, const Number<i> &,
288  const Number<j> &, const Number<k> &, const Number<l> &,
289  const Number<m> &, const Number<n> &) const {
290  if constexpr (a != nb - 1)
292  Number<k>(), Number<l>(), Number<m>(), Number<n>()) +
295  else
298  }
299 
300  template <int a, int i, int j, int k, int l, int m, int n>
301  inline C eval_fdS2(const Number<1> &, const Number<a> &, const Number<i> &,
302  const Number<j> &, const Number<k> &, const Number<l> &,
303  const Number<m> &, const Number<n> &) const {
304  if constexpr (a != 0)
305  return term_fd2S(Number<a>(), Number<0>(), Number<i>(), Number<j>(),
307  else
308  return 0;
309  }
310 
311  template <int nb, int a, int i, int j, int k, int l, int m, int n>
312  inline C eval_SM(const Number<nb> &, const Number<a> &, const Number<i> &,
313  const Number<j> &, const Number<k> &, const Number<l> &,
314  const Number<m> &, const Number<n> &) const {
315  if constexpr (a != nb - 1)
318  Number<n>()) +
321 
322  else
325  }
326 
327  template <int a, int i, int j, int k, int l, int m, int n>
328  inline C eval_SM(const Number<1> &, const Number<a> &, const Number<i> &,
329  const Number<j> &, const Number<k> &, const Number<l> &,
330  const Number<m> &, const Number<n> &) const {
331  if constexpr (a != 0)
332  return term_SM(Number<a>(), Number<0>(), NumberNb(), Number<i>(),
334  Number<n>());
335  else
336  return 0;
337  }
338 
339  template <int nb, int a, int i, int j, int k, int l, int m, int n>
340  inline C eval(const Number<nb> &, const Number<a> &, const Number<i> &,
341  const Number<j> &, const Number<k> &, const Number<l> &,
342  const Number<m> &, const Number<n> &) const {
343  return
344 
345  (2 * e.fVal(a)) * eval_SM(NumberDim(), Number<a>(), Number<i>(),
346  Number<j>(), Number<k>(), Number<l>(),
347  Number<m>(), Number<n>())
348 
349  +
350 
353  }
354 };
355 
356 template <typename E, typename C> struct ReconstructMatImpl {
357  using Val = typename E::Val;
358  using Vec = typename E::Vec;
359  using Fun = typename E::Fun;
360 
361  template <int N> using Number = FTensor::Number<N>;
362 
364  E &e;
365 
366  template <int a, int i, int j> inline C term() const {
367  return e.aM[a](Number<i>(), Number<j>()) * e.fVal(a);
368  }
369 
370  template <int nb, int i, int j>
371  inline C eval(const Number<nb> &, const Number<i> &,
372  const Number<j> &) const {
373  return term<nb - 1, i, j>() +
375  }
376 
377  template <int i, int j>
378  inline C eval(const Number<1> &, const Number<i> &, const Number<j> &) const {
379  return term<0, i, j>();
380  }
381 };
382 
383 template <typename E, typename C> struct FirstMatrixDirectiveImpl {
384  using Val = typename E::Val;
385  using Vec = typename E::Vec;
386  using Fun = typename E::Fun;
387 
388  template <int N> using Number = FTensor::Number<N>;
389 
392  E &e;
393 
394  template <int a, int i, int j, int k, int l> inline C term() const {
395  return
396 
397  e.aMM[a][a](i, j, k, l) * e.dfVal(a)
398 
399  +
400 
401  r.eval(typename E::NumberDim(), Number<a>(), i, j, k, l) * 0.5;
402  }
403 
404  template <int nb, int i, int j, int k, int l>
405  inline C eval(const Number<nb> &, const Number<i> &, const Number<j> &,
406  const Number<k> &, const Number<l> &) const {
407  return term<nb - 1, i, j, k, l>() + eval(Number<nb - 1>(), Number<i>(),
408  Number<j>(), Number<k>(),
409  Number<l>());
410  }
411 
412  template <int i, int j, int k, int l>
413  inline C eval(const Number<1> &, const Number<i> &, const Number<j> &,
414  const Number<k> &, const Number<l> &) const {
415  return term<0, i, j, k, l>();
416  }
417 };
418 
419 template <typename E, typename C> struct SecondMatrixDirectiveImpl {
420  using Val = typename E::Val;
421  using Vec = typename E::Vec;
422  using Fun = typename E::Fun;
423 
424  using NumberDim = typename E::NumberDim;
425 
426  template <int N> using Number = FTensor::Number<N>;
427 
430  E &e;
431 
432  template <int a, int i, int j, int k, int l, int m, int n>
433  inline C term1() const {
434  return e.d2MType0[a][get_sym_index(Number<k>(), Number<l>(), NumberDim())](
436  };
437 
438  template <int a, int i, int j, int k, int l, int m, int n>
439  inline C term2() const {
440  return e.d2MType0[a][get_sym_index(Number<i>(), Number<j>(), NumberDim())](
442  }
443 
444  template <int a, int i, int j, int k, int l, int m, int n>
445  inline C term3() const {
446  return e.d2MType0[a][get_sym_index(Number<n>(), Number<m>(), NumberDim())](
448  }
449 
450  template <int a, int i, int j, int k, int l, int m, int n>
451  inline C term() const {
452 
453  return
454 
455  (term1<a, i, j, k, l, m, n>() + term2<a, i, j, k, l, m, n>() +
456  term3<a, i, j, k, l, m, n>()) *
457  0.5
458 
459  +
460 
461  (e.aMM[a][a](Number<i>(), Number<j>(), Number<k>(), Number<l>()) *
462  e.aM[a](Number<m>(), Number<n>())) *
463  e.ddfVal(a)
464 
465  +
466 
467  r.eval(typename E::NumberDim(), Number<a>(), Number<i>(), Number<j>(),
468  Number<k>(), Number<l>(), Number<m>(), Number<n>()) *
469  0.25;
470  }
471 
472  template <int nb, int i, int j, int k, int l, int m, int n>
473  inline C eval(const Number<nb> &, const Number<i> &, const Number<j> &,
474  const Number<k> &, const Number<l> &, const Number<m> &,
475  const Number<n> &) const {
476  return term<nb - 1, i, j, k, l, m, n>()
477 
478  +
479 
481  Number<l>(), Number<m>(), Number<n>());
482  }
483 
484  template <int i, int j, int k, int l, int m, int n>
485  inline C eval(const Number<1> &, const Number<i> &, const Number<j> &,
486  const Number<k> &, const Number<l> &, const Number<m> &,
487  const Number<n> &) const {
488  return term<0, i, j, k, l, m, n>();
489  }
490 };
491 
492 template <typename E, typename C, typename T> struct GetMatImpl {
493  using Val = typename E::Val;
494  using Vec = typename E::Vec;
495  using Fun = typename E::Fun;
496 
497  using NumberNb = typename E::NumberNb;
498  using NumberDim = typename E::NumberDim;
499 
500  template <int N> using Number = FTensor::Number<N>;
501 
502  GetMatImpl(E &e, T &t_a) : r(e), tA(t_a) {}
504  T &tA;
505 
506  template <int i, int j>
507  inline void set(const Number<i> &, const Number<j> &) {
508 
510 
512  r.eval(NumberDim(), Number<i - 1>(), Number<j - 1>());
513  }
514 
515  template <int i> inline void set(const Number<i> &, const Number<1> &) {
516 
518 
519  tA(Number<i - 1>(), Number<0>()) =
520  r.eval(NumberDim(), Number<i - 1>(), Number<0>());
521  }
522 
523  inline void set(const Number<1> &, const Number<1> &) {
524  tA(Number<0>(), Number<0>()) =
525  r.eval(NumberDim(), Number<0>(), Number<0>());
526  }
527 };
528 
529 template <typename E, typename C, typename T> struct GetDiffMatImpl {
530  using Val = typename E::Val;
531  using Vec = typename E::Vec;
532  using Fun = typename E::Fun;
533 
534  using NumberNb = typename E::NumberNb;
535  using NumberDim = typename E::NumberDim;
536 
537  template <int N> using Number = FTensor::Number<N>;
538 
539  GetDiffMatImpl(E &e, T &t_a) : r(e), tA(t_a) {}
541  T &tA;
542 
543  template <int i, int j, int k, int l>
544  inline void set(const Number<i> &, const Number<j> &, const Number<k> &,
545  const Number<l> &) {
546 
549  Number<l - 1>());
550 
552  }
553 
554  template <int j, int k, int l>
555  inline void set(const Number<1> &, const Number<j> &, const Number<k> &,
556  const Number<l> &) {
557 
560  Number<l - 1>());
561 
563  }
564 
565  template <int k, int l>
566  inline void set(const Number<1> &, const Number<1> &, const Number<k> &,
567  const Number<l> &) {
568 
570  r.eval(NumberDim(), Number<0>(), Number<0>(), Number<k - 1>(),
571  Number<l - 1>());
572 
574  }
575 
576  template <int l>
577  inline void set(const Number<1> &, const Number<1> &, const Number<1> &,
578  const Number<l> &) {
579 
581  r.eval(NumberDim(), Number<0>(), Number<0>(), Number<0>(),
582  Number<l - 1>());
583 
585  }
586 
587  inline void set(const Number<1> &, const Number<1> &, const Number<1> &,
588  const Number<1> &) {
589 
591  r.eval(NumberDim(), Number<0>(), Number<0>(), Number<0>(), Number<0>());
592  }
593 };
594 
595 template <typename E, typename C, typename T1, typename T2>
597  using Val = typename E::Val;
598  using Vec = typename E::Vec;
599  using Fun = typename E::Fun;
600 
601  using NumberNb = typename E::NumberNb;
602  using NumberDim = typename E::NumberDim;
603 
604  template <int N> using Number = FTensor::Number<N>;
605 
606  GetDiffDiffMatImpl(E &e, T1 &t_a, T2 &t_S) : r(e), e(e), tA(t_a), tS(t_S) {}
608  E &e;
609  T1 &tA;
610  T2 &tS;
611 
612  template <int I, int J, int K, int L, int M, int N>
613  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
614  const Number<L> &, const Number<M> &, const Number<N> &) {
615  if constexpr (N != M)
616  return (tS(M - 1, N - 1) + tS(N - 1, M - 1)) *
617  r.eval(NumberDim(), Number<M - 1>(), Number<N - 1>(),
619  Number<L - 1>())
620 
621  +
622 
624  Number<M>(), Number<N - 1>());
625  else
626  return tS(M - 1, N - 1) * r.eval(NumberDim(), Number<M - 1>(),
629  Number<L - 1>())
630 
631  +
632 
634  Number<M>(), Number<N - 1>());
635  }
636 
637  template <int I, int J, int K, int L, int M>
638  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
639  const Number<L> &, const Number<M> &, const Number<1> &) {
640  return (tS(M - 1, 0) + tS(0, M - 1)) *
641  r.eval(NumberDim(), Number<M - 1>(), Number<0>(),
643  Number<L - 1>())
644 
645  +
646 
649 
650  ;
651  }
652 
653  template <int I, int J, int K, int L>
654  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
655  const Number<L> &, const Number<1> &, const Number<1> &) {
656  return tS(0, 0) * r.eval(NumberDim(), Number<0>(), Number<0>(),
658  Number<L - 1>());
659  }
660 
661  template <int I, int J, int K, int L>
662  inline void set(const Number<I> &, const Number<J> &, const Number<K> &,
663  const Number<L> &) {
665  tA(I - 1, J - 1, K - 1, L - 1) = add(Number<I>(), Number<J>(), Number<K>(),
666  Number<L>(), NumberDim(), NumberDim());
667  // Major symmetry
668  if constexpr (K != I || L != J)
669  tA(K - 1, L - 1, I - 1, J - 1) = tA(I - 1, J - 1, K - 1, L - 1);
670  }
671 
672  template <int I, int J, int K>
673  inline void set(const Number<I> &, const Number<J> &, const Number<K> &,
674  const Number<0> &) {
676  }
677 
678  template <int I, int J>
679  inline void set(const Number<I> &, const Number<J> &, const Number<0> &,
680  const Number<0> &) {
682  }
683 
684  template <int I, int K>
685  inline void set(const Number<I> &, const Number<0> &, const Number<K> &,
686  const Number<0> &) {
688  }
689 
690  inline void set(const Number<0> &, const Number<0> &, const Number<0> &,
691  const Number<0> &) {}
692 };
693 
694 template <typename E, typename C, typename T1, typename VT2, int DimT2>
695 struct GetDiffDiffMatImpl<E, C, T1, FTensor::Tensor2_symmetric<VT2, DimT2>> {
696  using Val = typename E::Val;
697  using Vec = typename E::Vec;
698  using Fun = typename E::Fun;
699 
700  using NumberNb = typename E::NumberNb;
701  using NumberDim = typename E::NumberDim;
702 
703  template <int N> using Number = FTensor::Number<N>;
704 
706  : r(e), e(e), tA(t_a), tS(t_S) {}
708  E &e;
709  T1 &tA;
711 
712  template <int I, int J, int K, int L, int M, int N>
713  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
714  const Number<L> &, const Number<M> &, const Number<N> &) {
715 
716  if constexpr (N != M)
717  return (2 * tS(Number<M - 1>(), Number<N - 1>())) *
718  r.eval(NumberDim(), Number<M - 1>(), Number<N - 1>(),
720  Number<L - 1>())
721 
722  +
723 
725  Number<M>(), Number<N - 1>());
726  else
727  return tS(Number<M - 1>(), Number<N - 1>()) *
728  r.eval(NumberDim(), Number<M - 1>(), Number<N - 1>(),
730  Number<L - 1>())
731 
732  +
733 
735  Number<M>(), Number<N - 1>());
736  }
737 
738  template <int I, int J, int K, int L, int M>
739  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
740  const Number<L> &, const Number<M> &, const Number<1> &) {
741  return (2 * tS(Number<M - 1>(), Number<0>())) *
742  r.eval(NumberDim(), Number<M - 1>(), Number<0>(),
744  Number<L - 1>())
745 
746  +
747 
750  }
751 
752  template <int I, int J, int K, int L>
753  inline auto add(const Number<I> &, const Number<J> &, const Number<K> &,
754  const Number<L> &, const Number<1> &, const Number<1> &) {
755  return tS(Number<0>(), Number<0>()) *
756  r.eval(NumberDim(), Number<0>(), Number<0>(), Number<I - 1>(),
758  }
759 
760  template <int I, int J, int K, int L>
761  inline void set(const Number<I> &, const Number<J> &, const Number<K> &,
762  const Number<L> &) {
763 
765 
768  NumberDim());
769 
770  // Major symmetry
771  if constexpr (K != I || L != J)
774  Number<L - 1>());
775  }
776 
777  template <int I, int J, int K>
778  inline void set(const Number<I> &, const Number<J> &, const Number<K> &,
779  const Number<0> &) {
781  }
782 
783  template <int I, int J>
784  inline void set(const Number<I> &, const Number<J> &, const Number<0> &,
785  const Number<0> &) {
787  }
788 
789  template <int I, int K>
790  inline void set(const Number<I> &, const Number<0> &, const Number<K> &,
791  const Number<0> &) {
793  }
794 
795  inline void set(const Number<0> &, const Number<0> &, const Number<0> &,
796  const Number<0> &) {}
797 };
798 
799 template <typename T1, typename T2, int NB, int Dim> struct EigenMatrixImp {
800 
803  using Fun = boost::function<double(const double)>;
804  using V = double; // typename FTensor::promote<T1, T2>::V;
805 
806  template <int N> using Number = FTensor::Number<N>;
807  template <char c> using I = typename FTensor::Index<c, Dim>;
808 
811 
812  EigenMatrixImp(Val &t_val, Vec &t_vec) : tVal(t_val), tVec(t_vec) {
813 
818 
819  for (auto aa = 0; aa != Dim; ++aa) {
820  auto &M = aM[aa];
821  for (auto ii = 0; ii != Dim; ++ii)
822  for (auto jj = 0; jj <= ii; ++jj)
823  M(ii, jj) = tVec(aa, ii) * tVec(aa, jj);
824  }
825 
826  for (auto aa = 0; aa != Dim; ++aa) {
827  for (auto bb = 0; bb != Dim; ++bb) {
828  auto &Ma = aM[aa];
829  auto &Mb = aM[bb];
830  auto &MM = aMM[aa][bb];
831  MM(i, j, k, l) = Ma(i, j) * Mb(k, l);
832  }
833  }
834 
835  for (auto aa = 0; aa != Dim; ++aa) {
836  for (auto bb = 0; bb != Dim; ++bb) {
837  if (aa != bb) {
838  auto &MM = aMM[aa][bb];
839  auto &G = aG[aa][bb];
840  G(i, j, k, l) = MM(i, k, j, l) || MM(i, l, j, k);
841  }
842  }
843  }
844 
845  for (auto aa = 0; aa != Dim; ++aa) {
846  for (auto bb = 0; bb != Dim; ++bb) {
847  if (aa != bb) {
848  auto &Gab = aG[aa][bb];
849  auto &Gba = aG[bb][aa];
850  auto &S = aS[get_sym_index(aa, bb, Dim)];
851  S(i, j, k, l) = Gab(i, j, k, l) + Gba(i, j, k, l);
852  }
853  }
854  }
855  }
856 
857  /**
858  * @brief Get matrix
859  *
860  * \f[
861  * \mathbf{B} = f(\mathbf{A})
862  * \f]
863  *
864  * \f[
865  * B_{ij} = \sum_{a}^3 f(\lambda^a) n^a_i n^a_j
866  * \f]
867  * where \f$a\f$ is eigen value number.
868  *
869  * @param t_val eiegn values vector
870  * @param t_vec eigen vectors matrix
871  * @param f function
872  * @return auto function symmetric tensor rank two
873  */
874  inline auto getMat(Fun f) {
875 
876  for (auto aa = 0; aa != Dim; ++aa)
877  fVal(aa) = f(tVal(aa));
878 
879  using T3 =
881  T3 t_A;
883  .set(NumberDim(), NumberDim());
884  return t_A;
885  }
886 
887  /**
888  * @brief Get derivative of matrix
889  *
890  * \f[
891  * P_{ijkl} = \frac{\partial B_{ij}}{\partial A_{kl}}
892  * \f]
893  *
894  * @param t_val eiegn values vector
895  * @param t_vec eiegn vectors matrix
896  * @param f function
897  * @param d_f directive of function
898  * @return auto derivatives, forth order tensor with minor simetries
899  */
900  inline auto getDiffMat(Fun f, Fun d_f) {
901 
902  for (auto aa = 0; aa != Dim; ++aa)
903  fVal(aa) = f(tVal(aa));
904 
905  for (auto aa = 0; aa != Dim; ++aa)
906  dfVal(aa) = d_f(tVal(aa));
907 
908  for (auto aa = 0; aa != Dim; ++aa)
909  for (auto bb = 0; bb != aa; ++bb) {
910  aF(aa, bb) = 1 / (tVal(aa) - tVal(bb));
911  aF(bb, aa) = -aF(aa, bb);
912  aF2(aa, bb) = aF(aa, bb) * aF(aa, bb);
913  }
914 
915  using T3 = FTensor::Ddg<V, Dim, Dim>;
916  T3 t_diff_A;
918  .set(NumberDim(), NumberDim(), NumberDim(), NumberDim());
919  return t_diff_A;
920  }
921 
922  /**
923  * @brief Get second directive of matrix
924  *
925  * \f[
926  * LS_{klmn} =
927  * S_{ij} \frac{\partial^2 B_{ij}}{\partial A_{kl} \partial A_{mn} }
928  * \f]
929  *
930  * @tparam T
931  * @param t_val eigen values vector
932  * @param t_vec eigen vectors matrix
933  * @param f function
934  * @param d_f derivative of function
935  * @param dd_f second derivative of function
936  * @param t_S second rank tensor S
937  * @return auto second derivatives, forth order tensor with minor symmetries
938  */
939  template <typename T>
940  inline auto getDiffDiffMat(Fun f, Fun d_f, Fun dd_f, T &t_S) {
941 
942  for (auto aa = 0; aa != Dim; ++aa)
943  fVal(aa) = f(tVal(aa));
944 
945  for (auto aa = 0; aa != Dim; ++aa)
946  dfVal(aa) = d_f(tVal(aa));
947 
948  for (auto aa = 0; aa != Dim; ++aa)
949  ddfVal(aa) = dd_f(tVal(aa));
950 
951  for (auto aa = 0; aa != Dim; ++aa)
952  for (auto bb = 0; bb != aa; ++bb) {
953  aF(aa, bb) = 1 / (tVal(aa) - tVal(bb));
954  aF(bb, aa) = -aF(aa, bb);
955  aF2(aa, bb) = aF(aa, bb) * aF(aa, bb);
956  }
957 
962 
963  for (auto aa = 0; aa != Dim; ++aa) {
964  for (auto bb = 0; bb != Dim; ++bb) {
965  if (aa != bb) {
966  const auto &S = aS[get_sym_index(aa, bb, Dim)];
967  const auto &M = aM[aa];
968  auto &SMmn = aSM[get_nodiag_index(aa, bb, Dim)];
969  for (auto mm = 0; mm != Dim; ++mm) {
970  for (auto nn = mm; nn != Dim; ++nn) {
971  SMmn[get_sym_index(mm, nn, Dim)](i, j, k, l) =
972  S(i, j, k, l) * M(mm, nn);
973  }
974  }
975  }
976  }
977  }
978 
979  for (auto aa = 0; aa != Dim; ++aa) {
980  for (auto mm = 0; mm != Dim; ++mm) {
981  for (auto nn = mm; nn != Dim; ++nn) {
982  d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) = 0;
983  }
984  }
985  }
986 
987  if constexpr (NB == 3)
988  for (auto aa = 0; aa != Dim; ++aa) {
989  for (auto bb = 0; bb != Dim; ++bb) {
990  if (aa != bb) {
991  const V v = dfVal(aa) * aF(aa, bb);
992  for (auto mm = 0; mm != Dim; ++mm) {
993  for (auto nn = mm; nn != Dim; ++nn) {
994  d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
995  v * aSM[get_nodiag_index(aa, bb, Dim)]
996  [get_sym_index(mm, nn, Dim)](i, j, k, l);
997  }
998  }
999  }
1000  }
1001  }
1002 
1003  if constexpr (NB == 2)
1004  for (auto aa = 0; aa != Dim; ++aa) {
1005  for (auto bb = 0; bb != Dim; ++bb) {
1006  if (aa != bb) {
1007  V v;
1008  if (aa == 1 || bb == 1)
1009  v = dfVal(aa) * aF(aa, bb);
1010  else
1011  v = ddfVal(aa) / 2;
1012  for (auto mm = 0; mm != Dim; ++mm) {
1013  for (auto nn = mm; nn != Dim; ++nn) {
1014  d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
1015  v * aSM[get_nodiag_index(aa, bb, Dim)]
1016  [get_sym_index(mm, nn, Dim)](i, j, k, l);
1017  }
1018  }
1019  }
1020  }
1021  }
1022 
1023  if constexpr (NB == 1)
1024  for (auto aa = 0; aa != Dim; ++aa) {
1025  for (auto bb = 0; bb != Dim; ++bb) {
1026  if (aa != bb) {
1027  const V v = ddfVal(aa) / 2;
1028  for (auto mm = 0; mm != Dim; ++mm) {
1029  for (auto nn = mm; nn != Dim; ++nn) {
1030  d2MType0[aa][get_sym_index(mm, nn, Dim)](i, j, k, l) +=
1031  v * aSM[get_nodiag_index(aa, bb, Dim)]
1032  [get_sym_index(mm, nn, Dim)](i, j, k, l);
1033  }
1034  }
1035  }
1036  }
1037  }
1038 
1039  for (auto aa = 0; aa != Dim; ++aa) {
1040  for (auto mm = 0; mm != Dim; ++mm) {
1041  for (auto nn = 0; nn != Dim; ++nn) {
1042  if (nn != mm)
1043  d2MType1[mm][get_sym_index(aa, nn, Dim)](i, j, k, l) = 0;
1044  }
1045  }
1046  }
1047 
1048  if constexpr (NB == 3)
1049  for (auto aa = 0; aa != Dim; ++aa) {
1050  for (auto bb = 0; bb != Dim; ++bb) {
1051  if (aa != bb) {
1052  const auto &S = aS[get_sym_index(aa, bb, Dim)];
1053  const auto v0 = aF(aa, bb);
1054  for (auto cc = 0; cc != Dim; ++cc) {
1055  for (auto dd = 0; dd != Dim; ++dd) {
1056  if (cc != dd) {
1057  const double v1 = fVal(cc) * aF(cc, dd);
1058  d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
1059  (v1 * v0) * S(i, j, k, l);
1060  }
1061  }
1062  }
1063  }
1064  }
1065  }
1066 
1067  if constexpr (NB == 2)
1068  for (auto aa = 0; aa != Dim; ++aa) {
1069  for (auto bb = 0; bb != Dim; ++bb) {
1070  if (aa != bb)
1071  for (auto cc = 0; cc != Dim; ++cc) {
1072  for (auto dd = 0; dd != Dim; ++dd)
1073  if (cc != dd) {
1074 
1075  V r;
1076 
1077  if ((cc == 1 || dd == 1) && (aa == 1 || bb == 1))
1078  r = fVal(cc) * aF(cc, dd) * aF(aa, bb);
1079  else if (cc != 1 && dd != 1 && aa != 1 && bb != 1) {
1080 
1081  if ((aa != bb && bb != dd) && (aa != dd && bb != cc))
1082  r = ddfVal(cc) / 4;
1083  else
1084  r = 0;
1085 
1086  } else if ((cc != 1 && dd != 1) && (aa == 1 || bb == 1))
1087  r = dfVal(cc) * aF(aa, bb) / 2;
1088  else if ((cc == 1 || dd == 1) && (aa != 1 && bb != 1)) {
1089 
1090  if ((cc == 2 && dd == 1) || (cc == 2 && dd == 1))
1091  r = (
1092 
1093  dfVal(cc)
1094 
1095  - (fVal(cc) - fVal(dd)) * aF(cc, dd)
1096 
1097  ) *
1098  aF(cc, dd);
1099 
1100  else
1101  r = 0;
1102 
1103  } else
1104  r = 0;
1105 
1106  if (r)
1107  d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
1108  r * aS[get_sym_index(aa, bb, Dim)](i, j, k, l);
1109  }
1110  }
1111  }
1112  }
1113 
1114  if constexpr (NB == 1)
1115  for (auto aa = 0; aa != Dim; ++aa) {
1116  for (auto bb = 0; bb != Dim; ++bb) {
1117  if (aa != bb) {
1118  for (auto cc = 0; cc != Dim; ++cc) {
1119  for (auto dd = 0; dd != Dim; ++dd) {
1120  if (cc != dd) {
1121  if ((bb != dd) && (aa != dd && bb != cc)) {
1122  const double r = ddfVal(cc) / 4;
1123  d2MType1[dd][get_sym_index(aa, cc, Dim)](i, j, k, l) +=
1124  r * aS[get_sym_index(aa, bb, Dim)](i, j, k, l);
1125  }
1126  }
1127  }
1128  }
1129  }
1130  }
1131  }
1132 
1133  using THIS = EigenMatrixImp<T1, T2, NB, Dim>;
1134  using T3 = FTensor::Ddg<V, Dim, Dim>;
1135 
1136  T3 t_diff_A;
1137  GetDiffDiffMatImpl<THIS, V, T3, T>(*this, t_diff_A, t_S)
1139  return t_diff_A;
1140  }
1141 
1142 private:
1148  FTensor::Ddg<V, Dim, Dim> aS[(Dim * (Dim + 1)) / 2];
1149  FTensor::Ddg<V, Dim, Dim> aSM[(Dim - 1) * Dim][(Dim * (Dim + 1)) / 2];
1150  FTensor::Ddg<V, Dim, Dim> d2MType0[Dim][(Dim * (Dim + 1)) / 2];
1151  FTensor::Ddg<V, Dim, Dim> d2MType1[Dim][(Dim * (Dim + 1)) / 2];
1157 
1158  template <typename E, typename C> friend struct d2MCoefficients;
1159  template <typename E, typename C, typename G> friend struct d2MImpl;
1160  template <typename E, typename C> friend struct Fdd4MImpl;
1161  template <typename E, typename C> friend struct ReconstructMatImpl;
1162  template <typename E, typename C> friend struct FirstMatrixDirectiveImpl;
1163  template <typename E, typename C> friend struct SecondMatrixDirectiveImpl;
1164  template <typename E, typename C, typename T> friend struct GetDiffMatImpl;
1165  template <typename E, typename C, typename T3, typename T4>
1166  friend struct GetDiffDiffMatImpl;
1167 
1168 }; // namespace EigenMatrix
1169 } // namespace EigenMatrix
EigenMatrix::GetDiffMatImpl::GetDiffMatImpl
GetDiffMatImpl(E &e, T &t_a)
Definition: MatrixFunctionTemplate.hpp:539
EigenMatrix::SecondMatrixDirectiveImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:422
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
EigenMatrix::FirstMatrixDirectiveImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:386
EigenMatrix::d2MImpl::g
G g
Definition: MatrixFunctionTemplate.hpp:139
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< 1 > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:753
EigenMatrix::EigenMatrixImp::ddfVal
FTensor::Tensor1< V, Dim > ddfVal
Definition: MatrixFunctionTemplate.hpp:1156
EigenMatrix::GetDiffMatImpl::set
void set(const Number< 1 > &, const Number< 1 > &, const Number< 1 > &, const Number< l > &)
Definition: MatrixFunctionTemplate.hpp:577
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::e
E & e
Definition: MatrixFunctionTemplate.hpp:708
FTensor
JSON compatible output.
Definition: Christof_constructor.hpp:6
EigenMatrix::SecondMatrixDirectiveImpl::term1
C term1() const
Definition: MatrixFunctionTemplate.hpp:433
EigenMatrix::EigenMatrixImp::getMat
auto getMat(Fun f)
Get matrix.
Definition: MatrixFunctionTemplate.hpp:874
EigenMatrix::FirstMatrixDirectiveImpl::eval
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &) const
Definition: MatrixFunctionTemplate.hpp:413
EigenMatrix::d2MImpl< E, C, EigenMatrix::d2MCoefficients< E, C > >::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:134
EigenMatrix::Fdd4MImpl::eval_fdS2
C eval_fdS2(const Number< 1 > &, const Number< a > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:301
EigenMatrix::ReconstructMatImpl::eval
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &) const
Definition: MatrixFunctionTemplate.hpp:378
EigenMatrix::SecondMatrixDirectiveImpl::eval
C eval(const Number< nb > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:473
EigenMatrix::d2MCoefficients::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:93
EigenMatrix::Fdd4MImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:167
EigenMatrix::EigenMatrixImp::V
double V
Definition: MatrixFunctionTemplate.hpp:804
EigenMatrix::GetDiffMatImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:531
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EigenMatrix::Fdd4MImpl::eval_SM
C eval_SM(const Number< nb > &, const Number< a > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:312
EigenMatrix::ReconstructMatImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:364
EigenMatrix::GetMatImpl::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:497
EigenMatrix::d2MImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:140
EigenMatrix::d2MCoefficients
Definition: MatrixFunctionTemplate.hpp:90
EigenMatrix::GetDiffDiffMatImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:598
EigenMatrix::SecondMatrixDirectiveImpl::term3
C term3() const
Definition: MatrixFunctionTemplate.hpp:445
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::r
SecondMatrixDirectiveImpl< E, C > r
Definition: MatrixFunctionTemplate.hpp:707
EigenMatrix::d2MCoefficients::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:91
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:701
EigenMatrix::GetMatImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:495
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< N > &)
Definition: MatrixFunctionTemplate.hpp:713
EigenMatrix::SecondMatrixDirectiveImpl::eval
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:485
EigenMatrix::EigenMatrixImp::tVec
Vec & tVec
Definition: MatrixFunctionTemplate.hpp:1144
EigenMatrix::FirstMatrixDirectiveImpl::eval
C eval(const Number< nb > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &) const
Definition: MatrixFunctionTemplate.hpp:405
E
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
EigenMatrix::GetDiffMatImpl::set
void set(const Number< 1 > &, const Number< 1 > &, const Number< k > &, const Number< l > &)
Definition: MatrixFunctionTemplate.hpp:566
EigenMatrix::GetMatImpl::r
ReconstructMatImpl< E, C > r
Definition: MatrixFunctionTemplate.hpp:503
EigenMatrix::d2MImpl< E, C, EigenMatrix::d2MCoefficients< E, C > >::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:133
EigenMatrix::get_sym_index
auto get_sym_index(const Number< N1 > &, const Number< N2 > &, const Number< Dim > &)
Definition: MatrixFunctionTemplate.hpp:58
EigenMatrix::Fdd4MImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:168
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::tA
T1 & tA
Definition: MatrixFunctionTemplate.hpp:709
EigenMatrix::FirstMatrixDirectiveImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:392
EigenMatrix::GetDiffMatImpl::r
FirstMatrixDirectiveImpl< E, C > r
Definition: MatrixFunctionTemplate.hpp:540
EigenMatrix::SecondMatrixDirectiveImpl
Definition: MatrixFunctionTemplate.hpp:419
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:698
EigenMatrix::Fdd4MImpl::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:173
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
EigenMatrix::GetDiffDiffMatImpl::set
void set(const Number< 0 > &, const Number< 0 > &, const Number< 0 > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:690
EigenMatrix::GetMatImpl::set
void set(const Number< i > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:515
EigenMatrix::GetDiffMatImpl::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:534
FTensor::Tensor2_symmetric< VT2, DimT2 >
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::set
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:778
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::set
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &)
Definition: MatrixFunctionTemplate.hpp:761
EigenMatrix::d2MImpl::d2MImpl
d2MImpl(E &e)
Definition: MatrixFunctionTemplate.hpp:138
EigenMatrix::Fdd4MImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:166
EigenMatrix::Fdd4MImpl::term_fd2S
auto term_fd2S(const Number< a > &, const Number< b > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:184
EigenMatrix::Fdd4MImpl::eval_SM
C eval_SM(const Number< 1 > &, const Number< a > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:328
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
EigenMatrix::EigenMatrixImp::aM
FTensor::Tensor2_symmetric< V, Dim > aM[Dim]
Definition: MatrixFunctionTemplate.hpp:1145
EigenMatrix::SecondMatrixDirectiveImpl::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:424
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::set
void set(const Number< 0 > &, const Number< 0 > &, const Number< 0 > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:795
EigenMatrix::GetDiffMatImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:532
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
EigenMatrix::EigenMatrixImp::Fun
boost::function< double(const double)> Fun
Definition: MatrixFunctionTemplate.hpp:803
EigenMatrix::ReconstructMatImpl
Definition: MatrixFunctionTemplate.hpp:356
EigenMatrix::d2MImpl::eval
C eval(const Number< 1 > &, const Number< a > &, const int i, const int j, const int k, const int l) const
Definition: MatrixFunctionTemplate.hpp:159
EigenMatrix::SecondMatrixDirectiveImpl::term2
C term2() const
Definition: MatrixFunctionTemplate.hpp:439
EigenMatrix::GetDiffDiffMatImpl::GetDiffDiffMatImpl
GetDiffDiffMatImpl(E &e, T1 &t_a, T2 &t_S)
Definition: MatrixFunctionTemplate.hpp:606
EigenMatrix::GetDiffDiffMatImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:597
EigenMatrix::FirstMatrixDirectiveImpl::term
C term() const
Definition: MatrixFunctionTemplate.hpp:394
EigenMatrix::EigenMatrixImp::aSM
FTensor::Ddg< V, Dim, Dim > aSM[(Dim - 1) *Dim][(Dim *(Dim+1))/2]
Definition: MatrixFunctionTemplate.hpp:1149
EigenMatrix::GetDiffMatImpl::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:535
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
EigenMatrix::GetDiffMatImpl
Definition: MatrixFunctionTemplate.hpp:529
FTensor::Number
Definition: Number.hpp:11
EigenMatrix::ReconstructMatImpl::eval
C eval(const Number< nb > &, const Number< i > &, const Number< j > &) const
Definition: MatrixFunctionTemplate.hpp:371
EigenMatrix::d2MCoefficients::get
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 2 > &) const
Definition: MatrixFunctionTemplate.hpp:113
EigenMatrix::GetDiffDiffMatImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:608
EigenMatrix::EigenMatrixImp::aMM
FTensor::Ddg< V, Dim, Dim > aMM[Dim][Dim]
Definition: MatrixFunctionTemplate.hpp:1146
EigenMatrix::GetDiffDiffMatImpl
Definition: MatrixFunctionTemplate.hpp:596
EigenMatrix::GetDiffMatImpl::set
void set(const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &)
Definition: MatrixFunctionTemplate.hpp:544
a
constexpr double a
Definition: approx_sphere.cpp:30
EigenMatrix::SecondMatrixDirectiveImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:430
EigenMatrix::EigenMatrixImp::getDiffMat
auto getDiffMat(Fun f, Fun d_f)
Get derivative of matrix.
Definition: MatrixFunctionTemplate.hpp:900
EigenMatrix::SecondMatrixDirectiveImpl::SecondMatrixDirectiveImpl
SecondMatrixDirectiveImpl(E &e)
Definition: MatrixFunctionTemplate.hpp:428
EigenMatrix::Fdd4MImpl::e
E & e
Definition: MatrixFunctionTemplate.hpp:171
double
EigenMatrix
Definition: MatrixFunction.cpp:8
EigenMatrix::ReconstructMatImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:358
EigenMatrix::GetDiffDiffMatImpl::r
SecondMatrixDirectiveImpl< E, C > r
Definition: MatrixFunctionTemplate.hpp:607
EigenMatrix::ReconstructMatImpl::term
C term() const
Definition: MatrixFunctionTemplate.hpp:366
EigenMatrix::ReconstructMatImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:357
EigenMatrix::FirstMatrixDirectiveImpl
Definition: MatrixFunctionTemplate.hpp:383
EigenMatrix::GetDiffDiffMatImpl::tS
T2 & tS
Definition: MatrixFunctionTemplate.hpp:610
EigenMatrix::ReconstructMatImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:359
EigenMatrix::d2MCoefficients::get
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 3 > &) const
Definition: MatrixFunctionTemplate.hpp:104
EigenMatrix::GetDiffDiffMatImpl::set
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:673
EigenMatrix::Val
const FTensor::Tensor1< T, Dim > Val
Definition: MatrixFunction.hpp:65
EigenMatrix::GetMatImpl
Definition: MatrixFunctionTemplate.hpp:492
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
EigenMatrix::FirstMatrixDirectiveImpl::FirstMatrixDirectiveImpl
FirstMatrixDirectiveImpl(E &e)
Definition: MatrixFunctionTemplate.hpp:390
EigenMatrix::GetDiffMatImpl::set
void set(const Number< 1 > &, const Number< 1 > &, const Number< 1 > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:587
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EigenMatrix::Fdd4MImpl::eval
C eval(const Number< nb > &, const Number< a > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:340
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::tS
FTensor::Tensor2_symmetric< VT2, DimT2 > & tS
Definition: MatrixFunctionTemplate.hpp:710
EigenMatrix::GetDiffMatImpl::tA
T & tA
Definition: MatrixFunctionTemplate.hpp:541
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:697
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::set
void set(const Number< I > &, const Number< J > &, const Number< 0 > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:784
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EigenMatrix::EigenMatrixImp::getDiffDiffMat
auto getDiffDiffMat(Fun f, Fun d_f, Fun dd_f, T &t_S)
Get second directive of matrix.
Definition: MatrixFunctionTemplate.hpp:940
EigenMatrix::GetDiffDiffMatImpl::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< 1 > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:654
FTensor::Index
Definition: Index.hpp:23
EigenMatrix::d2MCoefficients::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:92
EigenMatrix::d2MImpl::eval
C eval(const Number< nb > &, const Number< a > &, const int i, const int j, const int k, const int l) const
Definition: MatrixFunctionTemplate.hpp:152
EigenMatrix::EigenMatrixImp::aF
FTensor::Tensor2< V, Dim, Dim > aF
Definition: MatrixFunctionTemplate.hpp:1153
EigenMatrix::EigenMatrixImp::d2MType1
FTensor::Ddg< V, Dim, Dim > d2MType1[Dim][(Dim *(Dim+1))/2]
Definition: MatrixFunctionTemplate.hpp:1151
convert.n
n
Definition: convert.py:82
EigenMatrix::d2MCoefficients::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:95
EigenMatrix::EigenMatrixImp::d2MType0
FTensor::Ddg< V, Dim, Dim > d2MType0[Dim][(Dim *(Dim+1))/2]
Definition: MatrixFunctionTemplate.hpp:1150
EigenMatrix::d2MImpl::term
C term(const int i, const int j, const int k, const int l) const
Definition: MatrixFunctionTemplate.hpp:143
EigenMatrix::EigenMatrixImp::EigenMatrixImp
EigenMatrixImp(Val &t_val, Vec &t_vec)
Definition: MatrixFunctionTemplate.hpp:812
EigenMatrix::d2MCoefficients::d2MCoefficients
d2MCoefficients(E &e)
Definition: MatrixFunctionTemplate.hpp:100
EigenMatrix::GetDiffDiffMatImpl::set
void set(const Number< I > &, const Number< J > &, const Number< 0 > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:679
N
const int N
Definition: speed_test.cpp:3
EigenMatrix::d2MImpl< E, C, EigenMatrix::d2MCoefficients< E, C > >::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:135
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:739
EigenMatrix::GetMatImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:493
EigenMatrix::GetDiffDiffMatImpl::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:601
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
EigenMatrix::Fdd4MImpl::eval_fdS2
C eval_fdS2(const Number< nb > &, const Number< a > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:287
EigenMatrix::GetMatImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:494
EigenMatrix::Fdd4MImpl::fd2M
auto fd2M() const
Definition: MatrixFunctionTemplate.hpp:178
EigenMatrix::get_nodiag_index
auto get_nodiag_index(const Number< N1 > &, const Number< N2 > &, const Number< Dim > &)
Definition: MatrixFunctionTemplate.hpp:74
EigenMatrix::EigenMatrixImp::aS
FTensor::Ddg< V, Dim, Dim > aS[(Dim *(Dim+1))/2]
Definition: MatrixFunctionTemplate.hpp:1148
EigenMatrix::SecondMatrixDirectiveImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:420
EigenMatrix::Fdd4MImpl::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:174
EigenMatrix::ReconstructMatImpl::ReconstructMatImpl
ReconstructMatImpl(E &e)
Definition: MatrixFunctionTemplate.hpp:363
EigenMatrix::GetMatImpl::set
void set(const Number< 1 > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:523
EigenMatrix::GetDiffMatImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:530
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
EigenMatrix::Fdd4MImpl::Fdd4MImpl
Fdd4MImpl(E &e)
Definition: MatrixFunctionTemplate.hpp:170
EigenMatrix::GetDiffDiffMatImpl::set
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &)
Definition: MatrixFunctionTemplate.hpp:662
EigenMatrix::GetMatImpl::set
void set(const Number< i > &, const Number< j > &)
Definition: MatrixFunctionTemplate.hpp:507
EigenMatrix::FirstMatrixDirectiveImpl::r
d2MImpl< E, C, d2MCoefficients< E, C > > r
Definition: MatrixFunctionTemplate.hpp:391
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
EigenMatrix::EigenMatrixImp::aF2
FTensor::Tensor2_symmetric< V, Dim > aF2
Definition: MatrixFunctionTemplate.hpp:1152
EigenMatrix::d2MCoefficients::get
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 1 >) const
Definition: MatrixFunctionTemplate.hpp:123
FTensor::Ddg
Definition: Ddg_value.hpp:7
EigenMatrix::GetDiffDiffMatImpl::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< N > &)
Definition: MatrixFunctionTemplate.hpp:613
G
constexpr IntegrationType G
Definition: level_set.cpp:33
EigenMatrix::GetDiffMatImpl::set
void set(const Number< 1 > &, const Number< j > &, const Number< k > &, const Number< l > &)
Definition: MatrixFunctionTemplate.hpp:555
EigenMatrix::EigenMatrixImp
Definition: MatrixFunctionTemplate.hpp:799
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:696
EigenMatrix::FirstMatrixDirectiveImpl::Val
typename E::Val Val
Definition: MatrixFunctionTemplate.hpp:384
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
EigenMatrix::GetDiffDiffMatImpl::add
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< 1 > &)
Definition: MatrixFunctionTemplate.hpp:638
EigenMatrix::EigenMatrixImp::fVal
FTensor::Tensor1< V, Dim > fVal
Definition: MatrixFunctionTemplate.hpp:1154
EigenMatrix::EigenMatrixImp::dfVal
FTensor::Tensor1< V, Dim > dfVal
Definition: MatrixFunctionTemplate.hpp:1155
EigenMatrix::SecondMatrixDirectiveImpl::r
Fdd4MImpl< E, C > r
Definition: MatrixFunctionTemplate.hpp:429
EigenMatrix::Fun
boost::function< T(const T)> Fun
Definition: MatrixFunction.hpp:67
EigenMatrix::SecondMatrixDirectiveImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:421
EigenMatrix::Fdd4MImpl
Definition: MatrixFunctionTemplate.hpp:165
EigenMatrix::Fdd4MImpl::term_SM
C term_SM(const Number< a > &, const Number< b > &, const Number< NB > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
Definition: MatrixFunctionTemplate.hpp:245
EigenMatrix::d2MImpl
Definition: MatrixFunctionTemplate.hpp:132
EigenMatrix::EigenMatrixImp::I
typename FTensor::Index< c, Dim > I
Definition: MatrixFunctionTemplate.hpp:807
EigenMatrix::d2MCoefficients::e
E & e
Definition: MatrixFunctionTemplate.hpp:101
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
EigenMatrix::GetMatImpl::tA
T & tA
Definition: MatrixFunctionTemplate.hpp:504
EigenMatrix::EigenMatrixImp::aG
FTensor::Ddg< V, Dim, Dim > aG[Dim][Dim]
Definition: MatrixFunctionTemplate.hpp:1147
EigenMatrix::d2MCoefficients::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:96
EigenMatrix::EigenMatrixImp::NumberDim
Number< Dim > NumberDim
Definition: MatrixFunctionTemplate.hpp:810
EigenMatrix::GetDiffDiffMatImpl::set
void set(const Number< I > &, const Number< 0 > &, const Number< K > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:685
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::NumberNb
typename E::NumberNb NumberNb
Definition: MatrixFunctionTemplate.hpp:700
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::set
void set(const Number< I > &, const Number< 0 > &, const Number< K > &, const Number< 0 > &)
Definition: MatrixFunctionTemplate.hpp:790
EigenMatrix::GetMatImpl::GetMatImpl
GetMatImpl(E &e, T &t_a)
Definition: MatrixFunctionTemplate.hpp:502
EigenMatrix::SecondMatrixDirectiveImpl::term
C term() const
Definition: MatrixFunctionTemplate.hpp:451
EigenMatrix::EigenMatrixImp::tVal
Val & tVal
Definition: MatrixFunctionTemplate.hpp:1143
EigenMatrix::GetDiffDiffMatImpl::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:602
EigenMatrix::GetMatImpl::NumberDim
typename E::NumberDim NumberDim
Definition: MatrixFunctionTemplate.hpp:498
EigenMatrix::GetDiffDiffMatImpl< E, C, T1, FTensor::Tensor2_symmetric< VT2, DimT2 > >::GetDiffDiffMatImpl
GetDiffDiffMatImpl(E &e, T1 &t_a, FTensor::Tensor2_symmetric< VT2, DimT2 > &t_S)
Definition: MatrixFunctionTemplate.hpp:705
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
EigenMatrix::GetDiffDiffMatImpl::tA
T1 & tA
Definition: MatrixFunctionTemplate.hpp:609
EigenMatrix::FirstMatrixDirectiveImpl::Vec
typename E::Vec Vec
Definition: MatrixFunctionTemplate.hpp:385
EigenMatrix::GetDiffDiffMatImpl::Fun
typename E::Fun Fun
Definition: MatrixFunctionTemplate.hpp:599