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