15#ifndef dealii_tensor_h
16#define dealii_tensor_h
27#include <Kokkos_Array.hpp>
29#ifdef DEAL_II_WITH_ADOLC
30# include <adolc/adouble.h>
41template <
typename ElementType,
typename MemorySpace>
44template <
int dim,
typename Number>
48template <
int rank_,
int dim, typename Number =
double>
50template <typename Number>
52template <typename number>
93template <
int dim,
typename Number>
97 static_assert(dim >= 0,
98 "Tensors must have a dimension greater than or equal to one.");
113 static constexpr unsigned int rank = 0;
158 template <
typename OtherNumber>
167 template <
typename OtherNumber>
171#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
214 template <
typename OtherNumber>
218#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
231#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
245 template <
typename OtherNumber>
254 template <
typename OtherNumber>
261 template <
typename OtherNumber>
268 template <
typename OtherNumber>
277 template <
typename OtherNumber>
286 template <
typename OtherNumber>
295 template <
typename OtherNumber>
304 template <
typename OtherNumber>
355 template <
class Iterator>
357 unroll(
const Iterator begin,
const Iterator end)
const;
364 template <
class Archive>
381 template <
int,
int,
typename>
460template <
int rank_,
int dim,
typename Number>
464 static_assert(rank_ >= 1,
465 "Tensors must have a rank greater than or equal to one.");
466 static_assert(dim >= 0,
467 "Tensors must have a dimension greater than or equal to zero.");
481 static constexpr unsigned int rank = rank_;
506 std::conditional_t<rank_ == 1, Number,
Tensor<rank_ - 1, dim, Number>>;
521 Number[(dim != 0) ? dim : 1],
551 template <
typename ElementType,
typename MemorySpace>
562 template <
typename OtherNumber>
569 template <
typename OtherNumber>
576 template <
typename OtherNumber>
580#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
611 constexpr const Number &
655 template <
typename OtherNumber>
676#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
693 template <
typename OtherNumber>
700 template <
typename OtherNumber>
709 template <
typename OtherNumber>
718 template <
typename OtherNumber>
728 template <
typename OtherNumber>
737 template <
typename OtherNumber>
795 template <
class Iterator>
797 unroll(
const Iterator begin,
const Iterator end)
const;
818 static constexpr std::size_t
826 template <
class Archive>
842#if KOKKOS_VERSION >= 30700
843 std::conditional_t<rank_ == 1,
844 Kokkos::Array<Number, dim>,
845 Kokkos::Array<
Tensor<rank_ - 1, dim, Number>, dim>>
847 std::conditional_t<rank_ == 1,
848 std::array<Number, dim>,
849 std::array<
Tensor<rank_ - 1, dim, Number>, dim>>
859 template <
typename ArrayLike, std::size_t... Indices>
861 Tensor(
const ArrayLike &initializer, std::index_sequence<Indices...>);
864 template <
int,
int,
typename>
869 friend class Point<dim, Number>;
880 template <
int rank,
int dim,
typename T,
typename U>
881 struct ProductTypeImpl<
Tensor<rank, dim, T>,
std::complex<U>>
887 template <
int rank,
int dim,
typename T,
typename U>
888 struct ProductTypeImpl<
Tensor<rank, dim,
std::complex<T>>, std::complex<U>>
894 template <
typename T,
int rank,
int dim,
typename U>
895 struct ProductTypeImpl<
std::complex<T>,
Tensor<rank, dim, U>>
901 template <
int rank,
int dim,
typename T,
typename U>
902 struct ProductTypeImpl<
std::complex<T>,
Tensor<rank, dim, std::complex<U>>>
913 template <
int rank,
int dim,
typename T>
914 struct NumberType<
Tensor<rank, dim, T>>
937template <
int dim,
typename Number>
947template <
int dim,
typename Number>
948template <
typename OtherNumber>
956template <
int dim,
typename Number>
957template <
typename OtherNumber>
964# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
965template <
int dim,
typename Number>
973template <
int dim,
typename Number>
976 :
value{std::move(other.value)}
982template <
int dim,
typename Number>
987 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
992template <
int dim,
typename Number>
997 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1003template <
int dim,
typename Number>
1004template <
typename OtherNumber>
1013# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1014template <
int dim,
typename Number>
1023# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1024template <
int dim,
typename Number>
1028 value = std::move(other.value);
1035template <
int dim,
typename Number>
1036template <
typename OtherNumber>
1045template <
int dim,
typename Number>
1046template <
typename OtherNumber>
1047constexpr inline bool
1050# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1051 Assert(!(std::is_same_v<Number, adouble> ||
1052 std::is_same_v<OtherNumber, adouble>),
1054 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1055 "been extended to support advanced branching."));
1062template <
int dim,
typename Number>
1063template <
typename OtherNumber>
1067 return !((*this) == p);
1071template <
int dim,
typename Number>
1072template <
typename OtherNumber>
1081template <
int dim,
typename Number>
1082template <
typename OtherNumber>
1094 namespace ComplexWorkaround
1096 template <
typename Number,
typename OtherNumber>
1098 multiply_assign_scalar(Number &val,
const OtherNumber &s)
1103 template <
typename Number,
typename OtherNumber>
1105 multiply_assign_scalar(std::complex<Number> &val,
const OtherNumber &s)
1107# if KOKKOS_VERSION >= 30600
1108 KOKKOS_IF_ON_HOST((val *= s;))
1109 KOKKOS_IF_ON_DEVICE(({
1113 "This function is not implemented for std::complex<Number>!\n");
1116# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1122 "This function is not implemented for std::complex<Number>!\n");
1130template <
int dim,
typename Number>
1131template <
typename OtherNumber>
1135 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1141template <
int dim,
typename Number>
1142template <
typename OtherNumber>
1151template <
int dim,
typename Number>
1159template <
int dim,
typename Number>
1164 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1169template <
int dim,
typename Number>
1175 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1181template <
int dim,
typename Number>
1182constexpr inline void
1192template <
int dim,
typename Number>
1193template <
class Iterator>
1200 ExcMessage(
"Cannot unroll an object of type Tensor<0,0,Number>"));
1201 Assert(std::distance(begin, end) >= 1,
1202 ExcMessage(
"The provided iterator range must contain at least one "
1209template <
int dim,
typename Number>
1210template <
class Archive>
1218template <
int dim,
typename Number>
1224template <
int rank_,
int dim,
typename Number>
1225template <
typename ArrayLike, std::size_t... indices>
1228 std::index_sequence<indices...>)
1245 if constexpr (rank_ == 1)
1251 static_assert(
sizeof...(indices) == dim,
1252 "dim should match the number of indices");
1256# ifdef DEAL_II_HAVE_CXX20
1258template <
int rank_,
int dim,
typename Number>
1293 []<
std::size_t... I>(
1294 const
std::index_sequence<I...> &) constexpr -> decltype(
values) {
1295 if constexpr (dim == 0)
1299 else if constexpr (rank_ == 1)
1301 auto get_zero_and_ignore_argument = [](
int) {
1304 return {{(get_zero_and_ignore_argument(I))...}};
1308 auto get_zero_and_ignore_argument = [](
int) {
1309 return Tensor<rank_ - 1, dim, Number>();
1311 return {{(get_zero_and_ignore_argument(I))...}};
1313 }(std::make_index_sequence<dim>()))
1329 namespace TensorInitialization
1331 template <
int rank,
int dim,
typename Number, std::size_t... I>
1332# if KOKKOS_VERSION >= 30700
1333 constexpr Kokkos::Array<typename Tensor<rank, dim, Number>::value_type, dim>
1335 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1337 make_zero_array(
const std::index_sequence<I...> &)
1339 static_assert(
sizeof...(I) == dim,
"This is bad.");
1345 if constexpr (dim == 0)
1349 else if constexpr (rank == 1)
1351 auto get_zero_and_ignore_argument = [](
int) {
1354 return {{(get_zero_and_ignore_argument(I))...}};
1358 auto get_zero_and_ignore_argument = [](
int) {
1359 return Tensor<rank - 1, dim, Number>();
1361 return {{(get_zero_and_ignore_argument(I))...}};
1368template <
int rank_,
int dim,
typename Number>
1371 :
values(
internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1372 std::make_index_sequence<dim>()))
1379template <
int rank_,
int dim,
typename Number>
1382 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1387template <
int rank_,
int dim,
typename Number>
1388template <
typename ElementType,
typename MemorySpace>
1394 const int my_n_independent_components = n_independent_components;
1397 for (
unsigned int i = 0; i < my_n_independent_components; ++i)
1398 (*
this)[unrolled_to_component_indices(i)] = initializer[i];
1403template <
int rank_,
int dim,
typename Number>
1404template <
typename OtherNumber>
1408 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1413template <
int rank_,
int dim,
typename Number>
1414template <
typename OtherNumber>
1418 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1423template <
int rank_,
int dim,
typename Number>
1424template <
typename OtherNumber>
1426operator
Tensor<1, dim,
Tensor<rank_ - 1, dim, OtherNumber>>()
const
1429 std::copy(
values.begin(),
values.end(), x.values.begin());
1434# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1435template <
int rank_,
int dim,
typename Number>
1443template <
int rank_,
int dim,
typename Number>
1446 :
values(std::move(other.values))
1452template <
int rank_,
int dim,
typename Number>
1458 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1466template <
int rank_,
int dim,
typename Number>
1472 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1480template <
int rank_,
int dim,
typename Number>
1485 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1487 return TensorAccessors::extract<rank_>(*
this, indices);
1492template <
int rank_,
int dim,
typename Number>
1497 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1499 return TensorAccessors::extract<rank_>(*
this, indices);
1504template <
int rank_,
int dim,
typename Number>
1508 static_assert(rank_ == 1,
1509 "This function is only available for rank-1 tensors "
1510 "because higher-rank tensors may not store their elements "
1511 "in a contiguous array.");
1513 return std::addressof(
1514 this->
operator[](this->unrolled_to_component_indices(0)));
1519template <
int rank_,
int dim,
typename Number>
1520inline const Number *
1523 static_assert(rank_ == 1,
1524 "This function is only available for rank-1 tensors "
1525 "because higher-rank tensors may not store their elements "
1526 "in a contiguous array.");
1528 return std::addressof(
1529 this->
operator[](this->unrolled_to_component_indices(0)));
1534template <
int rank_,
int dim,
typename Number>
1538 static_assert(rank_ == 1,
1539 "This function is only available for rank-1 tensors "
1540 "because higher-rank tensors may not store their elements "
1541 "in a contiguous array.");
1543 return begin_raw() + n_independent_components;
1548template <
int rank_,
int dim,
typename Number>
1549inline const Number *
1552 static_assert(rank_ == 1,
1553 "This function is only available for rank-1 tensors "
1554 "because higher-rank tensors may not store their elements "
1555 "in a contiguous array.");
1557 return begin_raw() + n_independent_components;
1562template <
int rank_,
int dim,
typename Number>
1563template <
typename OtherNumber>
1569 for (
unsigned int i = 0; i < dim; ++i)
1576template <
int rank_,
int dim,
typename Number>
1584 for (
unsigned int i = 0; i < dim; ++i)
1590# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1591template <
int rank_,
int dim,
typename Number>
1595 for (
unsigned int i = 0; i < dim; ++i)
1596 values[i] = other.
values[i];
1602template <
int rank_,
int dim,
typename Number>
1607 for (
unsigned int i = 0; i < dim; ++i)
1608 values[i] = other.
values[i];
1614template <
int rank_,
int dim,
typename Number>
1615template <
typename OtherNumber>
1616constexpr inline bool
1620# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1621 Assert(!(std::is_same_v<Number, adouble> ||
1622 std::is_same_v<OtherNumber, adouble>),
1624 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1625 "been extended to support advanced branching."));
1628 for (
unsigned int i = 0; i < dim; ++i)
1642constexpr inline bool
1649template <
int rank_,
int dim,
typename Number>
1650template <
typename OtherNumber>
1655 return !((*this) == p);
1659template <
int rank_,
int dim,
typename Number>
1660template <
typename OtherNumber>
1666 for (
unsigned int i = 0; i < dim; ++i)
1667 values[i] += p.
values[i];
1672template <
int rank_,
int dim,
typename Number>
1673template <
typename OtherNumber>
1679 for (
unsigned int i = 0; i < dim; ++i)
1680 values[i] -= p.
values[i];
1685template <
int rank_,
int dim,
typename Number>
1686template <
typename OtherNumber>
1691 for (
unsigned int i = 0; i < dim; ++i)
1698template <
int rank_,
int dim,
typename Number>
1699template <
typename OtherNumber>
1704 if constexpr (std::is_integral<
1706 std::is_same_v<Number, Differentiation::SD::Expression>)
1709 for (
unsigned int d = 0;
d < dim; ++
d)
1716 const Number inverse_factor = Number(1.) / s;
1717 for (
unsigned int d = 0;
d < dim; ++
d)
1718 values[d] *= inverse_factor;
1725template <
int rank_,
int dim,
typename Number>
1732 for (
unsigned int i = 0; i < dim; ++i)
1733 tmp.
values[i] = -values[i];
1739template <
int rank_,
int dim,
typename Number>
1745 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1749 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1762 return sqrt(norm_square());
1767template <
int rank_,
int dim,
typename Number>
1772 if constexpr (dim == 0)
1775 else if constexpr (rank_ == 1)
1781 for (
unsigned int i = 1; i < dim; ++i)
1792 for (
unsigned int i = 1; i < dim; ++i)
1793 s += values[i].norm_square();
1801template <
int rank_,
int dim,
typename Number>
1802template <
class Iterator>
1805 const Iterator end)
const
1807 if constexpr (rank_ > 1)
1810 Iterator next =
begin;
1811 for (
unsigned int i = 0; i < dim; ++i)
1813 values[i].unroll(next, end);
1823 Assert(std::distance(begin, end) >= dim,
1825 "The provided iterator range must contain at least 'dim' "
1833template <
int rank_,
int dim,
typename Number>
1834constexpr inline unsigned int
1838 unsigned int index = 0;
1839 for (
int r = 0; r < rank_; ++r)
1840 index = index * dim + indices[r];
1847template <
int rank_,
int dim,
typename Number>
1852 unsigned int dummy = n_independent_components;
1856 if constexpr (dim == 0)
1860 "A tensor with dimension 0 does not store any elements. "
1861 "There is no indexing that can address its elements."));
1868 unsigned int remainder = i;
1869 for (
int r = rank_ - 1; r >= 0; --r)
1871 indices[r] = remainder % dim;
1872 remainder = remainder / dim;
1881template <
int rank_,
int dim,
typename Number>
1882constexpr inline void
1885 for (
unsigned int i = 0; i < dim; ++i)
1890template <
int rank_,
int dim,
typename Number>
1891constexpr std::size_t
1898template <
int rank_,
int dim,
typename Number>
1899template <
class Archive>
1903 for (
int i = 0; i < dim; ++i)
1910template <
int rank_,
int dim,
typename Number>
1929template <
int rank_,
int dim,
typename Number>
1930inline std::ostream &
1933 for (
unsigned int i = 0; i < dim; ++i)
1937 for (
unsigned int j = 0; j < rank_; ++j)
1951template <
int dim,
typename Number>
1952inline std::ostream &
1955 out << static_cast<const Number &>(p);
1978template <
int dim,
typename Number,
typename Other>
1983 return object *
static_cast<const Number &
>(t);
1998template <
int dim,
typename Number,
typename Other>
2003 return static_cast<const Number &
>(t) *
object;
2018template <
int dim,
typename Number,
typename OtherNumber>
2024 return static_cast<const Number &
>(src1) *
2025 static_cast<const OtherNumber &
>(src2);
2036template <
int dim,
typename Number,
typename OtherNumber>
2044 return static_cast<const Number &
>(t) / factor;
2055template <
int dim,
typename Number,
typename OtherNumber>
2061 return static_cast<const Number &
>(p) +
static_cast<const OtherNumber &
>(q);
2072template <
int dim,
typename Number,
typename OtherNumber>
2078 return static_cast<const Number &
>(p) -
static_cast<const OtherNumber &
>(q);
2094template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2120template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2143template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2166template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2187template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2204template <
int dim,
typename Number,
typename OtherNumber>
2233template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2241 for (
unsigned int i = 0; i < dim; ++i)
2292template <
int rank_1,
2296 typename OtherNumber,
2297 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2299 typename Tensor<rank_1 + rank_2 - 2,
2311 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2315 static_assert(dim > 0,
"Tensors cannot have dimension zero.");
2317 for (
unsigned int i = 1; i < dim; ++i)
2318 sum += src1[i] * src2[i];
2322 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2328 rank_1 + rank_2 - 2,
2331 for (
unsigned int i = 0; i < dim; ++i)
2332 result[i] += src1[i] * src2;
2341 rank_1 + rank_2 - 2,
2345 TensorAccessors::internal::
2346 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2347 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2348 TensorAccessors::contract<1, rank_1, rank_2, dim>(result,
2385template <
int index_1,
2391 typename OtherNumber>
2393 typename Tensor<rank_1 + rank_2 - 2,
2399 Assert(0 <= index_1 && index_1 < rank_1,
2401 "The specified index_1 must lie within the range [0,rank_1)"));
2402 Assert(0 <= index_2 && index_2 < rank_2,
2404 "The specified index_2 must lie within the range [0,rank_2)"));
2411 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2417 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2419 typename Tensor<rank_1 + rank_2 - 2,
2423 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2458template <
int index_1,
2466 typename OtherNumber>
2468 typename Tensor<rank_1 + rank_2 - 4,
2474 Assert(0 <= index_1 && index_1 < rank_1,
2476 "The specified index_1 must lie within the range [0,rank_1)"));
2477 Assert(0 <= index_3 && index_3 < rank_1,
2479 "The specified index_3 must lie within the range [0,rank_1)"));
2480 Assert(index_1 != index_3,
2481 ExcMessage(
"index_1 and index_3 must not be the same"));
2482 Assert(0 <= index_2 && index_2 < rank_2,
2484 "The specified index_2 must lie within the range [0,rank_2)"));
2485 Assert(0 <= index_4 && index_4 < rank_2,
2487 "The specified index_4 must lie within the range [0,rank_2)"));
2488 Assert(index_2 != index_4,
2489 ExcMessage(
"index_2 and index_4 must not be the same"));
2496 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2500 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2506 (index_3 < index_1 ? index_3 : index_3 - 1),
2518 (index_4 < index_2 ? index_4 : index_4 - 1),
2526 typename Tensor<rank_1 + rank_2 - 4,
2530 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2547template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2554 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2576template <
template <
int,
int,
typename>
class TensorT1,
2577 template <
int,
int,
typename>
2579 template <
int,
int,
typename>
2589 contract3(
const TensorT1<rank_1, dim, T1> &left,
2590 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2591 const TensorT3<rank_2, dim, T3> &right)
2595 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2611template <
int rank_1,
2615 typename OtherNumber>
2621 typename Tensor<rank_1 + rank_2,
2625 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2648template <
int dim,
typename Number>
2657 result[1] = -src[0];
2672template <
int dim,
typename Number1,
typename Number2>
2683 constexpr int s0 = 0 % dim;
2684 constexpr int s1 = 1 % dim;
2685 constexpr int s2 = 2 % dim;
2687 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2688 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2689 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2708template <
int dim,
typename Number>
2716 for (
unsigned int k = 0; k < dim; ++k)
2718 Tensor<2, dim - 1, Number> minor;
2719 for (
unsigned int i = 0; i < dim - 1; ++i)
2720 for (
unsigned int j = 0; j < dim - 1; ++j)
2721 minor[i][j] = t[i][j < k ? j : j + 1];
2723 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) *
determinant(minor);
2725 det += t[dim - 1][k] * cofactor;
2728 return ((dim % 2 == 0) ? 1. : -1.) * det;
2736template <
typename Number>
2748template <
typename Number>
2753 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2761template <
typename Number>
2772 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2782template <
int dim,
typename Number>
2787 for (
unsigned int i = 1; i < dim; ++i)
2801template <
int dim,
typename Number>
2805 Number return_tensor[dim][dim];
2818template <
typename Number>
2826 return return_tensor;
2830template <
typename Number>
2837 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2838 return_tensor[0][0] = t[1][1];
2839 return_tensor[0][1] = -t[0][1];
2840 return_tensor[1][0] = -t[1][0];
2841 return_tensor[1][1] = t[0][0];
2842 return_tensor *= inv_det_t;
2844 return return_tensor;
2847template <
typename Number>
2853 const auto value = [](
const auto &t) {
2857 return_tensor[0][0] =
value(t[1][1] * t[2][2]) -
value(t[1][2] * t[2][1]);
2858 return_tensor[0][1] =
value(t[0][2] * t[2][1]) -
value(t[0][1] * t[2][2]);
2859 return_tensor[0][2] =
value(t[0][1] * t[1][2]) -
value(t[0][2] * t[1][1]);
2860 return_tensor[1][0] =
value(t[1][2] * t[2][0]) -
value(t[1][0] * t[2][2]);
2861 return_tensor[1][1] =
value(t[0][0] * t[2][2]) -
value(t[0][2] * t[2][0]);
2862 return_tensor[1][2] =
value(t[0][2] * t[1][0]) -
value(t[0][0] * t[1][2]);
2863 return_tensor[2][0] =
value(t[1][0] * t[2][1]) -
value(t[1][1] * t[2][0]);
2864 return_tensor[2][1] =
value(t[0][1] * t[2][0]) -
value(t[0][0] * t[2][1]);
2865 return_tensor[2][2] =
value(t[0][0] * t[1][1]) -
value(t[0][1] * t[1][0]);
2867 const Number inv_det_t =
2868 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2869 t[0][2] * return_tensor[2][0]));
2870 return_tensor *= inv_det_t;
2872 return return_tensor;
2883template <
int dim,
typename Number>
2888 for (
unsigned int i = 0; i < dim; ++i)
2891 for (
unsigned int j = i + 1; j < dim; ++j)
2914template <
int dim,
typename Number>
2935template <
int dim,
typename Number>
3006template <
int dim,
typename Number>
3018template <
int dim,
typename Number>
3023 for (
unsigned int j = 0; j < dim; ++j)
3026 for (
unsigned int i = 0; i < dim; ++i)
3044template <
int dim,
typename Number>
3049 for (
unsigned int i = 0; i < dim; ++i)
3052 for (
unsigned int j = 0; j < dim; ++j)
3070# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3079 for (
unsigned int j = 0; j < dim; ++j)
3082 for (
unsigned int i = 0; i < dim; ++i)
3083 sum += fabs(t[i][j]);
3085 condassign(max, (sum > max), sum, max);
3097 for (
unsigned int i = 0; i < dim; ++i)
3100 for (
unsigned int j = 0; j < dim; ++j)
3101 sum +=
fabs(t[i][j]);
3103 condassign(max, (sum > max), sum, max);
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr Tensor & operator=(const OtherNumber &d) &&=delete
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor & operator/=(const OtherNumber &factor)
constexpr Tensor & operator-=(const Tensor< 0, dim, OtherNumber > &rhs)
constexpr Tensor(const Tensor< 0, dim, OtherNumber > &initializer)
constexpr Tensor(const OtherNumber &initializer)
constexpr real_type norm_square() const
constexpr bool operator!=(const Tensor< 0, dim, OtherNumber > &rhs) const
constexpr Tensor & operator=(const Tensor< 0, dim, OtherNumber > &rhs)
constexpr Tensor & operator+=(const Tensor< 0, dim, OtherNumber > &rhs)
void unroll(const Iterator begin, const Iterator end) const
constexpr bool operator==(const Tensor< 0, dim, OtherNumber > &rhs) const
typename numbers::NumberTraits< Number >::real_type real_type
constexpr Tensor & operator=(const OtherNumber &d) &
constexpr Tensor operator-() const
constexpr Tensor & operator/=(const OtherNumber &factor)
constexpr Tensor(const ArrayView< ElementType, MemorySpace > &initializer)
constexpr bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
constexpr Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
constexpr Tensor(const Tensor< 1, dim, Tensor< rank_ - 1, dim, OtherNumber > > &initializer)
constexpr const Number & operator[](const TableIndices< rank_ > &indices) const
constexpr Tensor< 0, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< 0, dim, Number > &t, const OtherNumber &factor)
static constexpr unsigned int rank
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
numbers::NumberTraits< Number >::real_type norm() const
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
void unroll(const Iterator begin, const Iterator end) const
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
constexpr ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
const Number * begin_raw() const
std::conditional_t< rank_==1, Number[(dim !=0) ? dim :1], typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1]> array_type
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
constexpr value_type & operator[](const unsigned int i)
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
static constexpr unsigned int dimension
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr Tensor & operator=(const Number &d) &&=delete
constexpr Tensor< rank, dim, typename ProductType< typename EnableIfScalar< Number >::type, OtherNumber >::type > operator*(const Number &factor, const Tensor< rank, dim, OtherNumber > &t)
static constexpr std::size_t memory_consumption()
constexpr Tensor & operator=(const Number &d) &
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
const Number * end_raw() const
constexpr Number & operator[](const TableIndices< rank_ > &indices)
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< rank, dim, Number > &src1, const Tensor< rank, dim, OtherNumber > &src2)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
constexpr Tensor(const ArrayLike &initializer, std::index_sequence< Indices... >)
constexpr Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr const value_type & operator[](const unsigned int i) const
constexpr Tensor(const array_type &initializer)
static constexpr unsigned int n_independent_components
constexpr Tensor operator-() const
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_CXX23_ASSUME(expr)
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
constexpr bool value_is_zero(const Number &value)
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
#define DEAL_II_HOST_DEVICE
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
static constexpr real_type abs_square(const number &x)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
std::ostream & operator<<(std::ostream &out, const Tensor< rank_, dim, Number > &p)
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Number l1_norm(const Tensor< 2, dim, Number > &t)