1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2021 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
#ifndef dealii_function_restriction_h
#define dealii_function_restriction_h
#include <deal.II/base/function.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
DEAL_II_NAMESPACE_OPEN
namespace Functions
{
/**
* This class takes a function in `dim + 1` dimensions and creates a new
* function in one dimension lower by restricting one of the coordinates to a
* given value. Mathematically this corresponds to taking a function
* $f = f(x, y, z)$,
* a fixed value, $Z$, and defining a new function (the restriction)
* $g = g(x, y) = f(x, y, Z)$.
* Using this class, this translates to
* @code
* Function<3> & function = ...
* double z = ...
* unsigned int restricted_direction = 2;
* CoordinateRestriction<2> restriction(function, restricted_direction, z);
* @endcode
*
* The `dim`-dimensional coordinates on the restriction are ordered starting
* from the restricted (`dim + 1`)-coordinate. In particular, this means that
* if the $y$-coordinate is locked to $Y$ in 3d, the coordinates are ordered
* as $(z, x)$ on the restriction:
* $g = g(z, x) = f(x, Y, z)$.
* This is the same convention as in BoundingBox::cross_section.
*/
template <int dim>
class CoordinateRestriction : public Function<dim>
{
public:
/**
* Constructor, takes the (`dim + 1`)-coordinate direction and the value
* that the incoming function should be restricted to.
*
* A pointer to the incoming function is stored internally, so the function
* must have a longer lifetime than the created restriction.
*/
CoordinateRestriction(const Function<dim + 1> &function,
const unsigned int direction,
const double coordinate_value);
double
value(const Point<dim> &point, const unsigned int component) const override;
Tensor<1, dim>
gradient(const Point<dim> &point,
const unsigned int component) const override;
SymmetricTensor<2, dim>
hessian(const Point<dim> &point,
const unsigned int component) const override;
private:
// The higher-dimensional function that has been restricted.
const ObserverPointer<const Function<dim + 1>> function;
// The (`dim + 1`)-coordinate direction that has been restricted.
const unsigned int restricted_direction;
// Value of the restricted coordinate.
const double coordinate_value;
};
/**
* This class creates a 1-dimensional function from a `dim + 1` dimensional
* function by restricting `dim` of the coordinate values to a given point.
* Mathematically this corresponds to taking a function, $f = f(x, y, z)$, and
* a point $(Y, Z)$, and defining a new function $g = g(x) = f(x, Y, Z)$.
* Using this class, this translates to
* @code
* Function<3> & function = ...
* Point<2> point(y, z);
* unsigned int open_direction = 0;
* PointRestriction<2> restriction(function, open_direction, point);
* @endcode
*
* The coordinates of the point will be expanded in the higher-dimensional
* functions coordinates starting from the open-direction (and wrapping
* around). In particular, if we restrict to a point $(Z, X)$ and choose to
* keep the y-direction open, the restriction that is created is the function
* $g(y) = f(X, y, Z)$.
* This is consistent with the convention in BoundingBox::cross_section.
*/
template <int dim>
class PointRestriction : public Function<1>
{
public:
/**
* Constructor, takes the point that the incoming function should be
* restricted to and which (`dim + 1`)-dimensional coordinate direction
* should be kept "open".
*
* A pointer to the incoming function is stored internally, so the function
* must have a longer lifetime than the created restriction.
*/
PointRestriction(const Function<dim + 1> &function,
const unsigned int open_direction,
const Point<dim> &point);
double
value(const Point<1> &point, const unsigned int component) const override;
Tensor<1, 1>
gradient(const Point<1> &point,
const unsigned int component) const override;
SymmetricTensor<2, 1>
hessian(const Point<1> &point, const unsigned int component) const override;
private:
// The higher-dimensional function that has been restricted.
const ObserverPointer<const Function<dim + 1>> function;
// The (`dim + 1`)-coordinate direction that is kept "open"
const unsigned int open_direction;
// The point that we have restricted the above function to.
const Point<dim> point;
};
} // namespace Functions
namespace internal
{
/**
* Creates a (`dim + 1`)-dimensional point by copying over the coordinates of
* the incoming `dim`-dimensional point and setting the "missing"
* (`dim + 1`)-dimensional component to the incoming coordinate value.
*
* For example, given the input
* $\{(x, y), 2, z \}$ this function creates the point $(x, y, z)$.
*
* The coordinates of the `dim`-dimensional point are written to the
* coordinates of the (`dim + 1`)-dimensional point in the order of the
* convention given by the function coordinate_to_one_dim_higher. Thus, the
* order of coordinates on the lower-dimensional point are not preserved:
* $\{(z, x), 1, y \}$ creates the point $(x, y, z)$.
*/
template <int dim>
Point<dim + 1>
create_higher_dim_point(const Point<dim> &point,
const unsigned int component_in_dim_plus_1,
const double coordinate_value);
} // namespace internal
DEAL_II_NAMESPACE_CLOSE
#endif /* dealii_function_restriction_h */
|
