Reference documentation for deal.II version 9.1.0-pre
notation.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_physics_notation_h
17 #define dealii_physics_notation_h
18 
19 #include <deal.II/base/exceptions.h>
20 #include <deal.II/base/numbers.h>
22 #include <deal.II/base/tensor.h>
23 
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/lac/vector.h>
26 
27 #include <type_traits>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 
32 namespace Physics
33 {
34  namespace Notation
35  {
285  namespace Kelvin
286  {
291  int,
292  int,
293  int,
294  << "The number of rows in the input matrix is " << arg1
295  << ", but needs to be either " << arg2 << " or " << arg3
296  << ".");
297 
298 
303  int,
304  int,
305  int,
306  int,
307  << "The number of rows in the input matrix is " << arg1
308  << ", but needs to be either " << arg2 << "," << arg3
309  << ", or " << arg4 << ".");
310 
311 
316  int,
317  int,
318  int,
319  << "The number of columns in the input matrix is " << arg1
320  << ", but needs to be either " << arg2 << " or " << arg3
321  << ".");
322 
323 
328  int,
329  int,
330  int,
331  int,
332  << "The number of columns in the input matrix is " << arg1
333  << ", but needs to be either " << arg2 << "," << arg3
334  << ", or " << arg4 << ".");
335 
336 
341 
347  template <typename Number>
348  Vector<Number>
349  to_vector(const Number &s);
350 
351 
357  template <int dim, typename Number>
358  Vector<Number>
360 
361 
367  template <int dim, typename Number>
368  Vector<Number>
370 
371 
377  template <int dim, typename Number>
378  Vector<Number>
380 
381 
388  template <int dim, typename Number>
389  Vector<Number>
391 
392 
398  template <typename Number>
400  to_matrix(const Number &s);
401 
402 
408  template <int dim, typename Number>
411 
412 
418  template <int dim, typename Number>
421 
422 
428  template <int dim, typename Number>
431 
432 
442  template <int dim, typename Number>
445 
478  template <int dim,
479  typename SubTensor1 = Tensor<2, dim>,
480  typename SubTensor2 = Tensor<1, dim>,
481  typename Number>
484 
485 
492  template <int dim, typename Number>
495 
496 
504  template <int dim, typename Number>
507 
509 
514 
518  template <typename Number>
519  void
520  to_tensor(const Vector<Number> &vec, Number &s);
521 
522 
526  template <int dim, typename Number>
527  void
528  to_tensor(const Vector<Number> &vec, Tensor<0, dim, Number> &s);
529 
530 
534  template <int dim, typename Number>
535  void
536  to_tensor(const Vector<Number> &vec, Tensor<1, dim, Number> &v);
537 
538 
542  template <int dim, typename Number>
543  void
544  to_tensor(const Vector<Number> &vec, Tensor<2, dim, Number> &t);
545 
546 
550  template <int dim, typename Number>
551  void
552  to_tensor(const Vector<Number> &vec, SymmetricTensor<2, dim, Number> &st);
553 
554 
558  template <typename Number>
559  void
560  to_tensor(const FullMatrix<Number> &mtrx, Number &s);
561 
562 
566  template <int dim, typename Number>
567  void
569 
570 
574  template <int dim, typename Number>
575  void
577 
578 
582  template <int dim, typename Number>
583  void
585 
586 
590  template <int dim, typename Number>
591  void
592  to_tensor(const FullMatrix<Number> & mtrx,
594 
595 
605  template <int dim, typename Number>
606  void
608 
609 
613  template <int dim, typename Number>
614  void
616 
617 
621  template <int dim, typename Number>
622  void
623  to_tensor(const FullMatrix<Number> & mtrx,
625 
626 
631  template <typename TensorType, typename Number>
632  TensorType
633  to_tensor(const Vector<Number> &vec);
634 
635 
640  template <typename TensorType, typename Number>
641  TensorType
642  to_tensor(const FullMatrix<Number> &vec);
644 
645  } // namespace Kelvin
646 
647  } // namespace Notation
648 } // namespace Physics
649 
650 
651 #ifndef DOXYGEN
652 
653 
654 // ------------------------- inline functions ------------------------
655 
656 
657 namespace Physics
658 {
659  namespace Notation
660  {
661  namespace Kelvin
662  {
663  namespace internal
664  {
672  template <int dim>
673  std::pair<unsigned int, unsigned int>
674  indices_from_component(const unsigned int component_n,
675  const bool symmetric);
676 
677 
678  template <int dim>
679  std::pair<unsigned int, unsigned int>
680  indices_from_component(const unsigned int component_n, const bool)
681  {
682  AssertThrow(false, ExcNotImplemented());
683  return std::pair<unsigned int, unsigned int>();
684  }
685 
686 
687  template <>
688  inline std::pair<unsigned int, unsigned int>
689  indices_from_component<1>(const unsigned int component_n, const bool)
690  {
691  Assert(component_n < 1, ExcIndexRange(component_n, 0, 1));
692 
693  return std::make_pair(0u, 0u);
694  }
695 
696 
697  template <>
698  inline std::pair<unsigned int, unsigned int>
699  indices_from_component<2>(const unsigned int component_n,
700  const bool symmetric)
701  {
702  if (symmetric == true)
703  {
704  Assert(
706  ExcIndexRange(component_n,
707  0,
709  }
710  else
711  {
713  ExcIndexRange(component_n,
714  0,
716  }
717 
718  static const unsigned int indices[4][2] = {{0, 0},
719  {1, 1},
720  {0, 1},
721  {1, 0}};
722  return std::make_pair(indices[component_n][0],
723  indices[component_n][1]);
724  }
725 
726 
727  template <>
728  inline std::pair<unsigned int, unsigned int>
729  indices_from_component<3>(const unsigned int component_n,
730  const bool symmetric)
731  {
732  if (symmetric == true)
733  {
734  Assert(
736  ExcIndexRange(component_n,
737  0,
739  }
740  else
741  {
743  ExcIndexRange(component_n,
744  0,
746  }
747 
748  static const unsigned int indices[9][2] = {{0, 0},
749  {1, 1},
750  {2, 2},
751  {1, 2},
752  {0, 2},
753  {0, 1},
754  {1, 0},
755  {2, 0},
756  {2, 1}};
757  return std::make_pair(indices[component_n][0],
758  indices[component_n][1]);
759  }
760 
761 
766  template <int dim>
767  double
768  vector_component_factor(const unsigned int component_i,
769  const bool symmetric)
770  {
771  if (symmetric == false)
772  return 1.0;
773 
774  if (component_i < dim)
775  return 1.0;
776  else
777  return numbers::SQRT2;
778  }
779 
780 
785  template <int dim>
786  double
787  matrix_component_factor(const unsigned int component_i,
788  const unsigned int component_j,
789  const bool symmetric)
790  {
791  if (symmetric == false)
792  return 1.0;
793 
794  // This case check returns equivalent of this result:
795  // internal::vector_component_factor<dim>(component_i,symmetric)*internal::vector_component_factor<dim>(component_j,symmetric);
796  if (component_i < dim && component_j < dim)
797  return 1.0;
798  else if (component_i >= dim && component_j >= dim)
799  return 2.0;
800  else // ((component_i >= dim && component_j < dim) || (component_i <
801  // dim && component_j >= dim))
802  return numbers::SQRT2;
803  }
804 
805  } // namespace internal
806 
807 
808  template <typename Number>
809  Vector<Number>
810  to_vector(const Number &s)
811  {
812  Vector<Number> out(1);
813  const unsigned int n_rows = out.size();
814  for (unsigned int r = 0; r < n_rows; ++r)
815  out(r) = s;
816  return out;
817  }
818 
819 
820  template <int dim, typename Number>
821  Vector<Number>
823  {
824  return to_vector(s.operator const Number &());
825  }
826 
827 
828  template <int dim, typename Number>
829  Vector<Number>
831  {
832  Vector<Number> out(v.n_independent_components);
833  const unsigned int n_rows = out.size();
834  for (unsigned int r = 0; r < n_rows; ++r)
835  {
836  const std::pair<unsigned int, unsigned int> indices =
837  internal::indices_from_component<dim>(r, false);
838  Assert(indices.first < dim, ExcInternalError());
839  const unsigned int i = indices.first;
840  out(r) = v[i];
841  }
842  return out;
843  }
844 
845 
846  template <int dim, typename Number>
847  Vector<Number>
849  {
850  Vector<Number> out(t.n_independent_components);
851  const unsigned int n_rows = out.size();
852  for (unsigned int r = 0; r < n_rows; ++r)
853  {
854  const std::pair<unsigned int, unsigned int> indices =
855  internal::indices_from_component<dim>(r, false);
856  Assert(indices.first < dim, ExcInternalError());
857  Assert(indices.second < dim, ExcInternalError());
858  const unsigned int i = indices.first;
859  const unsigned int j = indices.second;
860  out(r) = t[i][j];
861  }
862  return out;
863  }
864 
865 
866  template <int dim, typename Number>
867  Vector<Number>
869  {
870  Vector<Number> out(st.n_independent_components);
871  const unsigned int n_rows = out.size();
872  for (unsigned int r = 0; r < n_rows; ++r)
873  {
874  const std::pair<unsigned int, unsigned int> indices =
875  internal::indices_from_component<dim>(r, true);
876  Assert(indices.first < dim, ExcInternalError());
877  Assert(indices.second < dim, ExcInternalError());
878  Assert(indices.second >= indices.first, ExcInternalError());
879  const unsigned int i = indices.first;
880  const unsigned int j = indices.second;
881 
882  const double factor =
883  internal::vector_component_factor<dim>(r, true);
884 
885  out(r) = factor * st[i][j];
886  }
887  return out;
888  }
889 
890 
891  template <typename Number>
893  to_matrix(const Number &s)
894  {
895  FullMatrix<Number> out(1, 1);
896  out(0, 0) = s;
897  return out;
898  }
899 
900 
901  template <int dim, typename Number>
904  {
905  return to_matrix(s.operator const Number &());
906  }
907 
908 
909  template <int dim, typename Number>
912  {
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)
917  {
918  const std::pair<unsigned int, unsigned int> indices =
919  internal::indices_from_component<dim>(r, false);
920  Assert(indices.first < dim, ExcInternalError());
921  const unsigned int i = indices.first;
922 
923  for (unsigned int c = 0; c < n_cols; ++c)
924  {
925  Assert(c < 1, ExcInternalError());
926  out(r, c) = v[i];
927  }
928  }
929  return out;
930  }
931 
932 
933  template <int dim, typename Number>
936  {
937  FullMatrix<Number> out(dim, dim);
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)
941  {
942  const std::pair<unsigned int, unsigned int> indices_i =
943  internal::indices_from_component<dim>(r, false);
944  Assert(indices_i.first < dim, ExcInternalError());
945  Assert(indices_i.second < dim, ExcInternalError());
946  const unsigned int i = indices_i.first;
947 
948  for (unsigned int c = 0; c < n_cols; ++c)
949  {
950  const std::pair<unsigned int, unsigned int> indices_j =
951  internal::indices_from_component<dim>(c, false);
952  Assert(indices_j.first < dim, ExcInternalError());
953  Assert(indices_j.second < dim, ExcInternalError());
954  const unsigned int j = indices_j.second;
955 
956  out(r, c) = t[i][j];
957  }
958  }
959  return out;
960  }
961 
962 
963  template <int dim, typename Number>
966  {
967  return to_matrix(Tensor<2, dim, Number>(st));
968  }
969 
970 
971  namespace internal
972  {
973  template <typename TensorType>
974  struct is_rank_2_symmetric_tensor : std::false_type
975  {};
976 
977  template <int dim, typename Number>
978  struct is_rank_2_symmetric_tensor<SymmetricTensor<2, dim, Number>>
979  : std::true_type
980  {};
981  } // namespace internal
982 
983 
984  template <int dim,
985  typename SubTensor1,
986  typename SubTensor2,
987  typename Number>
990  {
991  static_assert(
992  (SubTensor1::dimension == dim && SubTensor2::dimension == dim),
993  "Sub-tensor spatial dimension is different from those of the input tensor.");
994 
995  static_assert(
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.");
999 
1000  FullMatrix<Number> out(SubTensor1::n_independent_components,
1001  SubTensor2::n_independent_components);
1002  const unsigned int n_rows = out.m();
1003  const unsigned int n_cols = out.n();
1004 
1005  if (SubTensor1::rank == 2 && SubTensor2::rank == 1)
1006  {
1007  const bool subtensor_is_rank_2_symmetric_tensor =
1008  internal::is_rank_2_symmetric_tensor<SubTensor1>::value;
1009 
1010  for (unsigned int r = 0; r < n_rows; ++r)
1011  {
1012  const std::pair<unsigned int, unsigned int> indices_ij =
1013  internal::indices_from_component<dim>(
1014  r, subtensor_is_rank_2_symmetric_tensor);
1015  Assert(indices_ij.first < dim, ExcInternalError());
1016  Assert(indices_ij.second < dim, ExcInternalError());
1017  if (subtensor_is_rank_2_symmetric_tensor)
1018  {
1019  Assert(indices_ij.second >= indices_ij.first,
1020  ExcInternalError());
1021  }
1022  const unsigned int i = indices_ij.first;
1023  const unsigned int j = indices_ij.second;
1024 
1025  const double factor = internal::vector_component_factor<dim>(
1026  r, subtensor_is_rank_2_symmetric_tensor);
1027 
1028  for (unsigned int c = 0; c < n_cols; ++c)
1029  {
1030  const std::pair<unsigned int, unsigned int> indices_k =
1031  internal::indices_from_component<dim>(c, false);
1032  Assert(indices_k.first < dim, ExcInternalError());
1033  const unsigned int k = indices_k.first;
1034 
1035  if (subtensor_is_rank_2_symmetric_tensor)
1036  out(r, c) = factor * t[i][j][k];
1037  else
1038  out(r, c) = t[i][j][k];
1039  }
1040  }
1041  }
1042  else if (SubTensor1::rank == 1 && SubTensor2::rank == 2)
1043  {
1044  const bool subtensor_is_rank_2_symmetric_tensor =
1045  internal::is_rank_2_symmetric_tensor<SubTensor2>::value;
1046 
1047  for (unsigned int r = 0; r < n_rows; ++r)
1048  {
1049  const std::pair<unsigned int, unsigned int> indices_k =
1050  internal::indices_from_component<dim>(r, false);
1051  Assert(indices_k.first < dim, ExcInternalError());
1052  const unsigned int k = indices_k.first;
1053 
1054  for (unsigned int c = 0; c < n_cols; ++c)
1055  {
1056  const std::pair<unsigned int, unsigned int> indices_ij =
1057  internal::indices_from_component<dim>(
1058  c, subtensor_is_rank_2_symmetric_tensor);
1059  Assert(indices_ij.first < dim, ExcInternalError());
1060  Assert(indices_ij.second < dim, ExcInternalError());
1061  if (subtensor_is_rank_2_symmetric_tensor)
1062  {
1063  Assert(indices_ij.second >= indices_ij.first,
1064  ExcInternalError());
1065  }
1066  const unsigned int i = indices_ij.first;
1067  const unsigned int j = indices_ij.second;
1068 
1069  if (subtensor_is_rank_2_symmetric_tensor)
1070  {
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];
1075  }
1076  else
1077  out(r, c) = t[k][i][j];
1078  }
1079  }
1080  }
1081  else
1082  {
1083  AssertThrow(false, ExcNotImplemented());
1084  }
1085 
1086  return out;
1087  }
1088 
1089 
1090  template <int dim, typename Number>
1093  {
1094  FullMatrix<Number> out(
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)
1100  {
1101  const std::pair<unsigned int, unsigned int> indices_ij =
1102  internal::indices_from_component<dim>(r, false);
1103  Assert(indices_ij.first < dim, ExcInternalError());
1104  Assert(indices_ij.second < dim, ExcInternalError());
1105  const unsigned int i = indices_ij.first;
1106  const unsigned int j = indices_ij.second;
1107 
1108  for (unsigned int c = 0; c < n_cols; ++c)
1109  {
1110  const std::pair<unsigned int, unsigned int> indices_kl =
1111  internal::indices_from_component<dim>(c, false);
1112  Assert(indices_kl.first < dim, ExcInternalError());
1113  Assert(indices_kl.second < dim, ExcInternalError());
1114  const unsigned int k = indices_kl.first;
1115  const unsigned int l = indices_kl.second;
1116 
1117  out(r, c) = t[i][j][k][l];
1118  }
1119  }
1120  return out;
1121  }
1122 
1123 
1124  template <int dim, typename Number>
1127  {
1128  FullMatrix<Number> out(
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)
1134  {
1135  const std::pair<unsigned int, unsigned int> indices_ij =
1136  internal::indices_from_component<dim>(r, true);
1137  Assert(indices_ij.first < dim, ExcInternalError());
1138  Assert(indices_ij.second < dim, ExcInternalError());
1139  Assert(indices_ij.second >= indices_ij.first, ExcInternalError());
1140  const unsigned int i = indices_ij.first;
1141  const unsigned int j = indices_ij.second;
1142 
1143  for (unsigned int c = 0; c < n_cols; ++c)
1144  {
1145  const std::pair<unsigned int, unsigned int> indices_kl =
1146  internal::indices_from_component<dim>(c, true);
1147  Assert(indices_kl.first < dim, ExcInternalError());
1148  Assert(indices_kl.second < dim, ExcInternalError());
1149  Assert(indices_kl.second >= indices_kl.first,
1150  ExcInternalError());
1151  const unsigned int k = indices_kl.first;
1152  const unsigned int l = indices_kl.second;
1153 
1154  const double factor =
1155  internal::matrix_component_factor<dim>(r, c, true);
1156 
1157  out(r, c) = factor * st[i][j][k][l];
1158  }
1159  }
1160  return out;
1161  }
1162 
1163 
1164  template <typename Number>
1165  void
1166  to_tensor(const Vector<Number> &vec, Number &s)
1167  {
1168  Assert(vec.size() == 1, ExcDimensionMismatch(vec.size(), 1));
1169  s = vec(0);
1170  }
1171 
1172 
1173  template <int dim, typename Number>
1174  void
1175  to_tensor(const Vector<Number> &vec, Tensor<0, dim, Number> &s)
1176  {
1177  return to_tensor(vec, s.operator Number &());
1178  }
1179 
1180 
1181  template <int dim, typename Number>
1182  void
1183  to_tensor(const Vector<Number> &vec, Tensor<1, dim, Number> &v)
1184  {
1185  Assert(vec.size() == v.n_independent_components,
1187  const unsigned int n_rows = vec.size();
1188  for (unsigned int r = 0; r < n_rows; ++r)
1189  {
1190  const std::pair<unsigned int, unsigned int> indices =
1191  internal::indices_from_component<dim>(r, false);
1192  Assert(indices.first < dim, ExcInternalError());
1193  const unsigned int i = indices.first;
1194  v[i] = vec(r);
1195  }
1196  }
1197 
1198 
1199  template <int dim, typename Number>
1200  void
1201  to_tensor(const Vector<Number> &vec, Tensor<2, dim, Number> &t)
1202  {
1203  Assert(vec.size() == t.n_independent_components,
1205  const unsigned int n_rows = vec.size();
1206  for (unsigned int r = 0; r < n_rows; ++r)
1207  {
1208  const std::pair<unsigned int, unsigned int> indices =
1209  internal::indices_from_component<dim>(r, false);
1210  Assert(indices.first < dim, ExcInternalError());
1211  Assert(indices.second < dim, ExcInternalError());
1212  const unsigned int i = indices.first;
1213  const unsigned int j = indices.second;
1214  t[i][j] = vec(r);
1215  }
1216  }
1217 
1218 
1219  template <int dim, typename Number>
1220  void
1221  to_tensor(const Vector<Number> &vec, SymmetricTensor<2, dim, Number> &st)
1222  {
1223  Assert(vec.size() == st.n_independent_components,
1225  const unsigned int n_rows = vec.size();
1226  for (unsigned int r = 0; r < n_rows; ++r)
1227  {
1228  const std::pair<unsigned int, unsigned int> indices =
1229  internal::indices_from_component<dim>(r, true);
1230  Assert(indices.first < dim, ExcInternalError());
1231  Assert(indices.second < dim, ExcInternalError());
1232  Assert(indices.second >= indices.first, ExcInternalError());
1233  const unsigned int i = indices.first;
1234  const unsigned int j = indices.second;
1235 
1236  const double inv_factor =
1237  1.0 / internal::vector_component_factor<dim>(r, true);
1238 
1239  st[i][j] = inv_factor * vec(r);
1240  }
1241  }
1242 
1243 
1244  template <typename Number>
1245  void
1246  to_tensor(const FullMatrix<Number> &mtrx, Number &s)
1247  {
1248  Assert(mtrx.m() == 1, ExcDimensionMismatch(mtrx.m(), 1));
1249  Assert(mtrx.n() == 1, ExcDimensionMismatch(mtrx.n(), 1));
1250  Assert(mtrx.n_elements() == 1,
1251  ExcDimensionMismatch(mtrx.n_elements(), 1));
1252  s = mtrx(0, 0);
1253  }
1254 
1255 
1256  template <int dim, typename Number>
1257  void
1259  {
1260  return to_tensor(mtrx, s.operator Number &());
1261  }
1262 
1263 
1264  template <int dim, typename Number>
1265  void
1267  {
1268  Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim));
1269  Assert(mtrx.n() == 1, ExcDimensionMismatch(mtrx.n(), 1));
1273 
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)
1277  {
1278  const std::pair<unsigned int, unsigned int> indices =
1279  internal::indices_from_component<dim>(r, false);
1280  Assert(indices.first < dim, ExcInternalError());
1281  Assert(indices.second == 0, ExcInternalError());
1282  const unsigned int i = indices.first;
1283 
1284  for (unsigned int c = 0; c < n_cols; ++c)
1285  {
1286  Assert(c < 1, ExcInternalError());
1287  v[i] = mtrx(r, c);
1288  }
1289  }
1290  }
1291 
1292 
1293  template <int dim, typename Number>
1294  void
1296  {
1297  Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim));
1298  Assert(mtrx.n() == dim, ExcDimensionMismatch(mtrx.n(), dim));
1302 
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)
1306  {
1307  const std::pair<unsigned int, unsigned int> indices_i =
1308  internal::indices_from_component<dim>(r, false);
1309  Assert(indices_i.first < dim, ExcInternalError());
1310  Assert(indices_i.second < dim, ExcInternalError());
1311  const unsigned int i = indices_i.first;
1312 
1313  for (unsigned int c = 0; c < n_cols; ++c)
1314  {
1315  const std::pair<unsigned int, unsigned int> indices_j =
1316  internal::indices_from_component<dim>(c, false);
1317  Assert(indices_j.first < dim, ExcInternalError());
1318  Assert(indices_j.second < dim, ExcInternalError());
1319  const unsigned int j = indices_j.second;
1320 
1321  t[i][j] = mtrx(r, c);
1322  }
1323  }
1324  }
1325 
1326 
1327  template <int dim, typename Number>
1328  void
1329  to_tensor(const FullMatrix<Number> & mtrx,
1331  {
1332  // Its impossible to fit the (dim^2 + dim)/2 entries into a square
1333  // matrix We therefore assume that its been converted to a standard
1334  // tensor format using to_matrix (SymmetricTensor<2,dim,Number>) at some
1335  // point...
1336  Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim));
1337  Assert(mtrx.n() == dim, ExcDimensionMismatch(mtrx.n(), dim));
1338  Assert((mtrx.n_elements() ==
1341  mtrx.n_elements(),
1343 
1345  to_tensor(mtrx, tmp);
1346  st = symmetrize(tmp);
1347  Assert((Tensor<2, dim, Number>(st) - tmp).norm() < 1e-12,
1348  ExcMessage(
1349  "The entries stored inside the matrix were not symmetric"));
1350  }
1351 
1352 
1353  template <int dim, typename Number>
1354  void
1356  {
1358  (mtrx.m() ==
1360  (mtrx.m() ==
1363  mtrx.m(),
1368  (mtrx.n() ==
1369  Tensor<2, dim, Number>::n_independent_components) ||
1370  (mtrx.n() ==
1373  mtrx.n(),
1375  Tensor<2, dim, Number>::n_independent_components,
1377 
1378  const unsigned int n_rows = mtrx.m();
1379  const unsigned int n_cols = mtrx.n();
1381  {
1382  Assert(
1383  (mtrx.m() == Tensor<2, dim, Number>::n_independent_components) ||
1384  (mtrx.m() ==
1387  mtrx.m(),
1388  Tensor<2, dim, Number>::n_independent_components,
1390 
1391  const bool subtensor_is_rank_2_symmetric_tensor =
1392  (mtrx.m() ==
1394 
1395  for (unsigned int r = 0; r < n_rows; ++r)
1396  {
1397  const std::pair<unsigned int, unsigned int> indices_ij =
1398  internal::indices_from_component<dim>(
1399  r, subtensor_is_rank_2_symmetric_tensor);
1400  Assert(indices_ij.first < dim, ExcInternalError());
1401  Assert(indices_ij.second < dim, ExcInternalError());
1402  if (subtensor_is_rank_2_symmetric_tensor)
1403  {
1404  Assert(indices_ij.second >= indices_ij.first,
1405  ExcInternalError());
1406  }
1407  const unsigned int i = indices_ij.first;
1408  const unsigned int j = indices_ij.second;
1409 
1410  const double inv_factor =
1411  1.0 / internal::vector_component_factor<dim>(
1412  r, subtensor_is_rank_2_symmetric_tensor);
1413 
1414  for (unsigned int c = 0; c < n_cols; ++c)
1415  {
1416  const std::pair<unsigned int, unsigned int> indices_k =
1417  internal::indices_from_component<dim>(c, false);
1418  Assert(indices_k.first < dim, ExcInternalError());
1419  const unsigned int k = indices_k.first;
1420 
1421  if (subtensor_is_rank_2_symmetric_tensor)
1422  {
1423  t[i][j][k] = inv_factor * mtrx(r, c);
1424  t[j][i][k] = t[i][j][k];
1425  }
1426  else
1427  t[i][j][k] = mtrx(r, c);
1428  }
1429  }
1430  }
1431  else
1432  {
1433  Assert(
1437  Assert(
1438  (mtrx.n() == Tensor<2, dim, Number>::n_independent_components) ||
1439  (mtrx.n() ==
1442  mtrx.n(),
1443  Tensor<2, dim, Number>::n_independent_components,
1445 
1446  const bool subtensor_is_rank_2_symmetric_tensor =
1447  (mtrx.n() ==
1449 
1450  for (unsigned int r = 0; r < n_rows; ++r)
1451  {
1452  const std::pair<unsigned int, unsigned int> indices_k =
1453  internal::indices_from_component<dim>(r, false);
1454  Assert(indices_k.first < dim, ExcInternalError());
1455  const unsigned int k = indices_k.first;
1456 
1457  for (unsigned int c = 0; c < n_cols; ++c)
1458  {
1459  const std::pair<unsigned int, unsigned int> indices_ij =
1460  internal::indices_from_component<dim>(
1461  c, subtensor_is_rank_2_symmetric_tensor);
1462  Assert(indices_ij.first < dim, ExcInternalError());
1463  Assert(indices_ij.second < dim, ExcInternalError());
1464  if (subtensor_is_rank_2_symmetric_tensor)
1465  {
1466  Assert(indices_ij.second >= indices_ij.first,
1467  ExcInternalError());
1468  }
1469  const unsigned int i = indices_ij.first;
1470  const unsigned int j = indices_ij.second;
1471 
1472  if (subtensor_is_rank_2_symmetric_tensor)
1473  {
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];
1479  }
1480  else
1481  t[k][i][j] = mtrx(r, c);
1482  }
1483  }
1484  }
1485  }
1486 
1487 
1488  template <int dim, typename Number>
1489  void
1491  {
1501 
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)
1505  {
1506  const std::pair<unsigned int, unsigned int> indices_ij =
1507  internal::indices_from_component<dim>(r, false);
1508  Assert(indices_ij.first < dim, ExcInternalError());
1509  Assert(indices_ij.second < dim, ExcInternalError());
1510  const unsigned int i = indices_ij.first;
1511  const unsigned int j = indices_ij.second;
1512 
1513  for (unsigned int c = 0; c < n_cols; ++c)
1514  {
1515  const std::pair<unsigned int, unsigned int> indices_kl =
1516  internal::indices_from_component<dim>(c, false);
1517  Assert(indices_kl.first < dim, ExcInternalError());
1518  Assert(indices_kl.second < dim, ExcInternalError());
1519  const unsigned int k = indices_kl.first;
1520  const unsigned int l = indices_kl.second;
1521 
1522  t[i][j][k][l] = mtrx(r, c);
1523  }
1524  }
1525  }
1526 
1527 
1528  template <int dim, typename Number>
1529  void
1530  to_tensor(const FullMatrix<Number> & mtrx,
1532  {
1533  Assert((mtrx.m() ==
1536  mtrx.m(),
1538  Assert((mtrx.n() ==
1541  mtrx.n(),
1546 
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)
1550  {
1551  const std::pair<unsigned int, unsigned int> indices_ij =
1552  internal::indices_from_component<dim>(r, false);
1553  Assert(indices_ij.first < dim, ExcInternalError());
1554  Assert(indices_ij.second < dim, ExcInternalError());
1555  const unsigned int i = indices_ij.first;
1556  const unsigned int j = indices_ij.second;
1557 
1558  for (unsigned int c = 0; c < n_cols; ++c)
1559  {
1560  const std::pair<unsigned int, unsigned int> indices_kl =
1561  internal::indices_from_component<dim>(c, false);
1562  Assert(indices_kl.first < dim, ExcInternalError());
1563  Assert(indices_kl.second < dim, ExcInternalError());
1564  const unsigned int k = indices_kl.first;
1565  const unsigned int l = indices_kl.second;
1566 
1567  const double inv_factor =
1568  1.0 / internal::matrix_component_factor<dim>(r, c, true);
1569 
1570  st[i][j][k][l] = inv_factor * mtrx(r, c);
1571  }
1572  }
1573  }
1574 
1575 
1576  template <typename TensorType, typename Number>
1577  inline TensorType
1578  to_tensor(const Vector<Number> &vec)
1579  {
1580  TensorType out;
1581  to_tensor(vec, out);
1582  return out;
1583  }
1584 
1585 
1586  template <typename TensorType, typename Number>
1587  inline TensorType
1588  to_tensor(const FullMatrix<Number> &mtrx)
1589  {
1590  TensorType out;
1591  to_tensor(mtrx, out);
1592  return out;
1593  }
1594 
1595  } // namespace Kelvin
1596  } // namespace Notation
1597 } // namespace Physics
1598 
1599 
1600 #endif // DOXYGEN
1601 
1602 
1603 DEAL_II_NAMESPACE_CLOSE
1604 
1605 #endif
static const double SQRT2
Definition: numbers.h:158
static constexpr unsigned int n_independent_components
void to_tensor(const Vector< Number > &vec, Number &s)
static::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize2(int arg1, int arg2, int arg3)
static const unsigned int n_independent_components
Definition: tensor.h:419
static::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize3(int arg1, int arg2, int arg3, int arg4)
Vector< Number > to_vector(const Number &s)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1227
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)
size_type m() const
size_type n_elements() const
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:444
static::ExceptionBase & ExcNotationExcFullMatrixToTensorRowSize2(int arg1, int arg2, int arg3)
static::ExceptionBase & ExcNotImplemented()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:432
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & ExcNotationExcFullMatrixToTensorColSize3(int arg1, int arg2, int arg3, int arg4)