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_scalar_operations.cc
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
17#include <deal.II/base/config.h>
18
19#ifdef DEAL_II_WITH_SYMENGINE
20
25
26# include <symengine/real_double.h>
27
29
30namespace Differentiation
31{
32 namespace SD
33 {
34 namespace SE = ::SymEngine;
35
36
37 /* ------------------------- Symbol creation -----------------------*/
38
39
40 Expression
41 make_symbol(const std::string &symbol)
42 {
43 return Expression(symbol);
44 }
45
46
47# ifndef DOXYGEN
48 Expression
49 make_symbolic_function(const std::string & symbol,
50 const SD::types::symbol_vector &arguments)
51 {
52 return Expression(symbol, arguments);
53 }
54
55
56 Expression
57 make_symbolic_function(const std::string & symbol,
58 const SD::types::substitution_map &arguments)
59 {
60 return make_symbolic_function(symbol,
62 }
63# endif
64
65
66 /* --------------------------- Differentiation ------------------------- */
67
68
69 Expression
70 differentiate(const Expression &func, const Expression &op)
71 {
72 return func.differentiate(op);
73 }
74
75
76 /* ---------------- Symbol map creation and manipulation --------------*/
77
78
79# ifndef DOXYGEN
80 namespace internal
81 {
82 bool
83 is_valid_substitution_symbol(const SE::Basic &entry)
84 {
85 // Allow substitution of a symbol
86 // It is pretty clear as to why this is wanted...
87 if (SE::is_a<SE::Symbol>(entry))
88 return true;
89
90 // Allow substitution of a function symbol
91 // If desired, we can transform general but undefined functional
92 // relationships to an explicit form that is concrete. This is
93 // required for a symbolic expression to be parsed by a Lambda or LLVM
94 // optimizer.
95 if (SE::is_a<SE::FunctionSymbol>(entry))
96 return true;
97
98 // Allow substitution of the explicit expression of the derivative of
99 // an implicitly defined symbol (i.e. the result of the derivative of
100 // a FunctionSymbol).
101 if (SE::is_a<SE::Derivative>(entry))
102 return true;
103
104 // When performing tensor differentiation, one may end up with a
105 // coefficient of one half due to symmetry operations, e.g.
106 // 0.5*Derivative(Qi_00(C_11, C_00, C_01), C_01)
107 // So we explicitly check for this exact case
108 if (SE::is_a<SE::Mul>(entry))
109 {
110 const SE::Mul &entry_mul = SE::down_cast<const SE::Mul &>(entry);
111 // Check that the factor is a half...
112 if (SE::eq(*(entry_mul.get_coef()), *SE::real_double(0.5)))
113 {
114 // ...and that there is only one entry and that its a
115 // Derivative type
116 const SE::map_basic_basic &entry_mul_dict =
117 entry_mul.get_dict();
118 if (entry_mul_dict.size() == 1 &&
119 SE::is_a<SE::Derivative>(*(entry_mul_dict.begin()->first)))
120 return true;
121 }
122 }
123
124 return false;
125 }
126
127
128 void
130 types::substitution_map & substitution_map,
131 const SymEngine::RCP<const SymEngine::Basic> &symbol,
132 const SymEngine::RCP<const SymEngine::Basic> &value)
133 {
134 Assert(
137 "Substitution with a number that does not represent a symbol or symbolic derivative"));
138
139 auto it_sym = substitution_map.find(Expression(symbol));
140 Assert(it_sym != substitution_map.end(),
141 ExcMessage("Did not find this symbol in the map."));
142
143 it_sym->second = Expression(value);
144 }
145
146 } // namespace internal
147
148
149 void
151 const Expression & symbol,
152 const Expression & value)
153 {
154 internal::set_value_in_symbol_map(substitution_map,
155 symbol.get_RCP(),
156 value.get_RCP());
157 }
158
159
160 void
162 const types::substitution_map &symbol_values)
163 {
164 for (const auto &entry : symbol_values)
166 }
167# endif
168
169
170 /* ------------------ Symbol substitution map creation ----------------*/
171
172
174 make_substitution_map(const Expression &symbol, const Expression &value)
175 {
176 types::substitution_map substitution_map;
177 add_to_substitution_map(substitution_map, symbol, value);
178 return substitution_map;
179 }
180
181
182# ifndef DOXYGEN
183 /* ---------------- Symbolic substitution map enlargement --------------*/
184
185
186 void
188 const types::substitution_map &symb_map_in)
189 {
190 // Do this by hand so that we can perform some sanity checks
191 for (const auto &entry : symb_map_in)
192 {
193 const typename types::substitution_map::const_iterator it_other =
194 symb_map_out.find(entry.first);
195 if (it_other == symb_map_out.end())
196 symb_map_out.insert(std::make_pair(entry.first, entry.second));
197 else
198 {
199 Assert(SE::eq(*(entry.second.get_RCP()),
200 *(it_other->second.get_RCP())),
201 ExcMessage("Key already in map, but values don't match"));
202 }
203 }
204 }
205
206
207 /* ---------------- Symbol substitution and evaluation --------------*/
208
209
212 const bool force_cyclic_dependency_resolution)
213 {
214 types::substitution_map symbol_values_resolved = symbol_values;
215 const std::size_t size = symbol_values.size();
216 (void)size;
217 for (auto &entry : symbol_values_resolved)
218 {
219 // Perform dictionary-based substitution to
220 // resolve all explicit relations in a map.
221 // Instead of checking by value (and thus having
222 // to store a temporary value), we check to see
223 // if the hash of the map entry changes.
224 Expression & out = entry.second;
225 SE::hash_t hash_old;
226 SE::hash_t hash_new = out.get_RCP()->hash();
227 unsigned int iter = 0;
228 do
229 {
230 // Write the substituted value straight back
231 // into the map.
232 if (force_cyclic_dependency_resolution)
233 {
234 // Here we substitute in symbol_values_resolved instead of
235 // symbol_values, in order to break any cyclic dependencies.
236 // The earlier entries in the dictionary are in this way
237 // guaranteed to be resolved before any subsequent entries,
238 // thereby breaking the dependency loop.
239 out = out.substitute(symbol_values_resolved);
240 }
241 else
242 {
243 out = out.substitute(symbol_values);
244 }
245
246 // Compute and store the hash of the new object
247 hash_old = hash_new;
248 hash_new = out.get_RCP()->hash();
250 iter < size,
252 "Unresolvable cyclic dependency detected in substitution map."));
253 ++iter;
254 }
255 while (hash_new != hash_old);
256 }
257
258 return symbol_values_resolved;
259 }
260
261
262 Expression
263 substitute(const Expression & expression,
264 const types::substitution_map &substitution_map)
265 {
266 return expression.substitute(substitution_map);
267 }
268# endif // DOXYGEN
269
270 } // namespace SD
271} // namespace Differentiation
272
274
275#endif // DEAL_II_WITH_SYMENGINE
Expression substitute(const types::substitution_map &substitution_values) const
Expression differentiate(const Expression &symbol) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
bool is_valid_substitution_symbol(const SymEngine::Basic &entry)
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::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
void merge_substitution_maps(types::substitution_map &substitution_map_out, const types::substitution_map &substitution_map_in)
Expression differentiate(const Expression &f, const Expression &x)
Expression make_symbolic_function(const std::string &symbol, const types::symbol_vector &arguments)
void set_value_in_symbol_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
void add_to_substitution_map(types::substitution_map &substitution_map, const Expression &symbol, const Expression &value)
types::substitution_map resolve_explicit_dependencies(const types::substitution_map &substitution_map, const bool force_cyclic_dependency_resolution=false)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
Expression make_symbol(const std::string &symbol)