Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
35namespace 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>
291
303 template <int rank, int dim>
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
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
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>
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
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,
637 const SymmetricTensor<rank, dim, ValueType> & value_tensor);
638
665 template <int rank, int dim>
668 const types::substitution_map & substitution_map);
669
689 template <int rank, int dim>
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
760 } // namespace SD
761} // namespace Differentiation
762
763
764/* -------------------- inline and template functions ------------------ */
765
766
767# ifndef DOXYGEN
768
769namespace 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 {
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>
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(
1181 types::substitution_map & substitution_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));
1190 set_value_in_symbol_map(substitution_map,
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(
1200 types::substitution_map & substitution_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);
1214 set_value_in_symbol_map(substitution_map,
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
1256 types::substitution_map & substitution_map,
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
1269 types::substitution_map & substitution_map,
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
1378 types::substitution_map & substitution_map,
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>(
1383 substitution_map,
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
1396 types::substitution_map & substitution_map,
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>(
1401 substitution_map,
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,
1419 const types::substitution_map & substitution_map)
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,
1436 const types::substitution_map & substitution_map)
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,
1447 const types::substitution_map & substitution_map)
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,
1474 const types::substitution_map & substitution_map)
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],
1483 substitution_map);
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,
1493 const types::substitution_map & substitution_map)
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,
1504 const types::substitution_map & substitution_map)
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],
1518 substitution_map);
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,
1528 const types::substitution_map & substitution_map)
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,
1537 const types::substitution_map & substitution_map)
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,
1547 const types::substitution_map & substitution_map)
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,
1558 const types::substitution_map & substitution_map)
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
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< rank, dim, Expression > make_tensor_of_symbols(const std::string &symbol)
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_symbols(const std::string &symbol)
Tensor< 1, dim, Expression > make_vector_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)
SymmetricTensor< rank, dim, Expression > make_symmetric_tensor_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
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)
Tensor< rank, dim, Expression > make_tensor_of_symbolic_functions(const std::string &symbol, const types::substitution_map &arguments)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)