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>
42template <
typename ElementType,
typename MemorySpace>
45template <
int dim,
typename Number>
49template <
int rank_,
int dim, typename Number =
double>
51template <typename Number>
53template <typename number>
94template <
int dim,
typename Number>
98 static_assert(dim >= 0,
99 "Tensors must have a dimension greater than or equal to one.");
114 static constexpr unsigned int rank = 0;
159 template <
typename OtherNumber>
168 template <
typename OtherNumber>
172#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
215 template <
typename OtherNumber>
219#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
232#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
246 template <
typename OtherNumber>
255 template <
typename OtherNumber>
262 template <
typename OtherNumber>
269 template <
typename OtherNumber>
278 template <
typename OtherNumber>
287 template <
typename OtherNumber>
296 template <
typename OtherNumber>
305 template <
typename OtherNumber>
356 template <
class Iterator>
358 unroll(
const Iterator begin,
const Iterator end)
const;
365 template <
class Archive>
382 template <
int,
int,
typename>
461template <
int rank_,
int dim,
typename Number>
465 static_assert(rank_ >= 1,
466 "Tensors must have a rank greater than or equal to one.");
467 static_assert(dim >= 0,
468 "Tensors must have a dimension greater than or equal to zero.");
482 static constexpr unsigned int rank = rank_;
507 std::conditional_t<rank_ == 1, Number,
Tensor<rank_ - 1, dim, Number>>;
522 Number[(dim != 0) ? dim : 1],
552 template <
typename ElementType,
typename MemorySpace>
563 template <
typename OtherNumber>
570 template <
typename OtherNumber>
577 template <
typename OtherNumber>
581#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
612 constexpr const Number &
656 template <
typename OtherNumber>
677#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
694 template <
typename OtherNumber>
701 template <
typename OtherNumber>
710 template <
typename OtherNumber>
719 template <
typename OtherNumber>
729 template <
typename OtherNumber>
738 template <
typename OtherNumber>
796 template <
class Iterator>
798 unroll(
const Iterator begin,
const Iterator end)
const;
819 static constexpr std::size_t
827 template <
class Archive>
843#if DEAL_II_KOKKOS_VERSION_GTE(3, 7, 0)
844 std::conditional_t<rank_ == 1,
845 Kokkos::Array<Number, dim>,
846 Kokkos::Array<
Tensor<rank_ - 1, dim, Number>, dim>>
848 std::conditional_t<rank_ == 1,
849 std::array<Number, dim>,
850 std::array<
Tensor<rank_ - 1, dim, Number>, dim>>
860 template <
typename ArrayLike, std::size_t... Indices>
862 Tensor(
const ArrayLike &initializer, std::index_sequence<Indices...>);
865 template <
int,
int,
typename>
870 friend class Point<dim, Number>;
881 template <
int rank,
int dim,
typename T,
typename U>
882 struct ProductTypeImpl<
Tensor<rank, dim, T>,
std::complex<U>>
888 template <
int rank,
int dim,
typename T,
typename U>
889 struct ProductTypeImpl<
Tensor<rank, dim,
std::complex<T>>, std::complex<U>>
895 template <
typename T,
int rank,
int dim,
typename U>
896 struct ProductTypeImpl<
std::complex<T>,
Tensor<rank, dim, U>>
902 template <
int rank,
int dim,
typename T,
typename U>
903 struct ProductTypeImpl<
std::complex<T>,
Tensor<rank, dim, std::complex<U>>>
914 template <
int rank,
int dim,
typename T>
915 struct NumberType<
Tensor<rank, dim,
T>>
938template <
int dim,
typename Number>
948template <
int dim,
typename Number>
949template <
typename OtherNumber>
957template <
int dim,
typename Number>
958template <
typename OtherNumber>
965# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
966template <
int dim,
typename Number>
974template <
int dim,
typename Number>
977 :
value{std::move(other.value)}
983template <
int dim,
typename Number>
988 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
993template <
int dim,
typename Number>
998 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1004template <
int dim,
typename Number>
1005template <
typename OtherNumber>
1014# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1015template <
int dim,
typename Number>
1024# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1025template <
int dim,
typename Number>
1029 value = std::move(other.value);
1036template <
int dim,
typename Number>
1037template <
typename OtherNumber>
1046template <
int dim,
typename Number>
1047template <
typename OtherNumber>
1048constexpr inline bool
1051# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1052 Assert(!(std::is_same_v<Number, adouble> ||
1053 std::is_same_v<OtherNumber, adouble>),
1055 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1056 "been extended to support advanced branching."));
1063template <
int dim,
typename Number>
1064template <
typename OtherNumber>
1068 return !((*this) == p);
1072template <
int dim,
typename Number>
1073template <
typename OtherNumber>
1082template <
int dim,
typename Number>
1083template <
typename OtherNumber>
1095 namespace ComplexWorkaround
1097 template <
typename Number,
typename OtherNumber>
1099 multiply_assign_scalar(Number &val,
const OtherNumber &s)
1104 template <
typename Number,
typename OtherNumber>
1106 multiply_assign_scalar(std::complex<Number> &val,
const OtherNumber &s)
1108# if DEAL_II_KOKKOS_VERSION_GTE(3, 6, 0)
1109 KOKKOS_IF_ON_HOST((val *= s;))
1110 KOKKOS_IF_ON_DEVICE(({
1114 "This function is not implemented for std::complex<Number>!\n");
1117# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1123 "This function is not implemented for std::complex<Number>!\n");
1131template <
int dim,
typename Number>
1132template <
typename OtherNumber>
1136 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1142template <
int dim,
typename Number>
1143template <
typename OtherNumber>
1152template <
int dim,
typename Number>
1160template <
int dim,
typename Number>
1165 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1170template <
int dim,
typename Number>
1176 ExcMessage(
"Cannot access an object of type Tensor<0,0,Number>"));
1182template <
int dim,
typename Number>
1183constexpr inline void
1193template <
int dim,
typename Number>
1194template <
class Iterator>
1201 ExcMessage(
"Cannot unroll an object of type Tensor<0,0,Number>"));
1202 Assert(std::distance(begin, end) >= 1,
1203 ExcMessage(
"The provided iterator range must contain at least one "
1210template <
int dim,
typename Number>
1211template <
class Archive>
1219template <
int dim,
typename Number>
1225template <
int rank_,
int dim,
typename Number>
1226template <
typename ArrayLike, std::size_t... indices>
1229 std::index_sequence<indices...>)
1246 if constexpr (rank_ == 1)
1252 static_assert(
sizeof...(indices) == dim,
1253 "dim should match the number of indices");
1257# if defined(DEAL_II_HAVE_CXX20) && !defined(__NVCC__)
1259template <
int rank_,
int dim,
typename Number>
1294 []<
std::size_t... I>(
1295 const
std::index_sequence<I...> &) constexpr -> decltype(
values) {
1296 if constexpr (dim == 0)
1300 else if constexpr (rank_ == 1)
1302 auto get_zero_and_ignore_argument = [](
int) {
1305 return {{(get_zero_and_ignore_argument(I))...}};
1309 auto get_zero_and_ignore_argument = [](
int) {
1310 return Tensor<rank_ - 1, dim, Number>();
1312 return {{(get_zero_and_ignore_argument(I))...}};
1314 }(std::make_index_sequence<dim>()))
1330 namespace TensorInitialization
1332 template <
int rank,
int dim,
typename Number, std::size_t... I>
1333# if DEAL_II_KOKKOS_VERSION_GTE(3, 7, 0)
1334 constexpr Kokkos::Array<typename Tensor<rank, dim, Number>::value_type, dim>
1336 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1338 make_zero_array(
const std::index_sequence<I...> &)
1340 static_assert(
sizeof...(I) == dim,
"This is bad.");
1346 if constexpr (dim == 0)
1350 else if constexpr (rank == 1)
1352 auto get_zero_and_ignore_argument = [](
int) {
1355 return {{(get_zero_and_ignore_argument(I))...}};
1359 auto get_zero_and_ignore_argument = [](
int) {
1360 return Tensor<rank - 1, dim, Number>();
1362 return {{(get_zero_and_ignore_argument(I))...}};
1369template <
int rank_,
int dim,
typename Number>
1372 :
values(
internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1373 std::make_index_sequence<dim>()))
1380template <
int rank_,
int dim,
typename Number>
1383 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1388template <
int rank_,
int dim,
typename Number>
1389template <
typename ElementType,
typename MemorySpace>
1395 const int my_n_independent_components = n_independent_components;
1398 for (
unsigned int i = 0; i < my_n_independent_components; ++i)
1399 (*
this)[unrolled_to_component_indices(i)] = initializer[i];
1404template <
int rank_,
int dim,
typename Number>
1405template <
typename OtherNumber>
1409 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1414template <
int rank_,
int dim,
typename Number>
1415template <
typename OtherNumber>
1419 :
Tensor(initializer,
std::make_index_sequence<dim>{})
1424template <
int rank_,
int dim,
typename Number>
1425template <
typename OtherNumber>
1427operator
Tensor<1, dim,
Tensor<rank_ - 1, dim, OtherNumber>>()
const
1435# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1436template <
int rank_,
int dim,
typename Number>
1444template <
int rank_,
int dim,
typename Number>
1447 :
values(std::move(other.values))
1453template <
int rank_,
int dim,
typename Number>
1459 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1467template <
int rank_,
int dim,
typename Number>
1473 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1481template <
int rank_,
int dim,
typename Number>
1486 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1488 return TensorAccessors::extract<rank_>(*
this, indices);
1493template <
int rank_,
int dim,
typename Number>
1498 ExcMessage(
"Cannot access an object of type Tensor<rank_,0,Number>"));
1500 return TensorAccessors::extract<rank_>(*
this, indices);
1505template <
int rank_,
int dim,
typename Number>
1509 static_assert(rank_ == 1,
1510 "This function is only available for rank-1 tensors "
1511 "because higher-rank tensors may not store their elements "
1512 "in a contiguous array.");
1514 return std::addressof(
1515 this->
operator[](this->unrolled_to_component_indices(0)));
1520template <
int rank_,
int dim,
typename Number>
1521inline const Number *
1524 static_assert(rank_ == 1,
1525 "This function is only available for rank-1 tensors "
1526 "because higher-rank tensors may not store their elements "
1527 "in a contiguous array.");
1529 return std::addressof(
1530 this->
operator[](this->unrolled_to_component_indices(0)));
1535template <
int rank_,
int dim,
typename Number>
1539 static_assert(rank_ == 1,
1540 "This function is only available for rank-1 tensors "
1541 "because higher-rank tensors may not store their elements "
1542 "in a contiguous array.");
1544 return begin_raw() + n_independent_components;
1549template <
int rank_,
int dim,
typename Number>
1550inline const Number *
1553 static_assert(rank_ == 1,
1554 "This function is only available for rank-1 tensors "
1555 "because higher-rank tensors may not store their elements "
1556 "in a contiguous array.");
1558 return begin_raw() + n_independent_components;
1563template <
int rank_,
int dim,
typename Number>
1564template <
typename OtherNumber>
1570 for (
unsigned int i = 0; i < dim; ++i)
1577template <
int rank_,
int dim,
typename Number>
1585 for (
unsigned int i = 0; i < dim; ++i)
1591# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1592template <
int rank_,
int dim,
typename Number>
1596 for (
unsigned int i = 0; i < dim; ++i)
1597 values[i] = other.
values[i];
1603template <
int rank_,
int dim,
typename Number>
1608 for (
unsigned int i = 0; i < dim; ++i)
1609 values[i] = other.
values[i];
1615template <
int rank_,
int dim,
typename Number>
1616template <
typename OtherNumber>
1617constexpr inline bool
1621# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1622 Assert(!(std::is_same_v<Number, adouble> ||
1623 std::is_same_v<OtherNumber, adouble>),
1625 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1626 "been extended to support advanced branching."));
1629 for (
unsigned int i = 0; i < dim; ++i)
1643constexpr inline bool
1650template <
int rank_,
int dim,
typename Number>
1651template <
typename OtherNumber>
1656 return !((*this) == p);
1660template <
int rank_,
int dim,
typename Number>
1661template <
typename OtherNumber>
1667 for (
unsigned int i = 0; i < dim; ++i)
1668 values[i] += p.
values[i];
1673template <
int rank_,
int dim,
typename Number>
1674template <
typename OtherNumber>
1680 for (
unsigned int i = 0; i < dim; ++i)
1681 values[i] -= p.
values[i];
1686template <
int rank_,
int dim,
typename Number>
1687template <
typename OtherNumber>
1692 for (
unsigned int i = 0; i < dim; ++i)
1699template <
int rank_,
int dim,
typename Number>
1700template <
typename OtherNumber>
1705 if constexpr (std::is_integral_v<
1707 std::is_same_v<Number, Differentiation::SD::Expression>)
1710 for (
unsigned int d = 0;
d < dim; ++
d)
1717 const Number inverse_factor = Number(1.) / s;
1718 for (
unsigned int d = 0;
d < dim; ++
d)
1719 values[d] *= inverse_factor;
1726template <
int rank_,
int dim,
typename Number>
1733 for (
unsigned int i = 0; i < dim; ++i)
1734 tmp.
values[i] = -values[i];
1740template <
int rank_,
int dim,
typename Number>
1746 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1750 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1763 return sqrt(norm_square());
1768template <
int rank_,
int dim,
typename Number>
1773 if constexpr (dim == 0)
1776 else if constexpr (rank_ == 1)
1782 for (
unsigned int i = 1; i < dim; ++i)
1793 for (
unsigned int i = 1; i < dim; ++i)
1794 s += values[i].norm_square();
1802template <
int rank_,
int dim,
typename Number>
1803template <
class Iterator>
1806 const Iterator end)
const
1808 if constexpr (rank_ > 1)
1811 Iterator next =
begin;
1812 for (
unsigned int i = 0; i < dim; ++i)
1814 values[i].unroll(next, end);
1824 Assert(std::distance(begin, end) >= dim,
1826 "The provided iterator range must contain at least 'dim' "
1834template <
int rank_,
int dim,
typename Number>
1835constexpr inline unsigned int
1839 unsigned int index = 0;
1840 for (
int r = 0; r < rank_; ++r)
1841 index = index * dim + indices[r];
1848template <
int rank_,
int dim,
typename Number>
1853 unsigned int dummy = n_independent_components;
1857 if constexpr (dim == 0)
1861 "A tensor with dimension 0 does not store any elements. "
1862 "There is no indexing that can address its elements."));
1869 unsigned int remainder = i;
1870 for (
int r = rank_ - 1; r >= 0; --r)
1872 indices[r] = remainder % dim;
1873 remainder = remainder / dim;
1882template <
int rank_,
int dim,
typename Number>
1883constexpr inline void
1886 for (
unsigned int i = 0; i < dim; ++i)
1891template <
int rank_,
int dim,
typename Number>
1892constexpr std::size_t
1899template <
int rank_,
int dim,
typename Number>
1900template <
class Archive>
1904 for (
int i = 0; i < dim; ++i)
1911template <
int rank_,
int dim,
typename Number>
1930template <
int rank_,
int dim,
typename Number>
1931inline std::ostream &
1934 for (
unsigned int i = 0; i < dim; ++i)
1938 for (
unsigned int j = 0; j < rank_; ++j)
1931inline std::ostream & {
…}
1952template <
int dim,
typename Number>
1953inline std::ostream &
1956 out << static_cast<const Number &>(p);
1953inline std::ostream & {
…}
1979template <
int dim,
typename Number,
typename Other>
1984 return object *
static_cast<const Number &
>(t);
1999template <
int dim,
typename Number,
typename Other>
2004 return static_cast<const Number &
>(t) *
object;
2019template <
int dim,
typename Number,
typename OtherNumber>
2025 return static_cast<const Number &
>(src1) *
2026 static_cast<const OtherNumber &
>(src2);
2037template <
int dim,
typename Number,
typename OtherNumber>
2045 return static_cast<const Number &
>(t) / factor;
2056template <
int dim,
typename Number,
typename OtherNumber>
2062 return static_cast<const Number &
>(p) +
static_cast<const OtherNumber &
>(q);
2073template <
int dim,
typename Number,
typename OtherNumber>
2079 return static_cast<const Number &
>(p) -
static_cast<const OtherNumber &
>(q);
2095template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2121template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2144template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2167template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2188template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2205template <
int dim,
typename Number,
typename OtherNumber>
2234template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2242 for (
unsigned int i = 0; i < dim; ++i)
2293template <
int rank_1,
2297 typename OtherNumber,
2298 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2300 typename Tensor<rank_1 + rank_2 - 2,
2312 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2316 static_assert(dim > 0,
"Tensors cannot have dimension zero.");
2318 for (
unsigned int i = 1; i < dim; ++i)
2319 sum += src1[i] * src2[i];
2323 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2329 rank_1 + rank_2 - 2,
2332 for (
unsigned int i = 0; i < dim; ++i)
2333 result[i] += src1[i] * src2;
2342 rank_1 + rank_2 - 2,
2346 TensorAccessors::internal::
2347 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2348 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2349 TensorAccessors::contract<1, rank_1, rank_2, dim>(result,
2386template <
int index_1,
2392 typename OtherNumber>
2394 typename Tensor<rank_1 + rank_2 - 2,
2400 Assert(0 <= index_1 && index_1 < rank_1,
2402 "The specified index_1 must lie within the range [0,rank_1)"));
2403 Assert(0 <= index_2 && index_2 < rank_2,
2405 "The specified index_2 must lie within the range [0,rank_2)"));
2412 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2418 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2420 typename Tensor<rank_1 + rank_2 - 2,
2424 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2459template <
int index_1,
2467 typename OtherNumber>
2469 typename Tensor<rank_1 + rank_2 - 4,
2475 Assert(0 <= index_1 && index_1 < rank_1,
2477 "The specified index_1 must lie within the range [0,rank_1)"));
2478 Assert(0 <= index_3 && index_3 < rank_1,
2480 "The specified index_3 must lie within the range [0,rank_1)"));
2481 Assert(index_1 != index_3,
2482 ExcMessage(
"index_1 and index_3 must not be the same"));
2483 Assert(0 <= index_2 && index_2 < rank_2,
2485 "The specified index_2 must lie within the range [0,rank_2)"));
2486 Assert(0 <= index_4 && index_4 < rank_2,
2488 "The specified index_4 must lie within the range [0,rank_2)"));
2489 Assert(index_2 != index_4,
2490 ExcMessage(
"index_2 and index_4 must not be the same"));
2497 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2501 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2507 (index_3 < index_1 ? index_3 : index_3 - 1),
2519 (index_4 < index_2 ? index_4 : index_4 - 1),
2527 typename Tensor<rank_1 + rank_2 - 4,
2531 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2548template <
int rank,
int dim,
typename Number,
typename OtherNumber>
2555 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2577template <
template <
int,
int,
typename>
class TensorT1,
2578 template <
int,
int,
typename>
2580 template <
int,
int,
typename>
2590 contract3(
const TensorT1<rank_1, dim, T1> &left,
2591 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2592 const TensorT3<rank_2, dim, T3> &right)
2596 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2612template <
int rank_1,
2616 typename OtherNumber>
2622 typename Tensor<rank_1 + rank_2,
2626 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2649template <
int dim,
typename Number>
2658 result[1] = -src[0];
2673template <
int dim,
typename Number1,
typename Number2>
2684 constexpr int s0 = 0 % dim;
2685 constexpr int s1 = 1 % dim;
2686 constexpr int s2 = 2 % dim;
2688 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2689 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2690 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2709template <
int dim,
typename Number>
2717 for (
unsigned int k = 0; k < dim; ++k)
2719 Tensor<2, dim - 1, Number> minor;
2720 for (
unsigned int i = 0; i < dim - 1; ++i)
2721 for (
unsigned int j = 0; j < dim - 1; ++j)
2722 minor[i][j] = t[i][j < k ? j : j + 1];
2724 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) *
determinant(minor);
2726 det += t[dim - 1][k] * cofactor;
2729 return ((dim % 2 == 0) ? 1. : -1.) * det;
2737template <
typename Number>
2749template <
typename Number>
2754 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2762template <
typename Number>
2773 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2783template <
int dim,
typename Number>
2788 for (
unsigned int i = 1; i < dim; ++i)
2802template <
int dim,
typename Number>
2806 Number return_tensor[dim][dim];
2819template <
typename Number>
2827 return return_tensor;
2831template <
typename Number>
2838 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2839 return_tensor[0][0] = t[1][1];
2840 return_tensor[0][1] = -t[0][1];
2841 return_tensor[1][0] = -t[1][0];
2842 return_tensor[1][1] = t[0][0];
2843 return_tensor *= inv_det_t;
2845 return return_tensor;
2848template <
typename Number>
2854 const auto value = [](
const auto &t) {
2858 return_tensor[0][0] =
value(t[1][1] * t[2][2]) -
value(t[1][2] * t[2][1]);
2859 return_tensor[0][1] =
value(t[0][2] * t[2][1]) -
value(t[0][1] * t[2][2]);
2860 return_tensor[0][2] =
value(t[0][1] * t[1][2]) -
value(t[0][2] * t[1][1]);
2861 return_tensor[1][0] =
value(t[1][2] * t[2][0]) -
value(t[1][0] * t[2][2]);
2862 return_tensor[1][1] =
value(t[0][0] * t[2][2]) -
value(t[0][2] * t[2][0]);
2863 return_tensor[1][2] =
value(t[0][2] * t[1][0]) -
value(t[0][0] * t[1][2]);
2864 return_tensor[2][0] =
value(t[1][0] * t[2][1]) -
value(t[1][1] * t[2][0]);
2865 return_tensor[2][1] =
value(t[0][1] * t[2][0]) -
value(t[0][0] * t[2][1]);
2866 return_tensor[2][2] =
value(t[0][0] * t[1][1]) -
value(t[0][1] * t[1][0]);
2868 const Number inv_det_t =
2869 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2870 t[0][2] * return_tensor[2][0]));
2871 return_tensor *= inv_det_t;
2873 return return_tensor;
2884template <
int dim,
typename Number>
2889 for (
unsigned int i = 0; i < dim; ++i)
2892 for (
unsigned int j = i + 1; j < dim; ++j)
2915template <
int dim,
typename Number>
2936template <
int dim,
typename Number>
3007template <
int dim,
typename Number>
3019template <
int dim,
typename Number>
3024 for (
unsigned int j = 0; j < dim; ++j)
3027 for (
unsigned int i = 0; i < dim; ++i)
3045template <
int dim,
typename Number>
3050 for (
unsigned int i = 0; i < dim; ++i)
3053 for (
unsigned int j = 0; j < dim; ++j)
3071# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3080 for (
unsigned int j = 0; j < dim; ++j)
3083 for (
unsigned int i = 0; i < dim; ++i)
3084 sum += fabs(t[i][j]);
3086 condassign(max, (sum > max), sum, max);
3098 for (
unsigned int i = 0; i < dim; ++i)
3101 for (
unsigned int j = 0; j < dim; ++j)
3102 sum +=
fabs(t[i][j]);
3104 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, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
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 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)
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< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
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_HOST_DEVICE
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
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 > &)
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr 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)