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/libmesh_device.h"
26#include "libmesh/type_vector.h"
27
28// C++ includes
29#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
30#include <cmath>
31#include <tuple>
32
33#ifdef LIBMESH_HAVE_METAPHYSICL
34#include "metaphysicl/dualnumber_decl.h"
35#endif
36
37namespace libMesh
38{
39
40// Forward declarations
41template <typename T> class TypeTensorColumn;
42template <typename T> class ConstTypeTensorColumn;
43template <unsigned int N, typename T> class TypeNTensor;
44
52template <typename T>
54{
55 template <typename T2>
56 friend class TypeTensor;
57
58public:
64
65protected:
72 explicit TypeTensor (const T & xx,
73 const T & xy=0,
74 const T & xz=0,
75 const T & yx=0,
76 const T & yy=0,
77 const T & yz=0,
78 const T & zx=0,
79 const T & zy=0,
80 const T & zz=0);
81
82
86 template <typename Scalar>
87 explicit TypeTensor (const Scalar & xx,
88 const Scalar & xy=0,
89 const Scalar & xz=0,
90 const Scalar & yx=0,
91 const Scalar & yy=0,
92 const Scalar & yz=0,
93 const Scalar & zx=0,
94 const Scalar & zy=0,
95 typename
96 std::enable_if<ScalarTraits<Scalar>::value,
97 const Scalar>::type & zz=0);
98
104 template<typename T2>
105 LIBMESH_DEVICE_INLINE
107
108 template<typename T2>
109 LIBMESH_DEVICE_INLINE
111 const TypeVector<T2> & vy);
112
113 template<typename T2>
114 LIBMESH_DEVICE_INLINE
116 const TypeVector<T2> & vy,
117 const TypeVector<T2> & vz);
118
119public:
120
124 typedef T value_type;
125
129 typedef std::tuple<unsigned int, unsigned int> index_type;
130
134 template<typename T2>
136
140 LIBMESH_DEVICE_INLINE
142
146 template<typename T2>
147 LIBMESH_DEVICE_INLINE
148 void assign (const TypeTensor<T2> &);
149
155 template <typename Scalar>
156 LIBMESH_DEVICE_INLINE
157 typename std::enable_if<
159 TypeTensor &>::type
160 operator = (const Scalar & libmesh_dbg_var(p))
161 { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }
162
166 const T & operator () (const unsigned int i, const unsigned int j) const;
167
171 T & operator () (const unsigned int i, const unsigned int j);
172
176 LIBMESH_DEVICE_INLINE
177 ConstTypeTensorColumn<T> slice (const unsigned int i) const;
178
182 LIBMESH_DEVICE_INLINE
183 TypeTensorColumn<T> slice (const unsigned int i);
184
188 TypeVector<T> row(const unsigned int r) const;
189
193 LIBMESH_DEVICE_INLINE
194 TypeVector<T> column(const unsigned int r) const;
195
201 template<typename T2>
204
210 template<typename T2>
212
216 template<typename T2>
217 void add (const TypeTensor<T2> &);
218
222 template <typename T2>
223 LIBMESH_DEVICE_INLINE
224 void add_scaled (const TypeTensor<T2> &, const T &);
225
231 template<typename T2>
234
240 template<typename T2>
242
246 template<typename T2>
247 void subtract (const TypeTensor<T2> &);
248
253 template <typename T2>
254 LIBMESH_DEVICE_INLINE
255 void subtract_scaled (const TypeTensor<T2> &, const T &);
256
261
267 template <typename Scalar>
268 auto
269 operator * (const Scalar & scalar) const -> typename std::enable_if<
271 TypeTensor<decltype(T() * scalar)>>::type;
272
278 template <typename Scalar, typename std::enable_if<
279 ScalarTraits<Scalar>::value, int>::type = 0>
280 LIBMESH_DEVICE_INLINE
281 const TypeTensor<T> & operator *= (const Scalar & factor)
282 {
283 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
284 _coords[i] *= factor;
285
286 return *this;
287 }
288
294 template <typename Scalar>
295 typename std::enable_if<
298 operator / (const Scalar &) const;
299
305 const TypeTensor<T> & operator /= (const T &);
306
313 template <typename T2>
316
322 template <typename T2>
324
334 template <typename T2>
335 LIBMESH_DEVICE_INLINE
337 contract (const TypeTensor<T2> &) const;
338
345 template <typename T2>
348
355 template <typename T2>
356 LIBMESH_DEVICE_INLINE
358 left_multiply (const TypeVector<T2> & p) const;
359
364
369
376 LIBMESH_DEVICE_INLINE
377 void solve(const TypeVector<T> & b, TypeVector<T> & x) const;
378
383 auto norm() const;
384
389 auto norm_sq() const;
390
394 LIBMESH_DEVICE_INLINE
395 bool is_zero() const;
396
403 T det() const;
404
408 T tr() const;
409
413 LIBMESH_DEVICE_INLINE
414 void zero();
415
419 LIBMESH_DEVICE_INLINE
420 bool operator == (const TypeTensor<T> & rhs) const;
421
426 bool operator < (const TypeTensor<T> & rhs) const;
427
431 bool operator > (const TypeTensor<T> & rhs) const;
432
436 void print(std::ostream & os = libMesh::out) const;
437
445 friend std::ostream & operator << (std::ostream & os, const TypeTensor<T> & t)
446 {
447 t.print(os);
448 return os;
449 }
450
455 void write_unformatted (std::ostream & out_stream,
456 const bool newline = true) const;
457
463 bool is_hpd(Real rel_tol = TOLERANCE*TOLERANCE) const;
464
465protected:
466
470 T _coords[LIBMESH_DIM*LIBMESH_DIM];
471};
472
473
474template <typename T>
476{
477public:
478
480 unsigned int j) :
481 _tensor(&tensor), _j(j) {}
482
487 T & operator () (const unsigned int i)
488 { return (*_tensor)(i,_j); }
489
490 T & slice (const unsigned int i)
491 { return (*_tensor)(i,_j); }
492
497 {
498 for (unsigned int i=0; i != LIBMESH_DIM; ++i)
499 (*this)(i) = rhs(i);
500 return *this;
501 }
502
503private:
505 const unsigned int _j;
506};
507
508
509template <typename T>
511{
512public:
513
515 unsigned int j) :
516 _tensor(&tensor), _j(j) {}
517
521 const T & operator () (const unsigned int i) const
522 { return (*_tensor)(i,_j); }
523
524 const T & slice (const unsigned int i) const
525 { return (*_tensor)(i,_j); }
526
527private:
529 const unsigned int _j;
530};
531
532//------------------------------------------------------
533// Inline functions
534template <typename T>
535LIBMESH_DEVICE_INLINE
537{
538 _coords[0] = {};
539
540#if LIBMESH_DIM > 1
541 _coords[1] = {};
542 _coords[2] = {};
543 _coords[3] = {};
544#endif
545
546#if LIBMESH_DIM > 2
547 _coords[4] = {};
548 _coords[5] = {};
549 _coords[6] = {};
550 _coords[7] = {};
551 _coords[8] = {};
552#endif
553}
554
555
556
557template <typename T>
558LIBMESH_DEVICE_INLINE
560 const T & xy,
561 const T & xz,
562 const T & yx,
563 const T & yy,
564 const T & yz,
565 const T & zx,
566 const T & zy,
567 const T & zz)
568{
569 _coords[0] = xx;
570
571#if LIBMESH_DIM == 2
572 _coords[1] = xy;
573 _coords[2] = yx;
574 _coords[3] = yy;
575#elif LIBMESH_DIM == 1
576 libmesh_assert_equal_to (xy, 0);
577 libmesh_assert_equal_to (yx, 0);
578 libmesh_assert_equal_to (yy, 0);
579 libmesh_ignore(xy, yx, yy);
580#endif
581
582#if LIBMESH_DIM == 3
583 _coords[1] = xy;
584 _coords[2] = xz;
585 _coords[3] = yx;
586 _coords[4] = yy;
587 _coords[5] = yz;
588 _coords[6] = zx;
589 _coords[7] = zy;
590 _coords[8] = zz;
591#else
592 libmesh_assert_equal_to (xz, 0);
593 libmesh_assert_equal_to (yz, 0);
594 libmesh_assert_equal_to (zx, 0);
595 libmesh_assert_equal_to (zy, 0);
596 libmesh_assert_equal_to (zz, 0);
597 libmesh_ignore(xz, yz, zx, zy, zz);
598#endif
599}
600
601
602template <typename T>
603template <typename Scalar>
604LIBMESH_DEVICE_INLINE
605TypeTensor<T>::TypeTensor (const Scalar & xx,
606 const Scalar & xy,
607 const Scalar & xz,
608 const Scalar & yx,
609 const Scalar & yy,
610 const Scalar & yz,
611 const Scalar & zx,
612 const Scalar & zy,
613 typename
614 std::enable_if<ScalarTraits<Scalar>::value,
615 const Scalar>::type & zz)
616{
617 _coords[0] = xx;
618
619#if LIBMESH_DIM == 2
620 _coords[1] = xy;
621 _coords[2] = yx;
622 _coords[3] = yy;
623#elif LIBMESH_DIM == 1
624 libmesh_assert_equal_to (xy, 0);
625 libmesh_assert_equal_to (yx, 0);
626 libmesh_assert_equal_to (yy, 0);
627 libmesh_ignore(xy, yx, yy);
628#endif
629
630#if LIBMESH_DIM == 3
631 _coords[1] = xy;
632 _coords[2] = xz;
633 _coords[3] = yx;
634 _coords[4] = yy;
635 _coords[5] = yz;
636 _coords[6] = zx;
637 _coords[7] = zy;
638 _coords[8] = zz;
639#else
640 libmesh_assert_equal_to (xz, 0);
641 libmesh_assert_equal_to (yz, 0);
642 libmesh_assert_equal_to (zx, 0);
643 libmesh_assert_equal_to (zy, 0);
644 libmesh_assert_equal_to (zz, 0);
645 libmesh_ignore(xz, yz, zx, zy, zz);
646#endif
647}
648
649
650
651template <typename T>
652template<typename T2>
653LIBMESH_DEVICE_INLINE
654TypeTensor<T>::TypeTensor (const TypeTensor <T2> & p)
655{
656 // copy the nodes from vector p to me
657 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
658 _coords[i] = p._coords[i];
659}
660
661
662template <typename T>
663template <typename T2>
664LIBMESH_DEVICE_INLINE
666{
667 libmesh_assert_equal_to (LIBMESH_DIM, 1);
668 _coords[0] = vx(0);
669}
670
671template <typename T>
672template <typename T2>
673LIBMESH_DEVICE_INLINE
675 const TypeVector<T2> & vy)
676{
677 libmesh_assert_equal_to (LIBMESH_DIM, 2);
678#if LIBMESH_DIM > 1
679 _coords[0] = vx(0);
680 _coords[1] = vx(1);
681 _coords[2] = vy(0);
682 _coords[3] = vy(1);
683#else
684 libmesh_ignore(vx, vy);
685#endif
686}
687
688template <typename T>
689template <typename T2>
690LIBMESH_DEVICE_INLINE
692 const TypeVector<T2> & vy,
693 const TypeVector<T2> & vz)
694{
695 libmesh_assert_equal_to (LIBMESH_DIM, 3);
696#if LIBMESH_DIM > 2
697 _coords[0] = vx(0);
698 _coords[1] = vx(1);
699 _coords[2] = vx(2);
700 _coords[3] = vy(0);
701 _coords[4] = vy(1);
702 _coords[5] = vy(2);
703 _coords[6] = vz(0);
704 _coords[7] = vz(1);
705 _coords[8] = vz(2);
706#else
707 libmesh_ignore(vx, vy, vz);
708#endif
709}
710
711
712
713
714template <typename T>
715LIBMESH_DEVICE_INLINE
719
720
721
722template <typename T>
723template<typename T2>
724LIBMESH_DEVICE_INLINE
726{
727 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
728 _coords[i] = p._coords[i];
729}
730
731
732
733template <typename T>
734LIBMESH_DEVICE_INLINE
735const T & TypeTensor<T>::operator () (const unsigned int i,
736 const unsigned int j) const
737{
738 libmesh_assert_less (i, 3);
739 libmesh_assert_less (j, 3);
740
741#if LIBMESH_DIM < 3
742 const static T my_zero = 0;
743 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
744 return my_zero;
745#endif
746
747 return _coords[i*LIBMESH_DIM+j];
748}
749
750
751
752template <typename T>
753LIBMESH_DEVICE_INLINE
754T & TypeTensor<T>::operator () (const unsigned int i,
755 const unsigned int j)
756{
757#if LIBMESH_DIM < 3
758
759 libmesh_error_msg_if(i >= LIBMESH_DIM || j >= LIBMESH_DIM,
760 "ERROR: You are assigning to a tensor component that is out of range for the compiled LIBMESH_DIM!");
761
762#endif
763
764 libmesh_assert_less (i, LIBMESH_DIM);
765 libmesh_assert_less (j, LIBMESH_DIM);
766
767 return _coords[i*LIBMESH_DIM+j];
768}
769
770
771template <typename T>
772LIBMESH_DEVICE_INLINE
774TypeTensor<T>::slice (const unsigned int i) const
775{
776 libmesh_assert_less (i, LIBMESH_DIM);
777 return ConstTypeTensorColumn<T>(*this, i);
778}
779
780
781template <typename T>
782LIBMESH_DEVICE_INLINE
784TypeTensor<T>::slice (const unsigned int i)
785{
786 libmesh_assert_less (i, LIBMESH_DIM);
787 return TypeTensorColumn<T>(*this, i);
788}
789
790
791template <typename T>
792LIBMESH_DEVICE_INLINE
794TypeTensor<T>::row(const unsigned int r) const
795{
796 TypeVector<T> return_vector;
797
798 for (unsigned int j=0; j<LIBMESH_DIM; j++)
799 return_vector._coords[j] = _coords[r*LIBMESH_DIM + j];
800
801 return return_vector;
802}
803
804
805template <typename T>
806LIBMESH_DEVICE_INLINE
808TypeTensor<T>::column(const unsigned int r) const
809{
810 TypeVector<T> return_vector;
811
812 for (unsigned int i=0; i<LIBMESH_DIM; i++)
813 return_vector._coords[i] = _coords[i*LIBMESH_DIM + r];
814
815 return return_vector;
816}
817
818
819template <typename T>
820template<typename T2>
821LIBMESH_DEVICE_INLINE
824{
825
826#if LIBMESH_DIM == 1
828#endif
829
830#if LIBMESH_DIM == 2
832 _coords[1] + p._coords[1],
833 0.,
834 _coords[2] + p._coords[2],
835 _coords[3] + p._coords[3]);
836#endif
837
838#if LIBMESH_DIM == 3
840 _coords[1] + p._coords[1],
841 _coords[2] + p._coords[2],
842 _coords[3] + p._coords[3],
843 _coords[4] + p._coords[4],
844 _coords[5] + p._coords[5],
845 _coords[6] + p._coords[6],
846 _coords[7] + p._coords[7],
847 _coords[8] + p._coords[8]);
848#endif
849
850}
851
852
853
854template <typename T>
855template<typename T2>
856LIBMESH_DEVICE_INLINE
858{
859 this->add (p);
860
861 return *this;
862}
863
864
865
866template <typename T>
867template<typename T2>
868LIBMESH_DEVICE_INLINE
870{
871 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
872 _coords[i] += p._coords[i];
873}
874
875
876
877template <typename T>
878template <typename T2>
879LIBMESH_DEVICE_INLINE
880void TypeTensor<T>::add_scaled (const TypeTensor<T2> & p, const T & factor)
881{
882 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
883 _coords[i] += factor*p._coords[i];
884
885}
886
887
888
889template <typename T>
890template<typename T2>
891LIBMESH_DEVICE_INLINE
894{
895
896#if LIBMESH_DIM == 1
898#endif
899
900#if LIBMESH_DIM == 2
902 _coords[1] - p._coords[1],
903 0.,
904 _coords[2] - p._coords[2],
905 _coords[3] - p._coords[3]);
906#endif
907
908#if LIBMESH_DIM == 3
910 _coords[1] - p._coords[1],
911 _coords[2] - p._coords[2],
912 _coords[3] - p._coords[3],
913 _coords[4] - p._coords[4],
914 _coords[5] - p._coords[5],
915 _coords[6] - p._coords[6],
916 _coords[7] - p._coords[7],
917 _coords[8] - p._coords[8]);
918#endif
919
920}
921
922
923
924template <typename T>
925template <typename T2>
926LIBMESH_DEVICE_INLINE
928{
929 this->subtract (p);
930
931 return *this;
932}
933
934
935
936template <typename T>
937template <typename T2>
938LIBMESH_DEVICE_INLINE
940{
941 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
942 _coords[i] -= p._coords[i];
943}
944
945
946
947template <typename T>
948template <typename T2>
949LIBMESH_DEVICE_INLINE
950void TypeTensor<T>::subtract_scaled (const TypeTensor<T2> & p, const T & factor)
951{
952 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
953 _coords[i] -= factor*p._coords[i];
954}
955
956
957
958template <typename T>
959LIBMESH_DEVICE_INLINE
961{
962
963#if LIBMESH_DIM == 1
964 return TypeTensor(-_coords[0]);
965#endif
966
967#if LIBMESH_DIM == 2
968 return TypeTensor(-_coords[0],
969 -_coords[1],
970 -_coords[2],
971 -_coords[3]);
972#endif
973
974#if LIBMESH_DIM == 3
975 return TypeTensor(-_coords[0],
976 -_coords[1],
977 -_coords[2],
978 -_coords[3],
979 -_coords[4],
980 -_coords[5],
981 -_coords[6],
982 -_coords[7],
983 -_coords[8]);
984#endif
985
986}
987
988
989
990template <typename T>
991template <typename Scalar>
992LIBMESH_DEVICE_INLINE
993auto
994TypeTensor<T>::operator * (const Scalar & factor) const -> typename std::enable_if<
996 TypeTensor<decltype(T() * factor)>>::type
997{
998 typedef decltype((*this)(0, 0) * factor) TS;
999
1000
1001#if LIBMESH_DIM == 1
1002 return TypeTensor<TS>(_coords[0]*factor);
1003#endif
1004
1005#if LIBMESH_DIM == 2
1006 return TypeTensor<TS>(_coords[0]*factor,
1007 _coords[1]*factor,
1008 _coords[2]*factor,
1009 _coords[3]*factor);
1010#endif
1011
1012#if LIBMESH_DIM == 3
1013 return TypeTensor<TS>(_coords[0]*factor,
1014 _coords[1]*factor,
1015 _coords[2]*factor,
1016 _coords[3]*factor,
1017 _coords[4]*factor,
1018 _coords[5]*factor,
1019 _coords[6]*factor,
1020 _coords[7]*factor,
1021 _coords[8]*factor);
1022#endif
1023}
1024
1025
1026
1027template <typename T, typename Scalar>
1028LIBMESH_DEVICE_INLINE
1029typename std::enable_if<
1032operator * (const Scalar & factor,
1033 const TypeTensor<T> & t)
1034{
1035 return t * factor;
1036}
1037
1038template <typename T>
1039template <typename Scalar>
1040LIBMESH_DEVICE_INLINE
1041typename std::enable_if<
1043 TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1044TypeTensor<T>::operator / (const Scalar & factor) const
1045{
1046 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1047
1048 typedef typename CompareTypes<T, Scalar>::supertype TS;
1049
1050#if LIBMESH_DIM == 1
1051 return TypeTensor<TS>(_coords[0]/factor);
1052#endif
1053
1054#if LIBMESH_DIM == 2
1055 return TypeTensor<TS>(_coords[0]/factor,
1056 _coords[1]/factor,
1057 _coords[2]/factor,
1058 _coords[3]/factor);
1059#endif
1060
1061#if LIBMESH_DIM == 3
1062 return TypeTensor<TS>(_coords[0]/factor,
1063 _coords[1]/factor,
1064 _coords[2]/factor,
1065 _coords[3]/factor,
1066 _coords[4]/factor,
1067 _coords[5]/factor,
1068 _coords[6]/factor,
1069 _coords[7]/factor,
1070 _coords[8]/factor);
1071#endif
1072
1073}
1074
1075
1076
1077template <typename T>
1078LIBMESH_DEVICE_INLINE
1080{
1081#if LIBMESH_DIM == 1
1082 return TypeTensor(_coords[0]);
1083#endif
1084
1085#if LIBMESH_DIM == 2
1086 return TypeTensor(_coords[0],
1087 _coords[2],
1088 _coords[1],
1089 _coords[3]);
1090#endif
1091
1092#if LIBMESH_DIM == 3
1093 return TypeTensor(_coords[0],
1094 _coords[3],
1095 _coords[6],
1096 _coords[1],
1097 _coords[4],
1098 _coords[7],
1099 _coords[2],
1100 _coords[5],
1101 _coords[8]);
1102#endif
1103}
1104
1105
1106
1107template <typename T>
1108LIBMESH_DEVICE_INLINE
1110{
1111#if LIBMESH_DIM == 1
1112 if (_coords[0] == static_cast<T>(0.))
1113 libmesh_convergence_failure();
1114 return TypeTensor(1. / _coords[0]);
1115#endif
1116
1117#if LIBMESH_DIM == 2
1118 // Get convenient reference to this.
1119 const TypeTensor<T> & A = *this;
1120
1121 // Use temporary variables, avoid multiple accesses
1122 T a = A(0,0), b = A(0,1),
1123 c = A(1,0), d = A(1,1);
1124
1125 // Make sure det = ad - bc is not zero
1126 T my_det = a*d - b*c;
1127
1128 if (my_det == static_cast<T>(0.))
1129 libmesh_convergence_failure();
1130
1131 return TypeTensor(d/my_det, -b/my_det, -c/my_det, a/my_det);
1132#endif
1133
1134#if LIBMESH_DIM == 3
1135 // Get convenient reference to this.
1136 const TypeTensor<T> & A = *this;
1137
1138 T a11 = A(0,0), a12 = A(0,1), a13 = A(0,2),
1139 a21 = A(1,0), a22 = A(1,1), a23 = A(1,2),
1140 a31 = A(2,0), a32 = A(2,1), a33 = A(2,2);
1141
1142 T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1143
1144 if (my_det == static_cast<T>(0.))
1145 libmesh_convergence_failure();
1146
1147 // Inline comment characters are for lining up columns.
1148 return TypeTensor( (a33*a22-a32*a23)/my_det, -(a33*a12-a32*a13)/my_det, (a23*a12-a22*a13)/my_det,
1149 -(a33*a21-a31*a23)/my_det, (a33*a11-a31*a13)/my_det, -(a23*a11-a21*a13)/my_det,
1150 (a32*a21-a31*a22)/my_det, -(a32*a11-a31*a12)/my_det, (a22*a11-a21*a12)/my_det);
1151#endif
1152}
1153
1154
1155
1156template <typename T>
1157LIBMESH_DEVICE_INLINE
1159{
1160#if LIBMESH_DIM == 1
1161 if (_coords[0] == static_cast<T>(0.))
1162 libmesh_convergence_failure();
1163 x(0) = b(0) / _coords[0];
1164#endif
1165
1166#if LIBMESH_DIM == 2
1167 T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1168
1169 if (my_det == static_cast<T>(0.))
1170 libmesh_convergence_failure();
1171
1172 T my_det_inv = 1./my_det;
1173
1174 x(0) = my_det_inv*( _coords[3]*b(0) - _coords[1]*b(1));
1175 x(1) = my_det_inv*(-_coords[2]*b(0) + _coords[0]*b(1));
1176#endif
1177
1178#if LIBMESH_DIM == 3
1179 T my_det =
1180 // a11*(a33 *a22 - a32 *a23)
1181 _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1182 // -a21*(a33 *a12 - a32 *a13)
1183 -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1184 // +a31*(a23 *a12 - a22 *a13)
1185 +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1186
1187 if (my_det == static_cast<T>(0.))
1188 libmesh_convergence_failure();
1189
1190 T my_det_inv = 1./my_det;
1191 x(0) = my_det_inv*((_coords[8]*_coords[4] - _coords[7]*_coords[5])*b(0) -
1192 (_coords[8]*_coords[1] - _coords[7]*_coords[2])*b(1) +
1193 (_coords[5]*_coords[1] - _coords[4]*_coords[2])*b(2));
1194
1195 x(1) = my_det_inv*((_coords[6]*_coords[5] - _coords[8]*_coords[3])*b(0) +
1196 (_coords[8]*_coords[0] - _coords[6]*_coords[2])*b(1) -
1197 (_coords[5]*_coords[0] - _coords[3]*_coords[2])*b(2));
1198
1199 x(2) = my_det_inv*((_coords[7]*_coords[3] - _coords[6]*_coords[4])*b(0) -
1200 (_coords[7]*_coords[0] - _coords[6]*_coords[1])*b(1) +
1201 (_coords[4]*_coords[0] - _coords[3]*_coords[1])*b(2));
1202#endif
1203}
1204
1205
1206
1207template <typename T>
1208LIBMESH_DEVICE_INLINE
1210{
1211 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1212
1213 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1214 _coords[i] /= factor;
1215
1216 return *this;
1217}
1218
1219
1220
1221
1222template <typename T>
1223template <typename T2>
1224LIBMESH_DEVICE_INLINE
1227{
1229 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1230 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1231 returnval(i) += (*this)(i,j)*p(j);
1232
1233 return returnval;
1234}
1235
1236template <typename T>
1237template <typename T2>
1238LIBMESH_DEVICE_INLINE
1241{
1243 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1244 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1245 returnval(i) += p(j)*(*this)(j,i);
1246
1247 return returnval;
1248}
1249
1250template <typename T, typename T2>
1251LIBMESH_DEVICE_INLINE
1254{
1255 return b.left_multiply(a);
1256}
1257
1258template <typename T>
1259template <typename T2>
1260LIBMESH_DEVICE_INLINE
1261TypeTensor<typename CompareTypes<T, T2>::supertype>
1263{
1265 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1266 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1267 for (unsigned int k=0; k<LIBMESH_DIM; k++)
1268 returnval(i,j) += (*this)(i,k)*p(k,j);
1269
1270 return returnval;
1271}
1272
1273template <typename T>
1274template <typename T2>
1275LIBMESH_DEVICE_INLINE
1277{
1278 TypeTensor<T> temp;
1279 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1280 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1281 for (unsigned int k=0; k<LIBMESH_DIM; k++)
1282 temp(i,j) += (*this)(i,k)*p(k,j);
1283
1284 this->assign(temp);
1285 return *this;
1286}
1287
1288
1293template <typename T>
1294template <typename T2>
1295LIBMESH_DEVICE_INLINE
1298{
1299 typename CompareTypes<T,T2>::supertype sum = 0.;
1300 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1301 sum += _coords[i]*t._coords[i];
1302 return sum;
1303}
1304
1305
1306
1307template <typename T>
1308LIBMESH_DEVICE_INLINE
1310{
1311 using std::sqrt;
1312 return sqrt(this->norm_sq());
1313}
1314
1315
1316template <typename T>
1317LIBMESH_DEVICE_INLINE
1319{
1320 for (const auto & val : _coords)
1321 if (val != T(0))
1322 return false;
1323 return true;
1324}
1325
1326template <typename T>
1327LIBMESH_DEVICE_INLINE
1329{
1330#if LIBMESH_DIM == 1
1331 return _coords[0];
1332#endif
1333
1334#if LIBMESH_DIM == 2
1335 return (_coords[0] * _coords[3]
1336 - _coords[1] * _coords[2]);
1337#endif
1338
1339#if LIBMESH_DIM == 3
1340 return (_coords[0] * _coords[4] * _coords[8]
1341 + _coords[1] * _coords[5] * _coords[6]
1342 + _coords[2] * _coords[3] * _coords[7]
1343 - _coords[0] * _coords[5] * _coords[7]
1344 - _coords[1] * _coords[3] * _coords[8]
1345 - _coords[2] * _coords[4] * _coords[6]);
1346#endif
1347}
1348
1349template <typename T>
1350LIBMESH_DEVICE_INLINE
1352{
1353#if LIBMESH_DIM == 1
1354 return _coords[0];
1355#endif
1356
1357#if LIBMESH_DIM == 2
1358 return _coords[0] + _coords[3];
1359#endif
1360
1361#if LIBMESH_DIM == 3
1362 return _coords[0] + _coords[4] + _coords[8];
1363#endif
1364}
1365
1366template <typename T>
1367LIBMESH_DEVICE_INLINE
1369{
1370 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1371 _coords[i] = 0.;
1372}
1373
1374
1375
1376template <typename T>
1377LIBMESH_DEVICE_INLINE
1379{
1380 Real sum = 0.;
1381 for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1382 sum += TensorTools::norm_sq(_coords[i]);
1383 return sum;
1384}
1385
1386
1387
1388template <typename T>
1389LIBMESH_DEVICE_INLINE
1391{
1392#if LIBMESH_DIM == 1
1393 return (std::abs(_coords[0] - rhs._coords[0])
1394 < TOLERANCE);
1395#endif
1396
1397#if LIBMESH_DIM == 2
1398 return ((std::abs(_coords[0] - rhs._coords[0]) +
1399 std::abs(_coords[1] - rhs._coords[1]) +
1400 std::abs(_coords[2] - rhs._coords[2]) +
1401 std::abs(_coords[3] - rhs._coords[3]))
1402 < 4.*TOLERANCE);
1403#endif
1404
1405#if LIBMESH_DIM == 3
1406 return ((std::abs(_coords[0] - rhs._coords[0]) +
1407 std::abs(_coords[1] - rhs._coords[1]) +
1408 std::abs(_coords[2] - rhs._coords[2]) +
1409 std::abs(_coords[3] - rhs._coords[3]) +
1410 std::abs(_coords[4] - rhs._coords[4]) +
1411 std::abs(_coords[5] - rhs._coords[5]) +
1412 std::abs(_coords[6] - rhs._coords[6]) +
1413 std::abs(_coords[7] - rhs._coords[7]) +
1414 std::abs(_coords[8] - rhs._coords[8]))
1415 < 9.*TOLERANCE);
1416#endif
1417
1418}
1419
1420template <typename T>
1421void TypeTensor<T>::print(std::ostream & os) const
1422{
1423#if LIBMESH_DIM == 1
1424
1425 os << "x=" << (*this)(0,0) << std::endl;
1426
1427#endif
1428#if LIBMESH_DIM == 2
1429
1430 os << "(xx,xy)=("
1431 << std::setw(8) << (*this)(0,0) << ", "
1432 << std::setw(8) << (*this)(0,1) << ")"
1433 << std::endl;
1434 os << "(yx,yy)=("
1435 << std::setw(8) << (*this)(1,0) << ", "
1436 << std::setw(8) << (*this)(1,1) << ")"
1437 << std::endl;
1438
1439#endif
1440#if LIBMESH_DIM == 3
1441
1442 os << "(xx,xy,xz)=("
1443 << std::setw(8) << (*this)(0,0) << ", "
1444 << std::setw(8) << (*this)(0,1) << ", "
1445 << std::setw(8) << (*this)(0,2) << ")"
1446 << std::endl;
1447 os << "(yx,yy,yz)=("
1448 << std::setw(8) << (*this)(1,0) << ", "
1449 << std::setw(8) << (*this)(1,1) << ", "
1450 << std::setw(8) << (*this)(1,2) << ")"
1451 << std::endl;
1452 os << "(zx,zy,zz)=("
1453 << std::setw(8) << (*this)(2,0) << ", "
1454 << std::setw(8) << (*this)(2,1) << ", "
1455 << std::setw(8) << (*this)(2,2) << ")"
1456 << std::endl;
1457#endif
1458}
1459
1460
1461// A matrix is finite iff every component is
1462template <typename T>
1463bool isfinite (const TypeTensor<T> & var)
1464{
1465 using std::isfinite;
1466 using libMesh::isfinite; // for T==complex
1467 for (auto i : make_range(LIBMESH_DIM))
1468 for (auto j : make_range(LIBMESH_DIM))
1469 if (!isfinite(var(i,j)))
1470 return false;
1471 return true;
1472}
1473
1474
1475// A matrix is infinite iff some component is infinite but no
1476// component is NaN.
1477//
1478// This is arguably inconsistent with our std::complex overload (and
1479// the C99 Annex G recommendations for _Complex, and C++ std::complex
1480// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1481// this is probably safer for users.
1482template <typename T>
1483bool isinf (const TypeTensor<T> & var)
1484{
1485 using std::isinf;
1486 using libMesh::isinf; // for T==complex
1487 using std::isnan;
1488 using libMesh::isnan;
1489 bool has_inf = false;
1490 for (auto i : make_range(LIBMESH_DIM))
1491 for (auto j : make_range(LIBMESH_DIM))
1492 {
1493 // NaN anywhere makes us NaN, not inf
1494 if (isnan(var(i,j)))
1495 return false;
1496 has_inf = has_inf || isinf(var(i,j));
1497 }
1498 return has_inf;
1499}
1500
1501
1502
1503// A matrix is NaN iff some component is NaN
1504//
1505// This is arguably inconsistent with our std::complex overload (and
1506// the C99 Annex G recommendations for _Complex, and C++ std::complex
1507// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
1508// this is probably safer for users.
1509template <typename T>
1510bool isnan (const TypeTensor<T> & var)
1511{
1512 using std::isnan;
1513 using libMesh::isnan; // for T==complex
1514 for (auto i : make_range(LIBMESH_DIM))
1515 for (auto j : make_range(LIBMESH_DIM))
1516 if (isnan(var(i,j)))
1517 return true;
1518 return false;
1519}
1520
1521
1522
1523template <typename T, typename T2>
1524LIBMESH_DEVICE_INLINE
1525TypeTensor<typename CompareTypes<T, T2>::supertype>
1527{
1529 for (unsigned int i=0; i<LIBMESH_DIM; i++)
1530 for (unsigned int j=0; j<LIBMESH_DIM; j++)
1531 ret(i,j) = a(i) * libmesh_conj(b(j));
1532
1533 return ret;
1534}
1535
1536
1537} // namespace libMesh
1538
1539#ifdef LIBMESH_HAVE_METAPHYSICL
1540namespace MetaPhysicL
1541{
1542template <typename T>
1543struct RawType<libMesh::TypeTensor<T>>
1544{
1546
1548 {
1549 value_type ret;
1550 for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1551 for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1552 ret(i,j) = raw_value(in(i,j));
1553
1554 return ret;
1555 }
1556};
1557
1558template <typename T, typename U>
1559struct ReplaceAlgebraicType<libMesh::TypeTensor<T>, U>
1560{
1561 typedef U type;
1562};
1563}
1564#endif
1565
1566#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:54
bool operator>(const TypeTensor< T > &rhs) const
const TypeTensor< T > & operator-=(const TypeTensor< T2 > &)
Subtract from this tensor.
LIBMESH_DEVICE_INLINE bool operator==(const TypeTensor< T > &rhs) const
bool operator<(const TypeTensor< T > &rhs) const
LIBMESH_DEVICE_INLINE ConstTypeTensorColumn< T > slice(const unsigned int i) const
friend class TypeTensor
Definition type_tensor.h:56
LIBMESH_DEVICE_INLINE TypeTensorColumn< T > slice(const unsigned int i)
LIBMESH_DEVICE_INLINE void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
TypeTensor< T > inverse() const
LIBMESH_DEVICE_INLINE TypeTensor(const TypeVector< T2 > &vx, const TypeVector< T2 > &vy, const TypeVector< T2 > &vz)
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
Add to this tensor.
void add(const TypeTensor< T2 > &)
Add to this tensor without creating a temporary.
LIBMESH_DEVICE_INLINE TypeTensor(const TypeVector< T2 > &vx, const TypeVector< T2 > &vy)
TypeTensor(const T &xx, const T &xy=0, const T &xz=0, const T &yx=0, const T &yy=0, const T &yz=0, const T &zx=0, const T &zy=0, const T &zz=0)
Constructor-from-T.
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.
LIBMESH_DEVICE_INLINE void assign(const TypeTensor< T2 > &)
Assign to this tensor without creating a temporary.
auto norm_sq() const
TypeTensor< T > transpose() const
TypeTensor(const TypeTensor< T2 > &p)
Copy-constructor.
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.
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.
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
LIBMESH_DEVICE_INLINE ~TypeTensor()
Destructor.
TypeTensor(const Scalar &xx, const Scalar &xy=0, const Scalar &xz=0, const Scalar &yx=0, const Scalar &yy=0, const Scalar &yz=0, const Scalar &zx=0, const Scalar &zy=0, typename std::enable_if< ScalarTraits< Scalar >::value, const Scalar >::type &zz=0)
Constructor-from-Scalar.
LIBMESH_DEVICE_INLINE TypeTensor(const TypeVector< T2 > &vx)
Constructor.
TypeTensor()
Empty constructor.
const T & operator()(const unsigned int i, const unsigned int j) const
LIBMESH_DEVICE_INLINE const TypeTensor< T > & operator*=(const Scalar &factor)
Multiply this tensor by a scalar value in place.
const TypeTensor< T > & operator/=(const T &)
Divide each entry of this tensor by a scalar value.
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.
LIBMESH_DEVICE_INLINE void subtract_scaled(const TypeTensor< T2 > &, const T &)
Subtract a scaled value from this tensor without creating a temporary.
friend std::ostream & operator<<(std::ostream &os, const TypeTensor< T > &t)
Formatted print as above but supports the syntax:
TypeTensor< T > operator-() const
LIBMESH_DEVICE_INLINE bool is_zero() const
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
LIBMESH_DEVICE_INLINE 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 ...
LIBMESH_DEVICE_INLINE CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
TypeVector< T > row(const unsigned int r) const
LIBMESH_DEVICE_INLINE TypeVector< typename CompareTypes< T, T2 >::supertype > left_multiply(const TypeVector< T2 > &p) const
Left-multiply this tensor by a vector, i.e.
LIBMESH_DEVICE_INLINE void zero()
Set all entries of the tensor to 0.
LIBMESH_DEVICE_INLINE TypeVector< T > column(const unsigned int r) const
LIBMESH_DEVICE_INLINE std::enable_if< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition type_vector.h:64
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
static const Real b
LIBMESH_DEVICE_INLINE auto norm_sq(const T &a)
The libMesh namespace provides an interface to certain functionality in the library.
LIBMESH_DEVICE_INLINE bool isfinite(std::complex< T > a)
LIBMESH_DEVICE_INLINE bool isnan(std::complex< T > a)
void libmesh_ignore(const Args &...)
OStreamProxy out
std::enable_if< ScalarTraits< Scalar >::value, TypeNTensor< N, typenameCompareTypes< Scalar, T >::supertype > >::type operator*(const Scalar &, const TypeNTensor< N, T > &)
LIBMESH_DEVICE_INLINE TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
LIBMESH_DEVICE_INLINE bool isinf(std::complex< T > a)
static constexpr Real TOLERANCE
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_DEVICE_INLINE T libmesh_conj(T a)
libMesh::TypeTensor< typename RawType< T >::value_type > value_type
static value_type value(const libMesh::TypeTensor< T > &in)
static const bool value