16 #ifndef dealii_differentiation_sd_symengine_tensor_operations_h
17 #define dealii_differentiation_sd_symengine_tensor_operations_h
21 #ifdef DEAL_II_WITH_SYMENGINE
85 template <
int rank,
int dim>
108 template <
int rank,
int dim>
156 template <
int rank,
int dim>
181 template <
int rank,
int dim>
184 const std::string &symbol,
203 template <
int rank,
int dim>
216 template <
int rank,
int dim>
230 template <
int rank,
int dim>
244 template <
int rank,
int dim>
258 template <
int rank,
int dim>
272 template <
int rank,
int dim>
287 template <
int rank,
int dim>
303 template <
int rank,
int dim>
318 template <
int rank_1,
int rank_2,
int dim>
335 template <
int rank_1,
int rank_2,
int dim>
350 template <
int rank_1,
int rank_2,
int dim>
365 template <
int rank_1,
int rank_2,
int dim>
401 template <
bool ignore_invalid_symbols =
false,
402 typename ValueType =
double,
405 typename SymbolicType>
434 template <
bool ignore_invalid_symbols =
false,
435 typename ValueType =
double,
438 typename SymbolicType>
464 template <
int rank,
int dim,
typename SymbolicType,
typename ValueType>
491 template <
int rank,
int dim,
typename SymbolicType,
typename ValueType>
526 template <
int rank,
int dim,
typename ExpressionType,
typename ValueType>
553 template <
int rank,
int dim,
typename ExpressionType,
typename ValueType>
592 template <
bool ignore_invalid_symbols =
false,
595 typename ExpressionType,
628 template <
bool ignore_invalid_symbols =
false,
631 typename ExpressionType,
665 template <
int rank,
int dim>
689 template <
int rank,
int dim>
720 template <
typename ValueType,
int rank,
int dim>
752 template <
typename ValueType,
int rank,
int dim>
780 make_rank_4_tensor_indices(
const unsigned int idx_i,
781 const unsigned int idx_j)
794 template <
int rank_1,
int rank_2>
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];
824 template <
int rank,
int dim,
typename ValueType>
833 template <
int rank,
int dim,
typename ValueType>
839 rank == 0 || rank == 2,
840 "Querying symmetric component for non rank-2 symmetric tensor index is not allowed.");
845 template <
int dim,
typename ValueType>
850 return table_indices[0] != table_indices[1];
855 typename ValueType = Expression,
856 template <
int,
int,
typename>
858 TensorType<0, dim, ValueType>
859 scalar_diff_tensor(
const ValueType &func,
860 const TensorType<0, dim, ValueType> &op)
868 typename ValueType = Expression,
869 template <
int,
int,
typename>
871 TensorType<rank, dim, ValueType>
872 scalar_diff_tensor(
const ValueType &func,
873 const TensorType<rank, dim, ValueType> &op)
875 TensorType<rank, dim, ValueType> out;
876 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
879 out.unrolled_to_component_indices(i));
882 if (is_symmetric_component(indices, op))
892 typename ValueType = Expression,
893 template <
int,
int,
typename>
895 TensorType<rank, dim, ValueType>
896 tensor_diff_tensor(
const TensorType<0, dim, ValueType> &func,
897 const TensorType<rank, dim, ValueType> &op)
899 TensorType<rank, dim, ValueType> out;
900 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
903 out.unrolled_to_component_indices(i));
906 if (is_symmetric_component(indices, op))
915 typename ValueType = Expression,
916 template <
int,
int,
typename>
918 TensorType<rank, dim, ValueType>
919 tensor_diff_scalar(
const TensorType<rank, dim, ValueType> &funcs,
922 TensorType<rank, dim, ValueType> out;
923 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
926 out.unrolled_to_component_indices(i));
936 typename ValueType = Expression,
937 template <
int,
int,
typename>
939 TensorType<rank, dim, ValueType>
940 tensor_diff_tensor(
const TensorType<rank, dim, ValueType> &funcs,
941 const TensorType<0, dim, ValueType> &op)
943 TensorType<rank, dim, ValueType> out;
944 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
947 out.unrolled_to_component_indices(i));
955 template <
int rank_1,
958 typename ValueType = Expression,
959 template <
int,
int,
typename>
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)
965 TensorType<rank_1 + rank_2, dim, ValueType> out;
966 for (
unsigned int i = 0; i < funcs.n_independent_components; ++i)
969 funcs.unrolled_to_component_indices(i));
970 for (
unsigned int j = 0; j < op.n_independent_components; ++j)
973 op.unrolled_to_component_indices(j));
975 concatenate_indices(indices_i, indices_j);
980 if (is_symmetric_component(indices_j, op))
981 out[indices_out] *= 0.5;
992 template <
int rank_1,
995 typename ValueType = Expression,
996 template <
int,
int,
typename>
998 template <
int,
int,
typename>
1001 tensor_diff_tensor(
const TensorType_1<rank_1, dim, ValueType> &funcs,
1002 const TensorType_2<rank_2, dim, ValueType> &op)
1005 for (
unsigned int i = 0; i < funcs.n_independent_components; ++i)
1008 funcs.unrolled_to_component_indices(i));
1009 for (
unsigned int j = 0; j < op.n_independent_components; ++j)
1012 op.unrolled_to_component_indices(j));
1014 concatenate_indices(indices_i, indices_j);
1019 if (is_symmetric_component(indices_j, op))
1020 out[indices_out] *= 0.5;
1024 TensorType_1<rank_1, dim, ValueType>,
1028 concatenate_indices(transpose_indices(indices_i),
1030 out[indices_out_t] = out[indices_out];
1032 else if (std::is_same_v<
1033 TensorType_2<rank_2, dim, ValueType>,
1038 concatenate_indices(indices_i,
1039 transpose_indices(indices_j));
1040 out[indices_out_t] = out[indices_out];
1047 "Expect mixed tensor differentiation to have at least "
1048 "one component stemming from a symmetric tensor."));
1058 template <
int rank,
int dim>
1063 return internal::scalar_diff_tensor(func, op);
1067 template <
int rank,
int dim>
1072 return internal::scalar_diff_tensor(func, op);
1076 template <
int rank,
int dim>
1081 return internal::scalar_diff_tensor(func, op);
1085 template <
int rank,
int dim>
1092 const Expression tmp = func;
1093 return internal::scalar_diff_tensor(tmp, op);
1097 template <
int rank,
int dim>
1100 const Expression &op)
1102 return internal::tensor_diff_scalar(symbol_tensor, op);
1106 template <
int rank,
int dim>
1111 return internal::tensor_diff_scalar(symbol_tensor, op);
1115 template <
int rank,
int dim>
1118 const Expression &op)
1120 return internal::tensor_diff_scalar(symbol_tensor, op);
1124 template <
int rank,
int dim>
1129 return internal::tensor_diff_scalar(symbol_tensor, op);
1133 template <
int rank_1,
int rank_2,
int dim>
1138 return internal::tensor_diff_tensor(symbol_tensor, op);
1142 template <
int rank_1,
int rank_2,
int dim>
1147 return internal::tensor_diff_tensor(symbol_tensor, op);
1151 template <
int rank_1,
int rank_2,
int dim>
1156 return internal::tensor_diff_tensor(symbol_tensor, op);
1160 template <
int rank_1,
int rank_2,
int dim>
1165 return internal::tensor_diff_tensor(symbol_tensor, op);
1174 template <
typename SymbolicType,
1178 template <
int,
int,
typename>
1181 set_tensor_value_in_symbol_map(
1183 const TensorType<rank, dim, SymbolicType> &symbol_tensor,
1184 const TensorType<rank, dim, ValueType> &value_tensor)
1186 TensorType<rank, dim, Expression> out;
1187 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
1190 out.unrolled_to_component_indices(i));
1192 symbol_tensor[indices],
1193 value_tensor[indices]);
1198 template <
typename SymbolicType,
typename ValueType,
int dim>
1200 set_tensor_value_in_symbol_map(
1206 for (
unsigned int i = 0;
1207 i < SymmetricTensor<2, dim>::n_independent_components;
1209 for (
unsigned int j = 0;
1210 j < SymmetricTensor<2, dim>::n_independent_components;
1214 make_rank_4_tensor_indices<dim>(i, j);
1216 symbol_tensor[indices],
1217 value_tensor[indices]);
1223 template <
bool ignore_invalid_symbols,
1227 typename SymbolicType>
1233 add_to_substitution_map<ignore_invalid_symbols>(
1238 template <
bool ignore_invalid_symbols,
1242 typename SymbolicType>
1249 add_to_substitution_map<ignore_invalid_symbols>(
1254 template <
int rank,
int dim,
typename SymbolicType,
typename ValueType>
1267 template <
int rank,
int dim,
typename SymbolicType,
typename ValueType>
1283 template <
int rank,
int dim,
typename ExpressionType,
typename ValueType>
1295 template <
int rank,
int dim,
typename ExpressionType,
typename ValueType>
1314 typename ExpressionType,
1316 template <
int,
int,
typename>
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)
1323 std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1324 for (
unsigned int i = 0; i < symbol_tensor.n_independent_components;
1328 symbol_tensor.unrolled_to_component_indices(i));
1329 symbol_values.push_back(
1330 std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1332 return symbol_values;
1336 template <
int dim,
typename ExpressionType,
typename ValueType>
1337 std::vector<std::pair<ExpressionType, ValueType>>
1338 make_tensor_entries_for_substitution_map(
1342 const ExpressionType &expression = symbol_tensor;
1343 const ValueType &
value = value_tensor;
1344 return {std::make_pair(expression, value)};
1348 template <
int dim,
typename ExpressionType,
typename ValueType>
1349 std::vector<std::pair<ExpressionType, ValueType>>
1350 make_tensor_entries_for_substitution_map(
1354 std::vector<std::pair<ExpressionType, ValueType>> symbol_values;
1355 for (
unsigned int i = 0;
1356 i < SymmetricTensor<2, dim>::n_independent_components;
1358 for (
unsigned int j = 0;
1359 j < SymmetricTensor<2, dim>::n_independent_components;
1363 make_rank_4_tensor_indices<dim>(i, j);
1364 symbol_values.push_back(
1365 std::make_pair(symbol_tensor[indices], value_tensor[indices]));
1367 return symbol_values;
1372 template <
bool ignore_invalid_symbols,
1375 typename ExpressionType,
1383 add_to_substitution_map<ignore_invalid_symbols>(
1385 internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1390 template <
bool ignore_invalid_symbols,
1393 typename ExpressionType,
1401 add_to_substitution_map<ignore_invalid_symbols>(
1403 internal::make_tensor_entries_for_substitution_map(symbol_tensor,
1415 template <
int,
int,
typename>
1417 TensorType<rank, dim, Expression>
1419 const TensorType<rank, dim, Expression> &expression_tensor,
1422 TensorType<rank, dim, Expression> out;
1423 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
1426 out.unrolled_to_component_indices(i));
1439 const Expression &expression = expression_tensor;
1451 for (
unsigned int i = 0;
1452 i < SymmetricTensor<2, dim>::n_independent_components;
1454 for (
unsigned int j = 0;
1455 j < SymmetricTensor<2, dim>::n_independent_components;
1459 make_rank_4_tensor_indices<dim>(i, j);
1467 template <
typename ValueType,
1470 template <
int,
int,
typename>
1472 TensorType<rank, dim, ValueType>
1473 substitute_and_evaluate_tensor(
1474 const TensorType<rank, dim, Expression> &expression_tensor,
1477 TensorType<rank, dim, ValueType> out;
1478 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
1481 out.unrolled_to_component_indices(i));
1483 substitute_and_evaluate<ValueType>(expression_tensor[indices],
1490 template <
typename ValueType,
int dim>
1492 substitute_and_evaluate_tensor(
1496 const Expression &expression = expression_tensor;
1501 template <
typename ValueType,
int dim>
1503 substitute_and_evaluate_tensor(
1508 for (
unsigned int i = 0;
1509 i < SymmetricTensor<2, dim>::n_independent_components;
1511 for (
unsigned int j = 0;
1512 j < SymmetricTensor<2, dim>::n_independent_components;
1516 make_rank_4_tensor_indices<dim>(i, j);
1518 substitute_and_evaluate<ValueType>(expression_tensor[indices],
1526 template <
int rank,
int dim>
1535 template <
int rank,
int dim>
1544 template <
typename ValueType,
int rank,
int dim>
1550 return internal::substitute_and_evaluate_tensor<ValueType>(
1555 template <
typename ValueType,
int rank,
int dim>
1561 return internal::substitute_and_evaluate_tensor<ValueType>(
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#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< 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)