Reference documentation for deal.II version GIT 3e4283dc79 2023-06-10 12:25:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symengine_tensor_operations.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2021 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_differentiation_sd_symengine_tensor_operations_h
17 #define dealii_differentiation_sd_symengine_tensor_operations_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_WITH_SYMENGINE
22 
24 # include <deal.II/base/tensor.h>
25 
29 
30 # include <utility>
31 # include <vector>
32 
34 
35 namespace Differentiation
36 {
37  namespace SD
38  {
62  template <int dim>
64  make_vector_of_symbols(const std::string &symbol);
65 
85  template <int rank, int dim>
87  make_tensor_of_symbols(const std::string &symbol);
88 
108  template <int rank, int dim>
110  make_symmetric_tensor_of_symbols(const std::string &symbol);
111 
131  template <int dim>
133  make_vector_of_symbolic_functions(const std::string & symbol,
134  const types::substitution_map &arguments);
135 
156  template <int rank, int dim>
158  make_tensor_of_symbolic_functions(const std::string & symbol,
159  const types::substitution_map &arguments);
160 
181  template <int rank, int dim>
184  const std::string & symbol,
185  const types::substitution_map &arguments);
186 
203  template <int rank, int dim>
206 
216  template <int rank, int dim>
220 
230  template <int rank, int dim>
234 
244  template <int rank, int dim>
248 
258  template <int rank, int dim>
261 
272  template <int rank, int dim>
275  const Expression & x);
276 
287  template <int rank, int dim>
290  const Tensor<0, dim, Expression> & x);
291 
303  template <int rank, int dim>
306  const Tensor<0, dim, Expression> & x);
307 
318  template <int rank_1, int rank_2, int dim>
322 
335  template <int rank_1, int rank_2, int dim>
339 
350  template <int rank_1, int rank_2, int dim>
354 
365  template <int rank_1, int rank_2, int dim>
369 
401  template <bool ignore_invalid_symbols = false,
402  typename ValueType = double,
403  int rank,
404  int dim,
405  typename SymbolicType>
406  void
408  const Tensor<rank, dim, SymbolicType> &symbol_tensor);
409 
434  template <bool ignore_invalid_symbols = false,
435  typename ValueType = double,
436  int rank,
437  int dim,
438  typename SymbolicType>
439  void
441  types::substitution_map & symbol_map,
442  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor);
443 
464  template <int rank, int dim, typename SymbolicType, typename ValueType>
465  void
468  const Tensor<rank, dim, SymbolicType> &symbol_tensor,
469  const Tensor<rank, dim, ValueType> & value_tensor);
470 
491  template <int rank, int dim, typename SymbolicType, typename ValueType>
492  void
495  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
496  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
497 
526  template <int rank, int dim, typename ExpressionType, typename ValueType>
529  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
530  const Tensor<rank, dim, ValueType> & value_tensor);
531 
553  template <int rank, int dim, typename ExpressionType, typename ValueType>
556  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
557  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
558 
592  template <bool ignore_invalid_symbols = false,
593  int rank,
594  int dim,
595  typename ExpressionType,
596  typename ValueType>
597  void
600  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
601  const Tensor<rank, dim, ValueType> & value_tensor);
602 
628  template <bool ignore_invalid_symbols = false,
629  int rank,
630  int dim,
631  typename ExpressionType,
632  typename ValueType>
633  void
636  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
637  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
638 
665  template <int rank, int dim>
667  substitute(const Tensor<rank, dim, Expression> &expression_tensor,
669 
689  template <int rank, int dim>
693 
720  template <typename ValueType, int rank, int dim>
723  const Tensor<rank, dim, Expression> &expression_tensor,
725 
752  template <typename ValueType, int rank, int dim>
755  const SymmetricTensor<rank, dim, Expression> &expression_tensor,
757 
760  } // namespace SD
761 } // namespace Differentiation
762 
763 
764 /* -------------------- inline and template functions ------------------ */
765 
766 
767 # ifndef DOXYGEN
768 
769 namespace Differentiation
770 {
771  namespace SD
772  {
773  /* ---------------- Symbolic differentiation --------------*/
774 
775 
776  namespace internal
777  {
778  template <int dim>
780  make_rank_4_tensor_indices(const unsigned int idx_i,
781  const unsigned int idx_j)
782  {
783  const TableIndices<2> indices_i(
785  const TableIndices<2> indices_j(
787  return TableIndices<4>(indices_i[0],
788  indices_i[1],
789  indices_j[0],
790  indices_j[1]);
791  }
792 
793 
794  template <int rank_1, int rank_2>
796  concatenate_indices(const TableIndices<rank_1> &indices_1,
797  const TableIndices<rank_2> &indices_2)
798  {
799  TableIndices<rank_1 + rank_2> indices_out;
800  for (unsigned int i = 0; i < rank_1; ++i)
801  indices_out[i] = indices_1[i];
802  for (unsigned int j = 0; j < rank_2; ++j)
803  indices_out[rank_1 + j] = indices_2[j];
804  return indices_out;
805  }
806 
807 
808  template <int rank>
810  transpose_indices(const TableIndices<rank> &indices)
811  {
812  return indices;
813  }
814 
815 
816  template <>
817  inline TableIndices<2>
818  transpose_indices(const TableIndices<2> &indices)
819  {
820  return TableIndices<2>(indices[1], indices[0]);
821  }
822 
823 
824  template <int rank, int dim, typename ValueType>
825  bool
826  is_symmetric_component(const TableIndices<rank> &,
828  {
829  return false;
830  }
831 
832 
833  template <int rank, int dim, typename ValueType>
834  bool
835  is_symmetric_component(const TableIndices<rank> &,
837  {
838  static_assert(
839  rank == 0 || rank == 2,
840  "Querying symmetric component for non rank-2 symmetric tensor index is not allowed.");
841  return false;
842  }
843 
844 
845  template <int dim, typename ValueType>
846  bool
847  is_symmetric_component(const TableIndices<2> &table_indices,
849  {
850  return table_indices[0] != table_indices[1];
851  }
852 
853 
854  template <int dim,
855  typename ValueType = Expression,
856  template <int, int, typename>
857  class TensorType>
858  TensorType<0, dim, ValueType>
859  scalar_diff_tensor(const ValueType & func,
860  const TensorType<0, dim, ValueType> &op)
861  {
862  return differentiate(func, op);
863  }
864 
865 
866  template <int rank,
867  int dim,
868  typename ValueType = Expression,
869  template <int, int, typename>
870  class TensorType>
871  TensorType<rank, dim, ValueType>
872  scalar_diff_tensor(const ValueType & func,
873  const TensorType<rank, dim, ValueType> &op)
874  {
875  TensorType<rank, dim, ValueType> out;
876  for (unsigned int i = 0; i < out.n_independent_components; ++i)
877  {
878  const TableIndices<rank> indices(
879  out.unrolled_to_component_indices(i));
880  out[indices] = differentiate(func, op[indices]);
881 
882  if (is_symmetric_component(indices, op))
883  out[indices] *= 0.5;
884  }
885  return out;
886  }
887 
888 
889  // Specialization for rank-0 tensor
890  template <int rank,
891  int dim,
892  typename ValueType = Expression,
893  template <int, int, typename>
894  class TensorType>
895  TensorType<rank, dim, ValueType>
896  tensor_diff_tensor(const TensorType<0, dim, ValueType> & func,
897  const TensorType<rank, dim, ValueType> &op)
898  {
899  TensorType<rank, dim, ValueType> out;
900  for (unsigned int i = 0; i < out.n_independent_components; ++i)
901  {
902  const TableIndices<rank> indices(
903  out.unrolled_to_component_indices(i));
904  out[indices] = differentiate(func, op[indices]);
905 
906  if (is_symmetric_component(indices, op))
907  out[indices] *= 0.5;
908  }
909  return out;
910  }
911 
912 
913  template <int rank,
914  int dim,
915  typename ValueType = Expression,
916  template <int, int, typename>
917  class TensorType>
918  TensorType<rank, dim, ValueType>
919  tensor_diff_scalar(const TensorType<rank, dim, ValueType> &funcs,
920  const ValueType & op)
921  {
922  TensorType<rank, dim, ValueType> out;
923  for (unsigned int i = 0; i < out.n_independent_components; ++i)
924  {
925  const TableIndices<rank> indices(
926  out.unrolled_to_component_indices(i));
927  out[indices] = differentiate(funcs[indices], op);
928  }
929  return out;
930  }
931 
932 
933  // Specialization for rank-0 tensor
934  template <int rank,
935  int dim,
936  typename ValueType = Expression,
937  template <int, int, typename>
938  class TensorType>
939  TensorType<rank, dim, ValueType>
940  tensor_diff_tensor(const TensorType<rank, dim, ValueType> &funcs,
941  const TensorType<0, dim, ValueType> & op)
942  {
943  TensorType<rank, dim, ValueType> out;
944  for (unsigned int i = 0; i < out.n_independent_components; ++i)
945  {
946  const TableIndices<rank> indices(
947  out.unrolled_to_component_indices(i));
948  out[indices] = differentiate(funcs[indices], op);
949  }
950  return out;
951  }
952 
953 
954  // For either symmetric or normal tensors
955  template <int rank_1,
956  int rank_2,
957  int dim,
958  typename ValueType = Expression,
959  template <int, int, typename>
960  class TensorType>
961  TensorType<rank_1 + rank_2, dim, ValueType>
962  tensor_diff_tensor(const TensorType<rank_1, dim, ValueType> &funcs,
963  const TensorType<rank_2, dim, ValueType> &op)
964  {
965  TensorType<rank_1 + rank_2, dim, ValueType> out;
966  for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
967  {
968  const TableIndices<rank_1> indices_i(
969  funcs.unrolled_to_component_indices(i));
970  for (unsigned int j = 0; j < op.n_independent_components; ++j)
971  {
972  const TableIndices<rank_2> indices_j(
973  op.unrolled_to_component_indices(j));
974  const TableIndices<rank_1 + rank_2> indices_out =
975  concatenate_indices(indices_i, indices_j);
976 
977  out[indices_out] =
978  differentiate(funcs[indices_i], op[indices_j]);
979 
980  if (is_symmetric_component(indices_j, op))
981  out[indices_out] *= 0.5;
982  }
983  }
984  return out;
985  }
986 
987 
988  // For mixed symmetric/standard tensors
989  // The return type is always a standard tensor, since we cannot be sure
990  // that any symmetries exist in either the function tensor or the
991  // differential operator.
992  template <int rank_1,
993  int rank_2,
994  int dim,
995  typename ValueType = Expression,
996  template <int, int, typename>
997  class TensorType_1,
998  template <int, int, typename>
999  class TensorType_2>
1001  tensor_diff_tensor(const TensorType_1<rank_1, dim, ValueType> &funcs,
1002  const TensorType_2<rank_2, dim, ValueType> &op)
1003  {
1005  for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
1006  {
1007  const TableIndices<rank_1> indices_i(
1008  funcs.unrolled_to_component_indices(i));
1009  for (unsigned int j = 0; j < op.n_independent_components; ++j)
1010  {
1011  const TableIndices<rank_2> indices_j(
1012  op.unrolled_to_component_indices(j));
1013  const TableIndices<rank_1 + rank_2> indices_out =
1014  concatenate_indices(indices_i, indices_j);
1015 
1016  out[indices_out] =
1017  differentiate(funcs[indices_i], op[indices_j]);
1018 
1019  if (is_symmetric_component(indices_j, op))
1020  out[indices_out] *= 0.5;
1021 
1022  // TODO: Implement for SymmetricTensor<4,dim,...>
1023  if (std::is_same<TensorType_1<rank_1, dim, ValueType>,
1025  value) // Symmetric function
1026  {
1027  const TableIndices<rank_1 + rank_2> indices_out_t =
1028  concatenate_indices(transpose_indices(indices_i),
1029  indices_j);
1030  out[indices_out_t] = out[indices_out];
1031  }
1032  else if (std::is_same<TensorType_2<rank_2, dim, ValueType>,
1034  value) // Symmetric operator
1035  {
1036  const TableIndices<rank_1 + rank_2> indices_out_t =
1037  concatenate_indices(indices_i,
1038  transpose_indices(indices_j));
1039  out[indices_out_t] = out[indices_out];
1040  }
1041  else
1042  {
1043  Assert(
1044  false,
1045  ExcMessage(
1046  "Expect mixed tensor differentiation to have at least "
1047  "one component stemming from a symmetric tensor."));
1048  }
1049  }
1050  }
1051  return out;
1052  }
1053 
1054  } // namespace internal
1055 
1056 
1057  template <int rank, int dim>
1059  differentiate(const Expression & func,
1061  {
1062  return internal::scalar_diff_tensor(func, op);
1063  }
1064 
1065 
1066  template <int rank, int dim>
1068  differentiate(const Expression & func,
1070  {
1071  return internal::scalar_diff_tensor(func, op);
1072  }
1073 
1074 
1075  template <int rank, int dim>
1079  {
1080  return internal::scalar_diff_tensor(func, op);
1081  }
1082 
1083 
1084  template <int rank, int dim>
1088  {
1089  // Ensure that this returns a symmetric tensor by
1090  // invoking the scalar value function
1091  const Expression tmp = func;
1092  return internal::scalar_diff_tensor(tmp, op);
1093  }
1094 
1095 
1096  template <int rank, int dim>
1098  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1099  const Expression & op)
1100  {
1101  return internal::tensor_diff_scalar(symbol_tensor, op);
1102  }
1103 
1104 
1105  template <int rank, int dim>
1107  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1108  const Tensor<0, dim, Expression> & op)
1109  {
1110  return internal::tensor_diff_scalar(symbol_tensor, op);
1111  }
1112 
1113 
1114  template <int rank, int dim>
1117  const Expression & op)
1118  {
1119  return internal::tensor_diff_scalar(symbol_tensor, op);
1120  }
1121 
1122 
1123  template <int rank, int dim>
1126  const Tensor<0, dim, Expression> & op)
1127  {
1128  return internal::tensor_diff_scalar(symbol_tensor, op);
1129  }
1130 
1131 
1132  template <int rank_1, int rank_2, int dim>
1134  differentiate(const Tensor<rank_1, dim, Expression> &symbol_tensor,
1136  {
1137  return internal::tensor_diff_tensor(symbol_tensor, op);
1138  }
1139 
1140 
1141  template <int rank_1, int rank_2, int dim>
1145  {
1146  return internal::tensor_diff_tensor(symbol_tensor, op);
1147  }
1148 
1149 
1150  template <int rank_1, int rank_2, int dim>
1152  differentiate(const Tensor<rank_1, dim, Expression> & symbol_tensor,
1154  {
1155  return internal::tensor_diff_tensor(symbol_tensor, op);
1156  }
1157 
1158 
1159  template <int rank_1, int rank_2, int dim>
1163  {
1164  return internal::tensor_diff_tensor(symbol_tensor, op);
1165  }
1166 
1167 
1168  /* ---------------- Symbol map creation and manipulation --------------*/
1169 
1170 
1171  namespace internal
1172  {
1173  template <typename SymbolicType,
1174  typename ValueType,
1175  int rank,
1176  int dim,
1177  template <int, int, typename>
1178  class TensorType>
1179  void
1180  set_tensor_value_in_symbol_map(
1182  const TensorType<rank, dim, SymbolicType> &symbol_tensor,
1183  const TensorType<rank, dim, ValueType> & value_tensor)
1184  {
1185  TensorType<rank, dim, Expression> out;
1186  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1187  {
1188  const TableIndices<rank> indices(
1189  out.unrolled_to_component_indices(i));
1191  symbol_tensor[indices],
1192  value_tensor[indices]);
1193  }
1194  }
1195 
1196 
1197  template <typename SymbolicType, typename ValueType, int dim>
1198  void
1199  set_tensor_value_in_symbol_map(
1201  const SymmetricTensor<4, dim, SymbolicType> &symbol_tensor,
1202  const SymmetricTensor<4, dim, ValueType> & value_tensor)
1203  {
1205  for (unsigned int i = 0;
1206  i < SymmetricTensor<2, dim>::n_independent_components;
1207  ++i)
1208  for (unsigned int j = 0;
1209  j < SymmetricTensor<2, dim>::n_independent_components;
1210  ++j)
1211  {
1212  const TableIndices<4> indices =
1213  make_rank_4_tensor_indices<dim>(i, j);
1215  symbol_tensor[indices],
1216  value_tensor[indices]);
1217  }
1218  }
1219  } // namespace internal
1220 
1221 
1222  template <bool ignore_invalid_symbols,
1223  typename ValueType,
1224  int rank,
1225  int dim,
1226  typename SymbolicType>
1227  void
1229  const Tensor<rank, dim, SymbolicType> &symbol_tensor)
1230  {
1231  // Call the above function
1232  add_to_substitution_map<ignore_invalid_symbols>(
1233  symbol_map, symbol_tensor, Tensor<rank, dim, ValueType>());
1234  }
1235 
1236 
1237  template <bool ignore_invalid_symbols,
1238  typename ValueType,
1239  int rank,
1240  int dim,
1241  typename SymbolicType>
1242  void
1244  types::substitution_map & symbol_map,
1245  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor)
1246  {
1247  // Call the above function
1248  add_to_substitution_map<ignore_invalid_symbols>(
1249  symbol_map, symbol_tensor, SymmetricTensor<rank, dim, ValueType>());
1250  }
1251 
1252 
1253  template <int rank, int dim, typename SymbolicType, typename ValueType>
1254  void
1257  const Tensor<rank, dim, SymbolicType> &symbol_tensor,
1258  const Tensor<rank, dim, ValueType> & value_tensor)
1259  {
1260  internal::set_tensor_value_in_symbol_map(substitution_map,
1261  symbol_tensor,
1262  value_tensor);
1263  }
1264 
1265 
1266  template <int rank, int dim, typename SymbolicType, typename ValueType>
1267  void
1270  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
1271  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1272  {
1273  internal::set_tensor_value_in_symbol_map(substitution_map,
1274  symbol_tensor,
1275  value_tensor);
1276  }
1277 
1278 
1279  /* ------------------ Symbol substitution map creation ----------------*/
1280 
1281 
1282  template <int rank, int dim, typename ExpressionType, typename ValueType>
1285  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1286  const Tensor<rank, dim, ValueType> & value_tensor)
1287  {
1289  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1290  return substitution_map;
1291  }
1292 
1293 
1294  template <int rank, int dim, typename ExpressionType, typename ValueType>
1297  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1298  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1299  {
1301  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1302  return substitution_map;
1303  }
1304 
1305 
1306  /* ---------------- Symbolic substitution map enlargement --------------*/
1307 
1308 
1309  namespace internal
1310  {
1311  template <int rank,
1312  int dim,
1313  typename ExpressionType,
1314  typename ValueType,
1315  template <int, int, typename>
1316  class TensorType>
1317  std::vector<std::pair<ExpressionType, ValueType>>
1318  make_tensor_entries_for_substitution_map(
1319  const TensorType<rank, dim, ExpressionType> &symbol_tensor,
1320  const TensorType<rank, dim, ValueType> & value_tensor)
1321  {
1322  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1323  for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1324  ++i)
1325  {
1326  const TableIndices<rank> indices(
1327  symbol_tensor.unrolled_to_component_indices(i));
1328  symbol_values.push_back(
1329  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1330  }
1331  return symbol_values;
1332  }
1333 
1334 
1335  template <int dim, typename ExpressionType, typename ValueType>
1336  std::vector<std::pair<ExpressionType, ValueType>>
1337  make_tensor_entries_for_substitution_map(
1338  const Tensor<0, dim, ExpressionType> &symbol_tensor,
1339  const Tensor<0, dim, ValueType> & value_tensor)
1340  {
1341  const ExpressionType &expression = symbol_tensor;
1342  const ValueType & value = value_tensor;
1343  return {std::make_pair(expression, value)};
1344  }
1345 
1346 
1347  template <int dim, typename ExpressionType, typename ValueType>
1348  std::vector<std::pair<ExpressionType, ValueType>>
1349  make_tensor_entries_for_substitution_map(
1350  const SymmetricTensor<4, dim, ExpressionType> &symbol_tensor,
1351  const SymmetricTensor<4, dim, ValueType> & value_tensor)
1352  {
1353  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1354  for (unsigned int i = 0;
1355  i < SymmetricTensor<2, dim>::n_independent_components;
1356  ++i)
1357  for (unsigned int j = 0;
1358  j < SymmetricTensor<2, dim>::n_independent_components;
1359  ++j)
1360  {
1361  const TableIndices<4> indices =
1362  make_rank_4_tensor_indices<dim>(i, j);
1363  symbol_values.push_back(
1364  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1365  }
1366  return symbol_values;
1367  }
1368  } // namespace internal
1369 
1370 
1371  template <bool ignore_invalid_symbols,
1372  int rank,
1373  int dim,
1374  typename ExpressionType,
1375  typename ValueType>
1376  void
1379  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1380  const Tensor<rank, dim, ValueType> & value_tensor)
1381  {
1382  add_to_substitution_map<ignore_invalid_symbols>(
1384  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1385  value_tensor));
1386  }
1387 
1388 
1389  template <bool ignore_invalid_symbols,
1390  int rank,
1391  int dim,
1392  typename ExpressionType,
1393  typename ValueType>
1394  void
1397  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1398  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1399  {
1400  add_to_substitution_map<ignore_invalid_symbols>(
1402  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1403  value_tensor));
1404  }
1405 
1406 
1407  /* ---------------- Symbol substitution and evaluation --------------*/
1408 
1409 
1410  namespace internal
1411  {
1412  template <int rank,
1413  int dim,
1414  template <int, int, typename>
1415  class TensorType>
1416  TensorType<rank, dim, Expression>
1417  substitute_tensor(
1418  const TensorType<rank, dim, Expression> &expression_tensor,
1420  {
1421  TensorType<rank, dim, Expression> out;
1422  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1423  {
1424  const TableIndices<rank> indices(
1425  out.unrolled_to_component_indices(i));
1426  out[indices] =
1427  substitute(expression_tensor[indices], substitution_map);
1428  }
1429  return out;
1430  }
1431 
1432 
1433  template <int dim>
1435  substitute_tensor(const Tensor<0, dim, Expression> &expression_tensor,
1437  {
1438  const Expression &expression = expression_tensor;
1439  return substitute(expression, substitution_map);
1440  }
1441 
1442 
1443  template <int dim>
1445  substitute_tensor(
1446  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1448  {
1450  for (unsigned int i = 0;
1451  i < SymmetricTensor<2, dim>::n_independent_components;
1452  ++i)
1453  for (unsigned int j = 0;
1454  j < SymmetricTensor<2, dim>::n_independent_components;
1455  ++j)
1456  {
1457  const TableIndices<4> indices =
1458  make_rank_4_tensor_indices<dim>(i, j);
1459  out[indices] =
1460  substitute(expression_tensor[indices], substitution_map);
1461  }
1462  return out;
1463  }
1464 
1465 
1466  template <typename ValueType,
1467  int rank,
1468  int dim,
1469  template <int, int, typename>
1470  class TensorType>
1471  TensorType<rank, dim, ValueType>
1472  substitute_and_evaluate_tensor(
1473  const TensorType<rank, dim, Expression> &expression_tensor,
1475  {
1476  TensorType<rank, dim, ValueType> out;
1477  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1478  {
1479  const TableIndices<rank> indices(
1480  out.unrolled_to_component_indices(i));
1481  out[indices] =
1482  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1484  }
1485  return out;
1486  }
1487 
1488 
1489  template <typename ValueType, int dim>
1491  substitute_and_evaluate_tensor(
1492  const Tensor<0, dim, Expression> &expression_tensor,
1494  {
1495  const Expression &expression = expression_tensor;
1496  return substitute_and_evaluate<ValueType>(expression, substitution_map);
1497  }
1498 
1499 
1500  template <typename ValueType, int dim>
1502  substitute_and_evaluate_tensor(
1503  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1505  {
1507  for (unsigned int i = 0;
1508  i < SymmetricTensor<2, dim>::n_independent_components;
1509  ++i)
1510  for (unsigned int j = 0;
1511  j < SymmetricTensor<2, dim>::n_independent_components;
1512  ++j)
1513  {
1514  const TableIndices<4> indices =
1515  make_rank_4_tensor_indices<dim>(i, j);
1516  out[indices] =
1517  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1519  }
1520  return out;
1521  }
1522  } // namespace internal
1523 
1524 
1525  template <int rank, int dim>
1527  substitute(const Tensor<rank, dim, Expression> &expression_tensor,
1529  {
1530  return internal::substitute_tensor(expression_tensor, substitution_map);
1531  }
1532 
1533 
1534  template <int rank, int dim>
1536  substitute(const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1538  {
1539  return internal::substitute_tensor(expression_tensor, substitution_map);
1540  }
1541 
1542 
1543  template <typename ValueType, int rank, int dim>
1546  const Tensor<rank, dim, Expression> &expression_tensor,
1548  {
1549  return internal::substitute_and_evaluate_tensor<ValueType>(
1550  expression_tensor, substitution_map);
1551  }
1552 
1553 
1554  template <typename ValueType, int rank, int dim>
1557  const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1559  {
1560  return internal::substitute_and_evaluate_tensor<ValueType>(
1561  expression_tensor, substitution_map);
1562  }
1563 
1564 
1565 
1566  } // namespace SD
1567 } // namespace Differentiation
1568 
1569 # endif // DOXYGEN
1570 
1572 
1573 #endif // DEAL_II_WITH_SYMENGINE
1574 
1575 #endif
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcMessage(std::string arg1)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const SymEngine::RCP< const SymEngine::Basic > &symbol, const SymEngine::RCP< const SymEngine::Basic > &value)
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
Tensor< 1, dim, Expression > make_vector_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
Tensor< rank, dim, Expression > make_tensor_of_symbols(const std::string &symbol)
Tensor< rank, dim, Expression > make_tensor_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
Expression differentiate(const Expression &f, const Expression &x)
ValueType substitute_and_evaluate(const Expression &expression, const types::substitution_map &substitution_map)
SymmetricTensor< rank, dim, Expression > make_symmetric_tensor_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
SymmetricTensor< rank, dim, Expression > make_symmetric_tensor_of_symbols(const std::string &symbol)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
void add_to_symbol_map(types::substitution_map &symbol_map, const Expression &symbol)
Tensor< 1, dim, Expression > make_vector_of_symbols(const std::string &symbol)
void add_to_substitution_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
static const char T