Reference documentation for deal.II version Git 20c3758ecc 2020-10-30 06:09:42 +0100
\(\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 
30 namespace 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 
48  make_symbolic_function(const std::string & symbol,
49  const SD::types::symbol_vector &arguments)
50  {
51  return Expression(symbol, arguments);
52  }
53 
54 
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 
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(
134  ExcMessage(
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  {
152  internal::set_value_in_symbol_map(substitution_map,
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  {
175  add_to_substitution_map(substitution_map, symbol, value);
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();
247  AssertThrow(
248  iter < size,
249  ExcMessage(
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,
262  const types::substitution_map &substitution_map)
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
Expression differentiate(const Expression &symbol) const
bool is_valid_substitution_symbol(const SymEngine::Basic &entry)
void merge_substitution_maps(types::substitution_map &substitution_map_out, const types::substitution_map &substitution_map_in)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
Expression make_symbolic_function(const std::string &symbol, const types::symbol_vector &arguments)
types::substitution_map resolve_explicit_dependencies(const types::substitution_map &substitution_map, const bool force_cyclic_dependency_resolution=false)
Expression differentiate(const Expression &f, const Expression &x)
Expression substitute(const types::substitution_map &substitution_values) const
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1423
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)
void add_to_substitution_map(types::substitution_map &substitution_map, const SymEngine::RCP< const SymEngine::Basic > &symbol, const SymEngine::RCP< const SymEngine::Basic > &value)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
std::vector< SD::Expression > symbol_vector
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
Expression make_symbol(const std::string &symbol)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map