15#ifndef dealii_physics_notation_h
16#define dealii_physics_notation_h
293 <<
"The number of rows in the input matrix is " << arg1
294 <<
", but needs to be either " << arg2 <<
" or " << arg3
306 <<
"The number of rows in the input matrix is " << arg1
307 <<
", but needs to be either " << arg2 <<
',' << arg3
308 <<
", or " << arg4 <<
'.');
318 <<
"The number of columns in the input matrix is " << arg1
319 <<
", but needs to be either " << arg2 <<
" or " << arg3
331 <<
"The number of columns in the input matrix is " << arg1
332 <<
", but needs to be either " << arg2 <<
',' << arg3
333 <<
", or " << arg4 <<
'.');
346 template <
typename Number>
356 template <
int dim,
typename Number>
366 template <
int dim,
typename Number>
376 template <
int dim,
typename Number>
387 template <
int dim,
typename Number>
397 template <
typename Number>
407 template <
int dim,
typename Number>
417 template <
int dim,
typename Number>
427 template <
int dim,
typename Number>
441 template <
int dim,
typename Number>
491 template <
int dim,
typename Number>
503 template <
int dim,
typename Number>
517 template <
typename Number>
525 template <
int dim,
typename Number>
533 template <
int dim,
typename Number>
541 template <
int dim,
typename Number>
549 template <
int dim,
typename Number>
557 template <
typename Number>
565 template <
int dim,
typename Number>
573 template <
int dim,
typename Number>
581 template <
int dim,
typename Number>
589 template <
int dim,
typename Number>
604 template <
int dim,
typename Number>
612 template <
int dim,
typename Number>
620 template <
int dim,
typename Number>
630 template <
typename TensorType,
typename Number>
639 template <
typename TensorType,
typename Number>
672 std::pair<unsigned int, unsigned int>
673 indices_from_component(
const unsigned int component_n,
674 const bool symmetric);
678 std::pair<unsigned int, unsigned int>
679 indices_from_component(
const unsigned int ,
const bool)
682 return std::make_pair(0u, 0u);
687 inline std::pair<unsigned int, unsigned int>
688 indices_from_component<1>(
const unsigned int component_n,
const bool)
692 return std::make_pair(0u, 0u);
697 inline std::pair<unsigned int, unsigned int>
698 indices_from_component<2>(
const unsigned int component_n,
699 const bool symmetric)
701 if (symmetric ==
true)
717 static const unsigned int indices[4][2] = {{0, 0},
721 return std::make_pair(indices[component_n][0],
722 indices[component_n][1]);
727 inline std::pair<unsigned int, unsigned int>
728 indices_from_component<3>(
const unsigned int component_n,
729 const bool symmetric)
731 if (symmetric ==
true)
747 static const unsigned int indices[9][2] = {{0, 0},
756 return std::make_pair(indices[component_n][0],
757 indices[component_n][1]);
767 vector_component_factor(
const unsigned int component_i,
768 const bool symmetric)
770 if (symmetric ==
false)
773 if (component_i < dim)
786 matrix_component_factor(
const unsigned int component_i,
787 const unsigned int component_j,
788 const bool symmetric)
790 if (symmetric ==
false)
795 if (component_i < dim && component_j < dim)
797 else if (component_i >= dim && component_j >= dim)
807 template <
typename Number>
812 const unsigned int n_rows =
out.size();
813 for (
unsigned int r = 0; r < n_rows; ++r)
819 template <
int dim,
typename Number>
823 return to_vector(s.operator
const Number &());
827 template <
int dim,
typename Number>
832 const unsigned int n_rows =
out.size();
833 for (
unsigned int r = 0; r < n_rows; ++r)
835 const std::pair<unsigned int, unsigned int> indices =
836 internal::indices_from_component<dim>(r,
false);
838 const unsigned int i = indices.first;
845 template <
int dim,
typename Number>
850 const unsigned int n_rows =
out.size();
851 for (
unsigned int r = 0; r < n_rows; ++r)
853 const std::pair<unsigned int, unsigned int> indices =
854 internal::indices_from_component<dim>(r,
false);
857 const unsigned int i = indices.first;
858 const unsigned int j = indices.second;
865 template <
int dim,
typename Number>
870 const unsigned int n_rows =
out.size();
871 for (
unsigned int r = 0; r < n_rows; ++r)
873 const std::pair<unsigned int, unsigned int> indices =
874 internal::indices_from_component<dim>(r,
true);
878 const unsigned int i = indices.first;
879 const unsigned int j = indices.second;
881 const double factor =
882 internal::vector_component_factor<dim>(r,
true);
884 out(r) = factor * st[i][j];
890 template <
typename Number>
900 template <
int dim,
typename Number>
904 return to_matrix(s.operator
const Number &());
908 template <
int dim,
typename Number>
913 const unsigned int n_rows =
out.m();
914 const unsigned int n_cols =
out.n();
915 for (
unsigned int r = 0; r < n_rows; ++r)
917 const std::pair<unsigned int, unsigned int> indices =
918 internal::indices_from_component<dim>(r,
false);
920 const unsigned int i = indices.first;
922 for (
unsigned int c = 0; c < n_cols; ++c)
932 template <
int dim,
typename Number>
937 const unsigned int n_rows =
out.m();
938 const unsigned int n_cols =
out.n();
939 for (
unsigned int r = 0; r < n_rows; ++r)
941 const std::pair<unsigned int, unsigned int> indices_i =
942 internal::indices_from_component<dim>(r,
false);
945 const unsigned int i = indices_i.first;
947 for (
unsigned int c = 0; c < n_cols; ++c)
949 const std::pair<unsigned int, unsigned int> indices_j =
950 internal::indices_from_component<dim>(c,
false);
953 const unsigned int j = indices_j.second;
962 template <
int dim,
typename Number>
972 template <
typename TensorType>
973 struct is_rank_2_symmetric_tensor : std::false_type
976 template <
int dim,
typename Number>
991 (SubTensor1::dimension == dim && SubTensor2::dimension == dim),
992 "Sub-tensor spatial dimension is different from those of the input tensor.");
995 (SubTensor1::rank == 2 && SubTensor2::rank == 1) ||
996 (SubTensor1::rank == 1 && SubTensor2::rank == 2),
997 "Cannot build a rank 3 tensor from the given combination of sub-tensors.");
1000 SubTensor2::n_independent_components);
1001 const unsigned int n_rows =
out.m();
1002 const unsigned int n_cols =
out.n();
1004 if (SubTensor1::rank == 2 && SubTensor2::rank == 1)
1006 const bool subtensor_is_rank_2_symmetric_tensor =
1007 internal::is_rank_2_symmetric_tensor<SubTensor1>::value;
1009 for (
unsigned int r = 0; r < n_rows; ++r)
1011 const std::pair<unsigned int, unsigned int> indices_ij =
1012 internal::indices_from_component<dim>(
1013 r, subtensor_is_rank_2_symmetric_tensor);
1016 if (subtensor_is_rank_2_symmetric_tensor)
1018 Assert(indices_ij.second >= indices_ij.first,
1021 const unsigned int i = indices_ij.first;
1022 const unsigned int j = indices_ij.second;
1024 const double factor = internal::vector_component_factor<dim>(
1025 r, subtensor_is_rank_2_symmetric_tensor);
1027 for (
unsigned int c = 0; c < n_cols; ++c)
1029 const std::pair<unsigned int, unsigned int> indices_k =
1030 internal::indices_from_component<dim>(c,
false);
1032 const unsigned int k = indices_k.first;
1034 if (subtensor_is_rank_2_symmetric_tensor)
1035 out(r, c) = factor * t[i][j][k];
1037 out(r, c) = t[i][j][k];
1041 else if (SubTensor1::rank == 1 && SubTensor2::rank == 2)
1043 const bool subtensor_is_rank_2_symmetric_tensor =
1044 internal::is_rank_2_symmetric_tensor<SubTensor2>::value;
1046 for (
unsigned int r = 0; r < n_rows; ++r)
1048 const std::pair<unsigned int, unsigned int> indices_k =
1049 internal::indices_from_component<dim>(r,
false);
1051 const unsigned int k = indices_k.first;
1053 for (
unsigned int c = 0; c < n_cols; ++c)
1055 const std::pair<unsigned int, unsigned int> indices_ij =
1056 internal::indices_from_component<dim>(
1057 c, subtensor_is_rank_2_symmetric_tensor);
1060 if (subtensor_is_rank_2_symmetric_tensor)
1062 Assert(indices_ij.second >= indices_ij.first,
1065 const unsigned int i = indices_ij.first;
1066 const unsigned int j = indices_ij.second;
1068 if (subtensor_is_rank_2_symmetric_tensor)
1070 const double factor =
1071 internal::vector_component_factor<dim>(
1072 c, subtensor_is_rank_2_symmetric_tensor);
1073 out(r, c) = factor * t[k][i][j];
1076 out(r, c) = t[k][i][j];
1089 template <
int dim,
typename Number>
1096 const unsigned int n_rows =
out.m();
1097 const unsigned int n_cols =
out.n();
1098 for (
unsigned int r = 0; r < n_rows; ++r)
1100 const std::pair<unsigned int, unsigned int> indices_ij =
1101 internal::indices_from_component<dim>(r,
false);
1104 const unsigned int i = indices_ij.first;
1105 const unsigned int j = indices_ij.second;
1107 for (
unsigned int c = 0; c < n_cols; ++c)
1109 const std::pair<unsigned int, unsigned int> indices_kl =
1110 internal::indices_from_component<dim>(c,
false);
1113 const unsigned int k = indices_kl.first;
1114 const unsigned int l = indices_kl.second;
1116 out(r, c) = t[i][j][k][
l];
1123 template <
int dim,
typename Number>
1130 const unsigned int n_rows =
out.m();
1131 const unsigned int n_cols =
out.n();
1132 for (
unsigned int r = 0; r < n_rows; ++r)
1134 const std::pair<unsigned int, unsigned int> indices_ij =
1135 internal::indices_from_component<dim>(r,
true);
1139 const unsigned int i = indices_ij.first;
1140 const unsigned int j = indices_ij.second;
1142 for (
unsigned int c = 0; c < n_cols; ++c)
1144 const std::pair<unsigned int, unsigned int> indices_kl =
1145 internal::indices_from_component<dim>(c,
true);
1148 Assert(indices_kl.second >= indices_kl.first,
1150 const unsigned int k = indices_kl.first;
1151 const unsigned int l = indices_kl.second;
1153 const double factor =
1154 internal::matrix_component_factor<dim>(r, c,
true);
1156 out(r, c) = factor * st[i][j][k][
l];
1163 template <
typename Number>
1172 template <
int dim,
typename Number>
1176 return to_tensor(vec, s.operator Number &());
1180 template <
int dim,
typename Number>
1186 const unsigned int n_rows = vec.
size();
1187 for (
unsigned int r = 0; r < n_rows; ++r)
1189 const std::pair<unsigned int, unsigned int> indices =
1190 internal::indices_from_component<dim>(r,
false);
1192 const unsigned int i = indices.first;
1198 template <
int dim,
typename Number>
1204 const unsigned int n_rows = vec.
size();
1205 for (
unsigned int r = 0; r < n_rows; ++r)
1207 const std::pair<unsigned int, unsigned int> indices =
1208 internal::indices_from_component<dim>(r,
false);
1211 const unsigned int i = indices.first;
1212 const unsigned int j = indices.second;
1218 template <
int dim,
typename Number>
1224 const unsigned int n_rows = vec.
size();
1225 for (
unsigned int r = 0; r < n_rows; ++r)
1227 const std::pair<unsigned int, unsigned int> indices =
1228 internal::indices_from_component<dim>(r,
true);
1232 const unsigned int i = indices.first;
1233 const unsigned int j = indices.second;
1235 const double inv_factor =
1236 1.0 / internal::vector_component_factor<dim>(r,
true);
1238 st[i][j] = inv_factor * vec(r);
1243 template <
typename Number>
1249 Assert(mtrx.n_elements() == 1,
1255 template <
int dim,
typename Number>
1259 return to_tensor(mtrx, s.operator Number &());
1263 template <
int dim,
typename Number>
1273 const unsigned int n_rows = mtrx.
m();
1274 const unsigned int n_cols = mtrx.
n();
1275 for (
unsigned int r = 0; r < n_rows; ++r)
1277 const std::pair<unsigned int, unsigned int> indices =
1278 internal::indices_from_component<dim>(r,
false);
1281 const unsigned int i = indices.first;
1283 for (
unsigned int c = 0; c < n_cols; ++c)
1292 template <
int dim,
typename Number>
1302 const unsigned int n_rows = mtrx.
m();
1303 const unsigned int n_cols = mtrx.
n();
1304 for (
unsigned int r = 0; r < n_rows; ++r)
1306 const std::pair<unsigned int, unsigned int> indices_i =
1307 internal::indices_from_component<dim>(r,
false);
1310 const unsigned int i = indices_i.first;
1312 for (
unsigned int c = 0; c < n_cols; ++c)
1314 const std::pair<unsigned int, unsigned int> indices_j =
1315 internal::indices_from_component<dim>(c,
false);
1318 const unsigned int j = indices_j.second;
1320 t[i][j] = mtrx(r, c);
1326 template <
int dim,
typename Number>
1337 Assert((mtrx.n_elements() ==
1348 "The entries stored inside the matrix were not symmetric"));
1352 template <
int dim,
typename Number>
1377 const unsigned int n_rows = mtrx.
m();
1378 const unsigned int n_cols = mtrx.
n();
1390 const bool subtensor_is_rank_2_symmetric_tensor =
1394 for (
unsigned int r = 0; r < n_rows; ++r)
1396 const std::pair<unsigned int, unsigned int> indices_ij =
1397 internal::indices_from_component<dim>(
1398 r, subtensor_is_rank_2_symmetric_tensor);
1401 if (subtensor_is_rank_2_symmetric_tensor)
1403 Assert(indices_ij.second >= indices_ij.first,
1406 const unsigned int i = indices_ij.first;
1407 const unsigned int j = indices_ij.second;
1409 const double inv_factor =
1410 1.0 / internal::vector_component_factor<dim>(
1411 r, subtensor_is_rank_2_symmetric_tensor);
1413 for (
unsigned int c = 0; c < n_cols; ++c)
1415 const std::pair<unsigned int, unsigned int> indices_k =
1416 internal::indices_from_component<dim>(c,
false);
1418 const unsigned int k = indices_k.first;
1420 if (subtensor_is_rank_2_symmetric_tensor)
1422 t[i][j][k] = inv_factor * mtrx(r, c);
1423 t[j][i][k] = t[i][j][k];
1426 t[i][j][k] = mtrx(r, c);
1445 const bool subtensor_is_rank_2_symmetric_tensor =
1449 for (
unsigned int r = 0; r < n_rows; ++r)
1451 const std::pair<unsigned int, unsigned int> indices_k =
1452 internal::indices_from_component<dim>(r,
false);
1454 const unsigned int k = indices_k.first;
1456 for (
unsigned int c = 0; c < n_cols; ++c)
1458 const std::pair<unsigned int, unsigned int> indices_ij =
1459 internal::indices_from_component<dim>(
1460 c, subtensor_is_rank_2_symmetric_tensor);
1463 if (subtensor_is_rank_2_symmetric_tensor)
1465 Assert(indices_ij.second >= indices_ij.first,
1468 const unsigned int i = indices_ij.first;
1469 const unsigned int j = indices_ij.second;
1471 if (subtensor_is_rank_2_symmetric_tensor)
1473 const double inv_factor =
1474 1.0 / internal::vector_component_factor<dim>(
1475 c, subtensor_is_rank_2_symmetric_tensor);
1476 t[k][i][j] = inv_factor * mtrx(r, c);
1477 t[k][j][i] = t[k][i][j];
1480 t[k][i][j] = mtrx(r, c);
1487 template <
int dim,
typename Number>
1501 const unsigned int n_rows = mtrx.
m();
1502 const unsigned int n_cols = mtrx.
n();
1503 for (
unsigned int r = 0; r < n_rows; ++r)
1505 const std::pair<unsigned int, unsigned int> indices_ij =
1506 internal::indices_from_component<dim>(r,
false);
1509 const unsigned int i = indices_ij.first;
1510 const unsigned int j = indices_ij.second;
1512 for (
unsigned int c = 0; c < n_cols; ++c)
1514 const std::pair<unsigned int, unsigned int> indices_kl =
1515 internal::indices_from_component<dim>(c,
false);
1518 const unsigned int k = indices_kl.first;
1519 const unsigned int l = indices_kl.second;
1521 t[i][j][k][
l] = mtrx(r, c);
1527 template <
int dim,
typename Number>
1546 const unsigned int n_rows = mtrx.
m();
1547 const unsigned int n_cols = mtrx.
n();
1548 for (
unsigned int r = 0; r < n_rows; ++r)
1550 const std::pair<unsigned int, unsigned int> indices_ij =
1551 internal::indices_from_component<dim>(r,
false);
1554 const unsigned int i = indices_ij.first;
1555 const unsigned int j = indices_ij.second;
1557 for (
unsigned int c = 0; c < n_cols; ++c)
1559 const std::pair<unsigned int, unsigned int> indices_kl =
1560 internal::indices_from_component<dim>(c,
false);
1563 const unsigned int k = indices_kl.first;
1564 const unsigned int l = indices_kl.second;
1566 const double inv_factor =
1567 1.0 / internal::matrix_component_factor<dim>(r, c,
true);
1569 st[i][j][k][
l] = inv_factor * mtrx(r, c);
1575 template <
typename TensorType,
typename Number>
1585 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
const InputIterator OutputIterator out
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)