15#ifndef dealii_physics_notation_h
16#define dealii_physics_notation_h
267 <<
"The number of rows in the input matrix is " << arg1
268 <<
", but needs to be either " << arg2 <<
" or " << arg3
280 <<
"The number of rows in the input matrix is " << arg1
281 <<
", but needs to be either " << arg2 <<
',' << arg3
282 <<
", or " << arg4 <<
'.');
292 <<
"The number of columns in the input matrix is " << arg1
293 <<
", but needs to be either " << arg2 <<
" or " << arg3
305 <<
"The number of columns in the input matrix is " << arg1
306 <<
", but needs to be either " << arg2 <<
',' << arg3
307 <<
", or " << arg4 <<
'.');
320 template <
typename Number>
330 template <
int dim,
typename Number>
340 template <
int dim,
typename Number>
350 template <
int dim,
typename Number>
361 template <
int dim,
typename Number>
371 template <
typename Number>
381 template <
int dim,
typename Number>
391 template <
int dim,
typename Number>
401 template <
int dim,
typename Number>
415 template <
int dim,
typename Number>
465 template <
int dim,
typename Number>
477 template <
int dim,
typename Number>
491 template <
typename Number>
499 template <
int dim,
typename Number>
507 template <
int dim,
typename Number>
515 template <
int dim,
typename Number>
523 template <
int dim,
typename Number>
531 template <
typename Number>
539 template <
int dim,
typename Number>
547 template <
int dim,
typename Number>
555 template <
int dim,
typename Number>
563 template <
int dim,
typename Number>
578 template <
int dim,
typename Number>
586 template <
int dim,
typename Number>
594 template <
int dim,
typename Number>
604 template <
typename TensorType,
typename Number>
613 template <
typename TensorType,
typename Number>
646 std::pair<unsigned int, unsigned int>
647 indices_from_component(
const unsigned int component_n,
648 const bool symmetric);
652 std::pair<unsigned int, unsigned int>
653 indices_from_component(
const unsigned int ,
const bool)
656 return std::make_pair(0u, 0u);
661 inline std::pair<unsigned int, unsigned int>
662 indices_from_component<1>(
const unsigned int component_n,
const bool)
666 return std::make_pair(0u, 0u);
671 inline std::pair<unsigned int, unsigned int>
672 indices_from_component<2>(
const unsigned int component_n,
673 const bool symmetric)
675 if (symmetric ==
true)
691 static const unsigned int indices[4][2] = {{0, 0},
695 return std::make_pair(indices[component_n][0],
696 indices[component_n][1]);
701 inline std::pair<unsigned int, unsigned int>
702 indices_from_component<3>(
const unsigned int component_n,
703 const bool symmetric)
705 if (symmetric ==
true)
721 static const unsigned int indices[9][2] = {{0, 0},
730 return std::make_pair(indices[component_n][0],
731 indices[component_n][1]);
741 vector_component_factor(
const unsigned int component_i,
742 const bool symmetric)
744 if (symmetric ==
false)
747 if (component_i < dim)
760 matrix_component_factor(
const unsigned int component_i,
761 const unsigned int component_j,
762 const bool symmetric)
764 if (symmetric ==
false)
769 if (component_i < dim && component_j < dim)
771 else if (component_i >= dim && component_j >= dim)
781 template <
typename Number>
786 const unsigned int n_rows = out.size();
787 for (
unsigned int r = 0; r < n_rows; ++r)
793 template <
int dim,
typename Number>
797 return to_vector(s.operator
const Number &());
801 template <
int dim,
typename Number>
806 const unsigned int n_rows = out.size();
807 for (
unsigned int r = 0; r < n_rows; ++r)
809 const std::pair<unsigned int, unsigned int> indices =
810 internal::indices_from_component<dim>(r,
false);
812 const unsigned int i = indices.first;
819 template <
int dim,
typename Number>
824 const unsigned int n_rows = out.size();
825 for (
unsigned int r = 0; r < n_rows; ++r)
827 const std::pair<unsigned int, unsigned int> indices =
828 internal::indices_from_component<dim>(r,
false);
831 const unsigned int i = indices.first;
832 const unsigned int j = indices.second;
839 template <
int dim,
typename Number>
844 const unsigned int n_rows = out.size();
845 for (
unsigned int r = 0; r < n_rows; ++r)
847 const std::pair<unsigned int, unsigned int> indices =
848 internal::indices_from_component<dim>(r,
true);
852 const unsigned int i = indices.first;
853 const unsigned int j = indices.second;
855 const double factor =
856 internal::vector_component_factor<dim>(r,
true);
858 out(r) = factor * st[i][j];
864 template <
typename Number>
874 template <
int dim,
typename Number>
878 return to_matrix(s.operator
const Number &());
882 template <
int dim,
typename Number>
887 const unsigned int n_rows = out.m();
888 const unsigned int n_cols = out.n();
889 for (
unsigned int r = 0; r < n_rows; ++r)
891 const std::pair<unsigned int, unsigned int> indices =
892 internal::indices_from_component<dim>(r,
false);
894 const unsigned int i = indices.first;
896 for (
unsigned int c = 0; c < n_cols; ++c)
906 template <
int dim,
typename Number>
911 const unsigned int n_rows = out.m();
912 const unsigned int n_cols = out.n();
913 for (
unsigned int r = 0; r < n_rows; ++r)
915 const std::pair<unsigned int, unsigned int> indices_i =
916 internal::indices_from_component<dim>(r,
false);
919 const unsigned int i = indices_i.first;
921 for (
unsigned int c = 0; c < n_cols; ++c)
923 const std::pair<unsigned int, unsigned int> indices_j =
924 internal::indices_from_component<dim>(c,
false);
927 const unsigned int j = indices_j.second;
936 template <
int dim,
typename Number>
946 template <
typename TensorType>
947 struct is_rank_2_symmetric_tensor : std::false_type
950 template <
int dim,
typename Number>
965 (SubTensor1::dimension == dim && SubTensor2::dimension == dim),
966 "Sub-tensor spatial dimension is different from those of the input tensor.");
969 (SubTensor1::rank == 2 && SubTensor2::rank == 1) ||
970 (SubTensor1::rank == 1 && SubTensor2::rank == 2),
971 "Cannot build a rank 3 tensor from the given combination of sub-tensors.");
974 SubTensor2::n_independent_components);
975 const unsigned int n_rows = out.m();
976 const unsigned int n_cols = out.n();
978 if (SubTensor1::rank == 2 && SubTensor2::rank == 1)
980 const bool subtensor_is_rank_2_symmetric_tensor =
981 internal::is_rank_2_symmetric_tensor<SubTensor1>::value;
983 for (
unsigned int r = 0; r < n_rows; ++r)
985 const std::pair<unsigned int, unsigned int> indices_ij =
986 internal::indices_from_component<dim>(
987 r, subtensor_is_rank_2_symmetric_tensor);
990 if (subtensor_is_rank_2_symmetric_tensor)
992 Assert(indices_ij.second >= indices_ij.first,
995 const unsigned int i = indices_ij.first;
996 const unsigned int j = indices_ij.second;
998 const double factor = internal::vector_component_factor<dim>(
999 r, subtensor_is_rank_2_symmetric_tensor);
1001 for (
unsigned int c = 0; c < n_cols; ++c)
1003 const std::pair<unsigned int, unsigned int> indices_k =
1004 internal::indices_from_component<dim>(c,
false);
1006 const unsigned int k = indices_k.first;
1008 if (subtensor_is_rank_2_symmetric_tensor)
1009 out(r, c) = factor * t[i][j][k];
1011 out(r, c) = t[i][j][k];
1015 else if (SubTensor1::rank == 1 && SubTensor2::rank == 2)
1017 const bool subtensor_is_rank_2_symmetric_tensor =
1018 internal::is_rank_2_symmetric_tensor<SubTensor2>::value;
1020 for (
unsigned int r = 0; r < n_rows; ++r)
1022 const std::pair<unsigned int, unsigned int> indices_k =
1023 internal::indices_from_component<dim>(r,
false);
1025 const unsigned int k = indices_k.first;
1027 for (
unsigned int c = 0; c < n_cols; ++c)
1029 const std::pair<unsigned int, unsigned int> indices_ij =
1030 internal::indices_from_component<dim>(
1031 c, subtensor_is_rank_2_symmetric_tensor);
1034 if (subtensor_is_rank_2_symmetric_tensor)
1036 Assert(indices_ij.second >= indices_ij.first,
1039 const unsigned int i = indices_ij.first;
1040 const unsigned int j = indices_ij.second;
1042 if (subtensor_is_rank_2_symmetric_tensor)
1044 const double factor =
1045 internal::vector_component_factor<dim>(
1046 c, subtensor_is_rank_2_symmetric_tensor);
1047 out(r, c) = factor * t[k][i][j];
1050 out(r, c) = t[k][i][j];
1063 template <
int dim,
typename Number>
1070 const unsigned int n_rows = out.m();
1071 const unsigned int n_cols = out.n();
1072 for (
unsigned int r = 0; r < n_rows; ++r)
1074 const std::pair<unsigned int, unsigned int> indices_ij =
1075 internal::indices_from_component<dim>(r,
false);
1078 const unsigned int i = indices_ij.first;
1079 const unsigned int j = indices_ij.second;
1081 for (
unsigned int c = 0; c < n_cols; ++c)
1083 const std::pair<unsigned int, unsigned int> indices_kl =
1084 internal::indices_from_component<dim>(c,
false);
1087 const unsigned int k = indices_kl.first;
1088 const unsigned int l = indices_kl.second;
1090 out(r, c) = t[i][j][k][
l];
1097 template <
int dim,
typename Number>
1104 const unsigned int n_rows = out.m();
1105 const unsigned int n_cols = out.n();
1106 for (
unsigned int r = 0; r < n_rows; ++r)
1108 const std::pair<unsigned int, unsigned int> indices_ij =
1109 internal::indices_from_component<dim>(r,
true);
1113 const unsigned int i = indices_ij.first;
1114 const unsigned int j = indices_ij.second;
1116 for (
unsigned int c = 0; c < n_cols; ++c)
1118 const std::pair<unsigned int, unsigned int> indices_kl =
1119 internal::indices_from_component<dim>(c,
true);
1122 Assert(indices_kl.second >= indices_kl.first,
1124 const unsigned int k = indices_kl.first;
1125 const unsigned int l = indices_kl.second;
1127 const double factor =
1128 internal::matrix_component_factor<dim>(r, c,
true);
1130 out(r, c) = factor * st[i][j][k][
l];
1137 template <
typename Number>
1146 template <
int dim,
typename Number>
1150 return to_tensor(vec, s.operator Number &());
1154 template <
int dim,
typename Number>
1160 const unsigned int n_rows = vec.
size();
1161 for (
unsigned int r = 0; r < n_rows; ++r)
1163 const std::pair<unsigned int, unsigned int> indices =
1164 internal::indices_from_component<dim>(r,
false);
1166 const unsigned int i = indices.first;
1172 template <
int dim,
typename Number>
1178 const unsigned int n_rows = vec.
size();
1179 for (
unsigned int r = 0; r < n_rows; ++r)
1181 const std::pair<unsigned int, unsigned int> indices =
1182 internal::indices_from_component<dim>(r,
false);
1185 const unsigned int i = indices.first;
1186 const unsigned int j = indices.second;
1192 template <
int dim,
typename Number>
1198 const unsigned int n_rows = vec.
size();
1199 for (
unsigned int r = 0; r < n_rows; ++r)
1201 const std::pair<unsigned int, unsigned int> indices =
1202 internal::indices_from_component<dim>(r,
true);
1206 const unsigned int i = indices.first;
1207 const unsigned int j = indices.second;
1209 const double inv_factor =
1210 1.0 / internal::vector_component_factor<dim>(r,
true);
1212 st[i][j] = inv_factor * vec(r);
1217 template <
typename Number>
1223 Assert(mtrx.n_elements() == 1,
1229 template <
int dim,
typename Number>
1233 return to_tensor(mtrx, s.operator Number &());
1237 template <
int dim,
typename Number>
1247 const unsigned int n_rows = mtrx.
m();
1248 const unsigned int n_cols = mtrx.
n();
1249 for (
unsigned int r = 0; r < n_rows; ++r)
1251 const std::pair<unsigned int, unsigned int> indices =
1252 internal::indices_from_component<dim>(r,
false);
1255 const unsigned int i = indices.first;
1257 for (
unsigned int c = 0; c < n_cols; ++c)
1266 template <
int dim,
typename Number>
1276 const unsigned int n_rows = mtrx.
m();
1277 const unsigned int n_cols = mtrx.
n();
1278 for (
unsigned int r = 0; r < n_rows; ++r)
1280 const std::pair<unsigned int, unsigned int> indices_i =
1281 internal::indices_from_component<dim>(r,
false);
1284 const unsigned int i = indices_i.first;
1286 for (
unsigned int c = 0; c < n_cols; ++c)
1288 const std::pair<unsigned int, unsigned int> indices_j =
1289 internal::indices_from_component<dim>(c,
false);
1292 const unsigned int j = indices_j.second;
1294 t[i][j] = mtrx(r, c);
1300 template <
int dim,
typename Number>
1311 Assert((mtrx.n_elements() ==
1322 "The entries stored inside the matrix were not symmetric"));
1326 template <
int dim,
typename Number>
1351 const unsigned int n_rows = mtrx.
m();
1352 const unsigned int n_cols = mtrx.
n();
1364 const bool subtensor_is_rank_2_symmetric_tensor =
1368 for (
unsigned int r = 0; r < n_rows; ++r)
1370 const std::pair<unsigned int, unsigned int> indices_ij =
1371 internal::indices_from_component<dim>(
1372 r, subtensor_is_rank_2_symmetric_tensor);
1375 if (subtensor_is_rank_2_symmetric_tensor)
1377 Assert(indices_ij.second >= indices_ij.first,
1380 const unsigned int i = indices_ij.first;
1381 const unsigned int j = indices_ij.second;
1383 const double inv_factor =
1384 1.0 / internal::vector_component_factor<dim>(
1385 r, subtensor_is_rank_2_symmetric_tensor);
1387 for (
unsigned int c = 0; c < n_cols; ++c)
1389 const std::pair<unsigned int, unsigned int> indices_k =
1390 internal::indices_from_component<dim>(c,
false);
1392 const unsigned int k = indices_k.first;
1394 if (subtensor_is_rank_2_symmetric_tensor)
1396 t[i][j][k] = inv_factor * mtrx(r, c);
1397 t[j][i][k] = t[i][j][k];
1400 t[i][j][k] = mtrx(r, c);
1419 const bool subtensor_is_rank_2_symmetric_tensor =
1423 for (
unsigned int r = 0; r < n_rows; ++r)
1425 const std::pair<unsigned int, unsigned int> indices_k =
1426 internal::indices_from_component<dim>(r,
false);
1428 const unsigned int k = indices_k.first;
1430 for (
unsigned int c = 0; c < n_cols; ++c)
1432 const std::pair<unsigned int, unsigned int> indices_ij =
1433 internal::indices_from_component<dim>(
1434 c, subtensor_is_rank_2_symmetric_tensor);
1437 if (subtensor_is_rank_2_symmetric_tensor)
1439 Assert(indices_ij.second >= indices_ij.first,
1442 const unsigned int i = indices_ij.first;
1443 const unsigned int j = indices_ij.second;
1445 if (subtensor_is_rank_2_symmetric_tensor)
1447 const double inv_factor =
1448 1.0 / internal::vector_component_factor<dim>(
1449 c, subtensor_is_rank_2_symmetric_tensor);
1450 t[k][i][j] = inv_factor * mtrx(r, c);
1451 t[k][j][i] = t[k][i][j];
1454 t[k][i][j] = mtrx(r, c);
1461 template <
int dim,
typename Number>
1475 const unsigned int n_rows = mtrx.
m();
1476 const unsigned int n_cols = mtrx.
n();
1477 for (
unsigned int r = 0; r < n_rows; ++r)
1479 const std::pair<unsigned int, unsigned int> indices_ij =
1480 internal::indices_from_component<dim>(r,
false);
1483 const unsigned int i = indices_ij.first;
1484 const unsigned int j = indices_ij.second;
1486 for (
unsigned int c = 0; c < n_cols; ++c)
1488 const std::pair<unsigned int, unsigned int> indices_kl =
1489 internal::indices_from_component<dim>(c,
false);
1492 const unsigned int k = indices_kl.first;
1493 const unsigned int l = indices_kl.second;
1495 t[i][j][k][
l] = mtrx(r, c);
1501 template <
int dim,
typename Number>
1520 const unsigned int n_rows = mtrx.
m();
1521 const unsigned int n_cols = mtrx.
n();
1522 for (
unsigned int r = 0; r < n_rows; ++r)
1524 const std::pair<unsigned int, unsigned int> indices_ij =
1525 internal::indices_from_component<dim>(r,
false);
1528 const unsigned int i = indices_ij.first;
1529 const unsigned int j = indices_ij.second;
1531 for (
unsigned int c = 0; c < n_cols; ++c)
1533 const std::pair<unsigned int, unsigned int> indices_kl =
1534 internal::indices_from_component<dim>(c,
false);
1537 const unsigned int k = indices_kl.first;
1538 const unsigned int l = indices_kl.second;
1540 const double inv_factor =
1541 1.0 / internal::matrix_component_factor<dim>(r, c,
true);
1543 st[i][j][k][
l] = inv_factor * mtrx(r, c);
1549 template <
typename TensorType,
typename Number>
1559 template <
typename TensorType,
typename Number>
static constexpr unsigned int n_independent_components
static constexpr unsigned int n_independent_components
virtual size_type size() const override
#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
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)