Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
function_derivative.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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/function_derivative.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 Point<dim> & dir,
28  const double h)
29  : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
30  , f(f)
31  , h(h)
32  , incr(1, h * dir)
33 {
34  set_formula();
35 }
36 
37 
38 
39 template <int dim>
41  const std::vector<Point<dim>> &dir,
42  const double h)
43  : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
44  , f(f)
45  , h(h)
46  , incr(dir.size())
47 {
48  for (unsigned int i = 0; i < incr.size(); ++i)
49  incr[i] = h * dir[i];
50  set_formula();
51 }
52 
53 
54 
55 template <int dim>
56 void
59 {
60  // go through all known formulas, reject ones we don't know about
61  // and don't handle in the member functions of this class
62  switch (form)
63  {
67  break;
68  default:
69  Assert(false,
70  ExcMessage("The argument passed to this function does not "
71  "match any known difference formula."));
72  }
73 
74  formula = form;
75 }
76 
77 
78 
79 template <int dim>
80 void
81 FunctionDerivative<dim>::set_h(const double new_h)
82 {
83  for (unsigned int i = 0; i < incr.size(); ++i)
84  incr[i] *= new_h / h;
85  h = new_h;
86 }
87 
88 
89 
90 template <int dim>
91 double
93  const unsigned int component) const
94 {
95  Assert(incr.size() == 1,
96  ExcMessage(
97  "FunctionDerivative was not initialized for constant direction"));
98 
99  switch (formula)
100  {
102  return (f.value(p + incr[0], component) -
103  f.value(p - incr[0], component)) /
104  (2 * h);
106  return (f.value(p, component) - f.value(p - incr[0], component)) / h;
108  return (-f.value(p + 2 * incr[0], component) +
109  8 * f.value(p + incr[0], component) -
110  8 * f.value(p - incr[0], component) +
111  f.value(p - 2 * incr[0], component)) /
112  (12 * h);
113  default:
114  Assert(false, ExcNotImplemented());
115  }
116  return 0.;
117 }
118 
119 
120 
121 template <int dim>
122 void
124  Vector<double> & result) const
125 {
126  Assert(incr.size() == 1,
127  ExcMessage(
128  "FunctionDerivative was not initialized for constant direction"));
129  Vector<double> aux(result.size());
130 
131  // Formulas are the same as in
132  // value, but here we have to use
133  // Vector arithmetic
134  switch (formula)
135  {
137  f.vector_value(p + incr[0], result);
138  f.vector_value(p - incr[0], aux);
139  result.sadd(1. / (2 * h), -1. / (2 * h), aux);
140  return;
142  f.vector_value(p, result);
143  f.vector_value(p - incr[0], aux);
144  result.sadd(1. / h, -1. / h, aux);
145  return;
147  f.vector_value(p - 2 * incr[0], result);
148  f.vector_value(p + 2 * incr[0], aux);
149  result.add(-1., aux);
150  f.vector_value(p - incr[0], aux);
151  result.add(-8., aux);
152  f.vector_value(p + incr[0], aux);
153  result.add(8., aux);
154  result /= (12. * h);
155  return;
156  default:
157  Assert(false, ExcNotImplemented());
158  }
159 }
160 
161 
162 
163 template <int dim>
164 void
166  std::vector<double> & values,
167  const unsigned int component) const
168 {
169  const unsigned int n = points.size();
170  const bool variable_direction = (incr.size() == 1) ? false : true;
171  if (variable_direction)
172  Assert(incr.size() == points.size(),
173  ExcDimensionMismatch(incr.size(), points.size()));
174 
175  // Vector of auxiliary values
176  std::vector<double> aux(n);
177  // Vector of auxiliary points
178  std::vector<Point<dim>> paux(n);
179  // Use the same formulas as in
180  // value, but with vector
181  // arithmetic
182  switch (formula)
183  {
185  for (unsigned int i = 0; i < n; ++i)
186  paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
187  f.value_list(paux, values, component);
188  for (unsigned int i = 0; i < n; ++i)
189  paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
190  f.value_list(paux, aux, component);
191  for (unsigned int i = 0; i < n; ++i)
192  values[i] = (values[i] - aux[i]) / (2 * h);
193  return;
195  f.value_list(points, values, component);
196  for (unsigned int i = 0; i < n; ++i)
197  paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
198  f.value_list(paux, aux, component);
199  for (unsigned int i = 0; i < n; ++i)
200  values[i] = (values[i] - aux[i]) / h;
201  return;
203  for (unsigned int i = 0; i < n; ++i)
204  paux[i] = points[i] - 2 * incr[(variable_direction) ? i : 0U];
205  f.value_list(paux, values, component);
206  for (unsigned int i = 0; i < n; ++i)
207  paux[i] = points[i] + 2 * incr[(variable_direction) ? i : 0U];
208  f.value_list(paux, aux, component);
209  for (unsigned int i = 0; i < n; ++i)
210  values[i] -= aux[i];
211  for (unsigned int i = 0; i < n; ++i)
212  paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
213  f.value_list(paux, aux, component);
214  for (unsigned int i = 0; i < n; ++i)
215  values[i] += 8. * aux[i];
216  for (unsigned int i = 0; i < n; ++i)
217  paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
218  f.value_list(paux, aux, component);
219  for (unsigned int i = 0; i < n; ++i)
220  values[i] = (values[i] - 8. * aux[i]) / (12 * h);
221  return;
222  default:
223  Assert(false, ExcNotImplemented());
224  }
225 }
226 
227 
228 
229 template <int dim>
230 std::size_t
232 {
233  // only simple data elements, so
234  // use sizeof operator
235  return sizeof(*this);
236 }
237 
238 
239 
240 // explicit instantiations
241 template class FunctionDerivative<1>;
242 template class FunctionDerivative<2>;
243 template class FunctionDerivative<3>;
244 
245 DEAL_II_NAMESPACE_CLOSE
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
void sadd(const Number s, const Vector< Number > &V)
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const override
void set_formula(typename AutoDerivativeFunction< dim >::DifferenceFormula formula=AutoDerivativeFunction< dim >::Euler)
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void vector_value(const Point< dim > &p, Vector< double > &value) const override
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::size_t memory_consumption() const
FunctionDerivative(const Function< dim > &f, const Point< dim > &direction, const double h=1.e-6)
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
std::vector< Tensor< 1, dim > > incr
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
void set_h(const double h)