Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 : #include <tuple>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Index.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "DataStructures/Variables.hpp"
14 : #include "Domain/Structure/Direction.hpp"
15 : #include "Domain/Structure/DirectionalId.hpp"
16 : #include "Domain/Structure/DirectionalIdMap.hpp"
17 : #include "Domain/Structure/ElementId.hpp"
18 : #include "Utilities/ContainerHelpers.hpp"
19 : #include "Utilities/ErrorHandling/Assert.hpp"
20 :
21 : namespace LinearSolver::Schwarz {
22 :
23 : /// Identifies a subdomain region that overlaps with another element
24 : template <size_t Dim>
25 1 : using OverlapId = DirectionalId<Dim>;
26 :
27 : /// Data structure that can store the `ValueType` on each possible overlap of an
28 : /// element-centered subdomain with its neighbors. Overlaps are identified by
29 : /// their `OverlapId`.
30 : template <size_t Dim, typename ValueType>
31 1 : using OverlapMap = DirectionalIdMap<Dim, ValueType>;
32 :
33 : /*!
34 : * \brief The number of points that an overlap extends into the `volume_extent`
35 : *
36 : * In a dimension where an element has `volume_extent` points, the overlap
37 : * extent is the largest number under these constraints:
38 : *
39 : * - It is at most `max_overlap`.
40 : * - It is smaller than the `volume_extent`.
41 : *
42 : * This means the overlap extent is always smaller than the `volume_extent`. The
43 : * reason for this constraint is that we define the _width_ of the overlap as
44 : * the element-logical coordinate distance from the face of the element to the
45 : * first collocation point _outside_ the overlap extent. Therefore, even an
46 : * overlap region that covers the full element in width does not include the
47 : * collocation point on the opposite side of the element.
48 : *
49 : * Here's a few notes on the definition of the overlap extent and width:
50 : *
51 : * - A typical smooth weighting function goes to zero at the overlap width, so
52 : * if the grid points located at the overlap width were included in the
53 : * subdomain, their solutions would not contribute to the weighted sum of
54 : * subdomain solutions.
55 : * - Defining the overlap width as the distance to the first point _outside_ the
56 : * overlap extent makes it non-zero even for a single point of overlap into a
57 : * Gauss-Lobatto grid (which has points located at the element face).
58 : * - Boundary contributions for many (but not all) discontinuous Galerkin
59 : * schemes on Gauss-Lobatto grids are limited to the grid points on the
60 : * element face, e.g. for a DG operator that is pre-multiplied by the mass
61 : * matrix, or one where boundary contributions are lifted using the diagonal
62 : * mass-matrix approximation. Not including the grid points facing away from
63 : * the subdomain in the overlap allows to ignore that face altogether in the
64 : * subdomain operator.
65 : */
66 1 : size_t overlap_extent(size_t volume_extent, size_t max_overlap);
67 :
68 : /*!
69 : * \brief Total number of grid points in an overlap region that extends
70 : * `overlap_extent` points into the `volume_extents` from either side in the
71 : * `overlap_dimension`
72 : *
73 : * The overlap region has `overlap_extent` points in the `overlap_dimension`,
74 : * and `volume_extents` points in the other dimensions. The number of grid
75 : * points returned by this function is the product of these extents.
76 : */
77 : template <size_t Dim>
78 1 : size_t overlap_num_points(const Index<Dim>& volume_extents,
79 : size_t overlap_extent, size_t overlap_dimension);
80 :
81 : /*!
82 : * \brief Width of an overlap extending `overlap_extent` points into the
83 : * `collocation_points` from either side.
84 : *
85 : * The "width" of an overlap is the element-logical coordinate distance from the
86 : * element boundary to the first collocation point outside the overlap region in
87 : * the overlap dimension, i.e. the dimension perpendicular to the element face.
88 : * See `LinearSolver::Schwarz::overlap_extent` for details.
89 : *
90 : * This function assumes the `collocation_points` are mirrored around 0.
91 : */
92 1 : double overlap_width(size_t overlap_extent,
93 : const DataVector& collocation_points);
94 :
95 : /*!
96 : * \brief Iterate over grid points in a region that extends partially into the
97 : * volume
98 : *
99 : * Here's an example how to use this iterator:
100 : *
101 : * \snippet Test_OverlapHelpers.cpp overlap_iterator
102 : */
103 1 : class OverlapIterator {
104 : public:
105 : template <size_t Dim>
106 0 : OverlapIterator(const Index<Dim>& volume_extents, size_t overlap_extent,
107 : const Direction<Dim>& direction);
108 :
109 0 : explicit operator bool() const;
110 :
111 0 : OverlapIterator& operator++();
112 :
113 : /// Offset into a DataVector that holds full volume data
114 1 : size_t volume_offset() const;
115 :
116 : /// Offset into a DataVector that holds data only on the overlap region
117 1 : size_t overlap_offset() const;
118 :
119 0 : void reset();
120 :
121 : private:
122 0 : size_t size_ = std::numeric_limits<size_t>::max();
123 0 : size_t num_slices_ = std::numeric_limits<size_t>::max();
124 0 : size_t stride_ = std::numeric_limits<size_t>::max();
125 0 : size_t stride_count_ = std::numeric_limits<size_t>::max();
126 0 : size_t jump_ = std::numeric_limits<size_t>::max();
127 0 : size_t initial_offset_ = std::numeric_limits<size_t>::max();
128 0 : size_t volume_offset_ = std::numeric_limits<size_t>::max();
129 0 : size_t overlap_offset_ = std::numeric_limits<size_t>::max();
130 : };
131 :
132 : /// @{
133 : /// The part of the tensor data that lies within the overlap region
134 : template <size_t Dim, typename DataType, typename... TensorStructure>
135 1 : void data_on_overlap(const gsl::not_null<Tensor<DataType, TensorStructure...>*>
136 : restricted_tensor,
137 : const Tensor<DataType, TensorStructure...>& tensor,
138 : const Index<Dim>& volume_extents,
139 : const size_t overlap_extent,
140 : const Direction<Dim>& direction) {
141 : for (OverlapIterator overlap_iterator{volume_extents, overlap_extent,
142 : direction};
143 : overlap_iterator; ++overlap_iterator) {
144 : for (size_t tensor_component = 0; tensor_component < tensor.size();
145 : ++tensor_component) {
146 : (*restricted_tensor)[tensor_component][overlap_iterator
147 : .overlap_offset()] =
148 : tensor[tensor_component][overlap_iterator.volume_offset()];
149 : }
150 : }
151 : }
152 :
153 : template <size_t Dim, typename DataType, typename... TensorStructure>
154 1 : Tensor<DataType, TensorStructure...> data_on_overlap(
155 : const Tensor<DataType, TensorStructure...>& tensor,
156 : const Index<Dim>& volume_extents, const size_t overlap_extent,
157 : const Direction<Dim>& direction) {
158 : Tensor<DataType, TensorStructure...> restricted_tensor{overlap_num_points(
159 : volume_extents, overlap_extent, direction.dimension())};
160 : data_on_overlap(make_not_null(&restricted_tensor), tensor, volume_extents,
161 : overlap_extent, direction);
162 : return restricted_tensor;
163 : }
164 :
165 : namespace detail {
166 : template <size_t Dim>
167 : void data_on_overlap_impl(double* overlap_data, const double* volume_data,
168 : size_t num_components,
169 : const Index<Dim>& volume_extents,
170 : size_t overlap_extent,
171 : const Direction<Dim>& direction);
172 : } // namespace detail
173 :
174 : template <size_t Dim, typename OverlapTagsList, typename VolumeTagsList>
175 1 : void data_on_overlap(
176 : const gsl::not_null<Variables<OverlapTagsList>*> overlap_data,
177 : const Variables<VolumeTagsList>& volume_data,
178 : const Index<Dim>& volume_extents, const size_t overlap_extent,
179 : const Direction<Dim>& direction) {
180 : constexpr size_t num_components =
181 : Variables<VolumeTagsList>::number_of_independent_components;
182 : ASSERT(volume_data.number_of_grid_points() == volume_extents.product(),
183 : "volume_data has wrong number of grid points. Expected "
184 : << volume_extents.product() << ", got "
185 : << volume_data.number_of_grid_points());
186 : ASSERT(overlap_data->number_of_grid_points() ==
187 : overlap_num_points(volume_extents, overlap_extent,
188 : direction.dimension()),
189 : "overlap_data has wrong number of grid points. Expected "
190 : << overlap_num_points(volume_extents, overlap_extent,
191 : direction.dimension())
192 : << ", got " << overlap_data->number_of_grid_points());
193 : detail::data_on_overlap_impl(overlap_data->data(), volume_data.data(),
194 : num_components, volume_extents, overlap_extent,
195 : direction);
196 : }
197 :
198 : template <size_t Dim, typename TagsList>
199 1 : Variables<TagsList> data_on_overlap(const Variables<TagsList>& volume_data,
200 : const Index<Dim>& volume_extents,
201 : const size_t overlap_extent,
202 : const Direction<Dim>& direction) {
203 : Variables<TagsList> overlap_data{overlap_num_points(
204 : volume_extents, overlap_extent, direction.dimension())};
205 : data_on_overlap(make_not_null(&overlap_data), volume_data, volume_extents,
206 : overlap_extent, direction);
207 : return overlap_data;
208 : }
209 : /// @}
210 :
211 : namespace detail {
212 : template <size_t Dim>
213 : void add_overlap_data_impl(double* volume_data, const double* overlap_data,
214 : size_t num_components,
215 : const Index<Dim>& volume_extents,
216 : size_t overlap_extent,
217 : const Direction<Dim>& direction);
218 : } // namespace detail
219 :
220 : /// Add the `overlap_data` to the `volume_data`
221 : template <size_t Dim, typename VolumeTagsList, typename OverlapTagsList>
222 1 : void add_overlap_data(
223 : const gsl::not_null<Variables<VolumeTagsList>*> volume_data,
224 : const Variables<OverlapTagsList>& overlap_data,
225 : const Index<Dim>& volume_extents, const size_t overlap_extent,
226 : const Direction<Dim>& direction) {
227 : constexpr size_t num_components =
228 : Variables<VolumeTagsList>::number_of_independent_components;
229 : ASSERT(volume_data->number_of_grid_points() == volume_extents.product(),
230 : "volume_data has wrong number of grid points. Expected "
231 : << volume_extents.product() << ", got "
232 : << volume_data->number_of_grid_points());
233 : ASSERT(overlap_data.number_of_grid_points() ==
234 : overlap_num_points(volume_extents, overlap_extent,
235 : direction.dimension()),
236 : "overlap_data has wrong number of grid points. Expected "
237 : << overlap_num_points(volume_extents, overlap_extent,
238 : direction.dimension())
239 : << ", got " << overlap_data.number_of_grid_points());
240 : detail::add_overlap_data_impl(volume_data->data(), overlap_data.data(),
241 : num_components, volume_extents, overlap_extent,
242 : direction);
243 : }
244 :
245 : /// @{
246 : /// Extend the overlap data to the full mesh by filling it with zeros outside
247 : /// the overlap region
248 : template <size_t Dim, typename ExtendedTagsList, typename OverlapTagsList>
249 1 : void extended_overlap_data(
250 : const gsl::not_null<Variables<ExtendedTagsList>*> extended_data,
251 : const Variables<OverlapTagsList>& overlap_data,
252 : const Index<Dim>& volume_extents, const size_t overlap_extent,
253 : const Direction<Dim>& direction) {
254 : *extended_data = Variables<ExtendedTagsList>{volume_extents.product(), 0.};
255 : add_overlap_data(extended_data, overlap_data, volume_extents, overlap_extent,
256 : direction);
257 : }
258 :
259 : template <size_t Dim, typename TagsList>
260 1 : Variables<TagsList> extended_overlap_data(
261 : const Variables<TagsList>& overlap_data, const Index<Dim>& volume_extents,
262 : const size_t overlap_extent, const Direction<Dim>& direction) {
263 : Variables<TagsList> extended_data{volume_extents.product()};
264 : extended_overlap_data(make_not_null(&extended_data), overlap_data,
265 : volume_extents, overlap_extent, direction);
266 : return extended_data;
267 : }
268 : /// @}
269 :
270 : } // namespace LinearSolver::Schwarz
|
