deal.II version GIT relicensing-2206-gaa53ff9447 2024-12-02 09:10:00+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\}}\)
Loading...
Searching...
No Matches
symengine_tensor_operations.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_differentiation_sd_symengine_tensor_operations_h
16#define dealii_differentiation_sd_symengine_tensor_operations_h
17
18#include <deal.II/base/config.h>
19
20#ifdef DEAL_II_WITH_SYMENGINE
21
23# include <deal.II/base/tensor.h>
24
28
29# include <utility>
30# include <vector>
31
33
34namespace Differentiation
35{
36 namespace SD
37 {
61 template <int dim>
63 make_vector_of_symbols(const std::string &symbol);
64
84 template <int rank, int dim>
86 make_tensor_of_symbols(const std::string &symbol);
87
107 template <int rank, int dim>
109 make_symmetric_tensor_of_symbols(const std::string &symbol);
110
130 template <int dim>
132 make_vector_of_symbolic_functions(const std::string &symbol,
133 const types::substitution_map &arguments);
134
155 template <int rank, int dim>
157 make_tensor_of_symbolic_functions(const std::string &symbol,
158 const types::substitution_map &arguments);
159
180 template <int rank, int dim>
183 const std::string &symbol,
184 const types::substitution_map &arguments);
185
202 template <int rank, int dim>
205
215 template <int rank, int dim>
219
229 template <int rank, int dim>
233
243 template <int rank, int dim>
247
257 template <int rank, int dim>
260
271 template <int rank, int dim>
274 const Expression &x);
275
286 template <int rank, int dim>
290
302 template <int rank, int dim>
306
317 template <int rank_1, int rank_2, int dim>
321
334 template <int rank_1, int rank_2, int dim>
338
349 template <int rank_1, int rank_2, int dim>
353
364 template <int rank_1, int rank_2, int dim>
368
400 template <bool ignore_invalid_symbols = false,
401 typename ValueType = double,
402 int rank,
403 int dim,
404 typename SymbolicType>
405 void
407 const Tensor<rank, dim, SymbolicType> &symbol_tensor);
408
433 template <bool ignore_invalid_symbols = false,
434 typename ValueType = double,
435 int rank,
436 int dim,
437 typename SymbolicType>
438 void
440 types::substitution_map &symbol_map,
441 const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor);
442
463 template <int rank, int dim, typename SymbolicType, typename ValueType>
464 void
466 types::substitution_map &substitution_map,
467 const Tensor<rank, dim, SymbolicType> &symbol_tensor,
468 const Tensor<rank, dim, ValueType> &value_tensor);
469
490 template <int rank, int dim, typename SymbolicType, typename ValueType>
491 void
493 types::substitution_map &substitution_map,
494 const SymmetricTensor<rank, dim, SymbolicType> &symbol_tensor,
495 const SymmetricTensor<rank, dim, ValueType> &value_tensor);
496
525 template <int rank, int dim, typename ExpressionType, typename ValueType>
528 const Tensor<rank, dim, ExpressionType> &symbol_tensor,
529 const Tensor<rank, dim, ValueType> &value_tensor);
530
552 template <int rank, int dim, typename ExpressionType, typename ValueType>
556 const SymmetricTensor<rank, dim, ValueType> &value_tensor);
557
591 template <bool ignore_invalid_symbols = false,
592 int rank,
593 int dim,
594 typename ExpressionType,
595 typename ValueType>
596 void
598 types::substitution_map &substitution_map,
599 const Tensor<rank, dim, ExpressionType> &symbol_tensor,
600 const Tensor<rank, dim, ValueType> &value_tensor);
601
627 template <bool ignore_invalid_symbols = false,
628 int rank,
629 int dim,
630 typename ExpressionType,
631 typename ValueType>
632 void
634 types::substitution_map &substitution_map,
636 const SymmetricTensor<rank, dim, ValueType> &value_tensor);
637
664 template <int rank, int dim>
667 const types::substitution_map &substitution_map);
668
688 template <int rank, int dim>
691 const types::substitution_map &substitution_map);
692
719 template <typename ValueType, int rank, int dim>
722 const Tensor<rank, dim, Expression> &expression_tensor,
723 const types::substitution_map &substitution_map);
724
751 template <typename ValueType, int rank, int dim>
754 const SymmetricTensor<rank, dim, Expression> &expression_tensor,
755 const types::substitution_map &substitution_map);
756
759 } // namespace SD
760} // namespace Differentiation
761
762
763/* -------------------- inline and template functions ------------------ */
764
765
766# ifndef DOXYGEN
767
768namespace Differentiation
769{
770 namespace SD
771 {
772 /* ---------------- Symbolic differentiation --------------*/
773
774
775 namespace internal
776 {
777 template <int dim>
779 make_rank_4_tensor_indices(const unsigned int idx_i,
780 const unsigned int idx_j)
781 {
782 const TableIndices<2> indices_i(
784 const TableIndices<2> indices_j(
786 return TableIndices<4>(indices_i[0],
787 indices_i[1],
788 indices_j[0],
789 indices_j[1]);
790 }
791
792
793 template <int rank_1, int rank_2>
795 concatenate_indices(const TableIndices<rank_1> &indices_1,
796 const TableIndices<rank_2> &indices_2)
797 {
799 for (unsigned int i = 0; i < rank_1; ++i)
800 indices_out[i] = indices_1[i];
801 for (unsigned int j = 0; j < rank_2; ++j)
802 indices_out[rank_1 + j] = indices_2[j];
803 return indices_out;
804 }
805
806
807 template <int rank>
809 transpose_indices(const TableIndices<rank> &indices)
810 {
811 return indices;
812 }
813
814
815 template <>
816 inline TableIndices<2>
817 transpose_indices(const TableIndices<2> &indices)
818 {
819 return TableIndices<2>(indices[1], indices[0]);
820 }
821
822
823 template <int rank, int dim, typename ValueType>
824 bool
825 is_symmetric_component(const TableIndices<rank> &,
827 {
828 return false;
829 }
830
831
832 template <int rank, int dim, typename ValueType>
833 bool
834 is_symmetric_component(const TableIndices<rank> &,
836 {
837 static_assert(
838 rank == 0 || rank == 2,
839 "Querying symmetric component for non rank-2 symmetric tensor index is not allowed.");
840 return false;
841 }
842
843
844 template <int dim, typename ValueType>
845 bool
846 is_symmetric_component(const TableIndices<2> &table_indices,
848 {
849 return table_indices[0] != table_indices[1];
850 }
851
852
853 template <int dim,
854 typename ValueType = Expression,
855 template <int, int, typename>
856 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>
869 class TensorType>
870 TensorType<rank, dim, ValueType>
871 scalar_diff_tensor(const ValueType &func,
872 const TensorType<rank, dim, ValueType> &op)
873 {
874 TensorType<rank, dim, ValueType> out;
875 for (unsigned int i = 0; i < out.n_independent_components; ++i)
876 {
877 const TableIndices<rank> indices(
878 out.unrolled_to_component_indices(i));
879 out[indices] = differentiate(func, op[indices]);
880
881 if (is_symmetric_component(indices, op))
882 out[indices] *= 0.5;
883 }
884 return out;
885 }
886
887
888 // Specialization for rank-0 tensor
889 template <int rank,
890 int dim,
891 typename ValueType = Expression,
892 template <int, int, typename>
893 class TensorType>
894 TensorType<rank, dim, ValueType>
895 tensor_diff_tensor(const TensorType<0, dim, ValueType> &func,
896 const TensorType<rank, dim, ValueType> &op)
897 {
898 TensorType<rank, dim, ValueType> out;
899 for (unsigned int i = 0; i < out.n_independent_components; ++i)
900 {
901 const TableIndices<rank> indices(
902 out.unrolled_to_component_indices(i));
903 out[indices] = differentiate(func, op[indices]);
904
905 if (is_symmetric_component(indices, op))
906 out[indices] *= 0.5;
907 }
908 return out;
909 }
910
911
912 template <int rank,
913 int dim,
914 typename ValueType = Expression,
915 template <int, int, typename>
916 class TensorType>
917 TensorType<rank, dim, ValueType>
918 tensor_diff_scalar(const TensorType<rank, dim, ValueType> &funcs,
919 const ValueType &op)
920 {
921 TensorType<rank, dim, ValueType> out;
922 for (unsigned int i = 0; i < out.n_independent_components; ++i)
923 {
924 const TableIndices<rank> indices(
925 out.unrolled_to_component_indices(i));
926 out[indices] = differentiate(funcs[indices], op);
927 }
928 return out;
929 }
930
931
932 // Specialization for rank-0 tensor
933 template <int rank,
934 int dim,
935 typename ValueType = Expression,
936 template <int, int, typename>
937 class TensorType>
938 TensorType<rank, dim, ValueType>
939 tensor_diff_tensor(const TensorType<rank, dim, ValueType> &funcs,
940 const TensorType<0, dim, ValueType> &op)
941 {
942 TensorType<rank, dim, ValueType> out;
943 for (unsigned int i = 0; i < out.n_independent_components; ++i)
944 {
945 const TableIndices<rank> indices(
946 out.unrolled_to_component_indices(i));
947 out[indices] = differentiate(funcs[indices], op);
948 }
949 return out;
950 }
951
952
953 // For either symmetric or normal tensors
954 template <int rank_1,
955 int rank_2,
956 int dim,
957 typename ValueType = Expression,
958 template <int, int, typename>
959 class TensorType>
960 TensorType<rank_1 + rank_2, dim, ValueType>
961 tensor_diff_tensor(const TensorType<rank_1, dim, ValueType> &funcs,
962 const TensorType<rank_2, dim, ValueType> &op)
963 {
964 TensorType<rank_1 + rank_2, dim, ValueType> out;
965 for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
966 {
967 const TableIndices<rank_1> indices_i(
968 funcs.unrolled_to_component_indices(i));
969 for (unsigned int j = 0; j < op.n_independent_components; ++j)
970 {
971 const TableIndices<rank_2> indices_j(
972 op.unrolled_to_component_indices(j));
973 const TableIndices<rank_1 + rank_2> indices_out =
974 concatenate_indices(indices_i, indices_j);
975
976 out[indices_out] =
977 differentiate(funcs[indices_i], op[indices_j]);
978
979 if (is_symmetric_component(indices_j, op))
980 out[indices_out] *= 0.5;
981 }
982 }
983 return out;
984 }
985
986
987 // For mixed symmetric/standard tensors
988 // The return type is always a standard tensor, since we cannot be sure
989 // that any symmetries exist in either the function tensor or the
990 // differential operator.
991 template <int rank_1,
992 int rank_2,
993 int dim,
994 typename ValueType = Expression,
995 template <int, int, typename>
996 class TensorType_1,
997 template <int, int, typename>
998 class TensorType_2>
1000 tensor_diff_tensor(const TensorType_1<rank_1, dim, ValueType> &funcs,
1001 const TensorType_2<rank_2, dim, ValueType> &op)
1002 {
1004 for (unsigned int i = 0; i < funcs.n_independent_components; ++i)
1005 {
1006 const TableIndices<rank_1> indices_i(
1007 funcs.unrolled_to_component_indices(i));
1008 for (unsigned int j = 0; j < op.n_independent_components; ++j)
1009 {
1010 const TableIndices<rank_2> indices_j(
1011 op.unrolled_to_component_indices(j));
1012 const TableIndices<rank_1 + rank_2> indices_out =
1013 concatenate_indices(indices_i, indices_j);
1014
1015 out[indices_out] =
1016 differentiate(funcs[indices_i], op[indices_j]);
1017
1018 if (is_symmetric_component(indices_j, op))
1019 out[indices_out] *= 0.5;
1020
1021 // TODO: Implement for SymmetricTensor<4,dim,...>
1022 if (std::is_same_v<
1023 TensorType_1<rank_1, dim, ValueType>,
1024 SymmetricTensor<2, dim, ValueType>>) // Symmetric function
1025 {
1026 const TableIndices<rank_1 + rank_2> indices_out_t =
1027 concatenate_indices(transpose_indices(indices_i),
1028 indices_j);
1029 out[indices_out_t] = out[indices_out];
1030 }
1031 else if (std::is_same_v<
1032 TensorType_2<rank_2, dim, ValueType>,
1034 // 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,
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>
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>
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#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)