Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
auto_derivative_function.cc
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
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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/auto_derivative_function.h>
17 #include <deal.II/base/point.h>
18 
19 #include <deal.II/lac/vector.h>
20 
21 #include <cmath>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 template <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);
36  set_formula();
37 }
38 
39 
40 
41 template <int dim>
42 void
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 
63 template <int dim>
64 void
66 {
67  h = hh;
68  for (unsigned int i = 0; i < dim; ++i)
69  ht[i][i] = h;
70 }
71 
72 
73 template <int dim>
76  const unsigned int comp) const
77 {
78  Tensor<1, dim> grad;
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];
98  grad[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  }
121  return grad;
122 }
123 
124 
125 template <int dim>
126 void
128  const Point<dim> & p,
129  std::vector<Tensor<1, dim>> &gradients) const
130 {
131  Assert(gradients.size() == this->n_components,
132  ExcDimensionMismatch(gradients.size(), this->n_components));
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)
189  gradients[comp][i] =
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 
202 template <int dim>
203 void
205  const std::vector<Point<dim>> &points,
206  std::vector<Tensor<1, dim>> & gradients,
207  const unsigned int comp) const
208 {
209  Assert(gradients.size() == points.size(),
210  ExcDimensionMismatch(gradients.size(), points.size()));
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];
221  gradients[p][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];
235  gradients[p][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];
251  gradients[p][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 
266 template <int dim>
267 void
269  const std::vector<Point<dim>> & points,
270  std::vector<std::vector<Tensor<1, dim>>> &gradients) const
271 {
272  Assert(gradients.size() == points.size(),
273  ExcDimensionMismatch(gradients.size(), points.size()));
274  for (unsigned int p = 0; p < points.size(); ++p)
275  Assert(gradients[p].size() == this->n_components,
276  ExcDimensionMismatch(gradients.size(), this->n_components));
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)
288  gradients[p][comp][i] =
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)
303  gradients[p][comp][i] =
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)
320  gradients[p][comp][i] =
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 
334 template <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 
355 template class AutoDerivativeFunction<1>;
356 template class AutoDerivativeFunction<2>;
357 template class AutoDerivativeFunction<3>;
358 
359 DEAL_II_NAMESPACE_CLOSE
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)
void set_formula(const DifferenceFormula formula=Euler)
virtual void vector_gradient_list(const std::vector< Point< dim >> &points, std::vector< std::vector< Tensor< 1, dim >>> &gradients) const override
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim >> &gradients) const override
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
AutoDerivativeFunction(const double h, const unsigned int n_components=1, const double initial_time=0.0)
static ::ExceptionBase & ExcNotImplemented()
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override