deal.II version GIT relicensing-2233-g9bd70d8389 2024-12-08 16:10:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 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 Point<dim> &dir,
27 const double h)
28 : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
29 , f(f)
30 , h(h)
31 , incr(1, h * dir)
32{
34}
35
36
37
38template <int dim>
40 const std::vector<Point<dim>> &dir,
41 const double h)
42 : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
43 , f(f)
44 , h(h)
45 , incr(dir.size())
46{
47 for (unsigned int i = 0; i < incr.size(); ++i)
48 incr[i] = h * dir[i];
50}
51
52
53
54template <int dim>
55void
58{
59 // go through all known formulas, reject ones we don't know about
60 // and don't handle in the member functions of this class
61 switch (form)
62 {
66 break;
67 default:
68 Assert(false,
69 ExcMessage("The argument passed to this function does not "
70 "match any known difference formula."));
71 }
72
73 formula = form;
74}
75
76
77
78template <int dim>
79void
81{
82 for (unsigned int i = 0; i < incr.size(); ++i)
83 incr[i] *= new_h / h;
84 h = new_h;
85}
86
87
88
89template <int dim>
90double
92 const unsigned int component) const
93{
94 Assert(incr.size() == 1,
96 "FunctionDerivative was not initialized for constant direction"));
97
98 switch (formula)
99 {
101 return (f.value(p + incr[0], component) -
102 f.value(p - incr[0], component)) /
103 (2 * h);
105 return (f.value(p, component) - f.value(p - incr[0], component)) / h;
107 return (-f.value(p + 2 * incr[0], component) +
108 8 * f.value(p + incr[0], component) -
109 8 * f.value(p - incr[0], component) +
110 f.value(p - 2 * incr[0], component)) /
111 (12 * h);
112 default:
114 }
115 return 0.;
116}
117
118
119
120template <int dim>
121void
123 Vector<double> &result) const
124{
125 Assert(incr.size() == 1,
127 "FunctionDerivative was not initialized for constant direction"));
128 Vector<double> aux(result.size());
129
130 // Formulas are the same as in
131 // value, but here we have to use
132 // Vector arithmetic
133 switch (formula)
134 {
136 f.vector_value(p + incr[0], result);
137 f.vector_value(p - incr[0], aux);
138 result.sadd(1. / (2 * h), -1. / (2 * h), aux);
139 return;
141 f.vector_value(p, result);
142 f.vector_value(p - incr[0], aux);
143 result.sadd(1. / h, -1. / h, aux);
144 return;
146 f.vector_value(p - 2 * incr[0], result);
147 f.vector_value(p + 2 * incr[0], aux);
148 result.add(-1., aux);
149 f.vector_value(p - incr[0], aux);
150 result.add(-8., aux);
151 f.vector_value(p + incr[0], aux);
152 result.add(8., aux);
153 result /= (12. * h);
154 return;
155 default:
157 }
158}
159
160
161
162template <int dim>
163void
165 std::vector<double> &values,
166 const unsigned int component) const
167{
168 const unsigned int n = points.size();
169 const bool variable_direction = (incr.size() == 1) ? false : true;
170 if (variable_direction)
171 Assert(incr.size() == points.size(),
172 ExcDimensionMismatch(incr.size(), points.size()));
173
174 // Vector of auxiliary values
175 std::vector<double> aux(n);
176 // Vector of auxiliary points
177 std::vector<Point<dim>> paux(n);
178 // Use the same formulas as in
179 // value, but with vector
180 // arithmetic
181 switch (formula)
182 {
184 for (unsigned int i = 0; i < n; ++i)
185 paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
186 f.value_list(paux, values, component);
187 for (unsigned int i = 0; i < n; ++i)
188 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
189 f.value_list(paux, aux, component);
190 for (unsigned int i = 0; i < n; ++i)
191 values[i] = (values[i] - aux[i]) / (2 * h);
192 return;
194 f.value_list(points, values, component);
195 for (unsigned int i = 0; i < n; ++i)
196 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
197 f.value_list(paux, aux, component);
198 for (unsigned int i = 0; i < n; ++i)
199 values[i] = (values[i] - aux[i]) / h;
200 return;
202 for (unsigned int i = 0; i < n; ++i)
203 paux[i] = points[i] - 2 * incr[(variable_direction) ? i : 0U];
204 f.value_list(paux, values, component);
205 for (unsigned int i = 0; i < n; ++i)
206 paux[i] = points[i] + 2 * incr[(variable_direction) ? i : 0U];
207 f.value_list(paux, aux, component);
208 for (unsigned int i = 0; i < n; ++i)
209 values[i] -= aux[i];
210 for (unsigned int i = 0; i < n; ++i)
211 paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
212 f.value_list(paux, aux, component);
213 for (unsigned int i = 0; i < n; ++i)
214 values[i] += 8. * aux[i];
215 for (unsigned int i = 0; i < n; ++i)
216 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
217 f.value_list(paux, aux, component);
218 for (unsigned int i = 0; i < n; ++i)
219 values[i] = (values[i] - 8. * aux[i]) / (12 * h);
220 return;
221 default:
223 }
224}
225
226
227
228template <int dim>
229std::size_t
231{
232 // only simple data elements, so
233 // use sizeof operator
234 return sizeof(*this);
235}
236
237
238
239// explicit instantiations
240template class FunctionDerivative<1>;
241template class FunctionDerivative<2>;
242template class FunctionDerivative<3>;
243
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:111
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
virtual size_type size() const override
void sadd(const Number s, const Vector< Number > &V)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#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()
std::size_t size
Definition mpi.cc:734