1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
|
// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2001 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
#ifndef dealii_auto_derivative_function_h
#define dealii_auto_derivative_function_h
#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
DEAL_II_NAMESPACE_OPEN
/**
* This class automatically computes the gradient of a function by employing
* numerical difference quotients. This only, if the user function does not
* provide the gradient function himself.
*
* The following example of an user defined function overloads and implements
* only the value() function but not the gradient() function. If the
* gradient() function is invoked then the gradient function implemented by
* the AutoDerivativeFunction is called, where the latter function employs
* numerical difference quotients.
*
* @code
* class UserFunction: public AutoDerivativeFunction
* {
* // access to one component at one point
* double value (const Point<dim> &p,
* const unsigned int component = 0) const override
* {
* // Implementation ....
* };
* };
*
* UserFunction user_function;
*
* // gradient by employing difference quotients.
* Tensor<1,dim> grad=user_function.gradient(some_point);
* @endcode
*
* If the user overloads and implements also the gradient function, then, of
* course, the users gradient function is called.
*
* Note, that the usage of the value() and gradient() functions explained
* above, also applies to the value_list() and gradient_list() functions as
* well as to the vector valued versions of these functions, see e.g.
* vector_value(), vector_gradient(), vector_value_list() and
* vector_gradient_list().
*
* The gradient() and gradient_list() functions make use of the
* Function::value() function. The vector_gradient() and
* vector_gradient_list() make use of the Function::vector_value() function.
* Make sure that the user defined function implements the value() function
* and the vector_value() function, respectively.
*
* Furthermore note, that an object of this class does <b>not</b> represent
* the derivative of a function, like FunctionDerivative, that gives a
* directional derivative by calling the value() function. In fact, this class
* (the AutoDerivativeFunction class) can substitute the Function class as
* base class for user defined classes. This class implements the gradient()
* functions for automatic computation of numerical difference quotients and
* serves as intermediate class between the base Function class and the user
* defined function class.
*
* @ingroup functions
*/
template <int dim>
class AutoDerivativeFunction : public Function<dim>
{
public:
/**
* Names of difference formulas.
*/
enum DifferenceFormula
{
/**
* The symmetric Euler formula of second order:
* @f[
* u'(t) \approx
* \frac{u(t+h) -
* u(t-h)}{2h}.
* @f]
*/
Euler,
/**
* The upwind Euler formula of first order:
* @f[
* u'(t) \approx
* \frac{u(t) -
* u(t-h)}{h}.
* @f]
*/
UpwindEuler,
/**
* The fourth order scheme
* @f[
* u'(t) \approx
* \frac{u(t-2h) - 8u(t-h)
* + 8u(t+h) - u(t+2h)}{12h}.
* @f]
*/
FourthOrder
};
/**
* Constructor. Takes the difference step size <tt>h</tt>. It's within the
* user's responsibility to choose an appropriate value here. <tt>h</tt>
* should be chosen taking into account the absolute value as well as the
* amount of local variation of the function. Setting <tt>h=1e-6</tt> might
* be a good choice for functions with an absolute value of about 1, that
* furthermore does not vary to much.
*
* <tt>h</tt> can be changed later using the set_h() function.
*
* Sets DifferenceFormula <tt>formula</tt> to the default <tt>Euler</tt>
* formula of the set_formula() function. Change this preset formula by
* calling the set_formula() function.
*/
AutoDerivativeFunction(const double h,
const unsigned int n_components = 1,
const double initial_time = 0.0);
/**
* Virtual destructor; absolutely necessary in this case.
*/
virtual ~AutoDerivativeFunction() override = default;
/**
* Choose the difference formula. See the enum #DifferenceFormula for
* available choices.
*/
void
set_formula(const DifferenceFormula formula = Euler);
/**
* Takes the difference step size <tt>h</tt>. It's within the user's
* responsibility to choose an appropriate value here. <tt>h</tt> should be
* chosen taking into account the absolute value of as well as the amount of
* local variation of the function. Setting <tt>h=1e-6</tt> might be a good
* choice for functions with an absolute value of about 1, that furthermore
* does not vary to much.
*/
void
set_h(const double h);
/**
* Return the gradient of the specified component of the function at the
* given point.
*
* Compute numerical difference quotients using the preset
* #DifferenceFormula.
*/
virtual Tensor<1, dim>
gradient(const Point<dim> &p,
const unsigned int component = 0) const override;
/**
* Return the gradient of all components of the function at the given point.
*
* Compute numerical difference quotients using the preset
* #DifferenceFormula.
*/
virtual void
vector_gradient(const Point<dim> &p,
std::vector<Tensor<1, dim>> &gradients) const override;
/**
* Set <tt>gradients</tt> to the gradients of the specified component of the
* function at the <tt>points</tt>. It is assumed that <tt>gradients</tt>
* already has the right size, i.e. the same size as the <tt>points</tt>
* array.
*
* Compute numerical difference quotients using the preset
* #DifferenceFormula.
*/
virtual void
gradient_list(const std::vector<Point<dim>> &points,
std::vector<Tensor<1, dim>> &gradients,
const unsigned int component = 0) const override;
/**
* Set <tt>gradients</tt> to the gradients of the function at the
* <tt>points</tt>, for all components. It is assumed that
* <tt>gradients</tt> already has the right size, i.e. the same size as the
* <tt>points</tt> array.
*
* The outer loop over <tt>gradients</tt> is over the points in the list,
* the inner loop over the different components of the function.
*
* Compute numerical difference quotients using the preset
* #DifferenceFormula.
*/
virtual void
vector_gradient_list(
const std::vector<Point<dim>> &points,
std::vector<std::vector<Tensor<1, dim>>> &gradients) const override;
/**
* Return a #DifferenceFormula of the order <tt>ord</tt> at minimum.
*/
static DifferenceFormula
get_formula_of_order(const unsigned int ord);
private:
/**
* Step size of the difference formula. Set by the set_h() function.
*/
double h;
/**
* Includes the unit vectors scaled by <tt>h</tt>.
*/
std::vector<Tensor<1, dim>> ht;
/**
* Difference formula. Set by the set_formula() function.
*/
DifferenceFormula formula;
};
DEAL_II_NAMESPACE_CLOSE
#endif
|
