Reference documentation for deal.II version GIT 8ad845fa53 2023-12-10 20:50: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 - 2022 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_v<
1024  TensorType_1<rank_1, dim, ValueType>,
1025  SymmetricTensor<2, dim, ValueType>>) // 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_v<
1033  TensorType_2<rank_2, dim, ValueType>,
1034  SymmetricTensor<2, dim, ValueType>>) // Symmetric
1035  // operator
1036  {
1037  const TableIndices<rank_1 + rank_2> indices_out_t =
1038  concatenate_indices(indices_i,
1039  transpose_indices(indices_j));
1040  out[indices_out_t] = out[indices_out];
1041  }
1042  else
1043  {
1044  Assert(
1045  false,
1046  ExcMessage(
1047  "Expect mixed tensor differentiation to have at least "
1048  "one component stemming from a symmetric tensor."));
1049  }
1050  }
1051  }
1052  return out;
1053  }
1054 
1055  } // namespace internal
1056 
1057 
1058  template <int rank, int dim>
1060  differentiate(const Expression &func,
1062  {
1063  return internal::scalar_diff_tensor(func, op);
1064  }
1065 
1066 
1067  template <int rank, int dim>
1069  differentiate(const Expression &func,
1071  {
1072  return internal::scalar_diff_tensor(func, op);
1073  }
1074 
1075 
1076  template <int rank, int dim>
1080  {
1081  return internal::scalar_diff_tensor(func, op);
1082  }
1083 
1084 
1085  template <int rank, int dim>
1089  {
1090  // Ensure that this returns a symmetric tensor by
1091  // invoking the scalar value function
1092  const Expression tmp = func;
1093  return internal::scalar_diff_tensor(tmp, op);
1094  }
1095 
1096 
1097  template <int rank, int dim>
1099  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1100  const Expression &op)
1101  {
1102  return internal::tensor_diff_scalar(symbol_tensor, op);
1103  }
1104 
1105 
1106  template <int rank, int dim>
1108  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1109  const Tensor<0, dim, Expression> &op)
1110  {
1111  return internal::tensor_diff_scalar(symbol_tensor, op);
1112  }
1113 
1114 
1115  template <int rank, int dim>
1118  const Expression &op)
1119  {
1120  return internal::tensor_diff_scalar(symbol_tensor, op);
1121  }
1122 
1123 
1124  template <int rank, int dim>
1127  const Tensor<0, dim, Expression> &op)
1128  {
1129  return internal::tensor_diff_scalar(symbol_tensor, op);
1130  }
1131 
1132 
1133  template <int rank_1, int rank_2, int dim>
1135  differentiate(const Tensor<rank_1, dim, Expression> &symbol_tensor,
1137  {
1138  return internal::tensor_diff_tensor(symbol_tensor, op);
1139  }
1140 
1141 
1142  template <int rank_1, int rank_2, int dim>
1146  {
1147  return internal::tensor_diff_tensor(symbol_tensor, op);
1148  }
1149 
1150 
1151  template <int rank_1, int rank_2, int dim>
1153  differentiate(const Tensor<rank_1, dim, Expression> &symbol_tensor,
1155  {
1156  return internal::tensor_diff_tensor(symbol_tensor, op);
1157  }
1158 
1159 
1160  template <int rank_1, int rank_2, int dim>
1164  {
1165  return internal::tensor_diff_tensor(symbol_tensor, op);
1166  }
1167 
1168 
1169  /* ---------------- Symbol map creation and manipulation --------------*/
1170 
1171 
1172  namespace internal
1173  {
1174  template <typename SymbolicType,
1175  typename ValueType,
1176  int rank,
1177  int dim,
1178  template <int, int, typename>
1179  class TensorType>
1180  void
1181  set_tensor_value_in_symbol_map(
1183  const TensorType<rank, dim, SymbolicType> &symbol_tensor,
1184  const TensorType<rank, dim, ValueType> &value_tensor)
1185  {
1186  TensorType<rank, dim, Expression> out;
1187  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1188  {
1189  const TableIndices<rank> indices(
1190  out.unrolled_to_component_indices(i));
1192  symbol_tensor[indices],
1193  value_tensor[indices]);
1194  }
1195  }
1196 
1197 
1198  template <typename SymbolicType, typename ValueType, int dim>
1199  void
1200  set_tensor_value_in_symbol_map(
1202  const SymmetricTensor<4, dim, SymbolicType> &symbol_tensor,
1203  const SymmetricTensor<4, dim, ValueType> &value_tensor)
1204  {
1206  for (unsigned int i = 0;
1207  i < SymmetricTensor<2, dim>::n_independent_components;
1208  ++i)
1209  for (unsigned int j = 0;
1210  j < SymmetricTensor<2, dim>::n_independent_components;
1211  ++j)
1212  {
1213  const TableIndices<4> indices =
1214  make_rank_4_tensor_indices<dim>(i, j);
1216  symbol_tensor[indices],
1217  value_tensor[indices]);
1218  }
1219  }
1220  } // namespace internal
1221 
1222 
1223  template <bool ignore_invalid_symbols,
1224  typename ValueType,
1225  int rank,
1226  int dim,
1227  typename SymbolicType>
1228  void
1230  const Tensor<rank, dim, SymbolicType> &symbol_tensor)
1231  {
1232  // Call the above function
1233  add_to_substitution_map<ignore_invalid_symbols>(
1234  symbol_map, symbol_tensor, Tensor<rank, dim, ValueType>());
1235  }
1236 
1237 
1238  template <bool ignore_invalid_symbols,
1239  typename ValueType,
1240  int rank,
1241  int dim,
1242  typename SymbolicType>
1243  void
1245  types::substitution_map &symbol_map,
1246  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor)
1247  {
1248  // Call the above function
1249  add_to_substitution_map<ignore_invalid_symbols>(
1250  symbol_map, symbol_tensor, SymmetricTensor<rank, dim, ValueType>());
1251  }
1252 
1253 
1254  template <int rank, int dim, typename SymbolicType, typename ValueType>
1255  void
1258  const Tensor<rank, dim, SymbolicType> &symbol_tensor,
1259  const Tensor<rank, dim, ValueType> &value_tensor)
1260  {
1261  internal::set_tensor_value_in_symbol_map(substitution_map,
1262  symbol_tensor,
1263  value_tensor);
1264  }
1265 
1266 
1267  template <int rank, int dim, typename SymbolicType, typename ValueType>
1268  void
1271  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
1272  const SymmetricTensor<rank, dim, ValueType> &value_tensor)
1273  {
1274  internal::set_tensor_value_in_symbol_map(substitution_map,
1275  symbol_tensor,
1276  value_tensor);
1277  }
1278 
1279 
1280  /* ------------------ Symbol substitution map creation ----------------*/
1281 
1282 
1283  template <int rank, int dim, typename ExpressionType, typename ValueType>
1286  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1287  const Tensor<rank, dim, ValueType> &value_tensor)
1288  {
1290  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1291  return substitution_map;
1292  }
1293 
1294 
1295  template <int rank, int dim, typename ExpressionType, typename ValueType>
1298  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1299  const SymmetricTensor<rank, dim, ValueType> &value_tensor)
1300  {
1302  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1303  return substitution_map;
1304  }
1305 
1306 
1307  /* ---------------- Symbolic substitution map enlargement --------------*/
1308 
1309 
1310  namespace internal
1311  {
1312  template <int rank,
1313  int dim,
1314  typename ExpressionType,
1315  typename ValueType,
1316  template <int, int, typename>
1317  class TensorType>
1318  std::vector<std::pair<ExpressionType, ValueType>>
1319  make_tensor_entries_for_substitution_map(
1320  const TensorType<rank, dim, ExpressionType> &symbol_tensor,
1321  const TensorType<rank, dim, ValueType> &value_tensor)
1322  {
1323  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1324  for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1325  ++i)
1326  {
1327  const TableIndices<rank> indices(
1328  symbol_tensor.unrolled_to_component_indices(i));
1329  symbol_values.push_back(
1330  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1331  }
1332  return symbol_values;
1333  }
1334 
1335 
1336  template <int dim, typename ExpressionType, typename ValueType>
1337  std::vector<std::pair<ExpressionType, ValueType>>
1338  make_tensor_entries_for_substitution_map(
1339  const Tensor<0, dim, ExpressionType> &symbol_tensor,
1340  const Tensor<0, dim, ValueType> &value_tensor)
1341  {
1342  const ExpressionType &expression = symbol_tensor;
1343  const ValueType &value = value_tensor;
1344  return {std::make_pair(expression, value)};
1345  }
1346 
1347 
1348  template <int dim, typename ExpressionType, typename ValueType>
1349  std::vector<std::pair<ExpressionType, ValueType>>
1350  make_tensor_entries_for_substitution_map(
1351  const SymmetricTensor<4, dim, ExpressionType> &symbol_tensor,
1352  const SymmetricTensor<4, dim, ValueType> &value_tensor)
1353  {
1354  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1355  for (unsigned int i = 0;
1356  i < SymmetricTensor<2, dim>::n_independent_components;
1357  ++i)
1358  for (unsigned int j = 0;
1359  j < SymmetricTensor<2, dim>::n_independent_components;
1360  ++j)
1361  {
1362  const TableIndices<4> indices =
1363  make_rank_4_tensor_indices<dim>(i, j);
1364  symbol_values.push_back(
1365  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1366  }
1367  return symbol_values;
1368  }
1369  } // namespace internal
1370 
1371 
1372  template <bool ignore_invalid_symbols,
1373  int rank,
1374  int dim,
1375  typename ExpressionType,
1376  typename ValueType>
1377  void
1380  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1381  const Tensor<rank, dim, ValueType> &value_tensor)
1382  {
1383  add_to_substitution_map<ignore_invalid_symbols>(
1385  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1386  value_tensor));
1387  }
1388 
1389 
1390  template <bool ignore_invalid_symbols,
1391  int rank,
1392  int dim,
1393  typename ExpressionType,
1394  typename ValueType>
1395  void
1398  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1399  const SymmetricTensor<rank, dim, ValueType> &value_tensor)
1400  {
1401  add_to_substitution_map<ignore_invalid_symbols>(
1403  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1404  value_tensor));
1405  }
1406 
1407 
1408  /* ---------------- Symbol substitution and evaluation --------------*/
1409 
1410 
1411  namespace internal
1412  {
1413  template <int rank,
1414  int dim,
1415  template <int, int, typename>
1416  class TensorType>
1417  TensorType<rank, dim, Expression>
1418  substitute_tensor(
1419  const TensorType<rank, dim, Expression> &expression_tensor,
1421  {
1422  TensorType<rank, dim, Expression> out;
1423  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1424  {
1425  const TableIndices<rank> indices(
1426  out.unrolled_to_component_indices(i));
1427  out[indices] =
1428  substitute(expression_tensor[indices], substitution_map);
1429  }
1430  return out;
1431  }
1432 
1433 
1434  template <int dim>
1436  substitute_tensor(const Tensor<0, dim, Expression> &expression_tensor,
1438  {
1439  const Expression &expression = expression_tensor;
1440  return substitute(expression, substitution_map);
1441  }
1442 
1443 
1444  template <int dim>
1446  substitute_tensor(
1447  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1449  {
1451  for (unsigned int i = 0;
1452  i < SymmetricTensor<2, dim>::n_independent_components;
1453  ++i)
1454  for (unsigned int j = 0;
1455  j < SymmetricTensor<2, dim>::n_independent_components;
1456  ++j)
1457  {
1458  const TableIndices<4> indices =
1459  make_rank_4_tensor_indices<dim>(i, j);
1460  out[indices] =
1461  substitute(expression_tensor[indices], substitution_map);
1462  }
1463  return out;
1464  }
1465 
1466 
1467  template <typename ValueType,
1468  int rank,
1469  int dim,
1470  template <int, int, typename>
1471  class TensorType>
1472  TensorType<rank, dim, ValueType>
1473  substitute_and_evaluate_tensor(
1474  const TensorType<rank, dim, Expression> &expression_tensor,
1476  {
1477  TensorType<rank, dim, ValueType> out;
1478  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1479  {
1480  const TableIndices<rank> indices(
1481  out.unrolled_to_component_indices(i));
1482  out[indices] =
1483  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1485  }
1486  return out;
1487  }
1488 
1489 
1490  template <typename ValueType, int dim>
1492  substitute_and_evaluate_tensor(
1493  const Tensor<0, dim, Expression> &expression_tensor,
1495  {
1496  const Expression &expression = expression_tensor;
1497  return substitute_and_evaluate<ValueType>(expression, substitution_map);
1498  }
1499 
1500 
1501  template <typename ValueType, int dim>
1503  substitute_and_evaluate_tensor(
1504  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1506  {
1508  for (unsigned int i = 0;
1509  i < SymmetricTensor<2, dim>::n_independent_components;
1510  ++i)
1511  for (unsigned int j = 0;
1512  j < SymmetricTensor<2, dim>::n_independent_components;
1513  ++j)
1514  {
1515  const TableIndices<4> indices =
1516  make_rank_4_tensor_indices<dim>(i, j);
1517  out[indices] =
1518  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1520  }
1521  return out;
1522  }
1523  } // namespace internal
1524 
1525 
1526  template <int rank, int dim>
1528  substitute(const Tensor<rank, dim, Expression> &expression_tensor,
1530  {
1531  return internal::substitute_tensor(expression_tensor, substitution_map);
1532  }
1533 
1534 
1535  template <int rank, int dim>
1537  substitute(const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1539  {
1540  return internal::substitute_tensor(expression_tensor, substitution_map);
1541  }
1542 
1543 
1544  template <typename ValueType, int rank, int dim>
1547  const Tensor<rank, dim, Expression> &expression_tensor,
1549  {
1550  return internal::substitute_and_evaluate_tensor<ValueType>(
1551  expression_tensor, substitution_map);
1552  }
1553 
1554 
1555  template <typename ValueType, int rank, int dim>
1558  const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1560  {
1561  return internal::substitute_and_evaluate_tensor<ValueType>(
1562  expression_tensor, substitution_map);
1563  }
1564 
1565 
1566 
1567  } // namespace SD
1568 } // namespace Differentiation
1569 
1570 # endif // DOXYGEN
1571 
1573 
1574 #endif // DEAL_II_WITH_SYMENGINE
1575 
1576 #endif
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1631
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