b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2tensor_calculus.H
Go to the documentation of this file.
1//------------------------------------------------------------------------
2// b2tensor_calculus.H --
3//
4//
5// written by Mathias Doreille
6// Thomas Blome <thomas.blome@dlr.de>
7//
8// Copyright (c) 2004-2012,2016,2017
9// SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// Copyright (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef __B2_TENSOR_CALCULUS_H__
22#define __B2_TENSOR_CALCULUS_H__
23
24#include <algorithm>
25#include <cfenv>
26#include <cmath>
27#include <complex>
28#include <limits>
29#include <type_traits>
30
31#include "b2csda.H"
32#include "b2util.H"
33
51namespace b2000 {
52
58template <typename T>
59inline T determinant_3x3(const std::array<std::array<T, 3>, 3>& mat_a) {
60 return mat_a[0][0] * (mat_a[1][1] * mat_a[2][2] - mat_a[1][2] * mat_a[2][1])
61 + mat_a[1][0] * (mat_a[0][2] * mat_a[2][1] - mat_a[0][1] * mat_a[2][2])
62 + mat_a[2][0] * (mat_a[0][1] * mat_a[1][2] - mat_a[0][2] * mat_a[1][1]);
63}
64
70template <typename T>
71inline T determinant_2x2(const std::array<std::array<T, 2>, 2>& mat_a) {
72 return mat_a[0][0] * mat_a[1][1] - mat_a[0][1] * mat_a[1][0];
73}
74
85template <typename T>
86inline T invert_3x3(
87 const std::array<std::array<T, 3>, 3>& mat_a, std::array<std::array<T, 3>, 3>& mat_b) {
88 const T determinant = determinant_3x3(mat_a);
89 if (determinant == T(0)) {
90 mat_b.fill({T(0), T(0), T(0)});
91 return T(0);
92 }
93 const T inv_det = T(1) / determinant;
94 mat_b[0][0] = (mat_a[1][1] * mat_a[2][2] - mat_a[1][2] * mat_a[2][1]) * inv_det;
95 mat_b[0][1] = (mat_a[0][2] * mat_a[2][1] - mat_a[0][1] * mat_a[2][2]) * inv_det;
96 mat_b[0][2] = (mat_a[0][1] * mat_a[1][2] - mat_a[0][2] * mat_a[1][1]) * inv_det;
97 mat_b[1][0] = (mat_a[1][2] * mat_a[2][0] - mat_a[1][0] * mat_a[2][2]) * inv_det;
98 mat_b[1][1] = (mat_a[0][0] * mat_a[2][2] - mat_a[0][2] * mat_a[2][0]) * inv_det;
99 mat_b[1][2] = (mat_a[0][2] * mat_a[1][0] - mat_a[0][0] * mat_a[1][2]) * inv_det;
100 mat_b[2][0] = (mat_a[1][0] * mat_a[2][1] - mat_a[1][1] * mat_a[2][0]) * inv_det;
101 mat_b[2][1] = (mat_a[0][1] * mat_a[2][0] - mat_a[0][0] * mat_a[2][1]) * inv_det;
102 mat_b[2][2] = (mat_a[0][0] * mat_a[1][1] - mat_a[0][1] * mat_a[1][0]) * inv_det;
103 return determinant;
104}
105
116template <typename T>
117inline T invert_2x2(
118 const std::array<std::array<T, 2>, 2>& mat_a, std::array<std::array<T, 2>, 2>& mat_b) {
119 T determinant = determinant_2x2(mat_a);
120 if (determinant == T(0)) {
121 mat_b.fill({T(0), T(0)});
122 return T(0);
123 }
124 const T inv_det = T(1) / determinant;
125 mat_b[0][0] = +mat_a[1][1] * inv_det;
126 mat_b[0][1] = -mat_a[0][1] * inv_det;
127 mat_b[1][0] = -mat_a[1][0] * inv_det;
128 mat_b[1][1] = +mat_a[0][0] * inv_det;
129 return determinant;
130}
131
141template <typename T>
142inline std::array<T, 2> eigenvalues_2(const std::array<std::array<T, 2>, 2>& matrix) {
143 const T trace = matrix[0][0] + matrix[1][1];
144 const T determinant = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
145 const T discriminant = trace * trace / 4 - determinant;
146 T s;
147
148 std::array<T, 2> eigenvalues;
149
150 // If the discriminant is of type floating point and close to zero
151 // but negativ, set s to zero. Throw an exception if the discriminant
152 // is of type floating point and much smaller than zero.
153 if constexpr (std::is_floating_point<T>::value) {
154 if (discriminant >= 0) {
155 s = std::sqrt(discriminant);
156 } else if (std::fabs(discriminant) <= 4096 * std::numeric_limits<T>::epsilon()) {
157 s = 0;
158 } else {
159 throw FE_INVALID; // Throw an exception if discriminant << 0.
160 }
161 } else {
162 s = std::sqrt(discriminant);
163 }
164
165 eigenvalues[0] = trace / 2 + s;
166 eigenvalues[1] = trace / 2 - s;
167
168 return eigenvalues;
169}
170
180template <typename T>
181inline void eigenvector_2(
182 const std::array<std::array<T, 2>, 2>& matrix, std::array<T, 2>& eigenvalues,
183 std::array<std::array<T, 2>, 2>& eigenvector) {
184 eigenvalues = eigenvalues_2(matrix);
185
186 if (std::abs(matrix[0][1]) > std::numeric_limits<T>::epsilon()) {
187 eigenvector[0][0] = eigenvalues[0] - matrix[1][1];
188 eigenvector[1][0] = eigenvalues[1] - matrix[1][1];
189 eigenvector[0][1] = eigenvector[1][1] = matrix[0][1];
190 } else if (std::abs(matrix[1][0]) > std::numeric_limits<T>::epsilon()) {
191 eigenvector[0][0] = eigenvector[1][0] = matrix[1][0];
192 eigenvector[0][1] = eigenvalues[0] - matrix[0][0];
193 eigenvector[1][1] = eigenvalues[1] - matrix[0][0];
194 } else {
195 eigenvector[0][0] = eigenvector[1][1] = 1;
196 eigenvector[1][0] = eigenvector[0][1] = 0;
197 }
198}
199
201template <typename T1>
202inline void set_zero_3(T1 a[3]) {
203 std::fill_n(a, 3, T1(0));
204}
205
207template <typename T1>
208inline void set_zero_2(T1 a[2]) {
209 std::fill_n(a, 2, T1(0));
210}
211
213template <typename T1>
214inline void copy_3(const T1 a[3], T1 b[3]) {
215 std::copy(a, a + 3, b);
216}
217
219template <typename T1, typename T2>
220inline void copy_3(const T1 a[3], T2 b[3]) {
221 for (int i = 0; i != 3; ++i) { b[i] = a[i]; }
222}
223
225template <typename T1>
226inline void copy_2(const T1 a[2], T1 b[2]) {
227 std::copy(a, a + 2, b);
228}
229
231template <typename T1, typename T2>
232inline void copy_2(const T1 a[2], T2 b[2]) {
233 for (int i = 0; i != 2; ++i) { b[i] = a[i]; }
234}
235
238template <typename T1, typename T2, typename T3>
239inline void sub_3(const T1 v1[3], const T2 v2[3], T3 dest[3]) {
240 dest[0] = v1[0] - v2[0];
241 dest[1] = v1[1] - v2[1];
242 dest[2] = v1[2] - v2[2];
243}
244
247template <typename T1, typename T2, typename T3>
248inline void sub_2(const T1 v1[2], const T2 v2[2], T3 dest[2]) {
249 dest[0] = v1[0] - v2[0];
250 dest[1] = v1[1] - v2[1];
251}
252
255template <typename T1, typename T2, typename T3>
256inline void add_3(const T1 v1[3], const T2 v2[3], T3 dest[3]) {
257 dest[0] = v1[0] + v2[0];
258 dest[1] = v1[1] + v2[1];
259 dest[2] = v1[2] + v2[2];
260}
261
264template <typename T1, typename T2, typename T3>
265inline void add_2(const T1 v1[2], const T2 v2[2], T3 dest[2]) {
266 dest[0] = v1[0] + v2[0];
267 dest[1] = v1[1] + v2[1];
268}
269
272template <typename T1, typename T2, typename T3, typename T4, typename T5>
273inline void add_scale_3(const T1 s1, const T2 v1[3], const T3 s2, const T4 v2[3], T5 dest[3]) {
274 dest[0] = s1 * v1[0] + s2 * v2[0];
275 dest[1] = s1 * v1[1] + s2 * v2[1];
276 dest[2] = s1 * v1[2] + s2 * v2[2];
277}
278
281template <typename T1, typename T2, typename T3, typename T4, typename T5>
282inline void add_scale_2(const T1 s1, const T2 v1[2], const T3 s2, const T4 v2[2], T5 dest[2]) {
283 dest[0] = s1 * v1[0] + s2 * v2[0];
284 dest[1] = s1 * v1[1] + s2 * v2[1];
285}
286
288template <typename T1, typename T2>
289inline void scale_3(T1 v[3], const T2 s) {
290 v[0] *= s;
291 v[1] *= s;
292 v[2] *= s;
293}
294
296template <typename T1, typename T2>
297inline void scale_2(T1 v[2], const T2 s) {
298 v[0] *= s;
299 v[1] *= s;
300}
301
303template <typename T1, typename T2>
304inline void add_3(T1 a[3], const T2 b[3]) {
305 for (int i = 0; i != 3; ++i) { a[i] += b[i]; }
306}
307
309template <typename T1, typename T2>
310inline void add_2(T1 a[2], const T2 b[2]) {
311 for (int i = 0; i != 2; ++i) { a[i] += b[i]; }
312}
313
315template <typename T1, typename T2, typename T3>
316inline void add_scale_3(T1 a[3], const T2 s, const T3 b[3]) {
317 for (int i = 0; i != 3; ++i) { a[i] += s * b[i]; }
318}
319
321template <typename T1, typename T2, typename T3>
322inline void add_scale_2(T1 a[2], const T2 s, const T3 b[2]) {
323 for (int i = 0; i != 2; ++i) { a[i] += s * b[i]; }
324}
325
327template <typename T>
328inline T dot_3(const T a[3], const T b[3]) {
329 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
330}
331
333template <typename T>
334inline std::complex<T> dot_3(const std::complex<T> a[3], const std::complex<T> b[3]) {
335 return a[0] * std::conj(b[0]) + a[1] * std::conj(b[1]) + a[2] * std::conj(b[2]);
336}
337
339template <typename T>
340inline b2000::csda<T> dot_3(const b2000::csda<T> a[3], const b2000::csda<T> b[3]) {
341 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
342}
343
345template <typename T>
346inline T dot_2(const T a[2], const T b[2]) {
347 return a[0] * b[0] + a[1] * b[1];
348}
349
351template <typename T>
352inline std::complex<T> dot_2(const std::complex<T> a[2], const std::complex<T> b[2]) {
353 return a[0] * std::conj(b[0]) + a[1] * std::conj(b[1]);
354}
355
357template <typename T>
358inline b2000::csda<T> dot_2(const b2000::csda<T> a[2], const b2000::csda<T> b[2]) {
359 return a[0] * b[0] + a[1] * b[1];
360}
361
363template <typename T>
364inline T norm_3(T a[3]) {
365 return std::sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
366}
367
369template <typename T>
370inline T norm_3(std::complex<T> a[3]) {
371 return std::sqrt(std::norm(a[0]) + std::norm(a[1]) + std::norm(a[2]));
372}
373
375template <typename T>
376inline T norm_2(T a[2]) {
377 return std::sqrt(a[0] * a[0] + a[1] * a[1]);
378}
379
381template <typename T>
382inline T norm_2(std::complex<T> a[2]) {
383 return std::sqrt(std::norm(a[0]) + std::norm(a[1]));
384}
385
388template <typename T1, typename T2, typename T3>
389inline void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3]) {
390 c[0] = a[1] * b[2] - a[2] * b[1];
391 c[1] = a[2] * b[0] - a[0] * b[2];
392 c[2] = a[0] * b[1] - a[1] * b[0];
393}
394
397template <typename T1, typename T2, typename T3>
398inline void outer_product_add_3(const T1 a[3], const T2 b[3], T3 c[3]) {
399 c[0] += a[1] * b[2] - a[2] * b[1];
400 c[1] += a[2] * b[0] - a[0] * b[2];
401 c[2] += a[0] * b[1] - a[1] * b[0];
402}
403
406template <typename T>
407inline T norm_outer_product_3(const T a[3], const T b[3]) {
408 T c[3];
409 c[0] = a[1] * b[2] - a[2] * b[1];
410 c[1] = a[2] * b[0] - a[0] * b[2];
411 c[2] = a[0] * b[1] - a[1] * b[0];
412 return norm_3(c);
413}
414
417template <typename T>
418inline T normalise_3(T a[3]) {
419 const T n = norm_3(a);
420 if (n != 0) {
421 T s = 1 / n;
422 a[0] *= s;
423 a[1] *= s;
424 a[2] *= s;
425 }
426 return n;
427}
428
431template <typename T>
432inline T normalise_3(std::complex<T> a[3]) {
433 const T n = norm_3(a);
434 if (n != 0) {
435 T s = T(1) / n;
436 a[0] *= s;
437 a[1] *= s;
438 a[2] *= s;
439 }
440 return n;
441}
442
445template <typename T>
446inline b2000::csda<T> normalise_3(b2000::csda<T> a[3]) {
447 const b2000::csda<T> n = norm_3(a);
448 if (n != 0) {
449 b2000::csda<T> s = b2000::csda<T>(1) / n;
450 a[0] *= s;
451 a[1] *= s;
452 a[2] *= s;
453 }
454 return n;
455}
456
459template <typename T>
460inline T normalise_2(T a[2]) {
461 const T n = norm_2(a);
462 if (n != 0) {
463 T s = T(1) / n;
464 a[0] *= s;
465 a[1] *= s;
466 }
467 return n;
468}
469
472template <typename T>
473inline T normalise_2(std::complex<T> a[2]) {
474 const T n = norm_2(a);
475 if (n != 0) {
476 T s = T(1) / n;
477 a[0] *= s;
478 a[1] *= s;
479 }
480 return n;
481}
482
485template <typename T>
486inline b2000::csda<T> normalise_2(b2000::csda<T> a[2]) {
487 const b2000::csda<T> n = norm_2(a);
488 if (n != 0) {
489 b2000::csda<T> s = b2000::csda<T>(1) / n;
490 a[0] *= s;
491 a[1] *= s;
492 }
493 return n;
494}
495
497template <typename T1>
498inline void set_zero_3_3(T1 a[3][3]) {
499 std::fill_n(&a[0][0], 9, T1(0));
500}
501
503template <typename T1>
504inline void set_zero_2_2(T1 a[2][2]) {
505 a[0][0] = a[0][1] = T1(0);
506 a[1][0] = a[1][1] = T1(0);
507}
508
510template <typename T1>
511inline void set_identity_3_3(T1 a[3][3]) {
512 std::fill_n(&a[0][0], 9, T1(0));
513 a[0][0] = a[1][1] = a[2][2] = T1(1);
514}
515
517template <typename T1>
518inline void set_identity_2_2(T1 a[2][2]) {
519 a[0][0] = a[1][1] = T1(1);
520 a[0][1] = a[1][0] = T1(0);
521}
522
524template <typename T1>
525inline void copy_3_3(const T1 a[3][3], T1 b[3][3]) {
526 std::copy(&a[0][0], &a[0][0] + 9, &b[0][0]);
527}
528
530template <typename T1, typename T2>
531inline void copy_3_3(const T1 a[3][3], T2 b[3][3]) {
532 for (int i = 0; i != 3; ++i) {
533 for (int j = 0; j != 3; ++j) { b[i][j] = a[i][j]; }
534 }
535}
536
538template <typename T1>
539inline void copy_2_2(const T1 a[2][2], T1 b[2][2]) {
540 std::copy(&a[0][0], &a[0][0] + 4, &b[0][0]);
541}
542
544template <typename T1, typename T2>
545inline void copy_2_2(const T1 a[2][2], T2 b[2][2]) {
546 for (int i = 0; i != 2; ++i) {
547 for (int j = 0; j != 2; ++j) { b[i][j] = a[i][j]; }
548 }
549}
550
552template <typename T>
553inline void scale_3_3(const T a, const T b[3][3], T c[3][3]) {
554 c[0][0] = a * b[0][0];
555 c[0][1] = a * b[0][1];
556 c[0][2] = a * b[0][2];
557 c[1][0] = a * b[1][0];
558 c[1][1] = a * b[1][1];
559 c[1][2] = a * b[1][2];
560 c[2][0] = a * b[2][0];
561 c[2][1] = a * b[2][1];
562 c[2][2] = a * b[2][2];
563}
564
566template <typename T>
567inline void scale_2_2(const T a, const T b[2][2], T c[2][2]) {
568 c[0][0] = a * b[0][0];
569 c[0][1] = a * b[0][1];
570 c[1][0] = a * b[1][0];
571 c[1][1] = a * b[1][1];
572}
573
575template <typename T1, typename T2>
576inline void add_3_3(T1 a[3][3], const T2 b[3][3]) {
577 T1* aa = a[0];
578 const T2* bb = b[0];
579 for (int i = 0; i != 9; ++i) { aa[i] += bb[i]; }
580}
581
583template <typename T1, typename T2>
584inline void add_2_2(T1 a[2][2], const T2 b[2][2]) {
585 T1* aa = a[0];
586 const T2* bb = b[0];
587 for (int i = 0; i != 4; ++i) { aa[i] += bb[i]; }
588}
589
591template <typename T>
592inline void add_scale_3_3(const T a, const T b[3][3], T c[3][3]) {
593 c[0][0] += a * b[0][0];
594 c[0][1] += a * b[0][1];
595 c[0][2] += a * b[0][2];
596 c[1][0] += a * b[1][0];
597 c[1][1] += a * b[1][1];
598 c[1][2] += a * b[1][2];
599 c[2][0] += a * b[2][0];
600 c[2][1] += a * b[2][1];
601 c[2][2] += a * b[2][2];
602}
603
605template <typename T>
606inline void add_scale_2_2(const T a, const T b[2][2], T c[2][2]) {
607 c[0][0] += a * b[0][0];
608 c[0][1] += a * b[0][1];
609 c[1][0] += a * b[1][0];
610 c[1][1] += a * b[1][1];
611}
612
614template <typename T1>
615inline void transpose_3_3(T1 a[3][3]) {
616 std::swap(a[0][1], a[1][0]);
617 std::swap(a[0][2], a[2][0]);
618 std::swap(a[1][2], a[2][1]);
619}
620
622template <typename T1>
623inline void transpose_2_2(T1 a[2][2]) {
624 std::swap(a[0][1], a[1][0]);
625}
626
628template <typename T1>
629inline void transpose_3_3(const T1 a[3][3], T1 b[3][3]) {
630 for (int i = 0; i != 3; ++i) {
631 for (int j = 0; j != 3; ++j) { b[i][j] = a[j][i]; }
632 }
633}
634
636template <typename T1>
637inline void transpose_2_2(const T1 a[2][2], T1 b[2][2]) {
638 for (int i = 0; i != 2; ++i) {
639 for (int j = 0; j != 2; ++j) { b[i][j] = a[j][i]; }
640 }
641}
642
644template <typename T1>
645inline T1 norm_inf_3_3(const T1 a[3][3]) {
646 const T1* aa = a[0];
647 T1 res = T1(0);
648 for (int i = 0; i < 9; i++) { res = std::max(res, aa[i]); }
649 return res;
650}
651
653template <typename T1>
654inline T1 norm_inf_2_2(const T1 a[3][3]) {
655 const T1* aa = a[0];
656 T1 res = T1(0);
657 for (int i = 0; i < 4; i++) { res = std::max(res, aa[i]); }
658 return res;
659}
660
668template <typename T>
669inline T determinant_3_3(const T a[3][3]) {
670 return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
671 + a[1][0] * (a[0][2] * a[2][1] - a[0][1] * a[2][2])
672 + a[2][0] * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
673}
674
682template <typename T>
683inline T determinant_2_2(const T a[2][2]) {
684 return a[0][0] * a[1][1] - a[0][1] * a[1][0];
685}
686
698template <typename T>
699inline T invert_3_3(const T a[3][3], T b[3][3]) {
700 const T det = determinant_3_3(a);
701 if (det == T(0)) {
702 std::fill(b[0], b[3], T(0));
703 return T(0);
704 }
705 const T inv_det = T(1) / det;
706 b[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) * inv_det;
707 b[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) * inv_det;
708 b[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * inv_det;
709 b[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) * inv_det;
710 b[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) * inv_det;
711 b[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * inv_det;
712 b[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) * inv_det;
713 b[2][1] = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) * inv_det;
714 b[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * inv_det;
715 return det;
716}
717
729template <typename T>
730inline T invert_2_2(const T a[2][2], T b[2][2]) {
731 T det = determinant_2_2(a);
732 const T inv_det = T(1) / det;
733 b[0][0] = +a[1][1] * inv_det;
734 b[0][1] = -a[0][1] * inv_det;
735 b[1][0] = -a[1][0] * inv_det;
736 b[1][1] = +a[0][0] * inv_det;
737 return det;
738}
739
741template <typename T>
742inline T invert_x_x(const T a[2][2], T b[2][2]) {
743 return invert_2_2(a, b);
744}
745
747template <typename T>
748inline T invert_x_x(const T a[3][3], T b[3][3]) {
749 return invert_3_3(a, b);
750}
751
763template <typename T>
764inline T transposed_invert_3_3(const T a[3][3], T b[3][3]) {
765 const T det = determinant_3_3(a);
766 if (det == T(0)) {
767 std::fill(b[0], b[3], T(0));
768 return T(0);
769 }
770 const T inv_det = T(1) / det;
771 b[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) * inv_det;
772 b[1][0] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) * inv_det;
773 b[2][0] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * inv_det;
774 b[0][1] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) * inv_det;
775 b[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) * inv_det;
776 b[2][1] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * inv_det;
777 b[0][2] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) * inv_det;
778 b[1][2] = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) * inv_det;
779 b[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * inv_det;
780 return det;
781}
782
794template <typename T>
795inline T transposed_invert_2_2(const T a[2][2], T b[2][2]) {
796 T det = determinant_2_2(a);
797 const T inv_det = T(1) / det;
798 b[0][0] = +a[1][1] * inv_det;
799 b[1][0] = -a[0][1] * inv_det;
800 b[0][1] = -a[1][0] * inv_det;
801 b[1][1] = +a[0][0] * inv_det;
802 return det;
803}
804
807template <typename T>
808inline void inner_product_3_3_NN(const T a[3][3], const T b[3][3], T c[3][3]) {
809 c[0][0] = a[0][0] * b[0][0] + a[1][0] * b[0][1] + a[2][0] * b[0][2];
810 c[0][1] = a[0][1] * b[0][0] + a[1][1] * b[0][1] + a[2][1] * b[0][2];
811 c[0][2] = a[0][2] * b[0][0] + a[1][2] * b[0][1] + a[2][2] * b[0][2];
812 c[1][0] = a[0][0] * b[1][0] + a[1][0] * b[1][1] + a[2][0] * b[1][2];
813 c[1][1] = a[0][1] * b[1][0] + a[1][1] * b[1][1] + a[2][1] * b[1][2];
814 c[1][2] = a[0][2] * b[1][0] + a[1][2] * b[1][1] + a[2][2] * b[1][2];
815 c[2][0] = a[0][0] * b[2][0] + a[1][0] * b[2][1] + a[2][0] * b[2][2];
816 c[2][1] = a[0][1] * b[2][0] + a[1][1] * b[2][1] + a[2][1] * b[2][2];
817 c[2][2] = a[0][2] * b[2][0] + a[1][2] * b[2][1] + a[2][2] * b[2][2];
818}
819
822template <typename T>
823inline void inner_product_3_3_TN(const T a[3][3], const T b[3][3], T c[3][3]) {
824 c[0][0] = a[0][0] * b[0][0] + a[0][1] * b[0][1] + a[0][2] * b[0][2];
825 c[0][1] = a[1][0] * b[0][0] + a[1][1] * b[0][1] + a[1][2] * b[0][2];
826 c[0][2] = a[2][0] * b[0][0] + a[2][1] * b[0][1] + a[2][2] * b[0][2];
827 c[1][0] = a[0][0] * b[1][0] + a[0][1] * b[1][1] + a[0][2] * b[1][2];
828 c[1][1] = a[1][0] * b[1][0] + a[1][1] * b[1][1] + a[1][2] * b[1][2];
829 c[1][2] = a[2][0] * b[1][0] + a[2][1] * b[1][1] + a[2][2] * b[1][2];
830 c[2][0] = a[0][0] * b[2][0] + a[0][1] * b[2][1] + a[0][2] * b[2][2];
831 c[2][1] = a[1][0] * b[2][0] + a[1][1] * b[2][1] + a[1][2] * b[2][2];
832 c[2][2] = a[2][0] * b[2][0] + a[2][1] * b[2][1] + a[2][2] * b[2][2];
833}
834
837template <typename T>
838inline void inner_product_2_2_NN(const T a[2][2], const T b[2][2], T c[2][2]) {
839 c[0][0] = a[0][0] * b[0][0] + a[1][0] * b[0][1];
840 c[0][1] = a[0][1] * b[0][0] + a[1][1] * b[0][1];
841 c[1][0] = a[0][0] * b[1][0] + a[1][0] * b[1][1];
842 c[1][1] = a[0][1] * b[1][0] + a[1][1] * b[1][1];
843}
844
847template <typename T>
848inline void inner_product_3_3_NT(const T a[3][3], const T b[3][3], T c[3][3]) {
849 c[0][0] = a[0][0] * b[0][0] + a[1][0] * b[1][0] + a[2][0] * b[2][0];
850 c[0][1] = a[0][1] * b[0][0] + a[1][1] * b[1][0] + a[2][1] * b[2][0];
851 c[0][2] = a[0][2] * b[0][0] + a[1][2] * b[1][0] + a[2][2] * b[2][0];
852 c[1][0] = a[0][0] * b[0][1] + a[1][0] * b[1][1] + a[2][0] * b[2][1];
853 c[1][1] = a[0][1] * b[0][1] + a[1][1] * b[1][1] + a[2][1] * b[2][1];
854 c[1][2] = a[0][2] * b[0][1] + a[1][2] * b[1][1] + a[2][2] * b[2][1];
855 c[2][0] = a[0][0] * b[0][2] + a[1][0] * b[1][2] + a[2][0] * b[2][2];
856 c[2][1] = a[0][1] * b[0][2] + a[1][1] * b[1][2] + a[2][1] * b[2][2];
857 c[2][2] = a[0][2] * b[0][2] + a[1][2] * b[1][2] + a[2][2] * b[2][2];
858}
859
862template <typename T>
863inline void inner_product_2_2_NT(const T a[2][2], const T b[2][2], T c[2][2]) {
864 c[0][0] = a[0][0] * b[0][0] + a[1][0] * b[1][0];
865 c[0][1] = a[0][1] * b[0][0] + a[1][1] * b[1][0];
866 c[1][0] = a[0][0] * b[0][1] + a[1][0] * b[1][1];
867 c[1][1] = a[0][1] * b[0][1] + a[1][1] * b[1][1];
868}
869
872template <typename T>
873inline void inner_product_3_3_TT(const T a[3][3], const T b[3][3], T c[3][3]) {
874 c[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0] + a[0][2] * b[2][0];
875 c[0][1] = a[1][0] * b[0][0] + a[1][1] * b[1][0] + a[1][2] * b[2][0];
876 c[0][2] = a[2][0] * b[0][0] + a[2][1] * b[1][0] + a[2][2] * b[2][0];
877 c[1][0] = a[0][0] * b[0][1] + a[0][1] * b[1][1] + a[0][2] * b[2][1];
878 c[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1] + a[1][2] * b[2][1];
879 c[1][2] = a[2][0] * b[0][1] + a[2][1] * b[1][1] + a[2][2] * b[2][1];
880 c[2][0] = a[0][0] * b[0][2] + a[0][1] * b[1][2] + a[0][2] * b[2][2];
881 c[2][1] = a[1][0] * b[0][2] + a[1][1] * b[1][2] + a[1][2] * b[2][2];
882 c[2][2] = a[2][0] * b[0][2] + a[2][1] * b[1][2] + a[2][2] * b[2][2];
883}
884
887template <typename T>
888inline void inner_product_2_2_TT(const T a[2][2], const T b[2][2], T c[2][2]) {
889 c[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0];
890 c[0][1] = a[1][0] * b[0][0] + a[1][1] * b[1][0];
891 c[1][0] = a[0][0] * b[0][1] + a[0][1] * b[1][1];
892 c[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1];
893}
894
898template <typename T>
899inline void inner_product_2_1_NN(const T a[2][2], const T b[2], T c[2]) {
900 c[0] = a[0][0] * b[0] + a[1][0] * b[1];
901 c[1] = a[0][1] * b[0] + a[1][1] * b[1];
902}
903
907template <typename T>
908inline void inner_product_2_1_TN(const T a[2][2], const T b[2], T c[2]) {
909 c[0] = a[0][0] * b[0] + a[0][1] * b[1];
910 c[1] = a[1][0] * b[0] + a[1][1] * b[1];
911}
912
915template <typename T>
916inline void inner_product_add_3_3_NN(const T a[3][3], const T b[3][3], T c[3][3]) {
917 c[0][0] += a[0][0] * b[0][0] + a[1][0] * b[0][1] + a[2][0] * b[0][2];
918 c[0][1] += a[0][1] * b[0][0] + a[1][1] * b[0][1] + a[2][1] * b[0][2];
919 c[0][2] += a[0][2] * b[0][0] + a[1][2] * b[0][1] + a[2][2] * b[0][2];
920 c[1][0] += a[0][0] * b[1][0] + a[1][0] * b[1][1] + a[2][0] * b[1][2];
921 c[1][1] += a[0][1] * b[1][0] + a[1][1] * b[1][1] + a[2][1] * b[1][2];
922 c[1][2] += a[0][2] * b[1][0] + a[1][2] * b[1][1] + a[2][2] * b[1][2];
923 c[2][0] += a[0][0] * b[2][0] + a[1][0] * b[2][1] + a[2][0] * b[2][2];
924 c[2][1] += a[0][1] * b[2][0] + a[1][1] * b[2][1] + a[2][1] * b[2][2];
925 c[2][2] += a[0][2] * b[2][0] + a[1][2] * b[2][1] + a[2][2] * b[2][2];
926}
927
930template <typename T>
931inline void inner_product_add_2_2_NN(const T a[2][2], const T b[2][2], T c[2][2]) {
932 c[0][0] += a[0][0] * b[0][0] + a[1][0] * b[0][1];
933 c[0][1] += a[0][1] * b[0][0] + a[1][1] * b[0][1];
934 c[1][0] += a[0][0] * b[1][0] + a[1][0] * b[1][1];
935 c[1][1] += a[0][1] * b[1][0] + a[1][1] * b[1][1];
936}
937
940template <typename T>
941inline void inner_product_add_3_3_TN(const T a[3][3], const T b[3][3], T c[3][3]) {
942 c[0][0] += a[0][0] * b[0][0] + a[0][1] * b[0][1] + a[0][2] * b[0][2];
943 c[0][1] += a[1][0] * b[0][0] + a[1][1] * b[0][1] + a[1][2] * b[0][2];
944 c[0][2] += a[2][0] * b[0][0] + a[2][1] * b[0][1] + a[2][2] * b[0][2];
945 c[1][0] += a[0][0] * b[1][0] + a[0][1] * b[1][1] + a[0][2] * b[1][2];
946 c[1][1] += a[1][0] * b[1][0] + a[1][1] * b[1][1] + a[1][2] * b[1][2];
947 c[1][2] += a[2][0] * b[1][0] + a[2][1] * b[1][1] + a[2][2] * b[1][2];
948 c[2][0] += a[0][0] * b[2][0] + a[0][1] * b[2][1] + a[0][2] * b[2][2];
949 c[2][1] += a[1][0] * b[2][0] + a[1][1] * b[2][1] + a[1][2] * b[2][2];
950 c[2][2] += a[2][0] * b[2][0] + a[2][1] * b[2][1] + a[2][2] * b[2][2];
951}
952
955template <typename T>
956inline void inner_product_add_2_2_TN(const T a[2][2], const T b[2][2], T c[2][2]) {
957 c[0][0] += a[0][0] * b[0][0] + a[0][1] * b[0][1];
958 c[0][1] += a[1][0] * b[0][0] + a[1][1] * b[0][1];
959 c[1][0] += a[0][0] * b[1][0] + a[0][1] * b[1][1];
960 c[1][1] += a[1][0] * b[1][0] + a[1][1] * b[1][1];
961}
962
965template <typename T>
966inline void inner_product_add_3_3_NT(const T a[3][3], const T b[3][3], T c[3][3]) {
967 c[0][0] += a[0][0] * b[0][0] + a[1][0] * b[1][0] + a[2][0] * b[2][0];
968 c[0][1] += a[0][1] * b[0][0] + a[1][1] * b[1][0] + a[2][1] * b[2][0];
969 c[0][2] += a[0][2] * b[0][0] + a[1][2] * b[1][0] + a[2][2] * b[2][0];
970 c[1][0] += a[0][0] * b[0][1] + a[1][0] * b[1][1] + a[2][0] * b[2][1];
971 c[1][1] += a[0][1] * b[0][1] + a[1][1] * b[1][1] + a[2][1] * b[2][1];
972 c[1][2] += a[0][2] * b[0][1] + a[1][2] * b[1][1] + a[2][2] * b[2][1];
973 c[2][0] += a[0][0] * b[0][2] + a[1][0] * b[1][2] + a[2][0] * b[2][2];
974 c[2][1] += a[0][1] * b[0][2] + a[1][1] * b[1][2] + a[2][1] * b[2][2];
975 c[2][2] += a[0][2] * b[0][2] + a[1][2] * b[1][2] + a[2][2] * b[2][2];
976}
977
980template <typename T>
981inline void inner_product_add_2_2_NT(const T a[2][2], const T b[2][2], T c[2][2]) {
982 c[0][0] += a[0][0] * b[0][0] + a[1][0] * b[1][0];
983 c[0][1] += a[0][1] * b[0][0] + a[1][1] * b[1][0];
984 c[1][0] += a[0][0] * b[0][1] + a[1][0] * b[1][1];
985 c[1][1] += a[0][1] * b[0][1] + a[1][1] * b[1][1];
986}
987
990template <typename T>
991inline void inner_product_2_2_TN(const T a[2][2], const T b[2][2], T c[2][2]) {
992 c[0][0] = a[0][0] * b[0][0] + a[0][1] * b[0][1];
993 c[0][1] = a[1][0] * b[0][0] + a[1][1] * b[0][1];
994 c[1][0] = a[0][0] * b[1][0] + a[0][1] * b[1][1];
995 c[1][1] = a[1][0] * b[1][0] + a[1][1] * b[1][1];
996}
997
1000template <typename T>
1001inline void inner_product_add_3_3_TT(const T a[3][3], const T b[3][3], T c[3][3]) {
1002 c[0][0] += a[0][0] * b[0][0] + a[0][1] * b[1][0] + a[0][2] * b[2][0];
1003 c[0][1] += a[1][0] * b[0][0] + a[1][1] * b[1][0] + a[1][2] * b[2][0];
1004 c[0][2] += a[2][0] * b[0][0] + a[2][1] * b[1][0] + a[2][2] * b[2][0];
1005 c[1][0] += a[0][0] * b[0][1] + a[0][1] * b[1][1] + a[0][2] * b[2][1];
1006 c[1][1] += a[1][0] * b[0][1] + a[1][1] * b[1][1] + a[1][2] * b[2][1];
1007 c[1][2] += a[2][0] * b[0][1] + a[2][1] * b[1][1] + a[2][2] * b[2][1];
1008 c[2][0] += a[0][0] * b[0][2] + a[0][1] * b[1][2] + a[0][2] * b[2][2];
1009 c[2][1] += a[1][0] * b[0][2] + a[1][1] * b[1][2] + a[1][2] * b[2][2];
1010 c[2][2] += a[2][0] * b[0][2] + a[2][1] * b[1][2] + a[2][2] * b[2][2];
1011}
1012
1015template <typename T>
1016inline void inner_product_add_2_2_TT(const T a[2][2], const T b[2][2], T c[2][2]) {
1017 c[0][0] += a[0][0] * b[0][0] + a[0][1] * b[1][0];
1018 c[0][1] += a[1][0] * b[0][0] + a[1][1] * b[1][0];
1019 c[1][0] += a[0][0] * b[0][1] + a[0][1] * b[1][1];
1020 c[1][1] += a[1][0] * b[0][1] + a[1][1] * b[1][1];
1021}
1022
1026template <typename T>
1027inline void inner_product_3_1_NN(const T a[3][3], const T b[3], T c[3]) {
1028 c[0] = a[0][0] * b[0] + a[1][0] * b[1] + a[2][0] * b[2];
1029 c[1] = a[0][1] * b[0] + a[1][1] * b[1] + a[2][1] * b[2];
1030 c[2] = a[0][2] * b[0] + a[1][2] * b[1] + a[2][2] * b[2];
1031}
1032
1036template <typename T>
1037inline void inner_product_3_1_TN(const T a[3][3], const T b[3], T c[3]) {
1038 c[0] = a[0][0] * b[0] + a[0][1] * b[1] + a[0][2] * b[2];
1039 c[1] = a[1][0] * b[0] + a[1][1] * b[1] + a[1][2] * b[2];
1040 c[2] = a[2][0] * b[0] + a[2][1] * b[1] + a[2][2] * b[2];
1041}
1042
1046template <typename T>
1047inline void inner_product_add_3_1_NN(const T a[3][3], const T b[3], T c[3]) {
1048 c[0] += a[0][0] * b[0] + a[1][0] * b[1] + a[2][0] * b[2];
1049 c[1] += a[0][1] * b[0] + a[1][1] * b[1] + a[2][1] * b[2];
1050 c[2] += a[0][2] * b[0] + a[1][2] * b[1] + a[2][2] * b[2];
1051}
1052
1055template <typename T>
1056inline void inner_product_add_3_1_TN(const T a[3][3], const T b[3], T c[3]) {
1057 c[0] += a[0][0] * b[0] + a[0][1] * b[1] + a[0][2] * b[2];
1058 c[1] += a[1][0] * b[0] + a[1][1] * b[1] + a[1][2] * b[2];
1059 c[2] += a[2][0] * b[0] + a[2][1] * b[1] + a[2][2] * b[2];
1060}
1061
1074template <typename T>
1075inline bool solve_2_2(const T a[2][2], const T y[2], T x[2]) {
1076 const T det = determinant_2_2(a);
1077 if (det == T(0)) {
1078 std::fill(x, x + 2, T(0));
1079 return false;
1080 }
1081 const T inv_det = T(1) / det;
1082 x[0] = (y[0] * a[1][1] - y[1] * a[1][0]) * inv_det;
1083 x[1] = (y[1] * a[0][0] - y[0] * a[0][1]) * inv_det;
1084 return true;
1085}
1086
1096template <typename T>
1097inline void eigenvector_2(const T m[2][2], T value[2], T vector[2][2] = 0) {
1098 const T t = m[0][0] + m[1][1];
1099 const T d = m[0][0] * m[1][1] - m[0][1] * m[1][0];
1100 const T s = std::sqrt(t * t / 4 - d);
1101 value[0] = t / 2 + s;
1102 value[1] = t / 2 - s;
1103 if (vector) {
1104 if (std::abs(m[0][1]) > std::numeric_limits<T>::epsilon()) {
1105 vector[0][0] = value[0] - m[1][1];
1106 vector[1][0] = value[1] - m[1][1];
1107 vector[0][1] = vector[1][1] = m[0][1];
1108 } else if (std::abs(m[1][0]) > std::numeric_limits<T>::epsilon()) {
1109 vector[0][0] = vector[1][0] = m[1][0];
1110 vector[0][1] = value[0] - m[0][0];
1111 vector[1][1] = value[1] - m[0][0];
1112 } else {
1113 vector[0][0] = vector[1][1] = 1;
1114 vector[1][0] = vector[0][1] = 0;
1115 }
1116 }
1117}
1118
1121template <typename T>
1122void stress_invariants(const T t[6], T& i1, T& i2, T& i3) {
1123 i1 = t[0] + t[1] + t[2];
1124 i2 = t[0] * t[1] + t[1] * t[2] + t[0] * t[2] - t[3] * t[3] - t[4] * t[4] - t[5] * t[5];
1125 i3 = t[0] * t[1] * t[2] + 2 * t[3] * t[4] * t[5] - t[0] * t[4] * t[4] - t[1] * t[5] * t[5]
1126 - t[2] * t[3] * t[3];
1127}
1128
1130template <typename T>
1131void stress_principals(const T i1, const T i2, const T i3, T& s1, T& s2, T& s3) {
1132 // TODO: a faster method exist using trigonometry (no need to solve a polynoms of order 3)
1133 // http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
1134 int nsol;
1135 solve_cubic_equation(T(1), -i1, i2, -i3, nsol, s1, s2, s3);
1136 if (nsol == T(1)) {
1137 s2 = s1;
1138 s3 = s1;
1139 } else {
1140 if (nsol == T(2)) { s3 = i1 - s1 - s2; }
1141 if (s1 < s2) { std::swap(s1, s2); }
1142 if (s1 < s3) { std::swap(s1, s3); }
1143 if (s2 < s3) { std::swap(s2, s3); }
1144 }
1145}
1146
1147} // namespace b2000
1148
1149#endif /* __B2_TENSOR_CALCULUS_H__ */
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::inner_product_3_3_TT
void inner_product_3_3_TT(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:873
b2000::set_zero_3_3
void set_zero_3_3(T1 a[3][3])
Definition b2tensor_calculus.H:498
b2000::add_2
void add_2(const T1 v1[2], const T2 v2[2], T3 dest[2])
Definition b2tensor_calculus.H:265
b2000::set_identity_3_3
void set_identity_3_3(T1 a[3][3])
Definition b2tensor_calculus.H:511
b2000::solve_2_2
bool solve_2_2(const T a[2][2], const T y[2], T x[2])
Definition b2tensor_calculus.H:1075
b2000::add_3
void add_3(const T1 v1[3], const T2 v2[3], T3 dest[3])
Definition b2tensor_calculus.H:256
b2000::inner_product_2_2_NT
void inner_product_2_2_NT(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:863
b2000::inner_product_2_2_TT
void inner_product_2_2_TT(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:888
b2000::inner_product_add_2_2_TN
void inner_product_add_2_2_TN(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:956
b2000::invert_2_2
T invert_2_2(const T a[2][2], T b[2][2])
Definition b2tensor_calculus.H:730
b2000::inner_product_add_2_2_NT
void inner_product_add_2_2_NT(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:981
b2000::sub_2
void sub_2(const T1 v1[2], const T2 v2[2], T3 dest[2])
Definition b2tensor_calculus.H:248
b2000::norm_outer_product_3
T norm_outer_product_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:407
b2000::add_scale_2_2
void add_scale_2_2(const T a, const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:606
b2000::norm_3
T norm_3(T a[3])
Definition b2tensor_calculus.H:364
b2000::copy_3
void copy_3(const T1 a[3], T1 b[3])
Definition b2tensor_calculus.H:214
b2000::inner_product_3_3_NN
void inner_product_3_3_NN(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:808
b2000::set_zero_2
void set_zero_2(T1 a[2])
Definition b2tensor_calculus.H:208
b2000::inner_product_2_1_TN
void inner_product_2_1_TN(const T a[2][2], const T b[2], T c[2])
Definition b2tensor_calculus.H:908
b2000::inner_product_add_3_3_TT
void inner_product_add_3_3_TT(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:1001
b2000::inner_product_add_3_3_NT
void inner_product_add_3_3_NT(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:966
b2000::determinant_3x3
T determinant_3x3(const std::array< std::array< T, 3 >, 3 > &mat_a)
Definition b2tensor_calculus.H:59
b2000::determinant_3_3
T determinant_3_3(const T a[3][3])
Definition b2tensor_calculus.H:669
b2000::inner_product_3_1_TN
void inner_product_3_1_TN(const T a[3][3], const T b[3], T c[3])
Definition b2tensor_calculus.H:1037
b2000::invert_x_x
T invert_x_x(const T a[2][2], T b[2][2])
Definition b2tensor_calculus.H:742
b2000::add_scale_3_3
void add_scale_3_3(const T a, const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:592
b2000::sub_3
void sub_3(const T1 v1[3], const T2 v2[3], T3 dest[3])
Definition b2tensor_calculus.H:239
b2000::set_identity_2_2
void set_identity_2_2(T1 a[2][2])
Definition b2tensor_calculus.H:518
b2000::add_scale_3
void add_scale_3(const T1 s1, const T2 v1[3], const T3 s2, const T4 v2[3], T5 dest[3])
Definition b2tensor_calculus.H:273
b2000::set_zero_2_2
void set_zero_2_2(T1 a[2][2])
Definition b2tensor_calculus.H:504
b2000::eigenvector_2
void eigenvector_2(const std::array< std::array< T, 2 >, 2 > &matrix, std::array< T, 2 > &eigenvalues, std::array< std::array< T, 2 >, 2 > &eigenvector)
Definition b2tensor_calculus.H:181
b2000::transposed_invert_3_3
T transposed_invert_3_3(const T a[3][3], T b[3][3])
Definition b2tensor_calculus.H:764
b2000::normalise_3
T normalise_3(T a[3])
Definition b2tensor_calculus.H:418
b2000::eigenvalues_2
std::array< T, 2 > eigenvalues_2(const std::array< std::array< T, 2 >, 2 > &matrix)
Definition b2tensor_calculus.H:142
b2000::add_scale_2
void add_scale_2(const T1 s1, const T2 v1[2], const T3 s2, const T4 v2[2], T5 dest[2])
Definition b2tensor_calculus.H:282
b2000::outer_product_3
void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:389
b2000::transpose_2_2
void transpose_2_2(T1 a[2][2])
Definition b2tensor_calculus.H:623
b2000::determinant_2x2
T determinant_2x2(const std::array< std::array< T, 2 >, 2 > &mat_a)
Definition b2tensor_calculus.H:71
b2000::copy_2_2
void copy_2_2(const T1 a[2][2], T1 b[2][2])
Definition b2tensor_calculus.H:539
b2000::inner_product_add_3_3_NN
void inner_product_add_3_3_NN(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:916
b2000::inner_product_add_2_2_NN
void inner_product_add_2_2_NN(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:931
b2000::solve_cubic_equation
void solve_cubic_equation(const T a, const T b, const T c, const T d, int &nsol, T &x1, T &x2, T &x3)
Definition b2util.H:770
b2000::inner_product_add_3_3_TN
void inner_product_add_3_3_TN(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:941
b2000::stress_invariants
void stress_invariants(const T t[6], T &i1, T &i2, T &i3)
Definition b2tensor_calculus.H:1122
b2000::inner_product_3_3_NT
void inner_product_3_3_NT(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:848
b2000::add_3_3
void add_3_3(T1 a[3][3], const T2 b[3][3])
Definition b2tensor_calculus.H:576
b2000::inner_product_add_3_1_TN
void inner_product_add_3_1_TN(const T a[3][3], const T b[3], T c[3])
Definition b2tensor_calculus.H:1056
b2000::inner_product_2_2_NN
void inner_product_2_2_NN(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:838
b2000::inner_product_3_1_NN
void inner_product_3_1_NN(const T a[3][3], const T b[3], T c[3])
Definition b2tensor_calculus.H:1027
b2000::add_2_2
void add_2_2(T1 a[2][2], const T2 b[2][2])
Definition b2tensor_calculus.H:584
b2000::copy_2
void copy_2(const T1 a[2], T1 b[2])
Definition b2tensor_calculus.H:226
b2000::set_zero_3
void set_zero_3(T1 a[3])
Definition b2tensor_calculus.H:202
b2000::transposed_invert_2_2
T transposed_invert_2_2(const T a[2][2], T b[2][2])
Definition b2tensor_calculus.H:795
b2000::invert_3_3
T invert_3_3(const T a[3][3], T b[3][3])
Definition b2tensor_calculus.H:699
b2000::norm_2
T norm_2(T a[2])
Definition b2tensor_calculus.H:376
b2000::norm_inf_3_3
T1 norm_inf_3_3(const T1 a[3][3])
Definition b2tensor_calculus.H:645
b2000::normalise_2
T normalise_2(T a[2])
Definition b2tensor_calculus.H:460
b2000::outer_product_add_3
void outer_product_add_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:398
b2000::copy_3_3
void copy_3_3(const T1 a[3][3], T1 b[3][3])
Definition b2tensor_calculus.H:525
b2000::norm_inf_2_2
T1 norm_inf_2_2(const T1 a[3][3])
Definition b2tensor_calculus.H:654
b2000::scale_2
void scale_2(T1 v[2], const T2 s)
Definition b2tensor_calculus.H:297
b2000::dot_3
T dot_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:328
b2000::inner_product_3_3_TN
void inner_product_3_3_TN(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:823
b2000::inner_product_add_2_2_TT
void inner_product_add_2_2_TT(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:1016
b2000::scale_3
void scale_3(T1 v[3], const T2 s)
Definition b2tensor_calculus.H:289
b2000::invert_3x3
T invert_3x3(const std::array< std::array< T, 3 >, 3 > &mat_a, std::array< std::array< T, 3 >, 3 > &mat_b)
Definition b2tensor_calculus.H:86
b2000::inner_product_2_2_TN
void inner_product_2_2_TN(const T a[2][2], const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:991
b2000::scale_2_2
void scale_2_2(const T a, const T b[2][2], T c[2][2])
Definition b2tensor_calculus.H:567
b2000::invert_2x2
T invert_2x2(const std::array< std::array< T, 2 >, 2 > &mat_a, std::array< std::array< T, 2 >, 2 > &mat_b)
Definition b2tensor_calculus.H:117
b2000::determinant_2_2
T determinant_2_2(const T a[2][2])
Definition b2tensor_calculus.H:683
b2000::dot_2
T dot_2(const T a[2], const T b[2])
Definition b2tensor_calculus.H:346
b2000::transpose_3_3
void transpose_3_3(T1 a[3][3])
Definition b2tensor_calculus.H:615
b2000::inner_product_add_3_1_NN
void inner_product_add_3_1_NN(const T a[3][3], const T b[3], T c[3])
Definition b2tensor_calculus.H:1047
b2000::inner_product_2_1_NN
void inner_product_2_1_NN(const T a[2][2], const T b[2], T c[2])
Definition b2tensor_calculus.H:899
b2000::scale_3_3
void scale_3_3(const T a, const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:553
b2000::stress_principals
void stress_principals(const T i1, const T i2, const T i3, T &s1, T &s2, T &s3)
Calculate the stress principals from the stress invariants.
Definition b2tensor_calculus.H:1131