16 #ifndef dealii_physics_notation_h 17 #define dealii_physics_notation_h 19 #include <deal.II/base/exceptions.h> 20 #include <deal.II/base/numbers.h> 21 #include <deal.II/base/tensor.h> 23 #include <deal.II/lac/full_matrix.h> 24 #include <deal.II/lac/vector.h> 26 #include <type_traits> 28 DEAL_II_NAMESPACE_OPEN
267 <<
"The number of rows in the input matrix is " << arg1
268 <<
", but needs to be either " << arg2
269 <<
" or " << arg3 <<
".");
276 <<
"The number of rows in the input matrix is " << arg1
277 <<
", but needs to be either " << arg2 <<
"," << arg3
278 <<
", or " << arg4 <<
".");
285 <<
"The number of columns in the input matrix is " << arg1
286 <<
", but needs to be either " << arg2
287 <<
" or " << arg3 <<
".");
294 <<
"The number of columns in the input matrix is " << arg1
295 <<
", but needs to be either " << arg2 <<
"," << arg3
296 <<
", or " << arg4 <<
".");
309 template<
typename Number>
319 template<
int dim,
typename Number>
329 template<
int dim,
typename Number>
339 template<
int dim,
typename Number>
350 template<
int dim,
typename Number>
360 template<
typename Number>
370 template<
int dim,
typename Number>
380 template<
int dim,
typename Number>
390 template<
int dim,
typename Number>
404 template<
int dim,
typename Number>
448 template<
int dim,
typename Number>
459 template<
int dim,
typename Number>
473 template<
typename Number>
482 template<
int dim,
typename Number>
491 template<
int dim,
typename Number>
500 template<
int dim,
typename Number>
509 template<
int dim,
typename Number>
518 template<
typename Number>
527 template<
int dim,
typename Number>
536 template<
int dim,
typename Number>
545 template<
int dim,
typename Number>
554 template<
int dim,
typename Number>
569 template<
int dim,
typename Number>
578 template<
int dim,
typename Number>
587 template<
int dim,
typename Number>
597 template<
typename TensorType,
typename Number>
606 template<
typename TensorType,
typename Number>
641 std::pair<unsigned int, unsigned int>
642 indices_from_component (
const unsigned int component_n,
643 const bool symmetric);
647 std::pair<unsigned int, unsigned int>
648 indices_from_component (
const unsigned int component_n,
652 return std::pair<unsigned int, unsigned int> ();
658 std::pair<unsigned int, unsigned int>
659 indices_from_component<1> (
const unsigned int component_n,
664 return std::make_pair(0u,0u);
670 std::pair<unsigned int, unsigned int>
671 indices_from_component<2> (
const unsigned int component_n,
674 if (symmetric ==
true)
685 static const unsigned int indices[4][2] =
690 return std::make_pair(indices[component_n][0], indices[component_n][1]);
696 std::pair<unsigned int, unsigned int>
697 indices_from_component<3> (
const unsigned int component_n,
700 if (symmetric ==
true)
711 static const unsigned int indices[9][2] =
717 return std::make_pair(indices[component_n][0], indices[component_n][1]);
727 vector_component_factor (
const unsigned int component_i,
728 const bool symmetric)
730 if (symmetric ==
false)
733 if (component_i < dim)
746 matrix_component_factor (
const unsigned int component_i,
747 const unsigned int component_j,
748 const bool symmetric)
750 if (symmetric ==
false)
755 if (component_i < dim && component_j < dim)
757 else if (component_i >= dim && component_j >= dim)
766 template<
typename Number>
770 Vector<Number> out (1);
771 const unsigned int n_rows = out.size();
772 for (
unsigned int r=0; r<n_rows; ++r)
778 template<
int dim,
typename Number>
782 return to_vector(s.operator
const Number &());
786 template<
int dim,
typename Number>
791 const unsigned int n_rows = out.size();
792 for (
unsigned int r=0; r<n_rows; ++r)
794 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
796 const unsigned int i = indices.first;
803 template<
int dim,
typename Number>
808 const unsigned int n_rows = out.size();
809 for (
unsigned int r=0; r<n_rows; ++r)
811 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
814 const unsigned int i = indices.first;
815 const unsigned int j = indices.second;
822 template<
int dim,
typename Number>
827 const unsigned int n_rows = out.size();
828 for (
unsigned int r=0; r<n_rows; ++r)
830 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
true);
834 const unsigned int i = indices.first;
835 const unsigned int j = indices.second;
837 const double factor = internal::vector_component_factor<dim>(r,
true);
839 out(r) = factor*st[i][j];
845 template<
typename Number>
855 template<
int dim,
typename Number>
859 return to_matrix(s.operator
const Number &());
863 template<
int dim,
typename Number>
868 const unsigned int n_rows = out.m();
869 const unsigned int n_cols = out.n();
870 for (
unsigned int r=0; r<n_rows; ++r)
872 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
874 const unsigned int i = indices.first;
876 for (
unsigned int c=0; c<n_cols; ++c)
886 template<
int dim,
typename Number>
891 const unsigned int n_rows = out.m();
892 const unsigned int n_cols = out.n();
893 for (
unsigned int r=0; r<n_rows; ++r)
895 const std::pair<unsigned int,unsigned int> indices_i = internal::indices_from_component<dim>(r,
false);
898 const unsigned int i = indices_i.first;
900 for (
unsigned int c=0; c<n_cols; ++c)
902 const std::pair<unsigned int,unsigned int> indices_j = internal::indices_from_component<dim>(c,
false);
905 const unsigned int j = indices_j.second;
914 template<
int dim,
typename Number>
926 template<
typename TensorType>
927 struct is_rank_2_symmetric_tensor : std::false_type
930 template<
int dim,
typename Number>
938 template<
int dim,
typename SubTensor1,
typename SubTensor2,
typename Number>
943 (SubTensor1::dimension == dim && SubTensor2::dimension == dim),
944 "Sub-tensor spatial dimension is different from those of the input tensor.");
947 (SubTensor1::rank == 2 && SubTensor2::rank == 1) ||
948 (SubTensor1::rank == 1 && SubTensor2::rank == 2),
949 "Cannot build a rank 3 tensor from the given combination of sub-tensors.");
951 FullMatrix<Number> out (SubTensor1::n_independent_components,SubTensor2::n_independent_components);
952 const unsigned int n_rows = out.m();
953 const unsigned int n_cols = out.n();
955 if (SubTensor1::rank == 2 && SubTensor2::rank == 1)
957 const bool subtensor_is_rank_2_symmetric_tensor = internal::is_rank_2_symmetric_tensor<SubTensor1>::value;
959 for (
unsigned int r=0; r<n_rows; ++r)
961 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,subtensor_is_rank_2_symmetric_tensor);
964 if (subtensor_is_rank_2_symmetric_tensor)
968 const unsigned int i = indices_ij.first;
969 const unsigned int j = indices_ij.second;
971 const double factor = internal::vector_component_factor<dim>(r,subtensor_is_rank_2_symmetric_tensor);
973 for (
unsigned int c=0; c<n_cols; ++c)
975 const std::pair<unsigned int,unsigned int> indices_k = internal::indices_from_component<dim>(c,
false);
977 const unsigned int k = indices_k.first;
979 if (subtensor_is_rank_2_symmetric_tensor)
980 out(r,c) = factor*t[i][j][k];
982 out(r,c) = t[i][j][k];
986 else if (SubTensor1::rank == 1 && SubTensor2::rank == 2)
988 const bool subtensor_is_rank_2_symmetric_tensor = internal::is_rank_2_symmetric_tensor<SubTensor2>::value;
990 for (
unsigned int r=0; r<n_rows; ++r)
992 const std::pair<unsigned int,unsigned int> indices_k = internal::indices_from_component<dim>(r,
false);
994 const unsigned int k = indices_k.first;
996 for (
unsigned int c=0; c<n_cols; ++c)
998 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(c,subtensor_is_rank_2_symmetric_tensor);
1001 if (subtensor_is_rank_2_symmetric_tensor)
1005 const unsigned int i = indices_ij.first;
1006 const unsigned int j = indices_ij.second;
1008 if (subtensor_is_rank_2_symmetric_tensor)
1010 const double factor = internal::vector_component_factor<dim>(c,subtensor_is_rank_2_symmetric_tensor);
1011 out(r,c) = factor*t[k][i][j];
1014 out(r,c) = t[k][i][j];
1027 template<
int dim,
typename Number>
1033 const unsigned int n_rows = out.m();
1034 const unsigned int n_cols = out.n();
1035 for (
unsigned int r=0; r<n_rows; ++r)
1037 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,
false);
1040 const unsigned int i = indices_ij.first;
1041 const unsigned int j = indices_ij.second;
1043 for (
unsigned int c=0; c<n_cols; ++c)
1045 const std::pair<unsigned int,unsigned int> indices_kl = internal::indices_from_component<dim>(c,
false);
1048 const unsigned int k = indices_kl.first;
1049 const unsigned int l = indices_kl.second;
1051 out(r,c) = t[i][j][k][
l];
1058 template<
int dim,
typename Number>
1064 const unsigned int n_rows = out.m();
1065 const unsigned int n_cols = out.n();
1066 for (
unsigned int r=0; r<n_rows; ++r)
1068 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,
true);
1072 const unsigned int i = indices_ij.first;
1073 const unsigned int j = indices_ij.second;
1075 for (
unsigned int c=0; c<n_cols; ++c)
1077 const std::pair<unsigned int,unsigned int> indices_kl = internal::indices_from_component<dim>(c,
true);
1081 const unsigned int k = indices_kl.first;
1082 const unsigned int l = indices_kl.second;
1084 const double factor = internal::matrix_component_factor<dim>(r,c,
true);
1086 out(r,c) = factor*st[i][j][k][
l];
1093 template<
typename Number>
1103 template<
int dim,
typename Number>
1108 return to_tensor(vec, s.operator Number &());
1112 template<
int dim,
typename Number>
1119 const unsigned int n_rows = vec.size();
1120 for (
unsigned int r=0; r<n_rows; ++r)
1122 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
1124 const unsigned int i = indices.first;
1130 template<
int dim,
typename Number>
1137 const unsigned int n_rows = vec.size();
1138 for (
unsigned int r=0; r<n_rows; ++r)
1140 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
1143 const unsigned int i = indices.first;
1144 const unsigned int j = indices.second;
1150 template<
int dim,
typename Number>
1157 const unsigned int n_rows = vec.size();
1158 for (
unsigned int r=0; r<n_rows; ++r)
1160 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
true);
1164 const unsigned int i = indices.first;
1165 const unsigned int j = indices.second;
1167 const double inv_factor = 1.0/internal::vector_component_factor<dim>(r,
true);
1169 st[i][j] = inv_factor*vec(r);
1174 template<
typename Number>
1186 template<
int dim,
typename Number>
1191 return to_tensor(mtrx, s.operator Number &());
1195 template<
int dim,
typename Number>
1205 const unsigned int n_rows = mtrx.
m();
1206 const unsigned int n_cols = mtrx.
n();
1207 for (
unsigned int r=0; r<n_rows; ++r)
1209 const std::pair<unsigned int,unsigned int> indices = internal::indices_from_component<dim>(r,
false);
1212 const unsigned int i = indices.first;
1214 for (
unsigned int c=0; c<n_cols; ++c)
1223 template<
int dim,
typename Number>
1233 const unsigned int n_rows = mtrx.
m();
1234 const unsigned int n_cols = mtrx.
n();
1235 for (
unsigned int r=0; r<n_rows; ++r)
1237 const std::pair<unsigned int,unsigned int> indices_i = internal::indices_from_component<dim>(r,
false);
1240 const unsigned int i = indices_i.first;
1242 for (
unsigned int c=0; c<n_cols; ++c)
1244 const std::pair<unsigned int,unsigned int> indices_j = internal::indices_from_component<dim>(c,
false);
1247 const unsigned int j = indices_j.second;
1249 t[i][j] = mtrx(r,c);
1256 template<
int dim,
typename Number>
1273 ExcMessage(
"The entries stored inside the matrix were not symmetric"));
1277 template<
int dim,
typename Number>
1299 const unsigned int n_rows = mtrx.
m();
1300 const unsigned int n_cols = mtrx.
n();
1311 const bool subtensor_is_rank_2_symmetric_tensor
1314 for (
unsigned int r=0; r<n_rows; ++r)
1316 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,subtensor_is_rank_2_symmetric_tensor);
1319 if (subtensor_is_rank_2_symmetric_tensor)
1323 const unsigned int i = indices_ij.first;
1324 const unsigned int j = indices_ij.second;
1326 const double inv_factor = 1.0/internal::vector_component_factor<dim>(r,subtensor_is_rank_2_symmetric_tensor);
1328 for (
unsigned int c=0; c<n_cols; ++c)
1330 const std::pair<unsigned int,unsigned int> indices_k = internal::indices_from_component<dim>(c,
false);
1332 const unsigned int k = indices_k.first;
1334 if (subtensor_is_rank_2_symmetric_tensor)
1336 t[i][j][k] = inv_factor*mtrx(r,c);
1337 t[j][i][k] = t[i][j][k];
1340 t[i][j][k] = mtrx(r,c);
1356 const bool subtensor_is_rank_2_symmetric_tensor
1359 for (
unsigned int r=0; r<n_rows; ++r)
1361 const std::pair<unsigned int,unsigned int> indices_k = internal::indices_from_component<dim>(r,
false);
1363 const unsigned int k = indices_k.first;
1365 for (
unsigned int c=0; c<n_cols; ++c)
1367 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(c,subtensor_is_rank_2_symmetric_tensor);
1370 if (subtensor_is_rank_2_symmetric_tensor)
1374 const unsigned int i = indices_ij.first;
1375 const unsigned int j = indices_ij.second;
1377 if (subtensor_is_rank_2_symmetric_tensor)
1379 const double inv_factor = 1.0/internal::vector_component_factor<dim>(c,subtensor_is_rank_2_symmetric_tensor);
1380 t[k][i][j] = inv_factor*mtrx(r,c);
1381 t[k][j][i] = t[k][i][j];
1384 t[k][i][j] = mtrx(r,c);
1391 template<
int dim,
typename Number>
1403 const unsigned int n_rows = mtrx.
m();
1404 const unsigned int n_cols = mtrx.
n();
1405 for (
unsigned int r=0; r<n_rows; ++r)
1407 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,
false);
1410 const unsigned int i = indices_ij.first;
1411 const unsigned int j = indices_ij.second;
1413 for (
unsigned int c=0; c<n_cols; ++c)
1415 const std::pair<unsigned int,unsigned int> indices_kl = internal::indices_from_component<dim>(c,
false);
1418 const unsigned int k = indices_kl.first;
1419 const unsigned int l = indices_kl.second;
1421 t[i][j][k][
l] = mtrx(r,c);
1427 template<
int dim,
typename Number>
1439 const unsigned int n_rows = mtrx.
m();
1440 const unsigned int n_cols = mtrx.
n();
1441 for (
unsigned int r=0; r<n_rows; ++r)
1443 const std::pair<unsigned int,unsigned int> indices_ij = internal::indices_from_component<dim>(r,
false);
1446 const unsigned int i = indices_ij.first;
1447 const unsigned int j = indices_ij.second;
1449 for (
unsigned int c=0; c<n_cols; ++c)
1451 const std::pair<unsigned int,unsigned int> indices_kl = internal::indices_from_component<dim>(c,
false);
1454 const unsigned int k = indices_kl.first;
1455 const unsigned int l = indices_kl.second;
1457 const double inv_factor = 1.0/internal::matrix_component_factor<dim>(r,c,
true);
1459 st[i][j][k][
l] = inv_factor*mtrx(r,c);
1465 template<
typename TensorType,
typename Number>
1475 template<
typename TensorType,
typename Number>
1492 DEAL_II_NAMESPACE_CLOSE
static const double SQRT2
void to_tensor(const Vector< Number > &vec, Number &s)
static const unsigned int n_independent_components
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize2(int arg1, int arg2, int arg3)
static const unsigned int n_independent_components
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize3(int arg1, int arg2, int arg3, int arg4)
Vector< Number > to_vector(const Number &s)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
size_type n_elements() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
FullMatrix< Number > to_matrix(const Number &s)
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize2(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcNotImplemented()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize3(int arg1, int arg2, int arg3, int arg4)