20#ifndef LIBMESH_TYPE_TENSOR_H
21#define LIBMESH_TYPE_TENSOR_H
24#include "libmesh/libmesh_common.h"
25#include "libmesh/libmesh_device.h"
26#include "libmesh/type_vector.h"
33#ifdef LIBMESH_HAVE_METAPHYSICL
34#include "metaphysicl/dualnumber_decl.h"
41template <
typename T>
class TypeTensorColumn;
42template <
typename T>
class ConstTypeTensorColumn;
43template <
unsigned int N,
typename T>
class TypeNTensor;
55 template <
typename T2>
86 template <
typename Scalar>
97 const Scalar>::type & zz=0);
104 template<
typename T2>
105 LIBMESH_DEVICE_INLINE
108 template<
typename T2>
109 LIBMESH_DEVICE_INLINE
113 template<
typename T2>
114 LIBMESH_DEVICE_INLINE
134 template<
typename T2>
140 LIBMESH_DEVICE_INLINE
146 template<
typename T2>
147 LIBMESH_DEVICE_INLINE
155 template <
typename Scalar>
156 LIBMESH_DEVICE_INLINE
157 typename std::enable_if<
161 { libmesh_assert_equal_to (p, Scalar(0)); this->
zero();
return *
this; }
166 const T &
operator () (
const unsigned int i,
const unsigned int j)
const;
176 LIBMESH_DEVICE_INLINE
182 LIBMESH_DEVICE_INLINE
193 LIBMESH_DEVICE_INLINE
201 template<
typename T2>
210 template<
typename T2>
216 template<
typename T2>
222 template <
typename T2>
223 LIBMESH_DEVICE_INLINE
231 template<
typename T2>
240 template<
typename T2>
246 template<
typename T2>
253 template <
typename T2>
254 LIBMESH_DEVICE_INLINE
267 template <
typename Scalar>
269 operator * (
const Scalar & scalar)
const ->
typename std::enable_if<
278 template <
typename Scalar,
typename std::enable_if<
280 LIBMESH_DEVICE_INLINE
283 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
294 template <
typename Scalar>
295 typename std::enable_if<
313 template <
typename T2>
322 template <
typename T2>
334 template <
typename T2>
335 LIBMESH_DEVICE_INLINE
345 template <
typename T2>
355 template <
typename T2>
356 LIBMESH_DEVICE_INLINE
376 LIBMESH_DEVICE_INLINE
394 LIBMESH_DEVICE_INLINE
413 LIBMESH_DEVICE_INLINE
419 LIBMESH_DEVICE_INLINE
456 const bool newline =
true)
const;
498 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
505 const unsigned int _j;
524 const T &
slice (
const unsigned int i)
const
529 const unsigned int _j;
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);
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);
603template <
typename Scalar>
615 const Scalar>::type & zz)
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);
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);
657 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
663template <
typename T2>
667 libmesh_assert_equal_to (LIBMESH_DIM, 1);
672template <
typename T2>
677 libmesh_assert_equal_to (LIBMESH_DIM, 2);
689template <
typename T2>
695 libmesh_assert_equal_to (LIBMESH_DIM, 3);
727 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
736 const unsigned int j)
const
738 libmesh_assert_less (i, 3);
739 libmesh_assert_less (j, 3);
742 const static T my_zero = 0;
743 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
747 return _coords[i*LIBMESH_DIM+j];
755 const unsigned int j)
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!");
764 libmesh_assert_less (i, LIBMESH_DIM);
765 libmesh_assert_less (j, LIBMESH_DIM);
767 return _coords[i*LIBMESH_DIM+j];
776 libmesh_assert_less (i, LIBMESH_DIM);
786 libmesh_assert_less (i, LIBMESH_DIM);
798 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
799 return_vector.
_coords[j] = _coords[r*LIBMESH_DIM + j];
801 return return_vector;
812 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
813 return_vector.
_coords[i] = _coords[i*LIBMESH_DIM + r];
815 return return_vector;
871 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
878template <
typename T2>
882 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
883 _coords[i] += factor*p.
_coords[i];
925template <
typename T2>
937template <
typename T2>
941 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
948template <
typename T2>
952 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
953 _coords[i] -= factor*p.
_coords[i];
991template <
typename Scalar>
998 typedef decltype((*this)(0, 0) * factor) TS;
1027template <
typename T,
typename Scalar>
1028LIBMESH_DEVICE_INLINE
1029typename std::enable_if<
1038template <
typename T>
1039template <
typename Scalar>
1040LIBMESH_DEVICE_INLINE
1041typename std::enable_if<
1043 TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1046 libmesh_assert_not_equal_to (factor,
static_cast<T
>(0.));
1077template <
typename T>
1078LIBMESH_DEVICE_INLINE
1107template <
typename T>
1108LIBMESH_DEVICE_INLINE
1112 if (_coords[0] ==
static_cast<T
>(0.))
1113 libmesh_convergence_failure();
1122 T a = A(0,0),
b = A(0,1),
1123 c = A(1,0), d = A(1,1);
1126 T my_det = a*d -
b*c;
1128 if (my_det ==
static_cast<T
>(0.))
1129 libmesh_convergence_failure();
1131 return TypeTensor(d/my_det, -
b/my_det, -c/my_det, a/my_det);
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);
1142 T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1144 if (my_det ==
static_cast<T
>(0.))
1145 libmesh_convergence_failure();
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);
1156template <
typename T>
1157LIBMESH_DEVICE_INLINE
1161 if (_coords[0] ==
static_cast<T
>(0.))
1162 libmesh_convergence_failure();
1163 x(0) =
b(0) / _coords[0];
1167 T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1169 if (my_det ==
static_cast<T
>(0.))
1170 libmesh_convergence_failure();
1172 T my_det_inv = 1./my_det;
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));
1181 _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1183 -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1185 +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1187 if (my_det ==
static_cast<T
>(0.))
1188 libmesh_convergence_failure();
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));
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));
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));
1207template <
typename T>
1208LIBMESH_DEVICE_INLINE
1211 libmesh_assert_not_equal_to (factor,
static_cast<T
>(0.));
1213 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1214 _coords[i] /= factor;
1222template <
typename T>
1223template <
typename T2>
1224LIBMESH_DEVICE_INLINE
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);
1236template <
typename T>
1237template <
typename T2>
1238LIBMESH_DEVICE_INLINE
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);
1250template <
typename T,
typename T2>
1251LIBMESH_DEVICE_INLINE
1255 return b.left_multiply(a);
1258template <
typename T>
1259template <
typename T2>
1260LIBMESH_DEVICE_INLINE
1261TypeTensor<typename CompareTypes<T, T2>::supertype>
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);
1273template <
typename T>
1274template <
typename T2>
1275LIBMESH_DEVICE_INLINE
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);
1293template <
typename T>
1294template <
typename T2>
1295LIBMESH_DEVICE_INLINE
1300 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1301 sum += _coords[i]*t.
_coords[i];
1307template <
typename T>
1308LIBMESH_DEVICE_INLINE
1312 return sqrt(this->norm_sq());
1316template <
typename T>
1317LIBMESH_DEVICE_INLINE
1320 for (
const auto & val : _coords)
1326template <
typename T>
1327LIBMESH_DEVICE_INLINE
1335 return (_coords[0] * _coords[3]
1336 - _coords[1] * _coords[2]);
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]);
1349template <
typename T>
1350LIBMESH_DEVICE_INLINE
1358 return _coords[0] + _coords[3];
1362 return _coords[0] + _coords[4] + _coords[8];
1366template <
typename T>
1367LIBMESH_DEVICE_INLINE
1370 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1376template <
typename T>
1377LIBMESH_DEVICE_INLINE
1381 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1388template <
typename T>
1389LIBMESH_DEVICE_INLINE
1393 return (std::abs(_coords[0] - rhs.
_coords[0])
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]))
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]))
1420template <
typename T>
1425 os <<
"x=" << (*this)(0,0) << std::endl;
1431 << std::setw(8) << (*this)(0,0) <<
", "
1432 << std::setw(8) << (*this)(0,1) <<
")"
1435 << std::setw(8) << (*this)(1,0) <<
", "
1436 << std::setw(8) << (*this)(1,1) <<
")"
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) <<
")"
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) <<
")"
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) <<
")"
1462template <
typename T>
1465 using std::isfinite;
1482template <
typename T>
1489 bool has_inf =
false;
1494 if (
isnan(var(i,j)))
1496 has_inf = has_inf ||
isinf(var(i,j));
1509template <
typename T>
1516 if (
isnan(var(i,j)))
1523template <
typename T,
typename T2>
1524LIBMESH_DEVICE_INLINE
1525TypeTensor<typename CompareTypes<T, T2>::supertype>
1529 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1530 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1539#ifdef LIBMESH_HAVE_METAPHYSICL
1542template <
typename T>
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));
1558template <
typename T,
typename U>
1559struct ReplaceAlgebraicType<
libMesh::TypeTensor<T>, U>
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)
TypeTensor< T > * _tensor
T & operator()(const unsigned int i)
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
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
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.
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.
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...
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.
T _coords[LIBMESH_DIM]
The coordinates of the TypeVector.
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 &...)
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...
LIBMESH_DEVICE_INLINE T libmesh_conj(T a)