16 #ifndef dealii_symmetric_tensor_h 17 #define dealii_symmetric_tensor_h 35 template <
int rank,
int dim,
typename Number =
double>
39 template <
int dim,
typename Number>
43 template <
int dim,
typename Number>
47 template <
int dim,
typename Number>
51 template <
int dim,
typename Number>
55 template <
int dim,
typename Number>
59 template <
int dim2,
typename Number>
63 template <
int dim,
typename Number>
67 template <
int dim,
typename Number>
79 template <
int rank,
int dim,
typename T,
typename U>
85 std::complex<typename ProductType<T, U>::type>>;
88 template <
int rank,
int dim,
typename T,
typename U>
95 std::complex<typename ProductType<T, U>::type>>;
98 template <
typename T,
int rank,
int dim,
typename U>
104 std::complex<typename ProductType<T, U>::type>>;
107 template <
int rank,
int dim,
typename T,
typename U>
114 std::complex<typename ProductType<T, U>::type>>;
122 namespace SymmetricTensorImplementation
128 template <
int rank,
int dim,
typename Number>
136 namespace SymmetricTensorAccessors
146 const unsigned int new_index,
147 const unsigned int position)
154 return {previous_indices[0], new_index};
167 const unsigned int new_index,
168 const unsigned int position)
178 numbers::invalid_unsigned_int};
180 return {previous_indices[0],
183 numbers::invalid_unsigned_int};
185 return {previous_indices[0],
188 numbers::invalid_unsigned_int};
190 return {previous_indices[0],
211 typename OtherNumber = Number>
226 template <
int dim,
typename Number,
typename OtherNumber>
246 template <
int rank,
int dim,
typename Number>
252 template <
int dim,
typename Number>
259 static const unsigned int n_independent_components =
260 (dim * dim + dim) / 2;
273 template <
int dim,
typename Number>
281 static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
286 static const unsigned int n_independent_components =
287 (n_rank2_components *
305 template <
int rank,
int dim,
bool constness,
typename Number>
314 template <
int rank,
int dim,
typename Number>
328 template <
int rank,
int dim,
typename Number>
369 template <
int rank,
int dim,
bool constness,
int P,
typename Number>
414 operator[](
const unsigned int i);
419 constexpr
Accessor<rank, dim, constness, P - 1, Number>
420 operator[](
const unsigned int i)
const;
431 template <
int,
int,
typename>
432 friend class ::SymmetricTensor;
433 template <
int,
int,
bool,
int,
typename>
435 friend class ::SymmetricTensor<rank, dim, Number>;
436 friend class Accessor<rank, dim, constness, P + 1, Number>;
448 template <
int rank,
int dim,
bool constness,
typename Number>
500 constexpr
reference operator[](
const unsigned int)
const;
511 template <
int,
int,
typename>
512 friend class ::SymmetricTensor;
513 template <
int,
int,
bool,
int,
typename>
515 friend class ::SymmetricTensor<rank, dim, Number>;
516 friend class SymmetricTensorAccessors::
596 template <int rank_, int dim, typename Number>
600 static_assert(rank_ % 2 == 0,
"A SymmetricTensor must have even rank!");
610 static const unsigned int dimension = dim;
615 static const unsigned int rank = rank_;
622 static constexpr
unsigned int n_independent_components =
624 n_independent_components;
645 template <
typename OtherNumber>
671 template <
typename OtherNumber>
706 template <
typename OtherNumber>
717 operator=(
const Number &
d);
740 template <
typename OtherNumber>
747 template <
typename OtherNumber>
755 template <
typename OtherNumber>
757 operator*=(
const OtherNumber &factor);
762 template <
typename OtherNumber>
764 operator/=(
const OtherNumber &factor);
798 template <
typename OtherNumber>
807 template <
typename OtherNumber>
828 constexpr internal::SymmetricTensorAccessors::
829 Accessor<rank_, dim,
true, rank_ - 1, Number>
830 operator[](
const unsigned int row)
const;
837 Accessor<rank_, dim,
false, rank_ - 1, Number>
838 operator[](
const unsigned int row);
861 access_raw_entry(
const unsigned int unrolled_index)
const;
870 access_raw_entry(
const unsigned int unrolled_index);
891 static constexpr
unsigned int 900 unrolled_to_component_indices(
const unsigned int i);
921 static constexpr std::size_t
929 template <
class Archive>
931 serialize(Archive &ar,
const unsigned int version);
951 template <
int,
int,
typename>
955 template <
int dim2,
typename Number2>
959 template <
int dim2,
typename Number2>
963 template <
int dim2,
typename Number2>
967 template <
int dim2,
typename Number2>
971 template <
int dim2,
typename Number2>
975 template <
int dim2,
typename Number2>
982 Inverse<2, dim, Number>;
985 Inverse<4, dim, Number>;
995 template <int rank, int dim, typename Number>
998 template <int rank_, int dim, typename Number>
999 constexpr unsigned
int 1004 namespace SymmetricTensorAccessors
1006 template <
int rank_,
int dim,
bool constness,
int P,
typename Number>
1008 Accessor<rank_, dim, constness, P, Number>::Accessor(
1009 tensor_type & tensor,
1012 , previous_indices(previous_indices)
1017 template <
int rank_,
int dim,
bool constness,
int P,
typename Number>
1019 Accessor<rank_, dim, constness, P - 1, Number>
1020 Accessor<rank_, dim, constness, P, Number>::
1021 operator[](
const unsigned int i)
1023 return Accessor<rank_, dim, constness, P - 1, Number>(
1024 tensor,
merge(previous_indices, i, rank_ - P));
1029 template <
int rank_,
int dim,
bool constness,
int P,
typename Number>
1031 Accessor<rank_, dim, constness, P - 1, Number>
1032 Accessor<rank_, dim, constness, P, Number>::
1033 operator[](
const unsigned int i)
const 1035 return Accessor<rank_, dim, constness, P - 1, Number>(
1036 tensor,
merge(previous_indices, i, rank_ - P));
1041 template <
int rank_,
int dim,
bool constness,
typename Number>
1043 Accessor<rank_, dim, constness, 1, Number>::Accessor(
1044 tensor_type & tensor,
1047 , previous_indices(previous_indices)
1052 template <
int rank_,
int dim,
bool constness,
typename Number>
1054 typename Accessor<rank_, dim, constness, 1, Number>::reference
1055 Accessor<rank_, dim, constness, 1, Number>::
1056 operator[](
const unsigned int i)
1058 return tensor(
merge(previous_indices, i, rank_ - 1));
1062 template <
int rank_,
int dim,
bool constness,
typename Number>
1064 typename Accessor<rank_, dim, constness, 1, Number>::reference
1065 Accessor<rank_, dim, constness, 1, Number>::
1066 operator[](
const unsigned int i)
const 1068 return tensor(
merge(previous_indices, i, rank_ - 1));
1075 template <
int rank_,
int dim,
typename Number>
1076 template <
typename OtherNumber>
1081 static_assert(rank == 2,
"This function is only implemented for rank==2");
1082 for (
unsigned int d = 0;
d < dim; ++
d)
1083 for (
unsigned int e = 0;
e <
d; ++
e)
1085 ExcMessage(
"The incoming Tensor must be exactly symmetric."));
1087 for (
unsigned int d = 0; d < dim; ++
d)
1090 for (
unsigned int d = 0, c = 0; d < dim; ++
d)
1091 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
1092 data[dim + c] = t[d][
e];
1097 template <
int rank_,
int dim,
typename Number>
1098 template <
typename OtherNumber>
1102 : data(initializer.
data)
1107 template <
int rank_,
int dim,
typename Number>
1110 const Number (&array)[n_independent_components])
1112 *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1115 Assert(
sizeof(
typename base_tensor_type::array_type) ==
sizeof(array),
1121 template <
int rank_,
int dim,
typename Number>
1122 template <
typename OtherNumber>
1134 template <
int rank_,
int dim,
typename Number>
1140 ExcMessage(
"Only assignment with zero is allowed"));
1151 namespace SymmetricTensorImplementation
1153 template <
int dim,
typename Number>
1155 ::Tensor<2, dim, Number>
1156 convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1161 for (
unsigned int d = 0; d < dim; ++
d)
1162 t[d][d] = s.access_raw_entry(d);
1165 for (
unsigned int d = 0, c = 0; d < dim; ++
d)
1166 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
1168 t[
d][
e] = s.access_raw_entry(dim + c);
1169 t[
e][
d] = s.access_raw_entry(dim + c);
1175 template <
int dim,
typename Number>
1176 DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1177 convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1184 for (
unsigned int i = 0; i < dim; ++i)
1185 for (
unsigned int j = i; j < dim; ++j)
1186 for (
unsigned int k = 0; k < dim; ++k)
1187 for (
unsigned int l = k;
l < dim; ++
l)
1197 template <
typename Number>
1198 struct Inverse<2, 1, Number>
1201 ::SymmetricTensor<2, 1, Number>
1202 value(const ::SymmetricTensor<2, 1, Number> &t)
1206 tmp[0][0] = 1.0 / t[0][0];
1213 template <
typename Number>
1214 struct Inverse<2, 2, Number>
1217 ::SymmetricTensor<2, 2, Number>
1218 value(const ::SymmetricTensor<2, 2, Number> &t)
1228 const Number inv_det_t =
1229 1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1230 tmp[idx_00] = t[idx_11];
1231 tmp[idx_01] = -t[idx_01];
1232 tmp[idx_11] = t[idx_00];
1240 template <
typename Number>
1241 struct Inverse<2, 3, Number>
1244 value(const ::SymmetricTensor<2, 3, Number> &t)
1288 const Number inv_det_t =
1289 1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1290 t[idx_00] * t[idx_12] * t[idx_12] -
1291 t[idx_01] * t[idx_01] * t[idx_22] +
1292 2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1293 t[idx_02] * t[idx_02] * t[idx_11]);
1294 tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1295 tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1296 tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1297 tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1298 tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1299 tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1307 template <
typename Number>
1308 struct Inverse<4, 1, Number>
1311 value(const ::SymmetricTensor<4, 1, Number> &t)
1314 tmp.
data[0][0] = 1.0 / t.data[0][0];
1320 template <
typename Number>
1321 struct Inverse<4, 2, Number>
1324 value(const ::SymmetricTensor<4, 2, Number> &t)
1350 const Number t4 = t.data[0][0] * t.data[1][1],
1351 t6 = t.data[0][0] * t.data[1][2],
1352 t8 = t.data[0][1] * t.data[1][0],
1353 t00 = t.data[0][2] * t.data[1][0],
1354 t01 = t.data[0][1] * t.data[2][0],
1355 t04 = t.data[0][2] * t.data[2][0],
1356 t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1357 t8 * t.data[2][2] + t00 * t.data[2][1] +
1358 t01 * t.data[1][2] - t04 * t.data[1][1]);
1360 (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1362 -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1364 -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1366 -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1367 tmp.
data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1368 tmp.
data[1][2] = -(t6 - t00) * t07;
1370 -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1371 tmp.
data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1372 tmp.
data[2][2] = (t4 - t8) * t07;
1376 tmp.
data[2][0] /= 2;
1377 tmp.
data[2][1] /= 2;
1378 tmp.
data[0][2] /= 2;
1379 tmp.
data[1][2] /= 2;
1380 tmp.
data[2][2] /= 4;
1387 template <
typename Number>
1388 struct Inverse<4, 3, Number>
1390 static ::SymmetricTensor<4, 3, Number>
1391 value(const ::SymmetricTensor<4, 3, Number> &t)
1401 const unsigned int N = 6;
1407 for (
unsigned int i = 0; i <
N; ++i)
1409 const Number typical_diagonal_element =
1410 diagonal_sum /
static_cast<double>(
N);
1411 (void)typical_diagonal_element;
1414 for (
unsigned int i = 0; i <
N; ++i)
1417 for (
unsigned int j = 0; j <
N; ++j)
1423 for (
unsigned int i = j + 1; i <
N; ++i)
1431 Assert(max > 1.
e-16 * typical_diagonal_element,
1432 ExcMessage(
"This tensor seems to be noninvertible"));
1437 for (
unsigned int k = 0; k <
N; ++k)
1444 const Number hr = 1. / tmp.
data[j][j];
1445 tmp.
data[j][j] = hr;
1446 for (
unsigned int k = 0; k <
N; ++k)
1450 for (
unsigned int i = 0; i <
N; ++i)
1454 tmp.
data[i][k] -= tmp.
data[i][j] * tmp.
data[j][k] * hr;
1457 for (
unsigned int i = 0; i <
N; ++i)
1459 tmp.
data[i][j] *= hr;
1460 tmp.
data[j][i] *= -hr;
1462 tmp.
data[j][j] = hr;
1467 for (
unsigned int i = 0; i <
N; ++i)
1469 for (
unsigned int k = 0; k <
N; ++k)
1470 hv[p[k]] = tmp.
data[i][k];
1471 for (
unsigned int k = 0; k <
N; ++k)
1472 tmp.
data[i][k] = hv[k];
1477 for (
unsigned int i = 3; i < 6; ++i)
1478 for (
unsigned int j = 0; j < 3; ++j)
1479 tmp.
data[i][j] /= 2;
1481 for (
unsigned int i = 0; i < 3; ++i)
1482 for (
unsigned int j = 3; j < 6; ++j)
1483 tmp.
data[i][j] /= 2;
1485 for (
unsigned int i = 3; i < 6; ++i)
1486 for (
unsigned int j = 3; j < 6; ++j)
1487 tmp.
data[i][j] /= 4;
1498 template <
int rank_,
int dim,
typename Number>
1502 return internal::SymmetricTensorImplementation::convert_to_tensor(*
this);
1507 template <
int rank_,
int dim,
typename Number>
1512 return data == t.
data;
1517 template <
int rank_,
int dim,
typename Number>
1522 return data != t.
data;
1527 template <
int rank_,
int dim,
typename Number>
1528 template <
typename OtherNumber>
1540 template <
int rank_,
int dim,
typename Number>
1541 template <
typename OtherNumber>
1553 template <
int rank_,
int dim,
typename Number>
1554 template <
typename OtherNumber>
1565 template <
int rank_,
int dim,
typename Number>
1566 template <
typename OtherNumber>
1577 template <
int rank_,
int dim,
typename Number>
1589 template <
int rank_,
int dim,
typename Number>
1598 template <
int rank_,
int dim,
typename Number>
1599 constexpr std::size_t
1611 template <
int dim,
typename Number,
typename OtherNumber = Number>
1615 perform_double_contraction(
1616 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1618 const typename SymmetricTensorAccessors::
1619 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1627 return data[0] * sdata[0];
1632 result_type
sum = data[dim] * sdata[dim];
1633 for (
unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++
d)
1634 sum += data[d] * sdata[d];
1638 for (
unsigned int d = 0; d < dim; ++
d)
1639 sum += data[d] * sdata[d];
1646 template <
int dim,
typename Number,
typename OtherNumber = Number>
1650 perform_double_contraction(
1651 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1653 const typename SymmetricTensorAccessors::
1654 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1661 const unsigned int data_dim = SymmetricTensorAccessors::
1662 StorageType<2, dim, value_type>::n_independent_components;
1663 value_type tmp[data_dim]{};
1664 for (
unsigned int i = 0; i < data_dim; ++i)
1666 perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1667 return result_type(tmp);
1672 template <
int dim,
typename Number,
typename OtherNumber = Number>
1674 typename SymmetricTensorAccessors::StorageType<
1680 perform_double_contraction(
1681 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1683 const typename SymmetricTensorAccessors::
1684 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1689 StorageType<2, dim, value_type>::base_tensor_type;
1692 for (
unsigned int i = 0; i < tmp.dimension; ++i)
1695 value_type
sum = data[dim] * sdata[dim][i];
1696 for (
unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++
d)
1697 sum += data[d] * sdata[d][i];
1701 for (
unsigned int d = 0; d < dim; ++
d)
1702 sum += data[d] * sdata[d][i];
1710 template <
int dim,
typename Number,
typename OtherNumber = Number>
1712 typename SymmetricTensorAccessors::StorageType<
1718 perform_double_contraction(
1719 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1721 const typename SymmetricTensorAccessors::
1722 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1727 StorageType<4, dim, value_type>::base_tensor_type;
1729 const unsigned int data_dim = SymmetricTensorAccessors::
1730 StorageType<2, dim, value_type>::n_independent_components;
1732 for (
unsigned int i = 0; i < data_dim; ++i)
1733 for (
unsigned int j = 0; j < data_dim; ++j)
1736 for (
unsigned int d = dim; d < (dim * (dim + 1) / 2); ++
d)
1737 tmp[i][j] += data[i][d] * sdata[d][j];
1738 tmp[i][j] += tmp[i][j];
1741 for (
unsigned int d = 0; d < dim; ++
d)
1742 tmp[i][j] += data[i][d] * sdata[d][j];
1751 template <
int rank_,
int dim,
typename Number>
1752 template <
typename OtherNumber>
1763 return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1769 template <
int rank_,
int dim,
typename Number>
1770 template <
typename OtherNumber>
1779 internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1802 template <
typename Type>
1803 struct Uninitialized
1808 template <
typename Type>
1811 template <
int dim,
typename Number>
1814 typename SymmetricTensorAccessors::
1815 StorageType<2, dim, Number>::base_tensor_type &data)
1823 if (indices[0] == indices[1])
1824 return data[indices[0]];
1831 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1832 ((indices[0] == 0) && (indices[1] == 1)),
1841 for (
unsigned int d = 0, c = 0; d < dim; ++
d)
1842 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
1843 if ((sorted_indices[0] == d) && (sorted_indices[1] ==
e))
1844 return data[dim + c];
1857 template <
int dim,
typename Number>
1860 const typename SymmetricTensorAccessors::
1861 StorageType<2, dim, Number>::base_tensor_type &data)
1869 if (indices[0] == indices[1])
1870 return data[indices[0]];
1877 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1878 ((indices[0] == 0) && (indices[1] == 1)),
1887 for (
unsigned int d = 0, c = 0; d < dim; ++
d)
1888 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
1889 if ((sorted_indices[0] == d) && (sorted_indices[1] ==
e))
1890 return data[dim + c];
1903 template <
int dim,
typename Number>
1906 typename SymmetricTensorAccessors::
1907 StorageType<4, dim, Number>::base_tensor_type &data)
1921 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1922 return data[base_index[indices[0]][indices[1]]]
1923 [base_index[indices[2]][indices[3]]];
1932 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1935 return data[base_index[indices[0]][indices[1]]]
1936 [base_index[indices[2]][indices[3]]];
1950 template <
int dim,
typename Number>
1953 const typename SymmetricTensorAccessors::
1954 StorageType<4, dim, Number>::base_tensor_type &data)
1968 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1969 return data[base_index[indices[0]][indices[1]]]
1970 [base_index[indices[2]][indices[3]]];
1979 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1982 return data[base_index[indices[0]][indices[1]]]
1983 [base_index[indices[2]][indices[3]]];
2000 template <
int rank_,
int dim,
typename Number>
2005 for (
unsigned int r = 0; r < rank; ++r)
2007 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2012 template <
int rank_,
int dim,
typename Number>
2017 for (
unsigned int r = 0; r < rank; ++r)
2019 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2026 namespace SymmetricTensorImplementation
2028 template <
int rank_>
2030 get_partially_filled_indices(
const unsigned int row,
2031 const std::integral_constant<int, 2> &)
2037 template <
int rank_>
2039 get_partially_filled_indices(
const unsigned int row,
2040 const std::integral_constant<int, 4> &)
2051 template <
int rank_,
int dim,
typename Number>
2053 Accessor<rank_, dim,
true, rank_ - 1, Number>
2057 return internal::SymmetricTensorAccessors::
2058 Accessor<rank_, dim,
true, rank_ - 1, Number>(
2060 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2061 rank_>(row, std::integral_constant<int, rank_>()));
2066 template <
int rank_,
int dim,
typename Number>
2068 SymmetricTensorAccessors::Accessor<rank_, dim,
false, rank_ - 1, Number>
2071 return internal::SymmetricTensorAccessors::
2072 Accessor<rank_, dim,
false, rank_ - 1, Number>(
2074 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2075 rank_>(row, std::integral_constant<int, rank_>()));
2080 template <
int rank_,
int dim,
typename Number>
2085 return operator()(indices);
2090 template <
int rank_,
int dim,
typename Number>
2095 return operator()(indices);
2100 template <
int rank_,
int dim,
typename Number>
2104 return std::addressof(this->access_raw_entry(0));
2109 template <
int rank_,
int dim,
typename Number>
2110 inline const Number *
2113 return std::addressof(this->access_raw_entry(0));
2118 template <
int rank_,
int dim,
typename Number>
2122 return begin_raw() + n_independent_components;
2127 template <
int rank_,
int dim,
typename Number>
2128 inline const Number *
2131 return begin_raw() + n_independent_components;
2138 namespace SymmetricTensorImplementation
2140 template <
int dim,
typename Number>
2141 constexpr
unsigned int 2142 entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2143 const unsigned int index)
2149 template <
int dim,
typename Number>
2150 constexpr ::TableIndices<2>
2151 entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2152 const unsigned int index)
2163 template <
int rank_,
int dim,
typename Number>
2166 const unsigned int index)
const 2169 return data[internal::SymmetricTensorImplementation::entry_to_indices(*
this,
2175 template <
int rank_,
int dim,
typename Number>
2180 return data[internal::SymmetricTensorImplementation::entry_to_indices(*
this,
2188 template <
int dim,
typename Number>
2190 compute_norm(
const typename SymmetricTensorAccessors::
2191 StorageType<2, dim, Number>::base_tensor_type &data)
2218 for (
unsigned int d = 0; d < dim; ++
d)
2221 for (
unsigned int d = dim; d < (dim * dim + dim) / 2; ++
d)
2225 return std::sqrt(return_value);
2232 template <
int dim,
typename Number>
2234 compute_norm(
const typename SymmetricTensorAccessors::
2235 StorageType<4, dim, Number>::base_tensor_type &data)
2247 const unsigned int n_independent_components = data.dimension;
2249 for (
unsigned int i = 0; i < dim; ++i)
2250 for (
unsigned int j = 0; j < dim; ++j)
2253 for (
unsigned int i = 0; i < dim; ++i)
2254 for (
unsigned int j = dim; j < n_independent_components; ++j)
2257 for (
unsigned int i = dim; i < n_independent_components; ++i)
2258 for (
unsigned int j = 0; j < dim; ++j)
2261 for (
unsigned int i = dim; i < n_independent_components; ++i)
2262 for (
unsigned int j = dim; j < n_independent_components; ++j)
2266 return std::sqrt(return_value);
2275 template <
int rank_,
int dim,
typename Number>
2279 return internal::compute_norm<dim, Number>(data);
2286 namespace SymmetricTensorImplementation
2309 constexpr
unsigned int table[2][2] = {{0, 2}, {2, 1}};
2310 return table[indices[0]][indices[1]];
2315 constexpr
unsigned int table[3][3] = {{0, 3, 4},
2318 return table[indices[0]][indices[1]];
2323 constexpr
unsigned int table[4][4] = {{0, 4, 5, 6},
2327 return table[indices[0]][indices[1]];
2333 if (indices[0] == indices[1])
2337 sorted_indices.
sort();
2339 for (
unsigned int d = 0, c = 0; d < dim; ++
d)
2340 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
2341 if ((sorted_indices[0] == d) && (sorted_indices[1] ==
e))
2357 template <
int dim,
int rank_>
2369 template <
int rank_,
int dim,
typename Number>
2370 constexpr
unsigned int 2374 return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2382 namespace SymmetricTensorImplementation
2393 unrolled_to_component_indices(
const unsigned int i,
2394 const std::integral_constant<int, 2> &)
2432 for (
unsigned int d = 0, c = dim; d < dim; ++
d)
2433 for (
unsigned int e = d + 1;
e < dim; ++
e, ++c)
2451 template <
int dim,
int rank_>
2453 typename std::enable_if<rank_ != 2, TableIndices<rank_>>
::type 2454 unrolled_to_component_indices(
const unsigned int i,
2455 const std::integral_constant<int, rank_> &)
2464 n_independent_components));
2472 template <
int rank_,
int dim,
typename Number>
2475 const unsigned int i)
2477 return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2478 dim>(i, std::integral_constant<int, rank_>());
2483 template <
int rank_,
int dim,
typename Number>
2484 template <
class Archive>
2509 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2534 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2554 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2571 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2588 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2605 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
2629 template <
int dim,
typename Number>
2645 return (tmp + tmp + t.
data[0] * t.
data[1] * t.
data[2] -
2669 template <
int dim,
typename Number>
2687 template <
int dim,
typename Number>
2691 Number t = d.
data[0];
2692 for (
unsigned int i = 1; i < dim; ++i)
2709 template <
int dim,
typename Number>
2728 template <
typename Number>
2755 template <
typename Number>
2759 return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2772 template <
typename Number>
2776 return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2777 t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2789 template <
typename Number>
2790 std::array<Number, 1>
2817 template <
typename Number>
2818 std::array<Number, 2>
2845 template <
typename Number>
2846 std::array<Number, 3>
2853 namespace SymmetricTensorImplementation
2892 template <
int dim,
typename Number>
2896 std::array<Number, dim> & d,
2897 std::array<Number, dim - 1> &
e);
2940 template <
int dim,
typename Number>
2941 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2942 ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2985 template <
int dim,
typename Number>
2986 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3004 template <
typename Number>
3005 std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3006 hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3042 template <
typename Number>
3043 std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3044 hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3050 template <
int dim,
typename Number>
3057 return lhs.first > rhs.first;
3160 template <
int dim,
typename Number>
3161 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3162 std::integral_constant<int, dim>::value>
3177 template <
int rank_,
int dim,
typename Number>
3196 template <
int dim,
typename Number>
3203 const Number tr =
trace(t) / dim;
3204 for (
unsigned int i = 0; i < dim; ++i)
3218 template <
int dim,
typename Number>
3238 for (
unsigned int d = 0; d < dim; ++
d)
3257 return unit_symmetric_tensor<dim, double>();
3290 template <
int dim,
typename Number>
3297 for (
unsigned int i = 0; i < dim; ++i)
3298 for (
unsigned int j = 0; j < dim; ++j)
3307 for (
unsigned int i = dim;
3308 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3329 return deviator_tensor<dim, double>();
3371 template <
int dim,
typename Number>
3378 for (
unsigned int i = 0; i < dim; ++i)
3386 for (
unsigned int i = dim;
3387 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3408 return identity_tensor<dim, double>();
3422 template <
int dim,
typename Number>
3442 template <
int dim,
typename Number>
3473 template <
int dim,
typename Number>
3481 for (
unsigned int i = 0; i < dim; ++i)
3482 for (
unsigned int j = i; j < dim; ++j)
3483 for (
unsigned int k = 0; k < dim; ++k)
3484 for (
unsigned int l = k;
l < dim; ++
l)
3485 tmp[i][j][k][
l] = t1[i][j] * t2[k][
l];
3499 template <
int dim,
typename Number>
3504 for (
unsigned int d = 0; d < dim; ++
d)
3505 result[d][d] = t[d][d];
3508 for (
unsigned int d = 0; d < dim; ++
d)
3509 for (
unsigned int e = d + 1;
e < dim; ++
e)
3510 result[d][
e] = (t[d][
e] + t[
e][d]) * half;
3523 template <
int rank_,
int dim,
typename Number>
3542 template <
int rank_,
int dim,
typename Number>
3577 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
3584 const OtherNumber & factor)
3607 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
3617 return (t * factor);
3627 template <
int rank_,
int dim,
typename Number,
typename OtherNumber>
3634 const OtherNumber & factor)
3650 template <
int rank_,
int dim>
3667 template <
int rank_,
int dim>
3683 template <
int rank_,
int dim>
3701 template <
int dim,
typename Number,
typename OtherNumber>
3723 template <
int dim,
typename Number,
typename OtherNumber>
3731 for (
unsigned int i = 0; i < dim; ++i)
3732 for (
unsigned int j = 0; j < dim; ++j)
3733 s += t1[i][j] * t2[i][j];
3751 template <
int dim,
typename Number,
typename OtherNumber>
3774 template <
typename Number,
typename OtherNumber>
3780 tmp[0][0] = t[0][0][0][0] * s[0][0];
3799 template <
typename Number,
typename OtherNumber>
3805 tmp[0][0] = t[0][0][0][0] * s[0][0];
3824 template <
typename Number,
typename OtherNumber>
3830 const unsigned int dim = 2;
3832 for (
unsigned int i = 0; i < dim; ++i)
3833 for (
unsigned int j = i; j < dim; ++j)
3834 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3835 2 * t[i][j][0][1] * s[0][1];
3854 template <
typename Number,
typename OtherNumber>
3860 const unsigned int dim = 2;
3862 for (
unsigned int i = 0; i < dim; ++i)
3863 for (
unsigned int j = i; j < dim; ++j)
3864 tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3865 2 * s[0][1] * t[0][1][i][j];
3884 template <
typename Number,
typename OtherNumber>
3890 const unsigned int dim = 3;
3892 for (
unsigned int i = 0; i < dim; ++i)
3893 for (
unsigned int j = i; j < dim; ++j)
3894 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3895 t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3896 2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3915 template <
typename Number,
typename OtherNumber>
3921 const unsigned int dim = 3;
3923 for (
unsigned int i = 0; i < dim; ++i)
3924 for (
unsigned int j = i; j < dim; ++j)
3925 tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3926 s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3927 2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3938 template <
int dim,
typename Number,
typename OtherNumber>
3945 for (
unsigned int i = 0; i < dim; ++i)
3946 for (
unsigned int j = 0; j < dim; ++j)
3947 dest[i] += src1[i][j] * src2[j];
3958 template <
int dim,
typename Number,
typename OtherNumber>
3988 template <
int rank_1,
3992 typename OtherNumber>
3994 typename Tensor<rank_1 + rank_2 - 2,
4024 template <
int rank_1,
4028 typename OtherNumber>
4030 typename Tensor<rank_1 + rank_2 - 2,
4050 template <
int dim,
typename Number>
4051 inline std::ostream &
4052 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4059 for (
unsigned int i = 0; i < dim; ++i)
4060 for (
unsigned int j = 0; j < dim; ++j)
4077 template <
int dim,
typename Number>
4078 inline std::ostream &
4079 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4086 for (
unsigned int i = 0; i < dim; ++i)
4087 for (
unsigned int j = 0; j < dim; ++j)
4088 for (
unsigned int k = 0; k < dim; ++k)
4089 for (
unsigned int l = 0;
l < dim; ++
l)
4090 tt[i][j][k][
l] = t[i][j][k][
l];
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
constexpr SymmetricTensor operator-() const
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
static constexpr const T & value(const T &t)
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
constexpr SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr numbers::NumberTraits< Number >::real_type norm() const
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
typename AccessorTypes< rank, dim, constness, Number >::reference reference
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
constexpr bool operator==(const SymmetricTensor &) const
static real_type abs(const number &x)
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
SymmetricTensorEigenvectorMethod
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
const TableIndices< rank > previous_indices
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
typename ProductType< Number, OtherNumber >::type type
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr bool operator!=(const SymmetricTensor &) const
constexpr SymmetricTensor()=default
typename base_tensor_descriptor::base_tensor_type base_tensor_type
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
std::pair< Number, Tensor< 1, dim, Number > > EigValsVecs
bool operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
constexpr TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
const TableIndices< rank > previous_indices
#define DEAL_II_NAMESPACE_CLOSE
typename ProductType< Number, OtherNumber >::type value_type
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
void serialize(Archive &ar, const unsigned int version)
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim - 1 > &e)
#define DEAL_II_ALWAYS_INLINE
typename AccessorTypes< rank, dim, constness, Number >::reference reference
constexpr void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
::SymmetricTensor< rank1+rank2 - 4, dim, value_type > type
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
static constexpr std::size_t memory_consumption()
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
#define DEAL_II_NAMESPACE_OPEN
T min(const T &t, const MPI_Comm &mpi_communicator)
constexpr bool value_is_zero(const Number &value)
decltype(std::declval< T >() *std::declval< U >()) type
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
static ::ExceptionBase & ExcNotImplemented()
constexpr Number & operator()(const TableIndices< rank_ > &indices)
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_CONSTEXPR
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)