Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
symengine_tensor_operations.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 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 
26 # include <deal.II/differentiation/sd/symengine_number_types.h>
27 # include <deal.II/differentiation/sd/symengine_scalar_operations.h>
28 # include <deal.II/differentiation/sd/symengine_types.h>
29 
30 # include <utility>
31 # include <vector>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 namespace Differentiation
36 {
37  namespace SD
38  {
43 
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 
188 
193 
203  template <int rank, int dim>
205  differentiate(const Expression &f, const Tensor<rank, dim, Expression> &T);
206 
216  template <int rank, int dim>
218  differentiate(const Expression & f,
220 
230  template <int rank, int dim>
234 
244  template <int rank, int dim>
248 
258  template <int rank, int dim>
260  differentiate(const Tensor<rank, dim, Expression> &T, const Expression &x);
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 
371 
376 
401  template <bool ignore_invalid_symbols = false,
402  typename ValueType = double,
403  int rank,
404  int dim,
405  typename SymbolicType>
406  void
407  add_to_symbol_map(types::substitution_map & symbol_map,
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
467  types::substitution_map & substitution_map,
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
494  types::substitution_map & substitution_map,
495  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
496  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
497 
499 
504 
526  template <int rank, int dim, typename ExpressionType, typename ValueType>
527  types::substitution_map
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>
554  types::substitution_map
556  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
557  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
558 
560 
565 
592  template <bool ignore_invalid_symbols = false,
593  int rank,
594  int dim,
595  typename ExpressionType,
596  typename ValueType>
597  void
599  types::substitution_map & substitution_map,
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
635  types::substitution_map & substitution_map,
636  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
637  const SymmetricTensor<rank, dim, ValueType> & value_tensor);
638 
640 
645 
665  template <int rank, int dim>
667  substitute(const Tensor<rank, dim, Expression> &expression_tensor,
668  const types::substitution_map & substitution_map);
669 
689  template <int rank, int dim>
691  substitute(const SymmetricTensor<rank, dim, Expression> &expression_tensor,
692  const types::substitution_map & substitution_map);
693 
720  template <typename ValueType, int rank, int dim>
723  const Tensor<rank, dim, Expression> &expression_tensor,
724  const types::substitution_map & substitution_map);
725 
752  template <typename ValueType, int rank, int dim>
755  const SymmetricTensor<rank, dim, Expression> &expression_tensor,
756  const types::substitution_map & substitution_map);
757 
759 
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> class TensorType>
857  TensorType<0, dim, ValueType>
858  scalar_diff_tensor(const ValueType & func,
859  const TensorType<0, dim, ValueType> &op)
860  {
861  return differentiate(func, op);
862  }
863 
864 
865  template <int rank,
866  int dim,
867  typename ValueType = Expression,
868  template <int, int, typename> class TensorType>
869  TensorType<rank, dim, ValueType>
870  scalar_diff_tensor(const ValueType & func,
871  const TensorType<rank, dim, ValueType> &op)
872  {
873  TensorType<rank, dim, ValueType> out;
874  for (unsigned int i = 0; i < out.n_independent_components; ++i)
875  {
876  const TableIndices<rank> indices(
877  out.unrolled_to_component_indices(i));
878  out[indices] = differentiate(func, op[indices]);
879 
880  if (is_symmetric_component(indices, op))
881  out[indices] *= 0.5;
882  }
883  return out;
884  }
885 
886 
887  // Specialization for rank-0 tensor
888  template <int rank,
889  int dim,
890  typename ValueType = Expression,
891  template <int, int, typename> class TensorType>
892  TensorType<rank, dim, ValueType>
893  tensor_diff_tensor(const TensorType<0, dim, ValueType> & func,
894  const TensorType<rank, dim, ValueType> &op)
895  {
896  TensorType<rank, dim, ValueType> out;
897  for (unsigned int i = 0; i < out.n_independent_components; ++i)
898  {
899  const TableIndices<rank> indices(
900  out.unrolled_to_component_indices(i));
901  out[indices] = differentiate(func, op[indices]);
902 
903  if (is_symmetric_component(indices, op))
904  out[indices] *= 0.5;
905  }
906  return out;
907  }
908 
909 
910  template <int rank,
911  int dim,
912  typename ValueType = Expression,
913  template <int, int, typename> class TensorType>
914  TensorType<rank, dim, ValueType>
915  tensor_diff_scalar(const TensorType<rank, dim, ValueType> &funcs,
916  const ValueType & op)
917  {
918  TensorType<rank, dim, ValueType> out;
919  for (unsigned int i = 0; i < out.n_independent_components; ++i)
920  {
921  const TableIndices<rank> indices(
922  out.unrolled_to_component_indices(i));
923  out[indices] = differentiate(funcs[indices], op);
924  }
925  return out;
926  }
927 
928 
929  // Specialization for rank-0 tensor
930  template <int rank,
931  int dim,
932  typename ValueType = Expression,
933  template <int, int, typename> class TensorType>
934  TensorType<rank, dim, ValueType>
935  tensor_diff_tensor(const TensorType<rank, dim, ValueType> &funcs,
936  const TensorType<0, dim, ValueType> & op)
937  {
938  TensorType<rank, dim, ValueType> out;
939  for (unsigned int i = 0; i < out.n_independent_components; ++i)
940  {
941  const TableIndices<rank> indices(
942  out.unrolled_to_component_indices(i));
943  out[indices] = differentiate(funcs[indices], op);
944  }
945  return out;
946  }
947 
948 
949  // For either symmetric or normal tensors
950  template <int rank_1,
951  int rank_2,
952  int dim,
953  typename ValueType = Expression,
954  template <int, int, typename> class TensorType>
955  TensorType<rank_1 + rank_2, dim, ValueType>
956  tensor_diff_tensor(const TensorType<rank_1, dim, ValueType> &funcs,
957  const TensorType<rank_2, dim, ValueType> &op)
958  {
959  TensorType<rank_1 + rank_2, dim, ValueType> out;
960  for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
961  {
962  const TableIndices<rank_1> indices_i(
963  funcs.unrolled_to_component_indices(i));
964  for (unsigned int j = 0; j < op.n_independent_components; ++j)
965  {
966  const TableIndices<rank_2> indices_j(
967  op.unrolled_to_component_indices(j));
968  const TableIndices<rank_1 + rank_2> indices_out =
969  concatenate_indices(indices_i, indices_j);
970 
971  out[indices_out] =
972  differentiate(funcs[indices_i], op[indices_j]);
973 
974  if (is_symmetric_component(indices_j, op))
975  out[indices_out] *= 0.5;
976  }
977  }
978  return out;
979  }
980 
981 
982  // For mixed symmetric/standard tensors
983  // The return type is always a standard tensor, since we cannot be sure
984  // that any symmetries exist in either the function tensor or the
985  // differential operator.
986  template <int rank_1,
987  int rank_2,
988  int dim,
989  typename ValueType = Expression,
990  template <int, int, typename> class TensorType_1,
991  template <int, int, typename> class TensorType_2>
993  tensor_diff_tensor(const TensorType_1<rank_1, dim, ValueType> &funcs,
994  const TensorType_2<rank_2, dim, ValueType> &op)
995  {
997  for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
998  {
999  const TableIndices<rank_1> indices_i(
1000  funcs.unrolled_to_component_indices(i));
1001  for (unsigned int j = 0; j < op.n_independent_components; ++j)
1002  {
1003  const TableIndices<rank_2> indices_j(
1004  op.unrolled_to_component_indices(j));
1005  const TableIndices<rank_1 + rank_2> indices_out =
1006  concatenate_indices(indices_i, indices_j);
1007 
1008  out[indices_out] =
1009  differentiate(funcs[indices_i], op[indices_j]);
1010 
1011  if (is_symmetric_component(indices_j, op))
1012  out[indices_out] *= 0.5;
1013 
1014  // TODO: Implement for SymmetricTensor<4,dim,...>
1015  if (std::is_same<TensorType_1<rank_1, dim, ValueType>,
1017  value) // Symmetric function
1018  {
1019  const TableIndices<rank_1 + rank_2> indices_out_t =
1020  concatenate_indices(transpose_indices(indices_i),
1021  indices_j);
1022  out[indices_out_t] = out[indices_out];
1023  }
1024  else if (std::is_same<TensorType_2<rank_2, dim, ValueType>,
1026  value) // Symmetric operator
1027  {
1028  const TableIndices<rank_1 + rank_2> indices_out_t =
1029  concatenate_indices(indices_i,
1030  transpose_indices(indices_j));
1031  out[indices_out_t] = out[indices_out];
1032  }
1033  else
1034  {
1035  Assert(
1036  false,
1037  ExcMessage(
1038  "Expect mixed tensor differentiation to have at least "
1039  "one component stemming from a symmetric tensor."));
1040  }
1041  }
1042  }
1043  return out;
1044  }
1045 
1046  } // namespace internal
1047 
1048 
1049  template <int rank, int dim>
1051  differentiate(const Expression & func,
1053  {
1054  return internal::scalar_diff_tensor(func, op);
1055  }
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>
1071  {
1072  return internal::scalar_diff_tensor(func, op);
1073  }
1074 
1075 
1076  template <int rank, int dim>
1080  {
1081  // Ensure that this returns a symmetric tensor by
1082  // invoking the scalar value function
1083  const Expression tmp = func;
1084  return internal::scalar_diff_tensor(tmp, op);
1085  }
1086 
1087 
1088  template <int rank, int dim>
1090  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1091  const Expression & op)
1092  {
1093  return internal::tensor_diff_scalar(symbol_tensor, op);
1094  }
1095 
1096 
1097  template <int rank, int dim>
1099  differentiate(const Tensor<rank, dim, Expression> &symbol_tensor,
1100  const Tensor<0, dim, Expression> & op)
1101  {
1102  return internal::tensor_diff_scalar(symbol_tensor, op);
1103  }
1104 
1105 
1106  template <int rank, int dim>
1109  const Expression & op)
1110  {
1111  return internal::tensor_diff_scalar(symbol_tensor, op);
1112  }
1113 
1114 
1115  template <int rank, int dim>
1118  const Tensor<0, dim, Expression> & op)
1119  {
1120  return internal::tensor_diff_scalar(symbol_tensor, op);
1121  }
1122 
1123 
1124  template <int rank_1, int rank_2, int dim>
1126  differentiate(const Tensor<rank_1, dim, Expression> &symbol_tensor,
1128  {
1129  return internal::tensor_diff_tensor(symbol_tensor, op);
1130  }
1131 
1132 
1133  template <int rank_1, int rank_2, int dim>
1137  {
1138  return internal::tensor_diff_tensor(symbol_tensor, op);
1139  }
1140 
1141 
1142  template <int rank_1, int rank_2, int dim>
1144  differentiate(const Tensor<rank_1, dim, Expression> & symbol_tensor,
1146  {
1147  return internal::tensor_diff_tensor(symbol_tensor, op);
1148  }
1149 
1150 
1151  template <int rank_1, int rank_2, int dim>
1155  {
1156  return internal::tensor_diff_tensor(symbol_tensor, op);
1157  }
1158 
1159 
1160  /* ---------------- Symbol map creation and manipulation --------------*/
1161 
1162 
1163  namespace internal
1164  {
1165  template <typename SymbolicType,
1166  typename ValueType,
1167  int rank,
1168  int dim,
1169  template <int, int, typename> class TensorType>
1170  void
1171  set_tensor_value_in_symbol_map(
1172  types::substitution_map & substitution_map,
1173  const TensorType<rank, dim, SymbolicType> &symbol_tensor,
1174  const TensorType<rank, dim, ValueType> & value_tensor)
1175  {
1176  TensorType<rank, dim, Expression> out;
1177  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1178  {
1179  const TableIndices<rank> indices(
1180  out.unrolled_to_component_indices(i));
1181  set_value_in_symbol_map(substitution_map,
1182  symbol_tensor[indices],
1183  value_tensor[indices]);
1184  }
1185  }
1186 
1187 
1188  template <typename SymbolicType, typename ValueType, int dim>
1189  void
1190  set_tensor_value_in_symbol_map(
1191  types::substitution_map & substitution_map,
1192  const SymmetricTensor<4, dim, SymbolicType> &symbol_tensor,
1193  const SymmetricTensor<4, dim, ValueType> & value_tensor)
1194  {
1196  for (unsigned int i = 0;
1197  i < SymmetricTensor<2, dim>::n_independent_components;
1198  ++i)
1199  for (unsigned int j = 0;
1200  j < SymmetricTensor<2, dim>::n_independent_components;
1201  ++j)
1202  {
1203  const TableIndices<4> indices =
1204  make_rank_4_tensor_indices<dim>(i, j);
1205  set_value_in_symbol_map(substitution_map,
1206  symbol_tensor[indices],
1207  value_tensor[indices]);
1208  }
1209  }
1210  } // namespace internal
1211 
1212 
1213  template <bool ignore_invalid_symbols,
1214  typename ValueType,
1215  int rank,
1216  int dim,
1217  typename SymbolicType>
1218  void
1219  add_to_symbol_map(types::substitution_map & symbol_map,
1220  const Tensor<rank, dim, SymbolicType> &symbol_tensor)
1221  {
1222  // Call the above function
1223  add_to_substitution_map<ignore_invalid_symbols>(
1224  symbol_map, symbol_tensor, Tensor<rank, dim, ValueType>());
1225  }
1226 
1227 
1228  template <bool ignore_invalid_symbols,
1229  typename ValueType,
1230  int rank,
1231  int dim,
1232  typename SymbolicType>
1233  void
1235  types::substitution_map & symbol_map,
1236  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor)
1237  {
1238  // Call the above function
1239  add_to_substitution_map<ignore_invalid_symbols>(
1240  symbol_map, symbol_tensor, SymmetricTensor<rank, dim, ValueType>());
1241  }
1242 
1243 
1244  template <int rank, int dim, typename SymbolicType, typename ValueType>
1245  void
1247  types::substitution_map & substitution_map,
1248  const Tensor<rank, dim, SymbolicType> &symbol_tensor,
1249  const Tensor<rank, dim, ValueType> & value_tensor)
1250  {
1251  internal::set_tensor_value_in_symbol_map(substitution_map,
1252  symbol_tensor,
1253  value_tensor);
1254  }
1255 
1256 
1257  template <int rank, int dim, typename SymbolicType, typename ValueType>
1258  void
1260  types::substitution_map & substitution_map,
1261  const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
1262  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1263  {
1264  internal::set_tensor_value_in_symbol_map(substitution_map,
1265  symbol_tensor,
1266  value_tensor);
1267  }
1268 
1269 
1270  /* ------------------ Symbol substitution map creation ----------------*/
1271 
1272 
1273  template <int rank, int dim, typename ExpressionType, typename ValueType>
1274  types::substitution_map
1276  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1277  const Tensor<rank, dim, ValueType> & value_tensor)
1278  {
1279  types::substitution_map substitution_map;
1280  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1281  return substitution_map;
1282  }
1283 
1284 
1285  template <int rank, int dim, typename ExpressionType, typename ValueType>
1286  types::substitution_map
1288  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1289  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1290  {
1291  types::substitution_map substitution_map;
1292  add_to_substitution_map(substitution_map, symbol_tensor, value_tensor);
1293  return substitution_map;
1294  }
1295 
1296 
1297  /* ---------------- Symbolic substitution map enlargement --------------*/
1298 
1299 
1300  namespace internal
1301  {
1302  template <int rank,
1303  int dim,
1304  typename ExpressionType,
1305  typename ValueType,
1306  template <int, int, typename> class TensorType>
1307  std::vector<std::pair<ExpressionType, ValueType>>
1308  make_tensor_entries_for_substitution_map(
1309  const TensorType<rank, dim, ExpressionType> &symbol_tensor,
1310  const TensorType<rank, dim, ValueType> & value_tensor)
1311  {
1312  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1313  for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1314  ++i)
1315  {
1316  const TableIndices<rank> indices(
1317  symbol_tensor.unrolled_to_component_indices(i));
1318  symbol_values.push_back(
1319  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1320  }
1321  return symbol_values;
1322  }
1323 
1324 
1325  template <int dim, typename ExpressionType, typename ValueType>
1326  std::vector<std::pair<ExpressionType, ValueType>>
1327  make_tensor_entries_for_substitution_map(
1328  const Tensor<0, dim, ExpressionType> &symbol_tensor,
1329  const Tensor<0, dim, ValueType> & value_tensor)
1330  {
1331  const ExpressionType &expression = symbol_tensor;
1332  const ValueType & value = value_tensor;
1333  return {std::make_pair(expression, value)};
1334  }
1335 
1336 
1337  template <int dim, typename ExpressionType, typename ValueType>
1338  std::vector<std::pair<ExpressionType, ValueType>>
1339  make_tensor_entries_for_substitution_map(
1340  const SymmetricTensor<4, dim, ExpressionType> &symbol_tensor,
1341  const SymmetricTensor<4, dim, ValueType> & value_tensor)
1342  {
1343  std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1344  for (unsigned int i = 0;
1345  i < SymmetricTensor<2, dim>::n_independent_components;
1346  ++i)
1347  for (unsigned int j = 0;
1348  j < SymmetricTensor<2, dim>::n_independent_components;
1349  ++j)
1350  {
1351  const TableIndices<4> indices =
1352  make_rank_4_tensor_indices<dim>(i, j);
1353  symbol_values.push_back(
1354  std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1355  }
1356  return symbol_values;
1357  }
1358  } // namespace internal
1359 
1360 
1361  template <bool ignore_invalid_symbols,
1362  int rank,
1363  int dim,
1364  typename ExpressionType,
1365  typename ValueType>
1366  void
1368  types::substitution_map & substitution_map,
1369  const Tensor<rank, dim, ExpressionType> &symbol_tensor,
1370  const Tensor<rank, dim, ValueType> & value_tensor)
1371  {
1372  add_to_substitution_map<ignore_invalid_symbols>(
1373  substitution_map,
1374  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1375  value_tensor));
1376  }
1377 
1378 
1379  template <bool ignore_invalid_symbols,
1380  int rank,
1381  int dim,
1382  typename ExpressionType,
1383  typename ValueType>
1384  void
1386  types::substitution_map & substitution_map,
1387  const SymmetricTensor<rank, dim, ExpressionType> &symbol_tensor,
1388  const SymmetricTensor<rank, dim, ValueType> & value_tensor)
1389  {
1390  add_to_substitution_map<ignore_invalid_symbols>(
1391  substitution_map,
1392  internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1393  value_tensor));
1394  }
1395 
1396 
1397  /* ---------------- Symbol substitution and evaluation --------------*/
1398 
1399 
1400  namespace internal
1401  {
1402  template <int rank,
1403  int dim,
1404  template <int, int, typename> class TensorType>
1405  TensorType<rank, dim, Expression>
1406  substitute_tensor(
1407  const TensorType<rank, dim, Expression> &expression_tensor,
1408  const types::substitution_map & substitution_map)
1409  {
1410  TensorType<rank, dim, Expression> out;
1411  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1412  {
1413  const TableIndices<rank> indices(
1414  out.unrolled_to_component_indices(i));
1415  out[indices] =
1416  substitute(expression_tensor[indices], substitution_map);
1417  }
1418  return out;
1419  }
1420 
1421 
1422  template <int dim>
1424  substitute_tensor(const Tensor<0, dim, Expression> &expression_tensor,
1425  const types::substitution_map & substitution_map)
1426  {
1427  const Expression &expression = expression_tensor;
1428  return substitute(expression, substitution_map);
1429  }
1430 
1431 
1432  template <int dim>
1434  substitute_tensor(
1435  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1436  const types::substitution_map & substitution_map)
1437  {
1439  for (unsigned int i = 0;
1440  i < SymmetricTensor<2, dim>::n_independent_components;
1441  ++i)
1442  for (unsigned int j = 0;
1443  j < SymmetricTensor<2, dim>::n_independent_components;
1444  ++j)
1445  {
1446  const TableIndices<4> indices =
1447  make_rank_4_tensor_indices<dim>(i, j);
1448  out[indices] =
1449  substitute(expression_tensor[indices], substitution_map);
1450  }
1451  return out;
1452  }
1453 
1454 
1455  template <typename ValueType,
1456  int rank,
1457  int dim,
1458  template <int, int, typename> class TensorType>
1459  TensorType<rank, dim, ValueType>
1460  substitute_and_evaluate_tensor(
1461  const TensorType<rank, dim, Expression> &expression_tensor,
1462  const types::substitution_map & substitution_map)
1463  {
1464  TensorType<rank, dim, ValueType> out;
1465  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1466  {
1467  const TableIndices<rank> indices(
1468  out.unrolled_to_component_indices(i));
1469  out[indices] =
1470  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1471  substitution_map);
1472  }
1473  return out;
1474  }
1475 
1476 
1477  template <typename ValueType, int dim>
1479  substitute_and_evaluate_tensor(
1480  const Tensor<0, dim, Expression> &expression_tensor,
1481  const types::substitution_map & substitution_map)
1482  {
1483  const Expression &expression = expression_tensor;
1484  return substitute_and_evaluate<ValueType>(expression, substitution_map);
1485  }
1486 
1487 
1488  template <typename ValueType, int dim>
1490  substitute_and_evaluate_tensor(
1491  const SymmetricTensor<4, dim, Expression> &expression_tensor,
1492  const types::substitution_map & substitution_map)
1493  {
1495  for (unsigned int i = 0;
1496  i < SymmetricTensor<2, dim>::n_independent_components;
1497  ++i)
1498  for (unsigned int j = 0;
1499  j < SymmetricTensor<2, dim>::n_independent_components;
1500  ++j)
1501  {
1502  const TableIndices<4> indices =
1503  make_rank_4_tensor_indices<dim>(i, j);
1504  out[indices] =
1505  substitute_and_evaluate<ValueType>(expression_tensor[indices],
1506  substitution_map);
1507  }
1508  return out;
1509  }
1510  } // namespace internal
1511 
1512 
1513  template <int rank, int dim>
1515  substitute(const Tensor<rank, dim, Expression> &expression_tensor,
1516  const types::substitution_map & substitution_map)
1517  {
1518  return internal::substitute_tensor(expression_tensor, substitution_map);
1519  }
1520 
1521 
1522  template <int rank, int dim>
1524  substitute(const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1525  const types::substitution_map & substitution_map)
1526  {
1527  return internal::substitute_tensor(expression_tensor, substitution_map);
1528  }
1529 
1530 
1531  template <typename ValueType, int rank, int dim>
1534  const Tensor<rank, dim, Expression> &expression_tensor,
1535  const types::substitution_map & substitution_map)
1536  {
1537  return internal::substitute_and_evaluate_tensor<ValueType>(
1538  expression_tensor, substitution_map);
1539  }
1540 
1541 
1542  template <typename ValueType, int rank, int dim>
1545  const SymmetricTensor<rank, dim, Expression> &expression_tensor,
1546  const types::substitution_map & substitution_map)
1547  {
1548  return internal::substitute_and_evaluate_tensor<ValueType>(
1549  expression_tensor, substitution_map);
1550  }
1551 
1552 
1553 
1554  } // namespace SD
1555 } // namespace Differentiation
1556 
1557 # endif // DOXYGEN
1558 
1559 DEAL_II_NAMESPACE_CLOSE
1560 
1561 #endif // DEAL_II_WITH_SYMENGINE
1562 
1563 #endif
SymmetricTensor< rank, dim, Expression > make_symmetric_tensor_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
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_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
Expression differentiate(const Expression &f, const Expression &x)
void add_to_substitution_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
ValueType substitute_and_evaluate(const Expression &expression, const types::substitution_map &substitution_map)
void add_to_symbol_map(types::substitution_map &symbol_map, const Expression &symbol)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
Definition: mpi.h:90
Tensor< 1, dim, Expression > make_vector_of_symbols(const std::string &symbol)
SymmetricTensor< rank, dim, Expression > make_symmetric_tensor_of_symbols(const std::string &symbol)
Tensor< rank, dim, Expression > make_tensor_of_symbols(const std::string &symbol)