deal.II version GIT relicensing-2865-ga6be64acbe 2025-03-19 21: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#include <vector>
22
23
25
26template <int dim>
28 const Point<dim> &dir,
29 const double h)
30 : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
31 , f(f)
32 , h(h)
33 , incr(1, h * dir)
34{
36}
37
38
39
40template <int dim>
42 const std::vector<Point<dim>> &dir,
43 const double h)
44 : AutoDerivativeFunction<dim>(h, f.n_components, f.get_time())
45 , f(f)
46 , h(h)
47 , incr(dir.size())
48{
49 for (unsigned int i = 0; i < incr.size(); ++i)
50 incr[i] = h * dir[i];
52}
53
54
55
56template <int dim>
57void
60{
61 // go through all known formulas, reject ones we don't know about
62 // and don't handle in the member functions of this class
63 switch (form)
64 {
68 break;
69 default:
70 Assert(false,
71 ExcMessage("The argument passed to this function does not "
72 "match any known difference formula."));
73 }
74
75 formula = form;
76}
77
78
79
80template <int dim>
81void
83{
84 for (unsigned int i = 0; i < incr.size(); ++i)
85 incr[i] *= new_h / h;
86 h = new_h;
87}
88
89
90
91template <int dim>
92double
94 const unsigned int component) const
95{
96 Assert(incr.size() == 1,
98 "FunctionDerivative was not initialized for constant direction"));
99
100 switch (formula)
101 {
103 return (f.value(p + incr[0], component) -
104 f.value(p - incr[0], component)) /
105 (2 * h);
107 return (f.value(p, component) - f.value(p - incr[0], component)) / h;
109 return (-f.value(p + 2 * incr[0], component) +
110 8 * f.value(p + incr[0], component) -
111 8 * f.value(p - incr[0], component) +
112 f.value(p - 2 * incr[0], component)) /
113 (12 * h);
114 default:
116 }
117 return 0.;
118}
119
120
121
122template <int dim>
123void
125 Vector<double> &result) const
126{
127 Assert(incr.size() == 1,
129 "FunctionDerivative was not initialized for constant direction"));
130 Vector<double> aux(result.size());
131
132 // Formulas are the same as in
133 // value, but here we have to use
134 // Vector arithmetic
135 switch (formula)
136 {
138 f.vector_value(p + incr[0], result);
139 f.vector_value(p - incr[0], aux);
140 result.sadd(1. / (2 * h), -1. / (2 * h), aux);
141 return;
143 f.vector_value(p, result);
144 f.vector_value(p - incr[0], aux);
145 result.sadd(1. / h, -1. / h, aux);
146 return;
148 f.vector_value(p - 2 * incr[0], result);
149 f.vector_value(p + 2 * incr[0], aux);
150 result.add(-1., aux);
151 f.vector_value(p - incr[0], aux);
152 result.add(-8., aux);
153 f.vector_value(p + incr[0], aux);
154 result.add(8., aux);
155 result /= (12. * h);
156 return;
157 default:
159 }
160}
161
162
163
164template <int dim>
165void
167 std::vector<double> &values,
168 const unsigned int component) const
169{
170 const unsigned int n = points.size();
171 const bool variable_direction = (incr.size() == 1) ? false : true;
172 if (variable_direction)
173 Assert(incr.size() == points.size(),
174 ExcDimensionMismatch(incr.size(), points.size()));
175
176 // Vector of auxiliary values
177 std::vector<double> aux(n);
178 // Vector of auxiliary points
179 std::vector<Point<dim>> paux(n);
180 // Use the same formulas as in
181 // value, but with vector
182 // arithmetic
183 switch (formula)
184 {
186 for (unsigned int i = 0; i < n; ++i)
187 paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
188 f.value_list(paux, values, component);
189 for (unsigned int i = 0; i < n; ++i)
190 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
191 f.value_list(paux, aux, component);
192 for (unsigned int i = 0; i < n; ++i)
193 values[i] = (values[i] - aux[i]) / (2 * h);
194 return;
196 f.value_list(points, values, component);
197 for (unsigned int i = 0; i < n; ++i)
198 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
199 f.value_list(paux, aux, component);
200 for (unsigned int i = 0; i < n; ++i)
201 values[i] = (values[i] - aux[i]) / h;
202 return;
204 for (unsigned int i = 0; i < n; ++i)
205 paux[i] = points[i] - 2 * incr[(variable_direction) ? i : 0U];
206 f.value_list(paux, values, component);
207 for (unsigned int i = 0; i < n; ++i)
208 paux[i] = points[i] + 2 * incr[(variable_direction) ? i : 0U];
209 f.value_list(paux, aux, component);
210 for (unsigned int i = 0; i < n; ++i)
211 values[i] -= aux[i];
212 for (unsigned int i = 0; i < n; ++i)
213 paux[i] = points[i] + incr[(variable_direction) ? i : 0U];
214 f.value_list(paux, aux, component);
215 for (unsigned int i = 0; i < n; ++i)
216 values[i] += 8. * aux[i];
217 for (unsigned int i = 0; i < n; ++i)
218 paux[i] = points[i] - incr[(variable_direction) ? i : 0U];
219 f.value_list(paux, aux, component);
220 for (unsigned int i = 0; i < n; ++i)
221 values[i] = (values[i] - 8. * aux[i]) / (12 * h);
222 return;
223 default:
225 }
226}
227
228
229
230template <int dim>
231std::size_t
233{
234 // only simple data elements, so
235 // use sizeof operator
236 return sizeof(*this);
237}
238
239
240
241// explicit instantiations
242template class FunctionDerivative<1>;
243template class FunctionDerivative<2>;
244template class FunctionDerivative<3>;
245
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:113
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:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::size_t size
Definition mpi.cc:745