libmesh解析
本工作只是尝试解析原libmesh的代码,供学习使用
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 
type_tensor.h
浏览该文件的文档.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2023 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 
36 namespace libMesh
37 {
38 
39 // Forward declarations
40 template <typename T> class TypeTensorColumn;
41 template <typename T> class ConstTypeTensorColumn;
42 template <unsigned int N, typename T> class TypeNTensor;
43 
51 template <typename T>
52 class TypeTensor
53 {
54  template <typename T2>
55  friend class TypeTensor;
56 
57 public:
62  TypeTensor ();
63 
64 protected:
69  explicit TypeTensor (const T & xx,
70  const T & xy=0,
71  const T & xz=0,
72  const T & yx=0,
73  const T & yy=0,
74  const T & yz=0,
75  const T & zx=0,
76  const T & zy=0,
77  const T & zz=0);
78 
79 
83  template <typename Scalar>
84  explicit TypeTensor (const Scalar & xx,
85  const Scalar & xy=0,
86  const Scalar & xz=0,
87  const Scalar & yx=0,
88  const Scalar & yy=0,
89  const Scalar & yz=0,
90  const Scalar & zx=0,
91  const Scalar & zy=0,
92  typename
93  boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
94  const Scalar>::type & zz=0);
95 
100  template<typename T2>
101  TypeTensor(const TypeVector<T2> & vx);
102 
103  template<typename T2>
104  TypeTensor(const TypeVector<T2> & vx,
105  const TypeVector<T2> & vy);
106 
107  template<typename T2>
108  TypeTensor(const TypeVector<T2> & vx,
109  const TypeVector<T2> & vy,
110  const TypeVector<T2> & vz);
111 
115  typedef T value_type;
116 
120  typedef std::tuple<unsigned int, unsigned int> index_type;
121 
128  template<typename T2>
129  TypeTensor(const TypeTensor<T2> &p);
130 
134  ~TypeTensor();
135 
143  template <typename Scalar>
144  typename boostcopy::enable_if_c<
145  ScalarTraits<Scalar>::value,
146  TypeTensor &>::type
147  operator=(const Scalar &libmesh_dbg_var(p))
148  {
149  libmesh_assert_equal_to(p, Scalar(0));
150  this->zero();
151  return *this;
152  }
153 
161  const T &operator()(const unsigned int i, const unsigned int j) const;
162 
170  T &operator()(const unsigned int i, const unsigned int j);
171 
178  ConstTypeTensorColumn<T> slice(const unsigned int i) const;
179 
186  TypeTensorColumn<T> slice(const unsigned int i);
187 
194  TypeVector<T> row(const unsigned int r) const;
195 
203  template<typename T2>
205  operator+(const TypeTensor<T2> &) const;
206 
214  template<typename T2>
215  const TypeTensor<T> &operator+=(const TypeTensor<T2> &);
216 
223  template<typename T2>
224  void add(const TypeTensor<T2> &);
225 
233  template <typename T2>
234  void add_scaled(const TypeTensor<T2> &, const T &scale);
235 
243  template<typename T2>
245  operator-(const TypeTensor<T2> &) const;
246 
254  template<typename T2>
255  const TypeTensor<T> &operator-=(const TypeTensor<T2> &);
256 
263  template<typename T2>
264  void subtract(const TypeTensor<T2> &);
265 
273  template<typename T2>
274  void subtract_scaled(const TypeTensor<T2> &, const T &);
275 
281  TypeTensor<T> operator-() const;
282 
290  template <typename Scalar>
291  auto operator*(const Scalar &scalar) const -> typename boostcopy::enable_if_c<
292  ScalarTraits<Scalar>::value,
294 
302  template <typename Scalar, typename boostcopy::enable_if_c<
303  ScalarTraits<Scalar>::value, int>::type = 0>
304  const TypeTensor<T> &operator*=(const Scalar &factor)
305  {
306  for (unsigned int i = 0; i < LIBMESH_DIM * LIBMESH_DIM; i++)
307  _coords[i] *= factor;
308 
309  return *this;
310  }
311 
319  template <typename Scalar>
320  typename boostcopy::enable_if_c<
321  ScalarTraits<Scalar>::value,
323  operator/(const Scalar &) const;
324 
331  const TypeTensor<T> &operator/=(const T &);
332 
341  template <typename T2>
343  operator*(const TypeTensor<T2> &) const;
344 
352  template <typename T2>
353  const TypeTensor<T> &operator*=(const TypeTensor<T2> &);
354 
363  template <typename T2>
364  typename CompareTypes<T, T2>::supertype
365  contract(const TypeTensor<T2> &) const;
366 
375  template <typename T2>
377  operator*(const TypeVector<T2> &) const;
378 
387  template <typename T2>
389  left_multiply(const TypeVector<T2> &p) const;
390 
396  TypeTensor<T> transpose() const;
397 
403  TypeTensor<T> inverse() const;
404 
411  void solve(const TypeVector<T> &b, TypeVector<T> &x) const;
412 
418  auto norm() const -> decltype(std::norm(T()));
419 
425  auto norm_sq() const -> decltype(std::norm(T()));
426 
432  bool is_zero() const;
433 
440  T det() const;
441 
447  T tr() const;
448 
452  void zero();
453 
460  bool operator==(const TypeTensor<T> &rhs) const;
461 
469  bool operator<(const TypeTensor<T> &rhs) const;
470 
477  bool operator>(const TypeTensor<T> &rhs) const;
478 
484  void print(std::ostream &os = libMesh::out) const;
485 
492  void write_unformatted(std::ostream &out_stream,
493  const bool newline = true) const;
494 
495 protected:
496 
500  T _coords[LIBMESH_DIM*LIBMESH_DIM];
501 };
502 
503 
511 template <typename T>
512 class TypeTensorColumn
513 {
514 public:
521  TypeTensorColumn(TypeTensor<T> &tensor, unsigned int j) :
522  _tensor(&tensor), _j(j) {}
523 
530  T &operator()(const unsigned int i)
531  {
532  return (*_tensor)(i, _j);
533  }
534 
541  T &slice(const unsigned int i)
542  {
543  return (*_tensor)(i, _j);
544  }
545 
553  {
554  for (unsigned int i = 0; i != LIBMESH_DIM; ++i)
555  (*this)(i) = rhs(i);
556  return *this;
557  }
558 
559 private:
561  const unsigned int _j;
562 };
563 
571 template <typename T>
573 {
574 public:
581  ConstTypeTensorColumn(const TypeTensor<T> &tensor, unsigned int j) :
582  _tensor(&tensor), _j(j) {}
583 
590  const T &operator()(const unsigned int i) const
591  {
592  return (*_tensor)(i, _j);
593  }
594 
601  const T &slice(const unsigned int i) const
602  {
603  return (*_tensor)(i, _j);
604  }
605 
606 private:
608  const unsigned int _j;
609 };
610 
611 //------------------------------------------------------
612 // Inline functions
613 template <typename T>
614 inline
616 {
617  _coords[0] = {};
618 
619 #if LIBMESH_DIM > 1
620  _coords[1] = {};
621  _coords[2] = {};
622  _coords[3] = {};
623 #endif
624 
625 #if LIBMESH_DIM > 2
626  _coords[4] = {};
627  _coords[5] = {};
628  _coords[6] = {};
629  _coords[7] = {};
630  _coords[8] = {};
631 #endif
632 }
633 
634 
635 
636 template <typename T>
637 inline
639  const T & xy,
640  const T & xz,
641  const T & yx,
642  const T & yy,
643  const T & yz,
644  const T & zx,
645  const T & zy,
646  const T & zz)
647 {
648  _coords[0] = xx;
649 
650 #if LIBMESH_DIM == 2
651  _coords[1] = xy;
652  _coords[2] = yx;
653  _coords[3] = yy;
654 #elif LIBMESH_DIM == 1
655  libmesh_assert_equal_to (xy, 0);
656  libmesh_assert_equal_to (yx, 0);
657  libmesh_assert_equal_to (yy, 0);
658  libmesh_ignore(xy, yx, yy);
659 #endif
660 
661 #if LIBMESH_DIM == 3
662  _coords[1] = xy;
663  _coords[2] = xz;
664  _coords[3] = yx;
665  _coords[4] = yy;
666  _coords[5] = yz;
667  _coords[6] = zx;
668  _coords[7] = zy;
669  _coords[8] = zz;
670 #else
671  libmesh_assert_equal_to (xz, 0);
672  libmesh_assert_equal_to (yz, 0);
673  libmesh_assert_equal_to (zx, 0);
674  libmesh_assert_equal_to (zy, 0);
675  libmesh_assert_equal_to (zz, 0);
676  libmesh_ignore(xz, yz, zx, zy, zz);
677 #endif
678 }
679 
680 
681 template <typename T>
682 template <typename Scalar>
683 inline
684 TypeTensor<T>::TypeTensor (const Scalar & xx,
685  const Scalar & xy,
686  const Scalar & xz,
687  const Scalar & yx,
688  const Scalar & yy,
689  const Scalar & yz,
690  const Scalar & zx,
691  const Scalar & zy,
692  typename
693  boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
694  const Scalar>::type & zz)
695 {
696  _coords[0] = xx;
697 
698 #if LIBMESH_DIM == 2
699  _coords[1] = xy;
700  _coords[2] = yx;
701  _coords[3] = yy;
702 #elif LIBMESH_DIM == 1
703  libmesh_assert_equal_to (xy, 0);
704  libmesh_assert_equal_to (yx, 0);
705  libmesh_assert_equal_to (yy, 0);
706  libmesh_ignore(xy, yx, yy);
707 #endif
708 
709 #if LIBMESH_DIM == 3
710  _coords[1] = xy;
711  _coords[2] = xz;
712  _coords[3] = yx;
713  _coords[4] = yy;
714  _coords[5] = yz;
715  _coords[6] = zx;
716  _coords[7] = zy;
717  _coords[8] = zz;
718 #else
719  libmesh_assert_equal_to (xz, 0);
720  libmesh_assert_equal_to (yz, 0);
721  libmesh_assert_equal_to (zx, 0);
722  libmesh_assert_equal_to (zy, 0);
723  libmesh_assert_equal_to (zz, 0);
724  libmesh_ignore(xz, yz, zx, zy, zz);
725 #endif
726 }
727 
728 
729 
730 template <typename T>
731 template<typename T2>
732 inline
734 {
735  // copy the nodes from vector p to me
736  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
737  _coords[i] = p._coords[i];
738 }
739 
740 
741 template <typename T>
742 template <typename T2>
744 {
745  libmesh_assert_equal_to (LIBMESH_DIM, 1);
746  _coords[0] = vx(0);
747 }
748 
749 template <typename T>
750 template <typename T2>
752  const TypeVector<T2> & vy)
753 {
754  libmesh_assert_equal_to (LIBMESH_DIM, 2);
755 #if LIBMESH_DIM > 2
756  _coords[0] = vx(0);
757  _coords[1] = vx(1);
758  _coords[2] = vy(0);
759  _coords[3] = vy(1);
760 #else
761  libmesh_ignore(vx, vy);
762 #endif
763 }
764 
765 template <typename T>
766 template <typename T2>
768  const TypeVector<T2> & vy,
769  const TypeVector<T2> & vz)
770 {
771  libmesh_assert_equal_to (LIBMESH_DIM, 3);
772 #if LIBMESH_DIM > 2
773  _coords[0] = vx(0);
774  _coords[1] = vx(1);
775  _coords[2] = vx(2);
776  _coords[3] = vy(0);
777  _coords[4] = vy(1);
778  _coords[5] = vy(2);
779  _coords[6] = vz(0);
780  _coords[7] = vz(1);
781  _coords[8] = vz(2);
782 #else
783  libmesh_ignore(vx, vy, vz);
784 #endif
785 }
786 
787 
788 
789 
790 template <typename T>
791 inline
793 {
794 }
795 
796 
797 
798 template <typename T>
799 template<typename T2>
800 inline
801 void TypeTensor<T>::assign (const TypeTensor<T2> & p)
802 {
803  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
804  _coords[i] = p._coords[i];
805 }
806 
807 
808 
809 template <typename T>
810 inline
811 const T & TypeTensor<T>::operator () (const unsigned int i,
812  const unsigned int j) const
813 {
814  libmesh_assert_less (i, 3);
815  libmesh_assert_less (j, 3);
816 
817 #if LIBMESH_DIM < 3
818  const static T my_zero = 0;
819  if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
820  return my_zero;
821 #endif
822 
823  return _coords[i*LIBMESH_DIM+j];
824 }
825 
826 
827 
828 template <typename T>
829 inline
830 T & TypeTensor<T>::operator () (const unsigned int i,
831  const unsigned int j)
832 {
833 #if LIBMESH_DIM < 3
834 
835  libmesh_error_msg_if(i >= LIBMESH_DIM || j >= LIBMESH_DIM,
836  "ERROR: You are assigning to a tensor component that is out of range for the compiled LIBMESH_DIM!");
837 
838 #endif
839 
840  libmesh_assert_less (i, LIBMESH_DIM);
841  libmesh_assert_less (j, LIBMESH_DIM);
842 
843  return _coords[i*LIBMESH_DIM+j];
844 }
845 
846 
847 template <typename T>
848 inline
850 TypeTensor<T>::slice (const unsigned int i) const
851 {
852  libmesh_assert_less (i, LIBMESH_DIM);
853  return ConstTypeTensorColumn<T>(*this, i);
854 }
855 
856 
857 template <typename T>
858 inline
860 TypeTensor<T>::slice (const unsigned int i)
861 {
862  libmesh_assert_less (i, LIBMESH_DIM);
863  return TypeTensorColumn<T>(*this, i);
864 }
865 
866 
867 template <typename T>
868 inline
870 TypeTensor<T>::row(const unsigned int r) const
871 {
872  TypeVector<T> return_vector;
873 
874  for (unsigned int j=0; j<LIBMESH_DIM; j++)
875  return_vector._coords[j] = _coords[r*LIBMESH_DIM + j];
876 
877  return return_vector;
878 }
879 
880 template <typename T>
881 template<typename T2>
882 inline
885 {
886 
887 #if LIBMESH_DIM == 1
889 #endif
890 
891 #if LIBMESH_DIM == 2
893  _coords[1] + p._coords[1],
894  0.,
895  _coords[2] + p._coords[2],
896  _coords[3] + p._coords[3]);
897 #endif
898 
899 #if LIBMESH_DIM == 3
901  _coords[1] + p._coords[1],
902  _coords[2] + p._coords[2],
903  _coords[3] + p._coords[3],
904  _coords[4] + p._coords[4],
905  _coords[5] + p._coords[5],
906  _coords[6] + p._coords[6],
907  _coords[7] + p._coords[7],
908  _coords[8] + p._coords[8]);
909 #endif
910 
911 }
912 
913 
914 
915 template <typename T>
916 template<typename T2>
917 inline
919 {
920  this->add (p);
921 
922  return *this;
923 }
924 
925 
926 
927 template <typename T>
928 template<typename T2>
929 inline
931 {
932  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
933  _coords[i] += p._coords[i];
934 }
935 
936 
937 
938 template <typename T>
939 template <typename T2>
940 inline
941 void TypeTensor<T>::add_scaled (const TypeTensor<T2> & p, const T & factor)
942 {
943  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
944  _coords[i] += factor*p._coords[i];
945 
946 }
947 
948 
949 
950 template <typename T>
951 template<typename T2>
952 inline
955 {
956 
957 #if LIBMESH_DIM == 1
959 #endif
960 
961 #if LIBMESH_DIM == 2
963  _coords[1] - p._coords[1],
964  0.,
965  _coords[2] - p._coords[2],
966  _coords[3] - p._coords[3]);
967 #endif
968 
969 #if LIBMESH_DIM == 3
971  _coords[1] - p._coords[1],
972  _coords[2] - p._coords[2],
973  _coords[3] - p._coords[3],
974  _coords[4] - p._coords[4],
975  _coords[5] - p._coords[5],
976  _coords[6] - p._coords[6],
977  _coords[7] - p._coords[7],
978  _coords[8] - p._coords[8]);
979 #endif
980 
981 }
982 
983 
984 
985 template <typename T>
986 template <typename T2>
987 inline
989 {
990  this->subtract (p);
991 
992  return *this;
993 }
994 
995 
996 
997 template <typename T>
998 template <typename T2>
999 inline
1001 {
1002  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1003  _coords[i] -= p._coords[i];
1004 }
1005 
1006 
1007 
1008 template <typename T>
1009 template <typename T2>
1010 inline
1011 void TypeTensor<T>::subtract_scaled (const TypeTensor<T2> & p, const T & factor)
1012 {
1013  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1014  _coords[i] -= factor*p._coords[i];
1015 }
1016 
1017 
1018 
1019 template <typename T>
1020 inline
1022 {
1023 
1024 #if LIBMESH_DIM == 1
1025  return TypeTensor(-_coords[0]);
1026 #endif
1027 
1028 #if LIBMESH_DIM == 2
1029  return TypeTensor(-_coords[0],
1030  -_coords[1],
1031  -_coords[2],
1032  -_coords[3]);
1033 #endif
1034 
1035 #if LIBMESH_DIM == 3
1036  return TypeTensor(-_coords[0],
1037  -_coords[1],
1038  -_coords[2],
1039  -_coords[3],
1040  -_coords[4],
1041  -_coords[5],
1042  -_coords[6],
1043  -_coords[7],
1044  -_coords[8]);
1045 #endif
1046 
1047 }
1048 
1049 
1050 
1051 template <typename T>
1052 template <typename Scalar>
1053 inline
1054 auto
1055 TypeTensor<T>::operator * (const Scalar & factor) const -> typename boostcopy::enable_if_c<
1056  ScalarTraits<Scalar>::value,
1058 {
1059  typedef decltype((*this)(0, 0) * factor) TS;
1060 
1061 
1062 #if LIBMESH_DIM == 1
1063  return TypeTensor<TS>(_coords[0]*factor);
1064 #endif
1065 
1066 #if LIBMESH_DIM == 2
1067  return TypeTensor<TS>(_coords[0]*factor,
1068  _coords[1]*factor,
1069  _coords[2]*factor,
1070  _coords[3]*factor);
1071 #endif
1072 
1073 #if LIBMESH_DIM == 3
1074  return TypeTensor<TS>(_coords[0]*factor,
1075  _coords[1]*factor,
1076  _coords[2]*factor,
1077  _coords[3]*factor,
1078  _coords[4]*factor,
1079  _coords[5]*factor,
1080  _coords[6]*factor,
1081  _coords[7]*factor,
1082  _coords[8]*factor);
1083 #endif
1084 }
1085 
1086 
1087 
1088 template <typename T, typename Scalar>
1089 inline
1090 typename boostcopy::enable_if_c<
1091  ScalarTraits<Scalar>::value,
1093 operator * (const Scalar & factor,
1094  const TypeTensor<T> & t)
1095 {
1096  return t * factor;
1097 }
1098 
1099 template <typename T>
1100 template <typename Scalar>
1101 inline
1102 typename boostcopy::enable_if_c<
1103  ScalarTraits<Scalar>::value,
1104  TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1105 TypeTensor<T>::operator / (const Scalar & factor) const
1106 {
1107  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1108 
1109  typedef typename CompareTypes<T, Scalar>::supertype TS;
1110 
1111 #if LIBMESH_DIM == 1
1112  return TypeTensor<TS>(_coords[0]/factor);
1113 #endif
1114 
1115 #if LIBMESH_DIM == 2
1116  return TypeTensor<TS>(_coords[0]/factor,
1117  _coords[1]/factor,
1118  _coords[2]/factor,
1119  _coords[3]/factor);
1120 #endif
1121 
1122 #if LIBMESH_DIM == 3
1123  return TypeTensor<TS>(_coords[0]/factor,
1124  _coords[1]/factor,
1125  _coords[2]/factor,
1126  _coords[3]/factor,
1127  _coords[4]/factor,
1128  _coords[5]/factor,
1129  _coords[6]/factor,
1130  _coords[7]/factor,
1131  _coords[8]/factor);
1132 #endif
1133 
1134 }
1135 
1136 
1137 
1138 template <typename T>
1139 inline
1141 {
1142 #if LIBMESH_DIM == 1
1143  return TypeTensor(_coords[0]);
1144 #endif
1145 
1146 #if LIBMESH_DIM == 2
1147  return TypeTensor(_coords[0],
1148  _coords[2],
1149  _coords[1],
1150  _coords[3]);
1151 #endif
1152 
1153 #if LIBMESH_DIM == 3
1154  return TypeTensor(_coords[0],
1155  _coords[3],
1156  _coords[6],
1157  _coords[1],
1158  _coords[4],
1159  _coords[7],
1160  _coords[2],
1161  _coords[5],
1162  _coords[8]);
1163 #endif
1164 }
1165 
1166 
1167 
1168 template <typename T>
1169 inline
1171 {
1172 #if LIBMESH_DIM == 1
1173  if (_coords[0] == static_cast<T>(0.))
1174  libmesh_convergence_failure();
1175  return TypeTensor(1. / _coords[0]);
1176 #endif
1177 
1178 #if LIBMESH_DIM == 2
1179  // Get convenient reference to this.
1180  const TypeTensor<T> & A = *this;
1181 
1182  // Use temporary variables, avoid multiple accesses
1183  T a = A(0,0), b = A(0,1),
1184  c = A(1,0), d = A(1,1);
1185 
1186  // Make sure det = ad - bc is not zero
1187  T my_det = a*d - b*c;
1188 
1189  if (my_det == static_cast<T>(0.))
1190  libmesh_convergence_failure();
1191 
1192  return TypeTensor(d/my_det, -b/my_det, -c/my_det, a/my_det);
1193 #endif
1194 
1195 #if LIBMESH_DIM == 3
1196  // Get convenient reference to this.
1197  const TypeTensor<T> & A = *this;
1198 
1199  T a11 = A(0,0), a12 = A(0,1), a13 = A(0,2),
1200  a21 = A(1,0), a22 = A(1,1), a23 = A(1,2),
1201  a31 = A(2,0), a32 = A(2,1), a33 = A(2,2);
1202 
1203  T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1204 
1205  if (my_det == static_cast<T>(0.))
1206  libmesh_convergence_failure();
1207 
1208  // Inline comment characters are for lining up columns.
1209  return TypeTensor( (a33*a22-a32*a23)/my_det, -(a33*a12-a32*a13)/my_det, (a23*a12-a22*a13)/my_det,
1210  -(a33*a21-a31*a23)/my_det, (a33*a11-a31*a13)/my_det, -(a23*a11-a21*a13)/my_det,
1211  (a32*a21-a31*a22)/my_det, -(a32*a11-a31*a12)/my_det, (a22*a11-a21*a12)/my_det);
1212 #endif
1213 }
1214 
1215 
1216 
1217 template <typename T>
1218 inline
1220 {
1221 #if LIBMESH_DIM == 1
1222  if (_coords[0] == static_cast<T>(0.))
1223  libmesh_convergence_failure();
1224  x(0) = b(0) / _coords[0];
1225 #endif
1226 
1227 #if LIBMESH_DIM == 2
1228  T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1229 
1230  if (my_det == static_cast<T>(0.))
1231  libmesh_convergence_failure();
1232 
1233  T my_det_inv = 1./my_det;
1234 
1235  x(0) = my_det_inv*( _coords[3]*b(0) - _coords[1]*b(1));
1236  x(1) = my_det_inv*(-_coords[2]*b(0) + _coords[0]*b(1));
1237 #endif
1238 
1239 #if LIBMESH_DIM == 3
1240  T my_det =
1241  // a11*(a33 *a22 - a32 *a23)
1242  _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1243  // -a21*(a33 *a12 - a32 *a13)
1244  -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1245  // +a31*(a23 *a12 - a22 *a13)
1246  +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1247 
1248  if (my_det == static_cast<T>(0.))
1249  libmesh_convergence_failure();
1250 
1251  T my_det_inv = 1./my_det;
1252  x(0) = my_det_inv*((_coords[8]*_coords[4] - _coords[7]*_coords[5])*b(0) -
1253  (_coords[8]*_coords[1] - _coords[7]*_coords[2])*b(1) +
1254  (_coords[5]*_coords[1] - _coords[4]*_coords[2])*b(2));
1255 
1256  x(1) = my_det_inv*((_coords[6]*_coords[5] - _coords[8]*_coords[3])*b(0) +
1257  (_coords[8]*_coords[0] - _coords[6]*_coords[2])*b(1) -
1258  (_coords[5]*_coords[0] - _coords[3]*_coords[2])*b(2));
1259 
1260  x(2) = my_det_inv*((_coords[7]*_coords[3] - _coords[6]*_coords[4])*b(0) -
1261  (_coords[7]*_coords[0] - _coords[6]*_coords[1])*b(1) +
1262  (_coords[4]*_coords[0] - _coords[3]*_coords[1])*b(2));
1263 #endif
1264 }
1265 
1266 
1267 
1268 template <typename T>
1269 inline
1271 {
1272  libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1273 
1274  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1275  _coords[i] /= factor;
1276 
1277  return *this;
1278 }
1279 
1280 
1281 
1282 
1283 template <typename T>
1284 template <typename T2>
1285 inline
1288 {
1290  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1291  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1292  returnval(i) += (*this)(i,j)*p(j);
1293 
1294  return returnval;
1295 }
1296 
1297 template <typename T>
1298 template <typename T2>
1299 inline
1302 {
1304  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1305  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1306  returnval(i) += p(j)*(*this)(j,i);
1307 
1308  return returnval;
1309 }
1310 
1311 template <typename T, typename T2>
1312 inline
1315 {
1316  return b.left_multiply(a);
1317 }
1318 
1319 template <typename T>
1320 template <typename T2>
1321 inline
1322 TypeTensor<typename CompareTypes<T, T2>::supertype>
1324 {
1326  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1327  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1328  for (unsigned int k=0; k<LIBMESH_DIM; k++)
1329  returnval(i,j) += (*this)(i,k)*p(k,j);
1330 
1331  return returnval;
1332 }
1333 
1334 template <typename T>
1335 template <typename T2>
1336 inline
1338 {
1339  TypeTensor<T> temp;
1340  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1341  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1342  for (unsigned int k=0; k<LIBMESH_DIM; k++)
1343  temp(i,j) += (*this)(i,k)*p(k,j);
1344 
1345  this->assign(temp);
1346  return *this;
1347 }
1348 
1349 
1354 template <typename T>
1355 template <typename T2>
1356 inline
1357 typename CompareTypes<T,T2>::supertype
1359 {
1360  typename CompareTypes<T,T2>::supertype sum = 0.;
1361  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1362  sum += _coords[i]*t._coords[i];
1363  return sum;
1364 }
1365 
1366 
1367 
1368 template <typename T>
1369 inline
1370 auto TypeTensor<T>::norm() const -> decltype(std::norm(T()))
1371 {
1372  return std::sqrt(this->norm_sq());
1373 }
1374 
1375 
1376 template <typename T>
1377 inline
1379 {
1380  for (const auto & val : _coords)
1381  if (val != T(0))
1382  return false;
1383  return true;
1384 }
1385 
1386 template <typename T>
1387 inline
1389 {
1390 #if LIBMESH_DIM == 1
1391  return _coords[0];
1392 #endif
1393 
1394 #if LIBMESH_DIM == 2
1395  return (_coords[0] * _coords[3]
1396  - _coords[1] * _coords[2]);
1397 #endif
1398 
1399 #if LIBMESH_DIM == 3
1400  return (_coords[0] * _coords[4] * _coords[8]
1401  + _coords[1] * _coords[5] * _coords[6]
1402  + _coords[2] * _coords[3] * _coords[7]
1403  - _coords[0] * _coords[5] * _coords[7]
1404  - _coords[1] * _coords[3] * _coords[8]
1405  - _coords[2] * _coords[4] * _coords[6]);
1406 #endif
1407 }
1408 
1409 template <typename T>
1410 inline
1412 {
1413 #if LIBMESH_DIM == 1
1414  return _coords[0];
1415 #endif
1416 
1417 #if LIBMESH_DIM == 2
1418  return _coords[0] + _coords[3];
1419 #endif
1420 
1421 #if LIBMESH_DIM == 3
1422  return _coords[0] + _coords[4] + _coords[8];
1423 #endif
1424 }
1425 
1426 template <typename T>
1427 inline
1429 {
1430  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1431  _coords[i] = 0.;
1432 }
1433 
1434 
1435 
1436 template <typename T>
1437 inline
1438 auto TypeTensor<T>::norm_sq () const -> decltype(std::norm(T()))
1439 {
1440  Real sum = 0.;
1441  for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1442  sum += TensorTools::norm_sq(_coords[i]);
1443  return sum;
1444 }
1445 
1446 
1447 
1448 template <typename T>
1449 inline
1451 {
1452 #if LIBMESH_DIM == 1
1453  return (std::abs(_coords[0] - rhs._coords[0])
1454  < TOLERANCE);
1455 #endif
1456 
1457 #if LIBMESH_DIM == 2
1458  return ((std::abs(_coords[0] - rhs._coords[0]) +
1459  std::abs(_coords[1] - rhs._coords[1]) +
1460  std::abs(_coords[2] - rhs._coords[2]) +
1461  std::abs(_coords[3] - rhs._coords[3]))
1462  < 4.*TOLERANCE);
1463 #endif
1464 
1465 #if LIBMESH_DIM == 3
1466  return ((std::abs(_coords[0] - rhs._coords[0]) +
1467  std::abs(_coords[1] - rhs._coords[1]) +
1468  std::abs(_coords[2] - rhs._coords[2]) +
1469  std::abs(_coords[3] - rhs._coords[3]) +
1470  std::abs(_coords[4] - rhs._coords[4]) +
1471  std::abs(_coords[5] - rhs._coords[5]) +
1472  std::abs(_coords[6] - rhs._coords[6]) +
1473  std::abs(_coords[7] - rhs._coords[7]) +
1474  std::abs(_coords[8] - rhs._coords[8]))
1475  < 9.*TOLERANCE);
1476 #endif
1477 
1478 }
1479 
1480 template <typename T, typename T2>
1481 inline
1484 {
1486  for (unsigned int i=0; i<LIBMESH_DIM; i++)
1487  for (unsigned int j=0; j<LIBMESH_DIM; j++)
1488  ret(i,j) = a(i) * libmesh_conj(b(j));
1489 
1490  return ret;
1491 }
1492 
1493 
1494 } // namespace libMesh
1495 
1496 #ifdef LIBMESH_HAVE_METAPHYSICL
1497 namespace MetaPhysicL
1498 {
1499 template <typename T>
1500 struct RawType<libMesh::TypeTensor<T>>
1501 {
1503 
1505  {
1506  value_type ret;
1507  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1508  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1509  ret(i,j) = raw_value(in(i,j));
1510 
1511  return ret;
1512  }
1513 };
1514 
1515 template <typename T, typename U>
1516 struct ReplaceAlgebraicType<libMesh::TypeTensor<T>, U>
1517 {
1518  typedef U type;
1519 };
1520 }
1521 #endif
1522 
1523 #endif // LIBMESH_TYPE_TENSOR_H
TypeTensorColumn< T > & operator=(const TypeVector< T > &rhs)
将该列的值设置为给定向量的值。
Definition: type_tensor.h:552
T _coords[LIBMESH_DIM]
TypeVector 的坐标。
Definition: type_vector.h:565
bool operator==(const TypeTensor< T > &rhs) const
检查两个 Tensor 是否相等。
Definition: type_tensor.h:1450
const unsigned int _j
Definition: type_tensor.h:561
const TypeTensor< T > & operator/=(const T &)
将该 Tensor 的每个元素除以标量值。
Definition: type_tensor.h:1270
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
Definition: type_tensor.h:1483
static value_type value(const libMesh::TypeTensor< T > &in)
Definition: type_tensor.h:1504
T libmesh_conj(T a)
TypeTensor< T > * _tensor
Definition: type_tensor.h:560
T & operator()(const unsigned int i)
返回该张量的第 行第 列元素的可写引用。
Definition: type_tensor.h:530
void add_scaled(const TypeTensor< T2 > &, const T &scale)
在不创建临时 Tensor 的情况下将一个缩放的 Tensor 加到该 Tensor 上。
Definition: type_tensor.h:941
T value_type
Helper typedef 用于 C++98 泛型编程。
Definition: type_tensor.h:115
static constexpr Real TOLERANCE
std::tuple< unsigned int, unsigned int > index_type
Helper typedef 用于泛型索引编程。
Definition: type_tensor.h:120
TypeTensor< typename CompareTypes< T, T2 >::supertype > operator+(const TypeTensor< T2 > &) const
将另一个 Tensor 加到该 Tensor 上。
Definition: type_tensor.h:884
void zero()
将 Tensor 的所有元素设置为零。
Definition: type_tensor.h:1428
void write_unformatted(std::ostream &out_stream, const bool newline=true) const
无格式输出到流,默认为打印元素并以空格和换行符分隔。
Definition: type_tensor.C:83
TypeTensor()
Empty constructor.
Definition: type_tensor.h:615
auto norm_sq() const -> decltype(std::norm(T()))
返回 Tensor 的 Frobenius 范数的平方,即元素模平方的和。
Definition: type_tensor.h:1438
const T & slice(const unsigned int i) const
返回该张量的第 行第 列元素的只读引用。
Definition: type_tensor.h:601
ConstTypeTensorColumn(const TypeTensor< T > &tensor, unsigned int j)
构造函数,初始化 ConstTypeTensorColumn。
Definition: type_tensor.h:581
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
计算自动微分实数向量的平方根。
Definition: type_vector.h:88
T det() const
返回 Tensor 的行列式。 由于这些是最多 3x3 的 Tensor,我们不像 DenseMatrix 那样做 LU 分解。
Definition: type_tensor.h:1388
T _coords[LIBMESH_DIM *LIBMESH_DIM]
TypeTensor 的坐标
Definition: type_tensor.h:500
const TypeTensor< T > * _tensor
Definition: type_tensor.h:607
void print(std::ostream &os=libMesh::out) const
格式化输出到流,默认为 libMesh::out。
Definition: type_tensor.C:39
libMesh::TypeTensor< typename RawType< T >::value_type > value_type
Definition: type_tensor.h:1502
ConstTypeTensorColumn< T > slice(const unsigned int i) const
返回 Tensor 的第 i 列的代理。
Definition: type_tensor.h:850
auto norm() const -> decltype(std::norm(T()))
返回 Tensor 的 Frobenius 范数,即元素平方和的平方根。
Definition: type_tensor.h:1370
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:36
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
计算自动微分实数向量的绝对值。
Definition: type_vector.h:112
T norm_sq(std::complex< T > a)
Definition: tensor_tools.h:74
TypeTensor< T > transpose() const
返回该 Tensor 的转置(对于复数不共轭)。
Definition: type_tensor.h:1140
表示 ConstTypeTensorColumn 类,用于访问 TypeTensor 的列并进行只读操作。
Definition: type_tensor.h:41
void libmesh_ignore(const Args &...)
TypeTensorColumn(TypeTensor< T > &tensor, unsigned int j)
构造函数,初始化 TypeTensorColumn。
Definition: type_tensor.h:521
void subtract_scaled(const TypeTensor< T2 > &, const T &)
在不创建临时 Tensor 的情况下从该 Tensor 减去一个缩放的 Tensor。
Definition: type_tensor.h:1011
auto operator*(const Scalar &scalar) const -> typename boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< decltype(T()*scalar)>>::type
将该 Tensor 与标量值相乘。
Definition: type_tensor.h:1055
TypeTensor< T > operator-() const
返回该 Tensor 的负值的副本。
Definition: type_tensor.h:1021
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< typename CompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
将该 Tensor 的每个元素除以标量值。
Definition: type_tensor.h:1105
void subtract(const TypeTensor< T2 > &)
在不创建临时 Tensor 的情况下从该 Tensor 减去另一个 Tensor。
Definition: type_tensor.h:1000
该类定义了一个在 LIBMESH_DIM 维度空间中类型为 T 的向量。
Definition: tensor_tools.h:34
OStreamProxy out
void add(const TypeTensor< T2 > &)
在不创建临时 Tensor 的情况下将另一个 Tensor 加到该 Tensor 上。
Definition: type_tensor.h:930
~TypeTensor()
析构函数。
Definition: type_tensor.h:792
TypeVector< T > row(const unsigned int r) const
返回 Tensor 的一行作为 TypeVector 的副本。
Definition: type_tensor.h:870
TypeTensor< T > inverse() const
返回该 Tensor 的逆作为独立对象。
Definition: type_tensor.h:1170
void solve(const TypeVector< T > &b, TypeVector< T > &x) const
解方程组 得到 ,其中 为该 Tensor。
Definition: type_tensor.h:1219
bool is_zero() const
如果 Tensor 中的所有值都为零,则返回 true。
Definition: type_tensor.h:1378
const T & operator()(const unsigned int i) const
返回该张量的第 行第 列元素的只读引用。
Definition: type_tensor.h:590
const T & operator()(const unsigned int i, const unsigned int j) const
返回 (i,j) 元素的常量引用。
Definition: type_tensor.h:811
该类最终将定义一个在类型为T的LIBMESH_DIM维空间中的N阶张量。
Definition: tensor_tools.h:38
TypeVector< typename CompareTypes< T, T2 >::supertype > left_multiply(const TypeVector< T2 > &p) const
将该 Tensor 左乘以一个向量,即矩阵-向量乘积。 这个 Tensor 和向量可以包含不同的数值类型。
Definition: type_tensor.h:1301
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
将两个 Tensor 相乘,返回一个标量,即 。 这些 Tensor 可能包含不同的数值类型。
Definition: type_tensor.h:1358
表示 TypeTensorColumn 类,用于访问 TypeTensor 的列并进行写操作。
Definition: type_tensor.h:40
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
在不创建临时 Tensor 的情况下对该 Tensor 进行赋值。
Definition: type_tensor.h:147
const TypeTensor< T > & operator-=(const TypeTensor< T2 > &)
从该 Tensor 减去另一个 Tensor。
Definition: type_tensor.h:988
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
在该 Tensor 上加上另一个 Tensor。
Definition: type_tensor.h:918
T & slice(const unsigned int i)
返回该张量的第 行第 列元素的可写引用。
Definition: type_tensor.h:541
const TypeTensor< T > & operator*=(const Scalar &factor)
在地方上将该 Tensor 乘以标量值。
Definition: type_tensor.h:304
T tr() const
返回 Tensor 的迹。
Definition: type_tensor.h:1411