Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
Searching...
No Matches
auto_derivative_function.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/point.h>
17
18#include <deal.II/lac/vector.h>
19
20#include <cmath>
21
23
24template <int dim>
26 const double hh,
27 const unsigned int n_components,
28 const double initial_time)
29 : Function<dim>(n_components, initial_time)
30 , h(1)
31 , ht(dim)
32 , formula(Euler)
33{
34 set_h(hh);
36}
37
38
39
40template <int dim>
41void
43{
44 // go through all known formulas, reject ones we don't know about
45 // and don't handle in the member functions of this class
46 switch (form)
47 {
48 case Euler:
49 case UpwindEuler:
50 case FourthOrder:
51 break;
52 default:
53 Assert(false,
54 ExcMessage("The argument passed to this function does not "
55 "match any known difference formula."));
56 }
57
58 formula = form;
59}
60
61
62template <int dim>
63void
65{
66 h = hh;
67 for (unsigned int i = 0; i < dim; ++i)
68 ht[i][i] = h;
69}
70
71
72template <int dim>
75 const unsigned int comp) const
76{
78 switch (formula)
79 {
80 case UpwindEuler:
81 {
82 Point<dim> q1;
83 for (unsigned int i = 0; i < dim; ++i)
84 {
85 q1 = p - ht[i];
86 grad[i] = (this->value(p, comp) - this->value(q1, comp)) / h;
87 }
88 break;
89 }
90 case Euler:
91 {
92 Point<dim> q1, q2;
93 for (unsigned int i = 0; i < dim; ++i)
94 {
95 q1 = p + ht[i];
96 q2 = p - ht[i];
98 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
99 }
100 break;
101 }
102 case FourthOrder:
103 {
104 Point<dim> q1, q2, q3, q4;
105 for (unsigned int i = 0; i < dim; ++i)
106 {
107 q2 = p + ht[i];
108 q1 = q2 + ht[i];
109 q3 = p - ht[i];
110 q4 = q3 - ht[i];
111 grad[i] = (-this->value(q1, comp) + 8 * this->value(q2, comp) -
112 8 * this->value(q3, comp) + this->value(q4, comp)) /
113 (12 * h);
114 }
115 break;
116 }
117 default:
119 }
121}
122
123
124template <int dim>
125void
127 const Point<dim> &p,
129{
132
133 switch (formula)
134 {
135 case UpwindEuler:
136 {
137 Point<dim> q1;
138 Vector<double> v(this->n_components), v1(this->n_components);
139 const double h_inv = 1. / h;
140 for (unsigned int i = 0; i < dim; ++i)
141 {
142 q1 = p - ht[i];
143 this->vector_value(p, v);
144 this->vector_value(q1, v1);
145
146 for (unsigned int comp = 0; comp < this->n_components; ++comp)
147 gradients[comp][i] = (v(comp) - v1(comp)) * h_inv;
148 }
149 break;
150 }
151
152 case Euler:
153 {
154 Point<dim> q1, q2;
155 Vector<double> v1(this->n_components), v2(this->n_components);
156 const double h_inv_2 = 1. / (2 * h);
157 for (unsigned int i = 0; i < dim; ++i)
158 {
159 q1 = p + ht[i];
160 q2 = p - ht[i];
161 this->vector_value(q1, v1);
162 this->vector_value(q2, v2);
163
164 for (unsigned int comp = 0; comp < this->n_components; ++comp)
165 gradients[comp][i] = (v1(comp) - v2(comp)) * h_inv_2;
166 }
167 break;
168 }
169
170 case FourthOrder:
171 {
172 Point<dim> q1, q2, q3, q4;
173 Vector<double> v1(this->n_components), v2(this->n_components),
174 v3(this->n_components), v4(this->n_components);
175 const double h_inv_12 = 1. / (12 * h);
176 for (unsigned int i = 0; i < dim; ++i)
177 {
178 q2 = p + ht[i];
179 q1 = q2 + ht[i];
180 q3 = p - ht[i];
181 q4 = q3 - ht[i];
182 this->vector_value(q1, v1);
183 this->vector_value(q2, v2);
184 this->vector_value(q3, v3);
185 this->vector_value(q4, v4);
186
187 for (unsigned int comp = 0; comp < this->n_components; ++comp)
189 (-v1(comp) + 8 * v2(comp) - 8 * v3(comp) + v4(comp)) *
190 h_inv_12;
191 }
192 break;
193 }
194
195 default:
197 }
198}
199
200
201template <int dim>
202void
204 const std::vector<Point<dim>> &points,
206 const unsigned int comp) const
207{
210
211 switch (formula)
212 {
213 case UpwindEuler:
214 {
215 Point<dim> q1;
216 for (unsigned int p = 0; p < points.size(); ++p)
217 for (unsigned int i = 0; i < dim; ++i)
218 {
219 q1 = points[p] - ht[i];
221 (this->value(points[p], comp) - this->value(q1, comp)) / h;
222 }
223 break;
224 }
225
226 case Euler:
227 {
228 Point<dim> q1, q2;
229 for (unsigned int p = 0; p < points.size(); ++p)
230 for (unsigned int i = 0; i < dim; ++i)
231 {
232 q1 = points[p] + ht[i];
233 q2 = points[p] - ht[i];
235 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
236 }
237 break;
238 }
239
240 case FourthOrder:
241 {
242 Point<dim> q1, q2, q3, q4;
243 for (unsigned int p = 0; p < points.size(); ++p)
244 for (unsigned int i = 0; i < dim; ++i)
245 {
246 q2 = points[p] + ht[i];
247 q1 = q2 + ht[i];
248 q3 = points[p] - ht[i];
249 q4 = q3 - ht[i];
251 (-this->value(q1, comp) + 8 * this->value(q2, comp) -
252 8 * this->value(q3, comp) + this->value(q4, comp)) /
253 (12 * h);
254 }
255 break;
256 }
257
258 default:
260 }
261}
262
263
264
265template <int dim>
266void
268 const std::vector<Point<dim>> &points,
270{
273 for (unsigned int p = 0; p < points.size(); ++p)
276
277 switch (formula)
278 {
279 case UpwindEuler:
280 {
281 Point<dim> q1;
282 for (unsigned int p = 0; p < points.size(); ++p)
283 for (unsigned int i = 0; i < dim; ++i)
284 {
285 q1 = points[p] - ht[i];
286 for (unsigned int comp = 0; comp < this->n_components; ++comp)
288 (this->value(points[p], comp) - this->value(q1, comp)) / h;
289 }
290 break;
291 }
292
293 case Euler:
294 {
295 Point<dim> q1, q2;
296 for (unsigned int p = 0; p < points.size(); ++p)
297 for (unsigned int i = 0; i < dim; ++i)
298 {
299 q1 = points[p] + ht[i];
300 q2 = points[p] - ht[i];
301 for (unsigned int comp = 0; comp < this->n_components; ++comp)
303 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
304 }
305 break;
306 }
307
308 case FourthOrder:
309 {
310 Point<dim> q1, q2, q3, q4;
311 for (unsigned int p = 0; p < points.size(); ++p)
312 for (unsigned int i = 0; i < dim; ++i)
313 {
314 q2 = points[p] + ht[i];
315 q1 = q2 + ht[i];
316 q3 = points[p] - ht[i];
317 q4 = q3 - ht[i];
318 for (unsigned int comp = 0; comp < this->n_components; ++comp)
320 (-this->value(q1, comp) + 8 * this->value(q2, comp) -
321 8 * this->value(q3, comp) + this->value(q4, comp)) /
322 (12 * h);
323 }
324 break;
325 }
326
327 default:
329 }
330}
331
332
333template <int dim>
336{
337 switch (ord)
338 {
339 case 0:
340 case 1:
341 return UpwindEuler;
342 case 2:
343 return Euler;
344 case 3:
345 case 4:
346 return FourthOrder;
347 default:
349 }
350 return Euler;
351}
352
353
354template class AutoDerivativeFunction<1>;
355template class AutoDerivativeFunction<2>;
356template class AutoDerivativeFunction<3>;
357
AutoDerivativeFunction(const double h, const unsigned int n_components=1, const double initial_time=0.0)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim > > &gradients) const override
virtual void vector_gradient_list(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const override
void set_formula(const DifferenceFormula formula=Euler)
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
static DifferenceFormula get_formula_of_order(const unsigned int ord)
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
const unsigned int v1
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()