Reference documentation for deal.II version 9.4.1
Searching...
No Matches
auto_derivative_function.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 2018 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
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.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17#include <deal.II/base/point.h>
18
19#include <deal.II/lac/vector.h>
20
21#include <cmath>
22
24
25template <int dim>
27 const double hh,
28 const unsigned int n_components,
29 const double initial_time)
30 : Function<dim>(n_components, initial_time)
31 , h(1)
32 , ht(dim)
33 , formula(Euler)
34{
35 set_h(hh);
37}
38
39
40
41template <int dim>
42void
44{
45 // go through all known formulas, reject ones we don't know about
46 // and don't handle in the member functions of this class
47 switch (form)
48 {
49 case Euler:
50 case UpwindEuler:
51 case FourthOrder:
52 break;
53 default:
54 Assert(false,
55 ExcMessage("The argument passed to this function does not "
56 "match any known difference formula."));
57 }
58
59 formula = form;
60}
61
62
63template <int dim>
64void
66{
67 h = hh;
68 for (unsigned int i = 0; i < dim; ++i)
69 ht[i][i] = h;
70}
71
72
73template <int dim>
76 const unsigned int comp) const
77{
79 switch (formula)
80 {
81 case UpwindEuler:
82 {
83 Point<dim> q1;
84 for (unsigned int i = 0; i < dim; ++i)
85 {
86 q1 = p - ht[i];
87 grad[i] = (this->value(p, comp) - this->value(q1, comp)) / h;
88 }
89 break;
90 }
91 case Euler:
92 {
93 Point<dim> q1, q2;
94 for (unsigned int i = 0; i < dim; ++i)
95 {
96 q1 = p + ht[i];
97 q2 = p - ht[i];
99 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
100 }
101 break;
102 }
103 case FourthOrder:
104 {
105 Point<dim> q1, q2, q3, q4;
106 for (unsigned int i = 0; i < dim; ++i)
107 {
108 q2 = p + ht[i];
109 q1 = q2 + ht[i];
110 q3 = p - ht[i];
111 q4 = q3 - ht[i];
112 grad[i] = (-this->value(q1, comp) + 8 * this->value(q2, comp) -
113 8 * this->value(q3, comp) + this->value(q4, comp)) /
114 (12 * h);
115 }
116 break;
117 }
118 default:
119 Assert(false, ExcNotImplemented());
120 }
122}
123
124
125template <int dim>
126void
128 const Point<dim> & p,
130{
133
134 switch (formula)
135 {
136 case UpwindEuler:
137 {
138 Point<dim> q1;
139 Vector<double> v(this->n_components), v1(this->n_components);
140 const double h_inv = 1. / h;
141 for (unsigned int i = 0; i < dim; ++i)
142 {
143 q1 = p - ht[i];
144 this->vector_value(p, v);
145 this->vector_value(q1, v1);
146
147 for (unsigned int comp = 0; comp < this->n_components; ++comp)
148 gradients[comp][i] = (v(comp) - v1(comp)) * h_inv;
149 }
150 break;
151 }
152
153 case Euler:
154 {
155 Point<dim> q1, q2;
156 Vector<double> v1(this->n_components), v2(this->n_components);
157 const double h_inv_2 = 1. / (2 * h);
158 for (unsigned int i = 0; i < dim; ++i)
159 {
160 q1 = p + ht[i];
161 q2 = p - ht[i];
162 this->vector_value(q1, v1);
163 this->vector_value(q2, v2);
164
165 for (unsigned int comp = 0; comp < this->n_components; ++comp)
166 gradients[comp][i] = (v1(comp) - v2(comp)) * h_inv_2;
167 }
168 break;
169 }
170
171 case FourthOrder:
172 {
173 Point<dim> q1, q2, q3, q4;
174 Vector<double> v1(this->n_components), v2(this->n_components),
175 v3(this->n_components), v4(this->n_components);
176 const double h_inv_12 = 1. / (12 * h);
177 for (unsigned int i = 0; i < dim; ++i)
178 {
179 q2 = p + ht[i];
180 q1 = q2 + ht[i];
181 q3 = p - ht[i];
182 q4 = q3 - ht[i];
183 this->vector_value(q1, v1);
184 this->vector_value(q2, v2);
185 this->vector_value(q3, v3);
186 this->vector_value(q4, v4);
187
188 for (unsigned int comp = 0; comp < this->n_components; ++comp)
190 (-v1(comp) + 8 * v2(comp) - 8 * v3(comp) + v4(comp)) *
191 h_inv_12;
192 }
193 break;
194 }
195
196 default:
197 Assert(false, ExcNotImplemented());
198 }
199}
200
201
202template <int dim>
203void
205 const std::vector<Point<dim>> &points,
207 const unsigned int comp) const
208{
211
212 switch (formula)
213 {
214 case UpwindEuler:
215 {
216 Point<dim> q1;
217 for (unsigned int p = 0; p < points.size(); ++p)
218 for (unsigned int i = 0; i < dim; ++i)
219 {
220 q1 = points[p] - ht[i];
222 (this->value(points[p], comp) - this->value(q1, comp)) / h;
223 }
224 break;
225 }
226
227 case Euler:
228 {
229 Point<dim> q1, q2;
230 for (unsigned int p = 0; p < points.size(); ++p)
231 for (unsigned int i = 0; i < dim; ++i)
232 {
233 q1 = points[p] + ht[i];
234 q2 = points[p] - ht[i];
236 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
237 }
238 break;
239 }
240
241 case FourthOrder:
242 {
243 Point<dim> q1, q2, q3, q4;
244 for (unsigned int p = 0; p < points.size(); ++p)
245 for (unsigned int i = 0; i < dim; ++i)
246 {
247 q2 = points[p] + ht[i];
248 q1 = q2 + ht[i];
249 q3 = points[p] - ht[i];
250 q4 = q3 - ht[i];
252 (-this->value(q1, comp) + 8 * this->value(q2, comp) -
253 8 * this->value(q3, comp) + this->value(q4, comp)) /
254 (12 * h);
255 }
256 break;
257 }
258
259 default:
260 Assert(false, ExcNotImplemented());
261 }
262}
263
264
265
266template <int dim>
267void
269 const std::vector<Point<dim>> & points,
271{
274 for (unsigned int p = 0; p < points.size(); ++p)
277
278 switch (formula)
279 {
280 case UpwindEuler:
281 {
282 Point<dim> q1;
283 for (unsigned int p = 0; p < points.size(); ++p)
284 for (unsigned int i = 0; i < dim; ++i)
285 {
286 q1 = points[p] - ht[i];
287 for (unsigned int comp = 0; comp < this->n_components; ++comp)
289 (this->value(points[p], comp) - this->value(q1, comp)) / h;
290 }
291 break;
292 }
293
294 case Euler:
295 {
296 Point<dim> q1, q2;
297 for (unsigned int p = 0; p < points.size(); ++p)
298 for (unsigned int i = 0; i < dim; ++i)
299 {
300 q1 = points[p] + ht[i];
301 q2 = points[p] - ht[i];
302 for (unsigned int comp = 0; comp < this->n_components; ++comp)
304 (this->value(q1, comp) - this->value(q2, comp)) / (2 * h);
305 }
306 break;
307 }
308
309 case FourthOrder:
310 {
311 Point<dim> q1, q2, q3, q4;
312 for (unsigned int p = 0; p < points.size(); ++p)
313 for (unsigned int i = 0; i < dim; ++i)
314 {
315 q2 = points[p] + ht[i];
316 q1 = q2 + ht[i];
317 q3 = points[p] - ht[i];
318 q4 = q3 - ht[i];
319 for (unsigned int comp = 0; comp < this->n_components; ++comp)
321 (-this->value(q1, comp) + 8 * this->value(q2, comp) -
322 8 * this->value(q3, comp) + this->value(q4, comp)) /
323 (12 * h);
324 }
325 break;
326 }
327
328 default:
329 Assert(false, ExcNotImplemented());
330 }
331}
332
333
334template <int dim>
337{
338 switch (ord)
339 {
340 case 0:
341 case 1:
342 return UpwindEuler;
343 case 2:
344 return Euler;
345 case 3:
346 case 4:
347 return FourthOrder;
348 default:
349 Assert(false, ExcNotImplemented());
350 }
351 return Euler;
352}
353
354
355template class AutoDerivativeFunction<1>;
356template class AutoDerivativeFunction<2>;
357template class AutoDerivativeFunction<3>;
358
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
Definition: tensor.h:503
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
const unsigned int v1
Definition: grid_tools.cc:1000
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)