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"
40template <
typename T>
class TypeTensorColumn;
41template <
typename T>
class ConstTypeTensorColumn;
42template <
unsigned int N,
typename T>
class TypeNTensor;
54 template <
typename T2>
85 template <
typename Scalar>
96 const Scalar>::type & zz=0);
103 template<
typename T2>
106 template<
typename T2>
110 template<
typename T2>
130 template<
typename T2>
141 template<
typename T2>
149 template <
typename Scalar>
150 typename std::enable_if<
154 { libmesh_assert_equal_to (p, Scalar(0)); this->
zero();
return *
this; }
159 const T &
operator () (
const unsigned int i,
const unsigned int j)
const;
164 T &
operator () (
const unsigned int i,
const unsigned int j);
191 template<
typename T2>
200 template<
typename T2>
206 template<
typename T2>
212 template <
typename T2>
220 template<
typename T2>
229 template<
typename T2>
235 template<
typename T2>
242 template <
typename T2>
255 template <
typename Scalar>
257 operator * (
const Scalar & scalar)
const ->
typename std::enable_if<
266 template <
typename Scalar,
typename std::enable_if<
270 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
281 template <
typename Scalar>
282 typename std::enable_if<
300 template <
typename T2>
309 template <
typename T2>
321 template <
typename T2>
331 template <
typename T2>
341 template <
typename T2>
437 const bool newline =
true)
const;
479 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
486 const unsigned int _j;
505 const T &
slice (
const unsigned int i)
const
510 const unsigned int _j;
556#elif LIBMESH_DIM == 1
557 libmesh_assert_equal_to (xy, 0);
558 libmesh_assert_equal_to (yx, 0);
559 libmesh_assert_equal_to (yy, 0);
573 libmesh_assert_equal_to (xz, 0);
574 libmesh_assert_equal_to (yz, 0);
575 libmesh_assert_equal_to (zx, 0);
576 libmesh_assert_equal_to (zy, 0);
577 libmesh_assert_equal_to (zz, 0);
584template <
typename Scalar>
596 const Scalar>::type & zz)
604#elif LIBMESH_DIM == 1
605 libmesh_assert_equal_to (xy, 0);
606 libmesh_assert_equal_to (yx, 0);
607 libmesh_assert_equal_to (yy, 0);
621 libmesh_assert_equal_to (xz, 0);
622 libmesh_assert_equal_to (yz, 0);
623 libmesh_assert_equal_to (zx, 0);
624 libmesh_assert_equal_to (zy, 0);
625 libmesh_assert_equal_to (zz, 0);
638 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
644template <
typename T2>
647 libmesh_assert_equal_to (LIBMESH_DIM, 1);
652template <
typename T2>
656 libmesh_assert_equal_to (LIBMESH_DIM, 2);
668template <
typename T2>
673 libmesh_assert_equal_to (LIBMESH_DIM, 3);
705 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
714 const unsigned int j)
const
716 libmesh_assert_less (i, 3);
717 libmesh_assert_less (j, 3);
720 const static T my_zero = 0;
721 if (i >= LIBMESH_DIM || j >= LIBMESH_DIM)
725 return _coords[i*LIBMESH_DIM+j];
733 const unsigned int j)
737 libmesh_error_msg_if(i >= LIBMESH_DIM || j >= LIBMESH_DIM,
738 "ERROR: You are assigning to a tensor component that is out of range for the compiled LIBMESH_DIM!");
742 libmesh_assert_less (i, LIBMESH_DIM);
743 libmesh_assert_less (j, LIBMESH_DIM);
745 return _coords[i*LIBMESH_DIM+j];
754 libmesh_assert_less (i, LIBMESH_DIM);
764 libmesh_assert_less (i, LIBMESH_DIM);
776 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
777 return_vector.
_coords[j] = _coords[r*LIBMESH_DIM + j];
779 return return_vector;
790 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
791 return_vector.
_coords[i] = _coords[i*LIBMESH_DIM + r];
793 return return_vector;
849 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
856template <
typename T2>
860 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
861 _coords[i] += factor*p.
_coords[i];
903template <
typename T2>
915template <
typename T2>
919 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
926template <
typename T2>
930 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
931 _coords[i] -= factor*p.
_coords[i];
969template <
typename Scalar>
976 typedef decltype((*this)(0, 0) * factor) TS;
1005template <
typename T,
typename Scalar>
1007typename std::enable_if<
1016template <
typename T>
1017template <
typename Scalar>
1019typename std::enable_if<
1021 TypeTensor<typename CompareTypes<T, Scalar>::supertype>>::type
1024 libmesh_assert_not_equal_to (factor,
static_cast<T
>(0.));
1055template <
typename T>
1085template <
typename T>
1090 if (_coords[0] ==
static_cast<T
>(0.))
1091 libmesh_convergence_failure();
1100 T a = A(0,0),
b = A(0,1),
1101 c = A(1,0), d = A(1,1);
1104 T my_det = a*d -
b*c;
1106 if (my_det ==
static_cast<T
>(0.))
1107 libmesh_convergence_failure();
1109 return TypeTensor(d/my_det, -
b/my_det, -c/my_det, a/my_det);
1116 T a11 = A(0,0), a12 = A(0,1), a13 = A(0,2),
1117 a21 = A(1,0), a22 = A(1,1), a23 = A(1,2),
1118 a31 = A(2,0), a32 = A(2,1), a33 = A(2,2);
1120 T my_det = a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
1122 if (my_det ==
static_cast<T
>(0.))
1123 libmesh_convergence_failure();
1126 return TypeTensor( (a33*a22-a32*a23)/my_det, -(a33*a12-a32*a13)/my_det, (a23*a12-a22*a13)/my_det,
1127 -(a33*a21-a31*a23)/my_det, (a33*a11-a31*a13)/my_det, -(a23*a11-a21*a13)/my_det,
1128 (a32*a21-a31*a22)/my_det, -(a32*a11-a31*a12)/my_det, (a22*a11-a21*a12)/my_det);
1134template <
typename T>
1139 if (_coords[0] ==
static_cast<T
>(0.))
1140 libmesh_convergence_failure();
1141 x(0) =
b(0) / _coords[0];
1145 T my_det = _coords[0]*_coords[3] - _coords[1]*_coords[2];
1147 if (my_det ==
static_cast<T
>(0.))
1148 libmesh_convergence_failure();
1150 T my_det_inv = 1./my_det;
1152 x(0) = my_det_inv*( _coords[3]*
b(0) - _coords[1]*
b(1));
1153 x(1) = my_det_inv*(-_coords[2]*
b(0) + _coords[0]*
b(1));
1159 _coords[0]*(_coords[8]*_coords[4] - _coords[7]*_coords[5])
1161 -_coords[3]*(_coords[8]*_coords[1] - _coords[7]*_coords[2]) +
1163 +_coords[6]*(_coords[5]*_coords[1] - _coords[4]*_coords[2]);
1165 if (my_det ==
static_cast<T
>(0.))
1166 libmesh_convergence_failure();
1168 T my_det_inv = 1./my_det;
1169 x(0) = my_det_inv*((_coords[8]*_coords[4] - _coords[7]*_coords[5])*
b(0) -
1170 (_coords[8]*_coords[1] - _coords[7]*_coords[2])*
b(1) +
1171 (_coords[5]*_coords[1] - _coords[4]*_coords[2])*
b(2));
1173 x(1) = my_det_inv*((_coords[6]*_coords[5] - _coords[8]*_coords[3])*
b(0) +
1174 (_coords[8]*_coords[0] - _coords[6]*_coords[2])*
b(1) -
1175 (_coords[5]*_coords[0] - _coords[3]*_coords[2])*
b(2));
1177 x(2) = my_det_inv*((_coords[7]*_coords[3] - _coords[6]*_coords[4])*
b(0) -
1178 (_coords[7]*_coords[0] - _coords[6]*_coords[1])*
b(1) +
1179 (_coords[4]*_coords[0] - _coords[3]*_coords[1])*
b(2));
1185template <
typename T>
1189 libmesh_assert_not_equal_to (factor,
static_cast<T
>(0.));
1191 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1192 _coords[i] /= factor;
1200template <
typename T>
1201template <
typename T2>
1207 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1208 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1209 returnval(i) += (*this)(i,j)*p(j);
1214template <
typename T>
1215template <
typename T2>
1221 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1222 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1223 returnval(i) += p(j)*(*this)(j,i);
1228template <
typename T,
typename T2>
1233 return b.left_multiply(a);
1236template <
typename T>
1237template <
typename T2>
1239TypeTensor<typename CompareTypes<T, T2>::supertype>
1243 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1244 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1245 for (
unsigned int k=0; k<LIBMESH_DIM; k++)
1246 returnval(i,j) += (*this)(i,k)*p(k,j);
1251template <
typename T>
1252template <
typename T2>
1257 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1258 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1259 for (
unsigned int k=0; k<LIBMESH_DIM; k++)
1260 temp(i,j) += (*this)(i,k)*p(k,j);
1271template <
typename T>
1272template <
typename T2>
1278 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1279 sum += _coords[i]*t.
_coords[i];
1285template <
typename T>
1290 return sqrt(this->norm_sq());
1294template <
typename T>
1298 for (
const auto & val : _coords)
1304template <
typename T>
1313 return (_coords[0] * _coords[3]
1314 - _coords[1] * _coords[2]);
1318 return (_coords[0] * _coords[4] * _coords[8]
1319 + _coords[1] * _coords[5] * _coords[6]
1320 + _coords[2] * _coords[3] * _coords[7]
1321 - _coords[0] * _coords[5] * _coords[7]
1322 - _coords[1] * _coords[3] * _coords[8]
1323 - _coords[2] * _coords[4] * _coords[6]);
1327template <
typename T>
1336 return _coords[0] + _coords[3];
1340 return _coords[0] + _coords[4] + _coords[8];
1344template <
typename T>
1348 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1354template <
typename T>
1359 for (
unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
1366template <
typename T>
1371 return (std::abs(_coords[0] - rhs.
_coords[0])
1376 return ((std::abs(_coords[0] - rhs.
_coords[0]) +
1377 std::abs(_coords[1] - rhs.
_coords[1]) +
1378 std::abs(_coords[2] - rhs.
_coords[2]) +
1379 std::abs(_coords[3] - rhs.
_coords[3]))
1384 return ((std::abs(_coords[0] - rhs.
_coords[0]) +
1385 std::abs(_coords[1] - rhs.
_coords[1]) +
1386 std::abs(_coords[2] - rhs.
_coords[2]) +
1387 std::abs(_coords[3] - rhs.
_coords[3]) +
1388 std::abs(_coords[4] - rhs.
_coords[4]) +
1389 std::abs(_coords[5] - rhs.
_coords[5]) +
1390 std::abs(_coords[6] - rhs.
_coords[6]) +
1391 std::abs(_coords[7] - rhs.
_coords[7]) +
1392 std::abs(_coords[8] - rhs.
_coords[8]))
1398template <
typename T>
1403 os <<
"x=" << (*this)(0,0) << std::endl;
1409 << std::setw(8) << (*this)(0,0) <<
", "
1410 << std::setw(8) << (*this)(0,1) <<
")"
1413 << std::setw(8) << (*this)(1,0) <<
", "
1414 << std::setw(8) << (*this)(1,1) <<
")"
1420 os <<
"(xx,xy,xz)=("
1421 << std::setw(8) << (*this)(0,0) <<
", "
1422 << std::setw(8) << (*this)(0,1) <<
", "
1423 << std::setw(8) << (*this)(0,2) <<
")"
1425 os <<
"(yx,yy,yz)=("
1426 << std::setw(8) << (*this)(1,0) <<
", "
1427 << std::setw(8) << (*this)(1,1) <<
", "
1428 << std::setw(8) << (*this)(1,2) <<
")"
1430 os <<
"(zx,zy,zz)=("
1431 << std::setw(8) << (*this)(2,0) <<
", "
1432 << std::setw(8) << (*this)(2,1) <<
", "
1433 << std::setw(8) << (*this)(2,2) <<
")"
1440template <
typename T>
1443 using std::isfinite;
1460template <
typename T>
1467 bool has_inf =
false;
1472 if (
isnan(var(i,j)))
1474 has_inf = has_inf ||
isinf(var(i,j));
1487template <
typename T>
1494 if (
isnan(var(i,j)))
1501template <
typename T,
typename T2>
1503TypeTensor<typename CompareTypes<T, T2>::supertype>
1507 for (
unsigned int i=0; i<LIBMESH_DIM; i++)
1508 for (
unsigned int j=0; j<LIBMESH_DIM; j++)
1517#ifdef LIBMESH_HAVE_METAPHYSICL
1520template <
typename T>
1528 for (
unsigned int i = 0; i < LIBMESH_DIM; ++i)
1529 for (
unsigned int j = 0; j < LIBMESH_DIM; ++j)
1530 ret(i,j) = raw_value(in(i,j));
1536template <
typename T,
typename U>
1537struct 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.
bool operator<(const TypeTensor< T > &rhs) const
std::enable_if< ScalarTraits< Scalar >::value, TypeTensor & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
TypeTensor< T > inverse() const
const T & operator()(const unsigned int i, const unsigned int j) const
void subtract_scaled(const TypeTensor< T2 > &, const T &)
Subtract a scaled value from this tensor without creating a temporary.
const TypeTensor< T > & operator+=(const TypeTensor< T2 > &)
Add to this tensor.
void add(const TypeTensor< T2 > &)
Add to this tensor without creating a temporary.
T value_type
Helper typedef for C++98 generic programming.
auto operator*(const Scalar &scalar) const -> typename std::enable_if< ScalarTraits< Scalar >::value, TypeTensor< decltype(T() *scalar)> >::type
Multiply this tensor by a scalar value.
void assign(const TypeTensor< T2 > &)
Assign to this tensor without creating a temporary.
const TypeTensor< T > & operator*=(const Scalar &factor)
Multiply this tensor by a scalar value in place.
std::enable_if< ScalarTraits< Scalar >::value, TypeTensor< typenameCompareTypes< T, Scalar >::supertype > >::type operator/(const Scalar &) const
Divide each entry of this tensor by a scalar value.
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
const TypeTensor< T > & operator/=(const T &)
Divide each entry of this tensor by a scalar value.
std::tuple< unsigned int, unsigned int > index_type
Helper typedef for generic index programming.
void subtract(const TypeTensor< T2 > &)
Subtract from this tensor without creating a temporary.
TypeVector< T > column(const unsigned int r) const
void print(std::ostream &os=libMesh::out) const
Formatted print to a stream which defaults to libMesh::out.
void write_unformatted(std::ostream &out_stream, const bool newline=true) const
Unformatted print to the stream out.
TypeVector< T > row(const unsigned int r) const
T _coords[LIBMESH_DIM *LIBMESH_DIM]
The coordinates of the TypeTensor.
TypeTensor< typename CompareTypes< T, T2 >::supertype > operator+(const TypeTensor< T2 > &) const
Add another tensor to this tensor.
friend std::ostream & operator<<(std::ostream &os, const TypeTensor< T > &t)
Formatted print as above but supports the syntax:
bool is_hpd(Real rel_tol=TOLERANCE *TOLERANCE) const
Returns true if the TypeTensor is Hermitian (or symmetric, when T==Real) and positive-definite to wit...
bool operator==(const TypeTensor< T > &rhs) const
void zero()
Set all entries of the tensor to 0.
TypeVector< typename CompareTypes< T, T2 >::supertype > left_multiply(const TypeVector< T2 > &p) const
Left-multiply this tensor by a vector, i.e.
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
TypeTensor< T > transpose() const
void solve(const TypeVector< T > &b, TypeVector< T > &x) const
Solve the 2x2 or 3x3 system of equations A*x = b for x by directly inverting A and multiplying it by ...
TypeTensor< T > operator-() const
ConstTypeTensorColumn< T > slice(const unsigned int i) const
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.
void libmesh_ignore(const Args &...)
bool isfinite(std::complex< T > a)
bool isnan(std::complex< T > a)
std::enable_if< ScalarTraits< Scalar >::value, TypeNTensor< N, typenameCompareTypes< Scalar, T >::supertype > >::type operator*(const Scalar &, const TypeNTensor< N, T > &)
TypeTensor< typename CompareTypes< T, T2 >::supertype > outer_product(const TypeVector< T > &a, const TypeVector< T2 > &b)
static constexpr Real TOLERANCE
bool isinf(std::complex< T > a)
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...