libMesh
Loading...
Searching...
No Matches
type_tensor.h
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19
20#ifndef LIBMESH_TYPE_TENSOR_H
21#define LIBMESH_TYPE_TENSOR_H
22
23// Local includes
24#include "libmesh/libmesh_common.h"
25#include "libmesh/type_vector.h"
26
27// C++ includes
28#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
29#include <cmath>
30#include <tuple>
31
32#ifdef LIBMESH_HAVE_METAPHYSICL
33#include "metaphysicl/dualnumber_decl.h"
34#endif
35
36namespace libMesh
37{
38
39// Forward declarations
40template <typename T> class TypeTensorColumn;
41template <typename T> class ConstTypeTensorColumn;
42template <unsigned int N, typename T> class TypeNTensor;
43
51template <typename T>
53{
54 template <typename T2>
55 friend class TypeTensor;
56
57public:
62 TypeTensor ();
63
64protected:
71 explicit TypeTensor (const T & xx,
72 const T & xy=0,
73 const T & xz=0,
74 const T & yx=0,
75 const T & yy=0,
76 const T & yz=0,
77 const T & zx=0,
78 const T & zy=0,
79 const T & zz=0);
80
81
85 template <typename Scalar>
86 explicit TypeTensor (const Scalar & xx,
87 const Scalar & xy=0,
88 const Scalar & xz=0,
89 const Scalar & yx=0,
90 const Scalar & yy=0,
91 const Scalar & yz=0,
92 const Scalar & zx=0,
93 const Scalar & zy=0,
94 typename
95 std::enable_if<ScalarTraits<Scalar>::value,
96 const Scalar>::type & zz=0);
97
103 template<typename T2>
104 TypeTensor(const TypeVector<T2> & vx);
105
106 template<typename T2>
107 TypeTensor(const TypeVector<T2> & vx,
108 const TypeVector<T2> & vy);
109
110 template<typename T2>
111 TypeTensor(const TypeVector<T2> & vx,
112 const TypeVector<T2> & vy,
113 const TypeVector<T2> & vz);
114
115public:
116
120 typedef T value_type;
121
125 typedef std::tuple<unsigned int, unsigned int> index_type;
126
130 template<typename T2>
131 TypeTensor(const TypeTensor<T2> & p);
132
136 ~TypeTensor();
137
141 template<typename T2>
142 void assign (const TypeTensor<T2> &);
143
149 template <typename Scalar>
150 typename std::enable_if<
152 TypeTensor &>::type
153 operator = (const Scalar & libmesh_dbg_var(p))
154 { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }
155
159 const T & operator () (const unsigned int i, const unsigned int j) const;
160
164 T & operator () (const unsigned int i, const unsigned int j);
165
169 ConstTypeTensorColumn<T> slice (const unsigned int i) const;
170
174 TypeTensorColumn<T> slice (const unsigned int i);
175
179 TypeVector<T> row(const unsigned int r) const;
180
184 TypeVector<T> column(const unsigned int r) const;
185
191 template<typename T2>
193 operator + (const TypeTensor<T2> &) const;
194
200 template<typename T2>
202
206 template<typename T2>
207 void add (const TypeTensor<T2> &);
208
212 template <typename T2>
213 void add_scaled (const TypeTensor<T2> &, const T &);
214
220 template<typename T2>
222 operator - (const TypeTensor<T2> &) const;
223
229 template<typename T2>
231
235 template<typename T2>
236 void subtract (const TypeTensor<T2> &);
237
242 template <typename T2>
243 void subtract_scaled (const TypeTensor<T2> &, const T &);
244
249
255 template <typename Scalar>
256 auto
257 operator * (const Scalar & scalar) const -> typename std::enable_if<
259 TypeTensor<decltype(T() * scalar)>>::type;
260
266 template <typename Scalar, typename std::enable_if<
267 ScalarTraits<Scalar>::value, int>::type = 0>
268 const TypeTensor<T> & operator *= (const Scalar & factor)
269 {
270 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
271 _coords[i] *= factor;
272
273 return *this;
274 }
275
281 template <typename Scalar>
282 typename std::enable_if<
285 operator / (const Scalar &) const;
286
292 const TypeTensor<T> & operator /= (const T &);
293
300 template <typename T2>
302 operator * (const TypeTensor<T2> &) const;
303
309 template <typename T2>
311
321 template <typename T2>
323 contract (const TypeTensor<T2> &) const;
324
331 template <typename T2>
333 operator * (const TypeVector<T2> &) const;
334
341 template <typename T2>
343 left_multiply (const TypeVector<T2> & p) const;
344
348 TypeTensor<T> transpose() const;
349
353 TypeTensor<T> inverse() const;
354
361 void solve(const TypeVector<T> & b, TypeVector<T> & x) const;
362
367 auto norm() const;
368
373 auto norm_sq() const;
374
378 bool is_zero() const;
379
386 T det() const;
387
391 T tr() const;
392
396 void zero();
397
401 bool operator == (const TypeTensor<T> & rhs) const;
402
407 bool operator < (const TypeTensor<T> & rhs) const;
408
412 bool operator > (const TypeTensor<T> & rhs) const;
413
417 void print(std::ostream & os = libMesh::out) const;
418
426 friend std::ostream & operator << (std::ostream & os, const TypeTensor<T> & t)
427 {
428 t.print(os);
429 return os;
430 }
431
436 void write_unformatted (std::ostream & out_stream,
437 const bool newline = true) const;
438
444 bool is_hpd(Real rel_tol = TOLERANCE*TOLERANCE) const;
445
446protected:
447
451 T _coords[LIBMESH_DIM*LIBMESH_DIM];
452};
453
454
455template <typename T>
457{
458public:
459
461 unsigned int j) :
462 _tensor(&tensor), _j(j) {}
463
468 T & operator () (const unsigned int i)
469 { return (*_tensor)(i,_j); }
470
471 T & slice (const unsigned int i)
472 { return (*_tensor)(i,_j); }
473
478 {
479 for (unsigned int i=0; i != LIBMESH_DIM; ++i)
480 (*this)(i) = rhs(i);
481 return *this;
482 }
483
484private:
486 const unsigned int _j;
487};
488
489
490template <typename T>
492{
493public:
494
496 unsigned int j) :
497 _tensor(&tensor), _j(j) {}
498
502 const T & operator () (const unsigned int i) const
503 { return (*_tensor)(i,_j); }
504
505 const T & slice (const unsigned int i) const
506 { return (*_tensor)(i,_j); }
507
508private:
510 const unsigned int _j;
511};
512
513//------------------------------------------------------
514// Inline functions
515template <typename T>
516inline
518{
519 _coords[0] = {};
520
521#if LIBMESH_DIM > 1
522 _coords[1] = {};
523 _coords[2] = {};
524 _coords[3] = {};
525#endif
526
527#if LIBMESH_DIM > 2
528 _coords[4] = {};
529 _coords[5] = {};
530 _coords[6] = {};
531 _coords[7] = {};
532 _coords[8] = {};
533#endif
534}
535
536
537
538template <typename T>
539inline
541 const T & xy,
542 const T & xz,
543 const T & yx,
544 const T & yy,
545 const T & yz,
546 const T & zx,
547 const T & zy,
548 const T & zz)
549{
550 _coords[0] = xx;
551
552#if LIBMESH_DIM == 2
553 _coords[1] = xy;
554 _coords[2] = yx;
555 _coords[3] = yy;
556#elif LIBMESH_DIM == 1
557 libmesh_assert_equal_to (xy, 0);
558 libmesh_assert_equal_to (yx, 0);
559 libmesh_assert_equal_to (yy, 0);
560 libmesh_ignore(xy, yx, yy);
561#endif
562
563#if LIBMESH_DIM == 3
564 _coords[1] = xy;
565 _coords[2] = xz;
566 _coords[3] = yx;
567 _coords[4] = yy;
568 _coords[5] = yz;
569 _coords[6] = zx;
570 _coords[7] = zy;
571 _coords[8] = zz;
572#else
573 libmesh_assert_equal_to (xz, 0);
574 libmesh_assert_equal_to (yz, 0);
575 libmesh_assert_equal_to (zx, 0);
576 libmesh_assert_equal_to (zy, 0);
577 libmesh_assert_equal_to (zz, 0);
578 libmesh_ignore(xz, yz, zx, zy, zz);
579#endif
580}
581
582
583template <typename T>
584template <typename Scalar>
585inline
586TypeTensor<T>::TypeTensor (const Scalar & xx,
587 const Scalar & xy,
588 const Scalar & xz,
589 const Scalar & yx,
590 const Scalar & yy,
591 const Scalar & yz,
592 const Scalar & zx,
593 const Scalar & zy,
594 typename
595 std::enable_if<ScalarTraits<Scalar>::value,
596 const Scalar>::type & zz)
597{
598 _coords[0] = xx;
599
600#if LIBMESH_DIM == 2
601 _coords[1] = xy;
602 _coords[2] = yx;
603 _coords[3] = yy;
604#elif LIBMESH_DIM == 1
605 libmesh_assert_equal_to (xy, 0);
606 libmesh_assert_equal_to (yx, 0);
607 libmesh_assert_equal_to (yy, 0);
608 libmesh_ignore(xy, yx, yy);
609#endif
610
611#if LIBMESH_DIM == 3
612 _coords[1] = xy;
613 _coords[2] = xz;
614 _coords[3] = yx;
615 _coords[4] = yy;
616 _coords[5] = yz;
617 _coords[6] = zx;
618 _coords[7] = zy;
619 _coords[8] = zz;
620#else
621 libmesh_assert_equal_to (xz, 0);
622 libmesh_assert_equal_to (yz, 0);
623 libmesh_assert_equal_to (zx, 0);
624 libmesh_assert_equal_to (zy, 0);
625 libmesh_assert_equal_to (zz, 0);
626 libmesh_ignore(xz, yz, zx, zy, zz);
627#endif
628}
629
630
631
632template <typename T>
633template<typename T2>
634inline
635TypeTensor<T>::TypeTensor (const TypeTensor <T2> & p)
636{
637 // copy the nodes from vector p to me
638 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
639 _coords[i] = p._coords[i];
640}
641
642
643template <typename T>
644template <typename T2>
646{
647 libmesh_assert_equal_to (LIBMESH_DIM, 1);
648 _coords[0] = vx(0);
649}
650
651template <typename T>
652template <typename T2>
654 const TypeVector<T2> & vy)
655{
656 libmesh_assert_equal_to (LIBMESH_DIM, 2);
657#if LIBMESH_DIM > 1
658 _coords[0] = vx(0);
659 _coords[1] = vx(1);
660 _coords[2] = vy(0);
661 _coords[3] = vy(1);
662#else
663 libmesh_ignore(vx, vy);
664#endif
665}
666
667template <typename T>
668template <typename T2>
670 const TypeVector<T2> & vy,
671 const TypeVector<T2> & vz)
672{
673 libmesh_assert_equal_to (LIBMESH_DIM, 3);
674#if LIBMESH_DIM > 2
675 _coords[0] = vx(0);
676 _coords[1] = vx(1);
677 _coords[2] = vx(2);
678 _coords[3] = vy(0);
679 _coords[4] = vy(1);
680 _coords[5] = vy(2);
681 _coords[6] = vz(0);
682 _coords[7] = vz(1);
683 _coords[8] = vz(2);
684#else
685 libmesh_ignore(vx, vy, vz);
686#endif
687}
688
689
690
691
692template <typename T>
693inline
697
698
699
700template <typename T>
701template<typename T2>
702inline
704{
705 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
706 _coords[i] = p._coords[i];
707}
708
709
710
711template <typename T>
712inline
713const T & TypeTensor<T>::operator () (const unsigned int i,
714 const unsigned int j) const
715{
716 libmesh_assert_less (i, 3);
717 libmesh_assert_less (j, 3);
718
719#if LIBMESH_DIM < 3
720 const static T my_zero = 0;
721 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
722 return my_zero;
723#endif
724
725 return _coords[i*LIBMESH_DIM+j];
726}
727
728
729
730template <typename T>
731inline
732T & TypeTensor<T>::operator () (const unsigned int i,
733 const unsigned int j)
734{
735#if LIBMESH_DIM < 3
736
737 libmesh_error_msg_if(i >= LIBMESH_DIM || j >= LIBMESH_DIM,
738 "ERROR: You are assigning to a tensor component that is out of range for the compiled LIBMESH_DIM!");
739
740#endif
741
742 libmesh_assert_less (i, LIBMESH_DIM);
743 libmesh_assert_less (j, LIBMESH_DIM);
744
745 return _coords[i*LIBMESH_DIM+j];
746}
747
748
749template <typename T>
750inline
752TypeTensor<T>::slice (const unsigned int i) const
753{
754 libmesh_assert_less (i, LIBMESH_DIM);
755 return ConstTypeTensorColumn<T>(*this, i);
756}
757
758
759template <typename T>
760inline
762TypeTensor<T>::slice (const unsigned int i)
763{
764 libmesh_assert_less (i, LIBMESH_DIM);
765 return TypeTensorColumn<T>(*this, i);
766}
767
768
769template <typename T>
770inline
772TypeTensor<T>::row(const unsigned int r) const
773{
774 TypeVector<T> return_vector;
775
776 for (unsigned int j=0; j<LIBMESH_DIM; j++)
777 return_vector._coords[j] = _coords[r*LIBMESH_DIM + j];
778
779 return return_vector;
780}
781
782
783template <typename T>
784inline
786TypeTensor<T>::column(const unsigned int r) const
787{
788 TypeVector<T> return_vector;
789
790 for (unsigned int i=0; i<LIBMESH_DIM; i++)
791 return_vector._coords[i] = _coords[i*LIBMESH_DIM + r];
792
793 return return_vector;
794}
795
796
797template <typename T>
798template<typename T2>
799inline
802{
803
804#if LIBMESH_DIM == 1
806#endif
807
808#if LIBMESH_DIM == 2
810 _coords[1] + p._coords[1],
811 0.,
812 _coords[2] + p._coords[2],
813 _coords[3] + p._coords[3]);
814#endif
815
816#if LIBMESH_DIM == 3
818 _coords[1] + p._coords[1],
819 _coords[2] + p._coords[2],
820 _coords[3] + p._coords[3],
821 _coords[4] + p._coords[4],
822 _coords[5] + p._coords[5],
823 _coords[6] + p._coords[6],
824 _coords[7] + p._coords[7],
825 _coords[8] + p._coords[8]);
826#endif
827
828}
829
830
831
832template <typename T>
833template<typename T2>
834inline
836{
837 this->add (p);
838
839 return *this;
840}
841
842
843
844template <typename T>
845template<typename T2>
846inline
848{
849 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
850 _coords[i] += p._coords[i];
851}
852
853
854
855template <typename T>
856template <typename T2>
857inline
858void TypeTensor<T>::add_scaled (const TypeTensor<T2> & p, const T & factor)
859{
860 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
861 _coords[i] += factor*p._coords[i];
862
863}
864
865
866
867template <typename T>
868template<typename T2>
869inline
872{
873
874#if LIBMESH_DIM == 1
876#endif
877
878#if LIBMESH_DIM == 2
880 _coords[1] - p._coords[1],
881 0.,
882 _coords[2] - p._coords[2],
883 _coords[3] - p._coords[3]);
884#endif
885
886#if LIBMESH_DIM == 3
888 _coords[1] - p._coords[1],
889 _coords[2] - p._coords[2],
890 _coords[3] - p._coords[3],
891 _coords[4] - p._coords[4],
892 _coords[5] - p._coords[5],
893 _coords[6] - p._coords[6],
894 _coords[7] - p._coords[7],
895 _coords[8] - p._coords[8]);
896#endif
897
898}
899
900
901
902template <typename T>
903template <typename T2>
904inline
906{
907 this->subtract (p);
908
909 return *this;
910}
911
912
913
914template <typename T>
915template <typename T2>
916inline
918{
919 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
920 _coords[i] -= p._coords[i];
921}
922
923
924
925template <typename T>
926template <typename T2>
927inline
928void TypeTensor<T>::subtract_scaled (const TypeTensor<T2> & p, const T & factor)
929{
930 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
931 _coords[i] -= factor*p._coords[i];
932}
933
934
935
936template <typename T>
937inline
939{
940
941#if LIBMESH_DIM == 1
942 return TypeTensor(-_coords[0]);
943#endif
944
945#if LIBMESH_DIM == 2
946 return TypeTensor(-_coords[0],
947 -_coords[1],
948 -_coords[2],
949 -_coords[3]);
950#endif
951
952#if LIBMESH_DIM == 3
953 return TypeTensor(-_coords[0],
954 -_coords[1],
955 -_coords[2],
956 -_coords[3],
957 -_coords[4],
958 -_coords[5],
959 -_coords[6],
960 -_coords[7],
961 -_coords[8]);
962#endif
963
964}
965
966
967
968template <typename T>
969template <typename Scalar>
970inline
971auto
972TypeTensor<T>::operator * (const Scalar & factor) const -> typename std::enable_if<
974 TypeTensor<decltype(T() * factor)>>::type
975{
976 typedef decltype((*this)(0, 0) * factor) TS;
977
978
979#if LIBMESH_DIM == 1
980 return TypeTensor<TS>(_coords[0]*factor);
981#endif
982
983#if LIBMESH_DIM == 2
984 return TypeTensor<TS>(_coords[0]*factor,
985 _coords[1]*factor,
986 _coords[2]*factor,
987 _coords[3]*factor);
988#endif
989
990#if LIBMESH_DIM == 3
991 return TypeTensor<TS>(_coords[0]*factor,
992 _coords[1]*factor,
993 _coords[2]*factor,
994 _coords[3]*factor,
995 _coords[4]*factor,
996 _coords[5]*factor,
997 _coords[6]*factor,
998 _coords[7]*factor,
999 _coords[8]*factor);
1000#endif
1001}
1002
1003
1004
1005template <typename T, typename Scalar>
1006inline
1007typename std::enable_if<
1010operator * (const Scalar & factor,
1011 const TypeTensor<T> & t)
1012{
1013 return t * factor;
1014}
1015
1016template <typename T>
1017template <typename Scalar>
1018inline
1019typename std::enable_if<
1021 TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1022TypeTensor<T>::operator / (const Scalar & factor) const
1023{
1024 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1025
1026 typedef typename CompareTypes<T, Scalar>::supertype TS;
1027
1028#if LIBMESH_DIM == 1
1029 return TypeTensor<TS>(_coords[0]/factor);
1030#endif
1031
1032#if LIBMESH_DIM == 2
1033 return TypeTensor<TS>(_coords[0]/factor,
1034 _coords[1]/factor,
1035 _coords[2]/factor,
1036 _coords[3]/factor);
1037#endif
1038
1039#if LIBMESH_DIM == 3
1040 return TypeTensor<TS>(_coords[0]/factor,
1041 _coords[1]/factor,
1042 _coords[2]/factor,
1043 _coords[3]/factor,
1044 _coords[4]/factor,
1045 _coords[5]/factor,
1046 _coords[6]/factor,
1047 _coords[7]/factor,
1048 _coords[8]/factor);
1049#endif
1050
1051}
1052
1053
1054
1055template <typename T>
1056inline
1058{
1059#if LIBMESH_DIM == 1
1060 return TypeTensor(_coords[0]);
1061#endif
1062
1063#if LIBMESH_DIM == 2
1064 return TypeTensor(_coords[0],
1065 _coords[2],
1066 _coords[1],
1067 _coords[3]);
1068#endif
1069
1070#if LIBMESH_DIM == 3
1071 return TypeTensor(_coords[0],
1072 _coords[3],
1073 _coords[6],
1074 _coords[1],
1075 _coords[4],
1076 _coords[7],
1077 _coords[2],
1078 _coords[5],
1079 _coords[8]);
1080#endif
1081}
1082
1083
1084
1085template <typename T>
1086inline
1088{
1089#if LIBMESH_DIM == 1
1090 if (_coords[0] == static_cast<T>(0.))
1091 libmesh_convergence_failure();
1092 return TypeTensor(1. / _coords[0]);
1093#endif
1094
1095#if LIBMESH_DIM == 2
1096 // Get convenient reference to this.
1097 const TypeTensor<T> & A = *this;
1098
1099 // Use temporary variables, avoid multiple accesses
1100 T a = A(0,0), b = A(0,1),
1101 c = A(1,0), d = A(1,1);
1102
1103 // Make sure det = ad - bc is not zero
1104 T my_det = a*d - b*c;
1105
1106 if (my_det == static_cast<T>(0.))
1107 libmesh_convergence_failure();
1108
1109 return TypeTensor(d/my_det, -b/my_det, -c/my_det, a/my_det);
1110#endif
1111
1112#if LIBMESH_DIM == 3
1113 // Get convenient reference to this.
1114 const TypeTensor<T> & A = *this;
1115
1116 T a11 = A(0,0), a12 = A(0,1), a13 = A(0,2),
1117 a21 = A(1,0), a22 = A(1,1), a23 = A(1,2),
1118 a31 = A(2,0), a32 = A(2,1), a33 = A(2,2);
1119
1120 T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1121
1122 if (my_det == static_cast<T>(0.))
1123 libmesh_convergence_failure();
1124
1125 // Inline comment characters are for lining up columns.
1126 return TypeTensor( (a33*a22-a32*a23)/my_det, -(a33*a12-a32*a13)/my_det, (a23*a12-a22*a13)/my_det,
1127 -(a33*a21-a31*a23)/my_det, (a33*a11-a31*a13)/my_det, -(a23*a11-a21*a13)/my_det,
1128 (a32*a21-a31*a22)/my_det, -(a32*a11-a31*a12)/my_det, (a22*a11-a21*a12)/my_det);
1129#endif
1130}
1131
1132
1133
1134template <typename T>
1135inline
1137{
1138#if LIBMESH_DIM == 1
1139 if (_coords[0] == static_cast<T>(0.))
1140 libmesh_convergence_failure();
1141 x(0) = b(0) / _coords[0];
1142#endif
1143
1144#if LIBMESH_DIM == 2
1145 T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1146
1147 if (my_det == static_cast<T>(0.))
1148 libmesh_convergence_failure();
1149
1150 T my_det_inv = 1./my_det;
1151
1152 x(0) = my_det_inv*( _coords[3]*b(0) - _coords[1]*b(1));
1153 x(1) = my_det_inv*(-_coords[2]*b(0) + _coords[0]*b(1));
1154#endif
1155
1156#if LIBMESH_DIM == 3
1157 T my_det =
1158 // a11*(a33 *a22 - a32 *a23)
1159 _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1160 // -a21*(a33 *a12 - a32 *a13)
1161 -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1162 // +a31*(a23 *a12 - a22 *a13)
1163 +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1164
1165 if (my_det == static_cast<T>(0.))
1166 libmesh_convergence_failure();
1167
1168 T my_det_inv = 1./my_det;
1169 x(0) = my_det_inv*((_coords[8]*_coords[4] - _coords[7]*_coords[5])*b(0) -
1170 (_coords[8]*_coords[1] - _coords[7]*_coords[2])*b(1) +
1171 (_coords[5]*_coords[1] - _coords[4]*_coords[2])*b(2));
1172
1173 x(1) = my_det_inv*((_coords[6]*_coords[5] - _coords[8]*_coords[3])*b(0) +
1174 (_coords[8]*_coords[0] - _coords[6]*_coords[2])*b(1) -
1175 (_coords[5]*_coords[0] - _coords[3]*_coords[2])*b(2));
1176
1177 x(2) = my_det_inv*((_coords[7]*_coords[3] - _coords[6]*_coords[4])*b(0) -
1178 (_coords[7]*_coords[0] - _coords[6]*_coords[1])*b(1) +
1179 (_coords[4]*_coords[0] - _coords[3]*_coords[1])*b(2));
1180#endif
1181}
1182
1183
1184
1185template <typename T>
1186inline
1188{
1189 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1190
1191 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1192 _coords[i] /= factor;
1193
1194 return *this;
1195}
1196
1197
1198
1199
1200template <typename T>
1201template <typename T2>
1202inline
1205{
1207 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1208 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1209 returnval(i) += (*this)(i,j)*p(j);
1210
1211 return returnval;
1212}
1213
1214template <typename T>
1215template <typename T2>
1216inline
1219{
1221 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1222 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1223 returnval(i) += p(j)*(*this)(j,i);
1224
1225 return returnval;
1226}
1227
1228template <typename T, typename T2>
1229inline
1232{
1233 return b.left_multiply(a);
1234}
1235
1236template <typename T>
1237template <typename T2>
1238inline
1239TypeTensor<typename CompareTypes<T, T2>::supertype>
1241{
1243 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1244 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1245 for (unsigned int k=0; k<LIBMESH_DIM; k++)
1246 returnval(i,j) += (*this)(i,k)*p(k,j);
1247
1248 return returnval;
1249}
1250
1251template <typename T>
1252template <typename T2>
1253inline
1255{
1256 TypeTensor<T> temp;
1257 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1258 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1259 for (unsigned int k=0; k<LIBMESH_DIM; k++)
1260 temp(i,j) += (*this)(i,k)*p(k,j);
1261
1262 this->assign(temp);
1263 return *this;
1264}
1265
1266
1271template <typename T>
1272template <typename T2>
1273inline
1276{
1277 typename CompareTypes<T,T2>::supertype sum = 0.;
1278 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1279 sum += _coords[i]*t._coords[i];
1280 return sum;
1281}
1282
1283
1284
1285template <typename T>
1286inline
1288{
1289 using std::sqrt;
1290 return sqrt(this->norm_sq());
1291}
1292
1293
1294template <typename T>
1295inline
1297{
1298 for (const auto & val : _coords)
1299 if (val != T(0))
1300 return false;
1301 return true;
1302}
1303
1304template <typename T>
1305inline
1307{
1308#if LIBMESH_DIM == 1
1309 return _coords[0];
1310#endif
1311
1312#if LIBMESH_DIM == 2
1313 return (_coords[0] * _coords[3]
1314 - _coords[1] * _coords[2]);
1315#endif
1316
1317#if LIBMESH_DIM == 3
1318 return (_coords[0] * _coords[4] * _coords[8]
1319 + _coords[1] * _coords[5] * _coords[6]
1320 + _coords[2] * _coords[3] * _coords[7]
1321 - _coords[0] * _coords[5] * _coords[7]
1322 - _coords[1] * _coords[3] * _coords[8]
1323 - _coords[2] * _coords[4] * _coords[6]);
1324#endif
1325}
1326
1327template <typename T>
1328inline
1330{
1331#if LIBMESH_DIM == 1
1332 return _coords[0];
1333#endif
1334
1335#if LIBMESH_DIM == 2
1336 return _coords[0] + _coords[3];
1337#endif
1338
1339#if LIBMESH_DIM == 3
1340 return _coords[0] + _coords[4] + _coords[8];
1341#endif
1342}
1343
1344template <typename T>
1345inline
1347{
1348 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1349 _coords[i] = 0.;
1350}
1351
1352
1353
1354template <typename T>
1355inline
1357{
1358 Real sum = 0.;
1359 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1360 sum += TensorTools::norm_sq(_coords[i]);
1361 return sum;
1362}
1363
1364
1365
1366template <typename T>
1367inline
1369{
1370#if LIBMESH_DIM == 1
1371 return (std::abs(_coords[0] - rhs._coords[0])
1372 < TOLERANCE);
1373#endif
1374
1375#if LIBMESH_DIM == 2
1376 return ((std::abs(_coords[0] - rhs._coords[0]) +
1377 std::abs(_coords[1] - rhs._coords[1]) +
1378 std::abs(_coords[2] - rhs._coords[2]) +
1379 std::abs(_coords[3] - rhs._coords[3]))
1380 < 4.*TOLERANCE);
1381#endif
1382
1383#if LIBMESH_DIM == 3
1384 return ((std::abs(_coords[0] - rhs._coords[0]) +
1385 std::abs(_coords[1] - rhs._coords[1]) +
1386 std::abs(_coords[2] - rhs._coords[2]) +
1387 std::abs(_coords[3] - rhs._coords[3]) +
1388 std::abs(_coords[4] - rhs._coords[4]) +
1389 std::abs(_coords[5] - rhs._coords[5]) +
1390 std::abs(_coords[6] - rhs._coords[6]) +
1391 std::abs(_coords[7] - rhs._coords[7]) +
1392 std::abs(_coords[8] - rhs._coords[8]))
1393 < 9.*TOLERANCE);
1394#endif
1395
1396}
1397
1398template <typename T>
1399void TypeTensor<T>::print(std::ostream & os) const
1400{
1401#if LIBMESH_DIM == 1
1402
1403 os << "x=" << (*this)(0,0) << std::endl;
1404
1405#endif
1406#if LIBMESH_DIM == 2
1407
1408 os << "(xx,xy)=("
1409 << std::setw(8) << (*this)(0,0) << ", "
1410 << std::setw(8) << (*this)(0,1) << ")"
1411 << std::endl;
1412 os << "(yx,yy)=("
1413 << std::setw(8) << (*this)(1,0) << ", "
1414 << std::setw(8) << (*this)(1,1) << ")"
1415 << std::endl;
1416
1417#endif
1418#if LIBMESH_DIM == 3
1419
1420 os << "(xx,xy,xz)=("
1421 << std::setw(8) << (*this)(0,0) << ", "
1422 << std::setw(8) << (*this)(0,1) << ", "
1423 << std::setw(8) << (*this)(0,2) << ")"
1424 << std::endl;
1425 os << "(yx,yy,yz)=("
1426 << std::setw(8) << (*this)(1,0) << ", "
1427 << std::setw(8) << (*this)(1,1) << ", "
1428 << std::setw(8) << (*this)(1,2) << ")"
1429 << std::endl;
1430 os << "(zx,zy,zz)=("
1431 << std::setw(8) << (*this)(2,0) << ", "
1432 << std::setw(8) << (*this)(2,1) << ", "
1433 << std::setw(8) << (*this)(2,2) << ")"
1434 << std::endl;
1435#endif
1436}
1437
1438
1439// A matrix is finite iff every component is
1440template <typename T>
1441bool isfinite (const TypeTensor<T> & var)
1442{
1443 using std::isfinite;
1444 using libMesh::isfinite; // for T==complex
1445 for (auto i : make_range(LIBMESH_DIM))
1446 for (auto j : make_range(LIBMESH_DIM))
1447 if (!isfinite(var(i,j)))
1448 return false;
1449 return true;
1450}
1451
1452
1453// A matrix is infinite iff some component is infinite but no
1454// component is NaN.
1455//
1456// This is arguably inconsistent with our std::complex overload (and
1457// the C99 Annex G recommendations for _Complex, and C++ std::complex
1458// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1459// this is probably safer for users.
1460template <typename T>
1461bool isinf (const TypeTensor<T> & var)
1462{
1463 using std::isinf;
1464 using libMesh::isinf; // for T==complex
1465 using std::isnan;
1466 using libMesh::isnan;
1467 bool has_inf = false;
1468 for (auto i : make_range(LIBMESH_DIM))
1469 for (auto j : make_range(LIBMESH_DIM))
1470 {
1471 // NaN anywhere makes us NaN, not inf
1472 if (isnan(var(i,j)))
1473 return false;
1474 has_inf = has_inf || isinf(var(i,j));
1475 }
1476 return has_inf;
1477}
1478
1479
1480
1481// A matrix is NaN iff some component is NaN
1482//
1483// This is arguably inconsistent with our std::complex overload (and
1484// the C99 Annex G recommendations for _Complex, and C++ std::complex
1485// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1486// this is probably safer for users.
1487template <typename T>
1488bool isnan (const TypeTensor<T> & var)
1489{
1490 using std::isnan;
1491 using libMesh::isnan; // for T==complex
1492 for (auto i : make_range(LIBMESH_DIM))
1493 for (auto j : make_range(LIBMESH_DIM))
1494 if (isnan(var(i,j)))
1495 return true;
1496 return false;
1497}
1498
1499
1500
1501template <typename T, typename T2>
1502inline
1503TypeTensor<typename CompareTypes<T, T2>::supertype>
1505{
1507 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1508 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1509 ret(i,j) = a(i) * libmesh_conj(b(j));
1510
1511 return ret;
1512}
1513
1514
1515} // namespace libMesh
1516
1517#ifdef LIBMESH_HAVE_METAPHYSICL
1518namespace MetaPhysicL
1519{
1520template <typename T>
1521struct RawType<libMesh::TypeTensor<T>>
1522{
1524
1526 {
1527 value_type ret;
1528 for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1529 for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1530 ret(i,j) = raw_value(in(i,j));
1531
1532 return ret;
1533 }
1534};
1535
1536template <typename T, typename U>
1537struct ReplaceAlgebraicType<libMesh::TypeTensor<T>, U>
1538{
1539 typedef U type;
1540};
1541}
1542#endif
1543
1544#endif // LIBMESH_TYPE_TENSOR_H
const TypeTensor< T > * _tensor
const T & slice(const unsigned int i) const
ConstTypeTensorColumn(const TypeTensor< T > &tensor, unsigned int j)
const T & operator()(const unsigned int i) const
Return the element of the tensor.
T & slice(const unsigned int i)
TypeTensorColumn< T > & operator=(const TypeVector< T > &rhs)
Assign values to this column of the tensor.
TypeTensorColumn(TypeTensor< T > &tensor, unsigned int j)
const unsigned int _j
TypeTensor< T > * _tensor
T & operator()(const unsigned int i)
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition type_tensor.h:53
bool operator>(const TypeTensor< T > &rhs) const
const TypeTensor< T > & operator-=(const TypeTensor< T2 > &)
Subtract from this tensor.
bool operator<(const TypeTensor< T > &rhs) const
friend class TypeTensor
Definition type_tensor.h:55
std::enable_if< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
TypeTensor< T > inverse() const
auto norm_sq() const
const T & operator()(const unsigned int i, const unsigned int j) const
void subtract_scaled(const TypeTensor< T2 > &, const T &)
Subtract a scaled value from this tensor without creating a temporary.
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
Add to this tensor.
void add(const TypeTensor< T2 > &)
Add to this tensor without creating a temporary.
T value_type
Helper typedef for C++98 generic programming.
auto operator*(const Scalar &scalar) const -> typename std::enable_if< ScalarTraits< Scalar >::value, TypeTensor< decltype(T() *scalar)> >::type
Multiply this tensor by a scalar value.
void assign(const TypeTensor< T2 > &)
Assign to this tensor without creating a temporary.
const TypeTensor< T > & operator*=(const Scalar &factor)
Multiply this tensor by a scalar value in place.
std::enable_if< ScalarTraits< Scalar >::value, TypeTensor< typenameCompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
Divide each entry of this tensor by a scalar value.
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
const TypeTensor< T > & operator/=(const T &)
Divide each entry of this tensor by a scalar value.
std::tuple< unsigned int, unsigned int > index_type
Helper typedef for generic index programming.
void subtract(const TypeTensor< T2 > &)
Subtract from this tensor without creating a temporary.
TypeVector< T > column(const unsigned int r) const
void print(std::ostream &os=libMesh::out) const
Formatted print to a stream which defaults to libMesh::out.
void write_unformatted(std::ostream &out_stream, const bool newline=true) const
Unformatted print to the stream out.
Definition type_tensor.C:39
TypeVector< T > row(const unsigned int r) const
T _coords[LIBMESH_DIM *LIBMESH_DIM]
The coordinates of the TypeTensor.
TypeTensor< typename CompareTypes< T, T2 >::supertype > operator+(const TypeTensor< T2 > &) const
Add another tensor to this tensor.
friend std::ostream & operator<<(std::ostream &os, const TypeTensor< T > &t)
Formatted print as above but supports the syntax:
bool is_hpd(Real rel_tol=TOLERANCE *TOLERANCE) const
Returns true if the TypeTensor is Hermitian (or symmetric, when T==Real) and positive-definite to wit...
Definition type_tensor.C:67
bool operator==(const TypeTensor< T > &rhs) const
void zero()
Set all entries of the tensor to 0.
TypeVector< typename CompareTypes< T, T2 >::supertype > left_multiply(const TypeVector< T2 > &p) const
Left-multiply this tensor by a vector, i.e.
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
TypeTensor< T > transpose() const
~TypeTensor()
Destructor.
void solve(const TypeVector< T > &b, TypeVector< T > &x) const
Solve the 2x2 or 3x3 system of equations A*x = b for x by directly inverting A and multiplying it by ...
TypeTensor< T > operator-() const
ConstTypeTensorColumn< T > slice(const unsigned int i) const
bool is_zero() const
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition type_vector.h:63
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
static const Real b
auto norm_sq(const T &a)
The libMesh namespace provides an interface to certain functionality in the library.
void libmesh_ignore(const Args &...)
bool isfinite(std::complex< T > a)
bool isnan(std::complex< T > a)
OStreamProxy out
std::enable_if< ScalarTraits< Scalar >::value, TypeNTensor< N, typenameCompareTypes< Scalar, T >::supertype > >::type operator*(const Scalar &, const TypeNTensor< N, T > &)
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
static constexpr Real TOLERANCE
bool isinf(std::complex< T > a)
T libmesh_conj(T a)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
libMesh::TypeTensor< typename RawType< T >::value_type > value_type
static value_type value(const libMesh::TypeTensor< T > &in)
static const bool value