16#ifndef dealii_physics_notation_h
17#define dealii_physics_notation_h
294 <<
"The number of rows in the input matrix is " << arg1
295 <<
", but needs to be either " << arg2 <<
" or " << arg3
307 <<
"The number of rows in the input matrix is " << arg1
308 <<
", but needs to be either " << arg2 <<
',' << arg3
309 <<
", or " << arg4 <<
'.');
319 <<
"The number of columns in the input matrix is " << arg1
320 <<
", but needs to be either " << arg2 <<
" or " << arg3
332 <<
"The number of columns in the input matrix is " << arg1
333 <<
", but needs to be either " << arg2 <<
',' << arg3
334 <<
", or " << arg4 <<
'.');
347 template <
typename Number>
357 template <
int dim,
typename Number>
367 template <
int dim,
typename Number>
377 template <
int dim,
typename Number>
388 template <
int dim,
typename Number>
398 template <
typename Number>
408 template <
int dim,
typename Number>
418 template <
int dim,
typename Number>
428 template <
int dim,
typename Number>
442 template <
int dim,
typename Number>
492 template <
int dim,
typename Number>
504 template <
int dim,
typename Number>
518 template <
typename Number>
526 template <
int dim,
typename Number>
534 template <
int dim,
typename Number>
542 template <
int dim,
typename Number>
550 template <
int dim,
typename Number>
558 template <
typename Number>
566 template <
int dim,
typename Number>
574 template <
int dim,
typename Number>
582 template <
int dim,
typename Number>
590 template <
int dim,
typename Number>
605 template <
int dim,
typename Number>
613 template <
int dim,
typename Number>
621 template <
int dim,
typename Number>
631 template <
typename TensorType,
typename Number>
640 template <
typename TensorType,
typename Number>
673 std::pair<unsigned int, unsigned int>
674 indices_from_component(
const unsigned int component_n,
675 const bool symmetric);
679 std::pair<unsigned int, unsigned int>
680 indices_from_component(
const unsigned int component_n,
const bool)
683 return std::make_pair(0u, 0u);
688 inline std::pair<unsigned int, unsigned int>
689 indices_from_component<1>(
const unsigned int component_n,
const bool)
693 return std::make_pair(0u, 0u);
698 inline std::pair<unsigned int, unsigned int>
699 indices_from_component<2>(
const unsigned int component_n,
700 const bool symmetric)
702 if (symmetric ==
true)
718 static const unsigned int indices[4][2] = {{0, 0},
722 return std::make_pair(indices[component_n][0],
723 indices[component_n][1]);
728 inline std::pair<unsigned int, unsigned int>
729 indices_from_component<3>(
const unsigned int component_n,
730 const bool symmetric)
732 if (symmetric ==
true)
748 static const unsigned int indices[9][2] = {{0, 0},
757 return std::make_pair(indices[component_n][0],
758 indices[component_n][1]);
768 vector_component_factor(
const unsigned int component_i,
769 const bool symmetric)
771 if (symmetric ==
false)
774 if (component_i < dim)
787 matrix_component_factor(
const unsigned int component_i,
788 const unsigned int component_j,
789 const bool symmetric)
791 if (symmetric ==
false)
796 if (component_i < dim && component_j < dim)
798 else if (component_i >= dim && component_j >= dim)
808 template <
typename Number>
813 const unsigned int n_rows = out.size();
814 for (
unsigned int r = 0; r < n_rows; ++r)
820 template <
int dim,
typename Number>
824 return to_vector(s.operator
const Number &());
828 template <
int dim,
typename Number>
833 const unsigned int n_rows = out.size();
834 for (
unsigned int r = 0; r < n_rows; ++r)
836 const std::pair<unsigned int, unsigned int> indices =
837 internal::indices_from_component<dim>(r,
false);
839 const unsigned int i = indices.first;
846 template <
int dim,
typename Number>
851 const unsigned int n_rows = out.size();
852 for (
unsigned int r = 0; r < n_rows; ++r)
854 const std::pair<unsigned int, unsigned int> indices =
855 internal::indices_from_component<dim>(r,
false);
858 const unsigned int i = indices.first;
859 const unsigned int j = indices.second;
866 template <
int dim,
typename Number>
871 const unsigned int n_rows = out.size();
872 for (
unsigned int r = 0; r < n_rows; ++r)
874 const std::pair<unsigned int, unsigned int> indices =
875 internal::indices_from_component<dim>(r,
true);
879 const unsigned int i = indices.first;
880 const unsigned int j = indices.second;
882 const double factor =
883 internal::vector_component_factor<dim>(r,
true);
885 out(r) = factor * st[i][j];
891 template <
typename Number>
901 template <
int dim,
typename Number>
905 return to_matrix(s.operator
const Number &());
909 template <
int dim,
typename Number>
914 const unsigned int n_rows = out.m();
915 const unsigned int n_cols = out.n();
916 for (
unsigned int r = 0; r < n_rows; ++r)
918 const std::pair<unsigned int, unsigned int> indices =
919 internal::indices_from_component<dim>(r,
false);
921 const unsigned int i = indices.first;
923 for (
unsigned int c = 0; c < n_cols; ++c)
933 template <
int dim,
typename Number>
938 const unsigned int n_rows = out.m();
939 const unsigned int n_cols = out.n();
940 for (
unsigned int r = 0; r < n_rows; ++r)
942 const std::pair<unsigned int, unsigned int> indices_i =
943 internal::indices_from_component<dim>(r,
false);
946 const unsigned int i = indices_i.first;
948 for (
unsigned int c = 0; c < n_cols; ++c)
950 const std::pair<unsigned int, unsigned int> indices_j =
951 internal::indices_from_component<dim>(c,
false);
954 const unsigned int j = indices_j.second;
963 template <
int dim,
typename Number>
973 template <
typename TensorType>
974 struct is_rank_2_symmetric_tensor : std::false_type
977 template <
int dim,
typename Number>
992 (SubTensor1::dimension == dim && SubTensor2::dimension == dim),
993 "Sub-tensor spatial dimension is different from those of the input tensor.");
996 (SubTensor1::rank == 2 && SubTensor2::rank == 1) ||
997 (SubTensor1::rank == 1 && SubTensor2::rank == 2),
998 "Cannot build a rank 3 tensor from the given combination of sub-tensors.");
1001 SubTensor2::n_independent_components);
1002 const unsigned int n_rows = out.m();
1003 const unsigned int n_cols = out.n();
1005 if (SubTensor1::rank == 2 && SubTensor2::rank == 1)
1007 const bool subtensor_is_rank_2_symmetric_tensor =
1008 internal::is_rank_2_symmetric_tensor<SubTensor1>::value;
1010 for (
unsigned int r = 0; r < n_rows; ++r)
1012 const std::pair<unsigned int, unsigned int> indices_ij =
1013 internal::indices_from_component<dim>(
1014 r, subtensor_is_rank_2_symmetric_tensor);
1017 if (subtensor_is_rank_2_symmetric_tensor)
1019 Assert(indices_ij.second >= indices_ij.first,
1022 const unsigned int i = indices_ij.first;
1023 const unsigned int j = indices_ij.second;
1025 const double factor = internal::vector_component_factor<dim>(
1026 r, subtensor_is_rank_2_symmetric_tensor);
1028 for (
unsigned int c = 0; c < n_cols; ++c)
1030 const std::pair<unsigned int, unsigned int> indices_k =
1031 internal::indices_from_component<dim>(c,
false);
1033 const unsigned int k = indices_k.first;
1035 if (subtensor_is_rank_2_symmetric_tensor)
1036 out(r, c) = factor * t[i][j][k];
1038 out(r, c) = t[i][j][k];
1042 else if (SubTensor1::rank == 1 && SubTensor2::rank == 2)
1044 const bool subtensor_is_rank_2_symmetric_tensor =
1045 internal::is_rank_2_symmetric_tensor<SubTensor2>::value;
1047 for (
unsigned int r = 0; r < n_rows; ++r)
1049 const std::pair<unsigned int, unsigned int> indices_k =
1050 internal::indices_from_component<dim>(r,
false);
1052 const unsigned int k = indices_k.first;
1054 for (
unsigned int c = 0; c < n_cols; ++c)
1056 const std::pair<unsigned int, unsigned int> indices_ij =
1057 internal::indices_from_component<dim>(
1058 c, subtensor_is_rank_2_symmetric_tensor);
1061 if (subtensor_is_rank_2_symmetric_tensor)
1063 Assert(indices_ij.second >= indices_ij.first,
1066 const unsigned int i = indices_ij.first;
1067 const unsigned int j = indices_ij.second;
1069 if (subtensor_is_rank_2_symmetric_tensor)
1071 const double factor =
1072 internal::vector_component_factor<dim>(
1073 c, subtensor_is_rank_2_symmetric_tensor);
1074 out(r, c) = factor * t[k][i][j];
1077 out(r, c) = t[k][i][j];
1090 template <
int dim,
typename Number>
1097 const unsigned int n_rows = out.m();
1098 const unsigned int n_cols = out.n();
1099 for (
unsigned int r = 0; r < n_rows; ++r)
1101 const std::pair<unsigned int, unsigned int> indices_ij =
1102 internal::indices_from_component<dim>(r,
false);
1105 const unsigned int i = indices_ij.first;
1106 const unsigned int j = indices_ij.second;
1108 for (
unsigned int c = 0; c < n_cols; ++c)
1110 const std::pair<unsigned int, unsigned int> indices_kl =
1111 internal::indices_from_component<dim>(c,
false);
1114 const unsigned int k = indices_kl.first;
1115 const unsigned int l = indices_kl.second;
1117 out(r, c) = t[i][j][k][
l];
1124 template <
int dim,
typename Number>
1131 const unsigned int n_rows = out.m();
1132 const unsigned int n_cols = out.n();
1133 for (
unsigned int r = 0; r < n_rows; ++r)
1135 const std::pair<unsigned int, unsigned int> indices_ij =
1136 internal::indices_from_component<dim>(r,
true);
1140 const unsigned int i = indices_ij.first;
1141 const unsigned int j = indices_ij.second;
1143 for (
unsigned int c = 0; c < n_cols; ++c)
1145 const std::pair<unsigned int, unsigned int> indices_kl =
1146 internal::indices_from_component<dim>(c,
true);
1149 Assert(indices_kl.second >= indices_kl.first,
1151 const unsigned int k = indices_kl.first;
1152 const unsigned int l = indices_kl.second;
1154 const double factor =
1155 internal::matrix_component_factor<dim>(r, c,
true);
1157 out(r, c) = factor * st[i][j][k][
l];
1164 template <
typename Number>
1173 template <
int dim,
typename Number>
1177 return to_tensor(vec, s.operator Number &());
1181 template <
int dim,
typename Number>
1187 const unsigned int n_rows = vec.
size();
1188 for (
unsigned int r = 0; r < n_rows; ++r)
1190 const std::pair<unsigned int, unsigned int> indices =
1191 internal::indices_from_component<dim>(r,
false);
1193 const unsigned int i = indices.first;
1199 template <
int dim,
typename Number>
1205 const unsigned int n_rows = vec.
size();
1206 for (
unsigned int r = 0; r < n_rows; ++r)
1208 const std::pair<unsigned int, unsigned int> indices =
1209 internal::indices_from_component<dim>(r,
false);
1212 const unsigned int i = indices.first;
1213 const unsigned int j = indices.second;
1219 template <
int dim,
typename Number>
1225 const unsigned int n_rows = vec.
size();
1226 for (
unsigned int r = 0; r < n_rows; ++r)
1228 const std::pair<unsigned int, unsigned int> indices =
1229 internal::indices_from_component<dim>(r,
true);
1233 const unsigned int i = indices.first;
1234 const unsigned int j = indices.second;
1236 const double inv_factor =
1237 1.0 / internal::vector_component_factor<dim>(r,
true);
1239 st[i][j] = inv_factor * vec(r);
1244 template <
typename Number>
1250 Assert(mtrx.n_elements() == 1,
1256 template <
int dim,
typename Number>
1260 return to_tensor(mtrx, s.operator Number &());
1264 template <
int dim,
typename Number>
1274 const unsigned int n_rows = mtrx.
m();
1275 const unsigned int n_cols = mtrx.
n();
1276 for (
unsigned int r = 0; r < n_rows; ++r)
1278 const std::pair<unsigned int, unsigned int> indices =
1279 internal::indices_from_component<dim>(r,
false);
1282 const unsigned int i = indices.first;
1284 for (
unsigned int c = 0; c < n_cols; ++c)
1293 template <
int dim,
typename Number>
1303 const unsigned int n_rows = mtrx.
m();
1304 const unsigned int n_cols = mtrx.
n();
1305 for (
unsigned int r = 0; r < n_rows; ++r)
1307 const std::pair<unsigned int, unsigned int> indices_i =
1308 internal::indices_from_component<dim>(r,
false);
1311 const unsigned int i = indices_i.first;
1313 for (
unsigned int c = 0; c < n_cols; ++c)
1315 const std::pair<unsigned int, unsigned int> indices_j =
1316 internal::indices_from_component<dim>(c,
false);
1319 const unsigned int j = indices_j.second;
1321 t[i][j] = mtrx(r, c);
1327 template <
int dim,
typename Number>
1338 Assert((mtrx.n_elements() ==
1349 "The entries stored inside the matrix were not symmetric"));
1353 template <
int dim,
typename Number>
1378 const unsigned int n_rows = mtrx.
m();
1379 const unsigned int n_cols = mtrx.
n();
1391 const bool subtensor_is_rank_2_symmetric_tensor =
1395 for (
unsigned int r = 0; r < n_rows; ++r)
1397 const std::pair<unsigned int, unsigned int> indices_ij =
1398 internal::indices_from_component<dim>(
1399 r, subtensor_is_rank_2_symmetric_tensor);
1402 if (subtensor_is_rank_2_symmetric_tensor)
1404 Assert(indices_ij.second >= indices_ij.first,
1407 const unsigned int i = indices_ij.first;
1408 const unsigned int j = indices_ij.second;
1410 const double inv_factor =
1411 1.0 / internal::vector_component_factor<dim>(
1412 r, subtensor_is_rank_2_symmetric_tensor);
1414 for (
unsigned int c = 0; c < n_cols; ++c)
1416 const std::pair<unsigned int, unsigned int> indices_k =
1417 internal::indices_from_component<dim>(c,
false);
1419 const unsigned int k = indices_k.first;
1421 if (subtensor_is_rank_2_symmetric_tensor)
1423 t[i][j][k] = inv_factor * mtrx(r, c);
1424 t[j][i][k] = t[i][j][k];
1427 t[i][j][k] = mtrx(r, c);
1446 const bool subtensor_is_rank_2_symmetric_tensor =
1450 for (
unsigned int r = 0; r < n_rows; ++r)
1452 const std::pair<unsigned int, unsigned int> indices_k =
1453 internal::indices_from_component<dim>(r,
false);
1455 const unsigned int k = indices_k.first;
1457 for (
unsigned int c = 0; c < n_cols; ++c)
1459 const std::pair<unsigned int, unsigned int> indices_ij =
1460 internal::indices_from_component<dim>(
1461 c, subtensor_is_rank_2_symmetric_tensor);
1464 if (subtensor_is_rank_2_symmetric_tensor)
1466 Assert(indices_ij.second >= indices_ij.first,
1469 const unsigned int i = indices_ij.first;
1470 const unsigned int j = indices_ij.second;
1472 if (subtensor_is_rank_2_symmetric_tensor)
1474 const double inv_factor =
1475 1.0 / internal::vector_component_factor<dim>(
1476 c, subtensor_is_rank_2_symmetric_tensor);
1477 t[k][i][j] = inv_factor * mtrx(r, c);
1478 t[k][j][i] = t[k][i][j];
1481 t[k][i][j] = mtrx(r, c);
1488 template <
int dim,
typename Number>
1502 const unsigned int n_rows = mtrx.
m();
1503 const unsigned int n_cols = mtrx.
n();
1504 for (
unsigned int r = 0; r < n_rows; ++r)
1506 const std::pair<unsigned int, unsigned int> indices_ij =
1507 internal::indices_from_component<dim>(r,
false);
1510 const unsigned int i = indices_ij.first;
1511 const unsigned int j = indices_ij.second;
1513 for (
unsigned int c = 0; c < n_cols; ++c)
1515 const std::pair<unsigned int, unsigned int> indices_kl =
1516 internal::indices_from_component<dim>(c,
false);
1519 const unsigned int k = indices_kl.first;
1520 const unsigned int l = indices_kl.second;
1522 t[i][j][k][
l] = mtrx(r, c);
1528 template <
int dim,
typename Number>
1547 const unsigned int n_rows = mtrx.
m();
1548 const unsigned int n_cols = mtrx.
n();
1549 for (
unsigned int r = 0; r < n_rows; ++r)
1551 const std::pair<unsigned int, unsigned int> indices_ij =
1552 internal::indices_from_component<dim>(r,
false);
1555 const unsigned int i = indices_ij.first;
1556 const unsigned int j = indices_ij.second;
1558 for (
unsigned int c = 0; c < n_cols; ++c)
1560 const std::pair<unsigned int, unsigned int> indices_kl =
1561 internal::indices_from_component<dim>(c,
false);
1564 const unsigned int k = indices_kl.first;
1565 const unsigned int l = indices_kl.second;
1567 const double inv_factor =
1568 1.0 / internal::matrix_component_factor<dim>(r, c,
true);
1570 st[i][j][k][
l] = inv_factor * mtrx(r, c);
1576 template <
typename TensorType,
typename Number>
1586 template <
typename TensorType,
typename Number>
static constexpr unsigned int n_independent_components
static constexpr unsigned int n_independent_components
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize2(int arg1, int arg2, int arg3)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize3(int arg1, int arg2, int arg3, int arg4)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize3(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize2(int arg1, int arg2, int arg3)
#define AssertThrow(cond, exc)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void to_tensor(const Vector< Number > &vec, Number &s)
FullMatrix< Number > to_matrix(const Number &s)
Vector< Number > to_vector(const Number &s)
static constexpr double SQRT2
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)