Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
function_derivative.cc
Go to the documentation of this file.
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
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 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{
35}
36
37
38
39template <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];
51}
52
53
54
55template <int dim>
56void
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
79template <int dim>
80void
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
90template <int dim>
91double
93 const unsigned int component) const
94{
95 Assert(incr.size() == 1,
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
121template <int dim>
122void
124 Vector<double> & result) const
125{
126 Assert(incr.size() == 1,
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
163template <int dim>
164void
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
229template <int dim>
230std::size_t
232{
233 // only simple data elements, so
234 // use sizeof operator
235 return sizeof(*this);
236}
237
238
239
240// explicit instantiations
241template class FunctionDerivative<1>;
242template class FunctionDerivative<2>;
243template class FunctionDerivative<3>;
244
virtual std::size_t memory_consumption() const override
virtual void vector_value(const Point< dim > &p, Vector< double > &value) const override
void set_h(const double h)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
FunctionDerivative(const Function< dim > &f, const Point< dim > &direction, const double h=1.e-6)
virtual double value(const Point< dim > &p, const unsigned int component=0) const override
std::vector< Tensor< 1, dim > > incr
void set_formula(typename AutoDerivativeFunction< dim >::DifferenceFormula formula=AutoDerivativeFunction< dim >::Euler)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Definition point.h:112
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
void sadd(const Number s, const Vector< Number > &V)
size_type size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)