20 #ifndef LIBMESH_TYPE_TENSOR_H
21 #define LIBMESH_TYPE_TENSOR_H
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/type_vector.h"
32 #ifdef LIBMESH_HAVE_METAPHYSICL
33 #include "metaphysicl/dualnumber_decl.h"
42 template <
unsigned int N,
typename T>
class TypeNTensor;
54 template <
typename T2>
83 template <
typename Scalar>
93 boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
94 const Scalar>::type & zz=0);
100 template<
typename T2>
103 template<
typename T2>
107 template<
typename T2>
128 template<
typename T2>
143 template <
typename Scalar>
144 typename boostcopy::enable_if_c<
145 ScalarTraits<Scalar>::value,
149 libmesh_assert_equal_to(p, Scalar(0));
161 const T &
operator()(
const unsigned int i,
const unsigned int j)
const;
170 T &
operator()(
const unsigned int i,
const unsigned int j);
203 template<
typename T2>
214 template<
typename T2>
223 template<
typename T2>
233 template <
typename T2>
243 template<
typename T2>
254 template<
typename T2>
263 template<
typename T2>
273 template<
typename T2>
290 template <
typename Scalar>
291 auto operator*(
const Scalar &scalar)
const ->
typename boostcopy::enable_if_c<
292 ScalarTraits<Scalar>::value,
302 template <
typename Scalar,
typename boostcopy::enable_if_c<
303 ScalarTraits<Scalar>::value,
int>::type = 0>
306 for (
unsigned int i = 0; i < LIBMESH_DIM * LIBMESH_DIM; i++)
319 template <
typename Scalar>
320 typename boostcopy::enable_if_c<
321 ScalarTraits<Scalar>::value,
341 template <
typename T2>
352 template <
typename T2>
363 template <
typename T2>
364 typename CompareTypes<T, T2>::supertype
375 template <
typename T2>
387 template <
typename T2>
418 auto norm() const -> decltype(std::
norm(T()));
460 bool operator==(const
TypeTensor<T> &rhs) const;
469 bool operator<(const
TypeTensor<T> &rhs) const;
477 bool operator>(const
TypeTensor<T> &rhs) const;
484 void print(std::ostream &os = libMesh::
out) const;
493 const
bool newline = true) const;
511 template <typename T>
522 _tensor(&tensor), _j(j) {}
532 return (*_tensor)(i, _j);
543 return (*_tensor)(i, _j);
554 for (
unsigned int i = 0; i != LIBMESH_DIM; ++i)
561 const unsigned int _j;
571 template <
typename T>
582 _tensor(&tensor), _j(j) {}
592 return (*_tensor)(i, _j);
601 const T &
slice(
const unsigned int i)
const
603 return (*_tensor)(i, _j);
608 const unsigned int _j;
613 template <
typename T>
636 template <
typename T>
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);
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);
681 template <
typename T>
682 template <
typename Scalar>
693 boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
694 const Scalar>::type & zz)
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);
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);
730 template <
typename T>
731 template<
typename T2>
736 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
741 template <
typename T>
742 template <
typename T2>
745 libmesh_assert_equal_to (LIBMESH_DIM, 1);
749 template <
typename T>
750 template <
typename T2>
754 libmesh_assert_equal_to (LIBMESH_DIM, 2);
765 template <
typename T>
766 template <
typename T2>
771 libmesh_assert_equal_to (LIBMESH_DIM, 3);
790 template <
typename T>
798 template <
typename T>
799 template<
typename T2>
803 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
809 template <
typename T>
812 const unsigned int j)
const
814 libmesh_assert_less (i, 3);
815 libmesh_assert_less (j, 3);
818 const static T my_zero = 0;
819 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
823 return _coords[i*LIBMESH_DIM+j];
828 template <
typename T>
831 const unsigned int j)
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!");
840 libmesh_assert_less (i, LIBMESH_DIM);
841 libmesh_assert_less (j, LIBMESH_DIM);
843 return _coords[i*LIBMESH_DIM+j];
847 template <
typename T>
852 libmesh_assert_less (i, LIBMESH_DIM);
857 template <
typename T>
862 libmesh_assert_less (i, LIBMESH_DIM);
867 template <
typename T>
874 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
877 return return_vector;
880 template <
typename T>
881 template<
typename T2>
915 template <
typename T>
916 template<
typename T2>
927 template <
typename T>
928 template<
typename T2>
932 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
938 template <
typename T>
939 template <
typename T2>
943 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
950 template <
typename T>
951 template<
typename T2>
985 template <
typename T>
986 template <
typename T2>
997 template <
typename T>
998 template <
typename T2>
1002 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1008 template <
typename T>
1009 template <
typename T2>
1013 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1019 template <
typename T>
1024 #if LIBMESH_DIM == 1
1028 #if LIBMESH_DIM == 2
1035 #if LIBMESH_DIM == 3
1051 template <
typename T>
1052 template <
typename Scalar>
1056 ScalarTraits<Scalar>::value,
1059 typedef decltype((*
this)(0, 0) * factor) TS;
1062 #if LIBMESH_DIM == 1
1066 #if LIBMESH_DIM == 2
1073 #if LIBMESH_DIM == 3
1088 template <
typename T,
typename Scalar>
1090 typename boostcopy::enable_if_c<
1091 ScalarTraits<Scalar>::value,
1099 template <
typename T>
1100 template <
typename Scalar>
1102 typename boostcopy::enable_if_c<
1103 ScalarTraits<Scalar>::value,
1104 TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1107 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1109 typedef typename CompareTypes<T, Scalar>::supertype TS;
1111 #if LIBMESH_DIM == 1
1115 #if LIBMESH_DIM == 2
1122 #if LIBMESH_DIM == 3
1138 template <
typename T>
1142 #if LIBMESH_DIM == 1
1146 #if LIBMESH_DIM == 2
1153 #if LIBMESH_DIM == 3
1168 template <
typename T>
1172 #if LIBMESH_DIM == 1
1173 if (
_coords[0] == static_cast<T>(0.))
1174 libmesh_convergence_failure();
1178 #if LIBMESH_DIM == 2
1183 T a = A(0,0), b = A(0,1),
1184 c = A(1,0), d = A(1,1);
1187 T my_det = a*d - b*c;
1189 if (my_det == static_cast<T>(0.))
1190 libmesh_convergence_failure();
1192 return TypeTensor(d/my_det, -b/my_det, -c/my_det, a/my_det);
1195 #if LIBMESH_DIM == 3
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);
1203 T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1205 if (my_det == static_cast<T>(0.))
1206 libmesh_convergence_failure();
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);
1217 template <
typename T>
1221 #if LIBMESH_DIM == 1
1222 if (
_coords[0] == static_cast<T>(0.))
1223 libmesh_convergence_failure();
1227 #if LIBMESH_DIM == 2
1230 if (my_det == static_cast<T>(0.))
1231 libmesh_convergence_failure();
1233 T my_det_inv = 1./my_det;
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));
1239 #if LIBMESH_DIM == 3
1242 _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1244 -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1246 +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1248 if (my_det == static_cast<T>(0.))
1249 libmesh_convergence_failure();
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));
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));
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));
1268 template <
typename T>
1272 libmesh_assert_not_equal_to (factor, static_cast<T>(0.));
1274 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1283 template <
typename T>
1284 template <
typename T2>
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);
1297 template <
typename T>
1298 template <
typename T2>
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);
1311 template <
typename T,
typename T2>
1319 template <
typename T>
1320 template <
typename T2>
1322 TypeTensor<typename CompareTypes<T, T2>::supertype>
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);
1334 template <
typename T>
1335 template <
typename T2>
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);
1354 template <
typename T>
1355 template <
typename T2>
1357 typename CompareTypes<T,T2>::supertype
1360 typename CompareTypes<T,T2>::supertype sum = 0.;
1361 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1368 template <
typename T>
1376 template <
typename T>
1380 for (
const auto & val :
_coords)
1386 template <
typename T>
1390 #if LIBMESH_DIM == 1
1394 #if LIBMESH_DIM == 2
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]);
1409 template <
typename T>
1413 #if LIBMESH_DIM == 1
1417 #if LIBMESH_DIM == 2
1421 #if LIBMESH_DIM == 3
1422 return _coords[0] + _coords[4] + _coords[8];
1426 template <
typename T>
1430 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1436 template <
typename T>
1441 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1448 template <
typename T>
1452 #if LIBMESH_DIM == 1
1457 #if LIBMESH_DIM == 2
1465 #if LIBMESH_DIM == 3
1480 template <
typename T,
typename T2>
1486 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1487 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1496 #ifdef LIBMESH_HAVE_METAPHYSICL
1497 namespace MetaPhysicL
1499 template <
typename T>
1500 struct RawType<libMesh::TypeTensor<T>>
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));
1515 template <
typename T,
typename U>
1516 struct ReplaceAlgebraicType<libMesh::TypeTensor<T>, U>
1523 #endif // LIBMESH_TYPE_TENSOR_H
TypeTensorColumn< T > & operator=(const TypeVector< T > &rhs)
将该列的值设置为给定向量的值。
T _coords[LIBMESH_DIM]
TypeVector 的坐标。
bool operator==(const TypeTensor< T > &rhs) const
检查两个 Tensor 是否相等。
const TypeTensor< T > & operator/=(const T &)
将该 Tensor 的每个元素除以标量值。
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
TypeTensor< T > * _tensor
T & operator()(const unsigned int i)
返回该张量的第 行第 列元素的可写引用。
void add_scaled(const TypeTensor< T2 > &, const T &scale)
在不创建临时 Tensor 的情况下将一个缩放的 Tensor 加到该 Tensor 上。
T value_type
Helper typedef 用于 C++98 泛型编程。
static constexpr Real TOLERANCE
std::tuple< unsigned int, unsigned int > index_type
Helper typedef 用于泛型索引编程。
TypeTensor< typename CompareTypes< T, T2 >::supertype > operator+(const TypeTensor< T2 > &) const
将另一个 Tensor 加到该 Tensor 上。
void zero()
将 Tensor 的所有元素设置为零。
void write_unformatted(std::ostream &out_stream, const bool newline=true) const
无格式输出到流,默认为打印元素并以空格和换行符分隔。
TypeTensor()
Empty constructor.
auto norm_sq() const -> decltype(std::norm(T()))
返回 Tensor 的 Frobenius 范数的平方,即元素模平方的和。
const T & slice(const unsigned int i) const
返回该张量的第 行第 列元素的只读引用。
ConstTypeTensorColumn(const TypeTensor< T > &tensor, unsigned int j)
构造函数,初始化 ConstTypeTensorColumn。
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
计算自动微分实数向量的平方根。
T det() const
返回 Tensor 的行列式。 由于这些是最多 3x3 的 Tensor,我们不像 DenseMatrix 那样做 LU 分解。
T _coords[LIBMESH_DIM *LIBMESH_DIM]
TypeTensor 的坐标
const TypeTensor< T > * _tensor
void print(std::ostream &os=libMesh::out) const
格式化输出到流,默认为 libMesh::out。
ConstTypeTensorColumn< T > slice(const unsigned int i) const
返回 Tensor 的第 i 列的代理。
auto norm() const -> decltype(std::norm(T()))
返回 Tensor 的 Frobenius 范数,即元素平方和的平方根。
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
计算自动微分实数向量的绝对值。
TypeTensor< T > transpose() const
返回该 Tensor 的转置(对于复数不共轭)。
表示 ConstTypeTensorColumn 类,用于访问 TypeTensor 的列并进行只读操作。
void libmesh_ignore(const Args &...)
TypeTensorColumn(TypeTensor< T > &tensor, unsigned int j)
构造函数,初始化 TypeTensorColumn。
void subtract_scaled(const TypeTensor< T2 > &, const T &)
在不创建临时 Tensor 的情况下从该 Tensor 减去一个缩放的 Tensor。
auto operator*(const Scalar &scalar) const -> typename boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< decltype(T()*scalar)>>::type
将该 Tensor 与标量值相乘。
TypeTensor< T > operator-() const
返回该 Tensor 的负值的副本。
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< typename CompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
将该 Tensor 的每个元素除以标量值。
void subtract(const TypeTensor< T2 > &)
在不创建临时 Tensor 的情况下从该 Tensor 减去另一个 Tensor。
该类定义了一个在 LIBMESH_DIM 维度空间中类型为 T 的向量。
void add(const TypeTensor< T2 > &)
在不创建临时 Tensor 的情况下将另一个 Tensor 加到该 Tensor 上。
TypeVector< T > row(const unsigned int r) const
返回 Tensor 的一行作为 TypeVector 的副本。
TypeTensor< T > inverse() const
返回该 Tensor 的逆作为独立对象。
void solve(const TypeVector< T > &b, TypeVector< T > &x) const
解方程组 得到 ,其中 为该 Tensor。
bool is_zero() const
如果 Tensor 中的所有值都为零,则返回 true。
const T & operator()(const unsigned int i) const
返回该张量的第 行第 列元素的只读引用。
const T & operator()(const unsigned int i, const unsigned int j) const
返回 (i,j) 元素的常量引用。
该类最终将定义一个在类型为T的LIBMESH_DIM维空间中的N阶张量。
TypeVector< typename CompareTypes< T, T2 >::supertype > left_multiply(const TypeVector< T2 > &p) const
将该 Tensor 左乘以一个向量,即矩阵-向量乘积。 这个 Tensor 和向量可以包含不同的数值类型。
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
将两个 Tensor 相乘,返回一个标量,即 。 这些 Tensor 可能包含不同的数值类型。
表示 TypeTensorColumn 类,用于访问 TypeTensor 的列并进行写操作。
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
在不创建临时 Tensor 的情况下对该 Tensor 进行赋值。
const TypeTensor< T > & operator-=(const TypeTensor< T2 > &)
从该 Tensor 减去另一个 Tensor。
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
在该 Tensor 上加上另一个 Tensor。
T & slice(const unsigned int i)
返回该张量的第 行第 列元素的可写引用。
const TypeTensor< T > & operator*=(const Scalar &factor)
在地方上将该 Tensor 乘以标量值。
T tr() const
返回 Tensor 的迹。