v0.14.0
Loading...
Searching...
No Matches
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
53namespace EigenMatrix {
54
55template <int N> using Number = FTensor::Number<N>;
56
57template <int N1, int N2, int Dim>
58inline 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
66inline 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
73template <int N1, int N2, int Dim>
74inline 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
83inline 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
90template <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
132template <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
165template <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>()) *
265
266 } else {
267 return 0;
268 }
269
270 } else {
271
272 return
273
274 e.aF2(Number<a>(), Number<b>()) *
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)
291 return term_fd2S(Number<a>(), Number<nb - 1>(), Number<i>(), Number<j>(),
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)
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)
316 return term_SM(Number<a>(), Number<nb - 1>(), NumberNb(), Number<i>(),
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>(),
347 Number<m>(), Number<n>())
348
349 +
350
353 }
354};
355
356template <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
383template <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
419template <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>(),
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
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
492template <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) {}
505
506 template <int i, int j>
507 inline void set(const Number<i> &, const Number<j> &) {
508
510
513 }
514
515 template <int i> inline void set(const Number<i> &, const Number<1> &) {
516
518
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
529template <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) {}
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
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
595template <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>(),
618 Number<I - 1>(), Number<J - 1>(), Number<K - 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>(),
630
631 +
632
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>(),
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>(),
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
694template <typename E, typename C, typename T1, typename VT2, int DimT2>
695struct 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>(),
719 Number<I - 1>(), Number<J - 1>(), Number<K - 1>(),
721
722 +
723
725 Number<M>(), Number<N - 1>());
726 else
727 return tS(Number<M - 1>(), Number<N - 1>()) *
731
732 +
733
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>(),
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>()) *
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)
772 tA(Number<K - 1>(), Number<L - 1>(), Number<I - 1>(), Number<J - 1>()) =
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
799template <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
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
1142private:
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
static Index< 'M', 3 > M
static Index< 'L', 3 > L
static Index< 'J', 3 > J
static Number< 2 > N2
static Number< 1 > N1
static Index< 'K', 3 > K
constexpr double a
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
auto get_sym_index(const Number< N1 > &, const Number< N2 > &, const Number< Dim > &)
auto get_nodiag_index(const Number< N1 > &, const Number< N2 > &, const Number< Dim > &)
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
constexpr IntegrationType I
const int N
Definition: speed_test.cpp:3
FTensor::Tensor1< V, Dim > fVal
FTensor::Ddg< V, Dim, Dim > d2MType0[Dim][(Dim *(Dim+1))/2]
FTensor::Ddg< V, Dim, Dim > aMM[Dim][Dim]
FTensor::Tensor2_symmetric< V, Dim > aF2
FTensor::Ddg< V, Dim, Dim > aS[(Dim *(Dim+1))/2]
boost::function< double(const double)> Fun
FTensor::Ddg< V, Dim, Dim > d2MType1[Dim][(Dim *(Dim+1))/2]
auto getDiffMat(Fun f, Fun d_f)
Get derivative of matrix.
EigenMatrixImp(Val &t_val, Vec &t_vec)
typename FTensor::Index< c, Dim > I
FTensor::Tensor1< V, Dim > ddfVal
FTensor::Tensor2_symmetric< V, Dim > aM[Dim]
FTensor::Ddg< V, Dim, Dim > aSM[(Dim - 1) *Dim][(Dim *(Dim+1))/2]
FTensor::Ddg< V, Dim, Dim > aG[Dim][Dim]
FTensor::Tensor2< V, Dim, Dim > aF
auto getDiffDiffMat(Fun f, Fun d_f, Fun dd_f, T &t_S)
Get second directive of matrix.
FTensor::Tensor1< V, Dim > dfVal
typename E::NumberDim NumberDim
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
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
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
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
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
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
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
C eval(const Number< nb > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &) const
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &) const
d2MImpl< E, C, d2MCoefficients< E, C > > r
void set(const Number< I > &, const Number< J > &, const Number< 0 > &, const Number< 0 > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< N > &)
void set(const Number< 0 > &, const Number< 0 > &, const Number< 0 > &, const Number< 0 > &)
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< 0 > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< 1 > &)
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< 1 > &, const Number< 1 > &)
GetDiffDiffMatImpl(E &e, T1 &t_a, FTensor::Tensor2_symmetric< VT2, DimT2 > &t_S)
void set(const Number< I > &, const Number< 0 > &, const Number< K > &, const Number< 0 > &)
GetDiffDiffMatImpl(E &e, T1 &t_a, T2 &t_S)
void set(const Number< I > &, const Number< J > &, const Number< 0 > &, const Number< 0 > &)
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< N > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< M > &, const Number< 1 > &)
void set(const Number< I > &, const Number< J > &, const Number< K > &, const Number< 0 > &)
void set(const Number< 0 > &, const Number< 0 > &, const Number< 0 > &, const Number< 0 > &)
auto add(const Number< I > &, const Number< J > &, const Number< K > &, const Number< L > &, const Number< 1 > &, const Number< 1 > &)
SecondMatrixDirectiveImpl< E, C > r
void set(const Number< I > &, const Number< 0 > &, const Number< K > &, const Number< 0 > &)
void set(const Number< 1 > &, const Number< 1 > &, const Number< k > &, const Number< l > &)
FirstMatrixDirectiveImpl< E, C > r
void set(const Number< 1 > &, const Number< 1 > &, const Number< 1 > &, const Number< l > &)
void set(const Number< 1 > &, const Number< j > &, const Number< k > &, const Number< l > &)
void set(const Number< 1 > &, const Number< 1 > &, const Number< 1 > &, const Number< 1 > &)
void set(const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &)
typename E::NumberDim NumberDim
void set(const Number< 1 > &, const Number< 1 > &)
ReconstructMatImpl< E, C > r
void set(const Number< i > &, const Number< j > &)
void set(const Number< i > &, const Number< 1 > &)
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &) const
C eval(const Number< nb > &, const Number< i > &, const Number< j > &) const
C eval(const Number< 1 > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
C eval(const Number< nb > &, const Number< i > &, const Number< j > &, const Number< k > &, const Number< l > &, const Number< m > &, const Number< n > &) const
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 2 > &) const
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 3 > &) const
auto get(const Number< a > &, const Number< b > &, const int i, const int j, const int k, const int l, const Number< 1 >) const
C term(const int i, const int j, const int k, const int l) const
C eval(const Number< nb > &, const Number< a > &, const int i, const int j, const int k, const int l) const
C eval(const Number< 1 > &, const Number< a > &, const int i, const int j, const int k, const int l) const