Reference documentation for deal.II version GIT relicensing-12-g89ba5bf9db 2024-02-24 00:00:03+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
bounding_box.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2023 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
17
18#include <limits>
19#include <numeric>
20
22
23template <int spacedim, typename Number>
24bool
26 const double tolerance) const
27{
28 for (unsigned int i = 0; i < spacedim; ++i)
29 {
30 // Bottom left-top right convention: the point is outside if it's smaller
31 // than the first or bigger than the second boundary point The bounding
32 // box is defined as a closed set
33 if ((p[i] <
34 this->boundary_points.first[i] - tolerance * side_length(i)) ||
35 (p[i] > this->boundary_points.second[i] + tolerance * side_length(i)))
36 return false;
37 }
38 return true;
39}
40
41
42
43template <int spacedim, typename Number>
44void
46 const BoundingBox<spacedim, Number> &other_bbox)
47{
48 for (unsigned int i = 0; i < spacedim; ++i)
49 {
50 this->boundary_points.first[i] =
51 std::min(this->boundary_points.first[i],
52 other_bbox.boundary_points.first[i]);
53 this->boundary_points.second[i] =
54 std::max(this->boundary_points.second[i],
55 other_bbox.boundary_points.second[i]);
56 }
57}
58
59template <int spacedim, typename Number>
60bool
62 const BoundingBox<spacedim, Number> &other_bbox,
63 const double tolerance) const
64{
65 for (unsigned int i = 0; i < spacedim; ++i)
66 {
67 // testing if the boxes are close enough to intersect
68 if ((other_bbox.boundary_points.second[i] <
69 this->boundary_points.first[i] - tolerance * side_length(i)) ||
70 (other_bbox.boundary_points.first[i] >
71 this->boundary_points.second[i] + tolerance * side_length(i)))
72 return false;
73 }
74 return true;
75}
76
77template <int spacedim, typename Number>
80 const BoundingBox<spacedim, Number> &other_bbox,
81 const double tolerance) const
82{
83 if (!has_overlap_with(other_bbox, tolerance))
85
86 if (spacedim == 1)
87 {
88 // In dimension 1 if the two bounding box are neighbors
89 // we can merge them
91 }
92 else
93 {
94 const std::array<Point<spacedim, Number>, 2> bbox1 = {
95 {this->get_boundary_points().first,
96 this->get_boundary_points().second}};
97 const std::array<Point<spacedim, Number>, 2> bbox2 = {
98 {other_bbox.get_boundary_points().first,
99 other_bbox.get_boundary_points().second}};
100
101 // The boxes intersect: we need to understand now how they intersect.
102 // We begin by computing the intersection:
103 std::array<double, spacedim> intersect_bbox_min;
104 std::array<double, spacedim> intersect_bbox_max;
105 for (unsigned int d = 0; d < spacedim; ++d)
106 {
107 intersect_bbox_min[d] = std::max(bbox1[0][d], bbox2[0][d]);
108 intersect_bbox_max[d] = std::min(bbox1[1][d], bbox2[1][d]);
109 }
110
111 // Finding the intersection's dimension
112 int intersect_dim = spacedim;
113 for (unsigned int d = 0; d < spacedim; ++d)
114 if (std::abs(intersect_bbox_min[d] - intersect_bbox_max[d]) <=
115 tolerance * (std::abs(intersect_bbox_min[d]) +
116 std::abs(intersect_bbox_max[d])))
117 --intersect_dim;
118
119 if (intersect_dim == 0 || intersect_dim == spacedim - 2)
121
122 // Checking the two mergeable cases: first if the boxes are aligned so
123 // that they can be merged
124 unsigned int not_align_1 = 0, not_align_2 = 0;
125 bool same_direction = true;
126 for (unsigned int d = 0; d < spacedim; ++d)
127 {
128 if (std::abs(bbox2[0][d] - bbox1[0][d]) >
129 tolerance * (std::abs(bbox2[0][d]) + std::abs(bbox1[0][d])))
130 ++not_align_1;
131 if (std::abs(bbox1[1][d] - bbox2[1][d]) >
132 tolerance * (std::abs(bbox1[1][d]) + std::abs(bbox2[1][d])))
133 ++not_align_2;
134 if (not_align_1 != not_align_2)
135 {
136 same_direction = false;
137 break;
138 }
139 }
140
141 if (not_align_1 <= 1 && not_align_2 <= 1 && same_direction)
143
144 // Second: one box is contained/equal to the other
145 if ((this->point_inside(bbox2[0]) && this->point_inside(bbox2[1])) ||
146 (other_bbox.point_inside(bbox1[0], tolerance) &&
147 other_bbox.point_inside(bbox1[1], tolerance)))
149
150 // Degenerate and mergeable cases have been found, it remains:
152 }
153}
154
155
156
157template <int spacedim, typename Number>
158double
160{
161 double vol = 1.0;
162 for (unsigned int i = 0; i < spacedim; ++i)
163 vol *= (this->boundary_points.second[i] - this->boundary_points.first[i]);
164 return vol;
165}
166
167
168
169template <int spacedim, typename Number>
170Number
171BoundingBox<spacedim, Number>::lower_bound(const unsigned int direction) const
172{
173 AssertIndexRange(direction, spacedim);
174
175 return boundary_points.first[direction];
176}
177
178
179
180template <int spacedim, typename Number>
181Number
182BoundingBox<spacedim, Number>::upper_bound(const unsigned int direction) const
183{
184 AssertIndexRange(direction, spacedim);
185
186 return boundary_points.second[direction];
187}
188
189
190
191template <int spacedim, typename Number>
194{
196 for (unsigned int i = 0; i < spacedim; ++i)
197 point[i] = .5 * (boundary_points.first[i] + boundary_points.second[i]);
198
199 return point;
200}
201
202
203
204template <int spacedim, typename Number>
206BoundingBox<spacedim, Number>::bounds(const unsigned int direction) const
208 AssertIndexRange(direction, spacedim);
209
210 std::pair<Point<1, Number>, Point<1, Number>> lower_upper_bounds;
211 lower_upper_bounds.first[0] = lower_bound(direction);
212 lower_upper_bounds.second[0] = upper_bound(direction);
213
214 return BoundingBox<1, Number>(lower_upper_bounds);
216
217
218
219template <int spacedim, typename Number>
220Number
221BoundingBox<spacedim, Number>::side_length(const unsigned int direction) const
222{
223 AssertIndexRange(direction, spacedim);
224
225 return boundary_points.second[direction] - boundary_points.first[direction];
226}
227
228
229
230template <int spacedim, typename Number>
232BoundingBox<spacedim, Number>::vertex(const unsigned int index) const
233{
235
236 const Point<spacedim> unit_cell_vertex =
238
240 for (unsigned int i = 0; i < spacedim; ++i)
241 point[i] = boundary_points.first[i] + side_length(i) * unit_cell_vertex[i];
242
243 return point;
244}
245
246
247
248template <int spacedim, typename Number>
250BoundingBox<spacedim, Number>::child(const unsigned int index) const
251{
253
254 // Vertex closest to child.
255 const Point<spacedim, Number> parent_vertex = vertex(index);
256 const Point<spacedim, Number> parent_center = center();
257
258 const Point<spacedim> upper_corner_unit_cell =
261
262 const Point<spacedim> lower_corner_unit_cell =
264
265 std::pair<Point<spacedim, Number>, Point<spacedim, Number>>
266 child_lower_upper_corner;
267 for (unsigned int i = 0; i < spacedim; ++i)
268 {
269 const double child_side_length = side_length(i) / 2;
270
271 const double child_center = (parent_center[i] + parent_vertex[i]) / 2;
272
273 child_lower_upper_corner.first[i] =
274 child_center + child_side_length * (lower_corner_unit_cell[i] - .5);
275 child_lower_upper_corner.second[i] =
276 child_center + child_side_length * (upper_corner_unit_cell[i] - .5);
277 }
279 return BoundingBox<spacedim, Number>(child_lower_upper_corner);
280}
281
282
283
284template <int spacedim, typename Number>
285BoundingBox<spacedim - 1, Number>
286BoundingBox<spacedim, Number>::cross_section(const unsigned int direction) const
287{
288 AssertIndexRange(direction, spacedim);
289
290 std::pair<Point<spacedim - 1, Number>, Point<spacedim - 1, Number>>
291 cross_section_lower_upper_corner;
292 for (unsigned int d = 0; d < spacedim - 1; ++d)
293 {
294 const int index_to_write_from =
295 internal::coordinate_to_one_dim_higher<spacedim - 1>(direction, d);
297 cross_section_lower_upper_corner.first[d] =
298 boundary_points.first[index_to_write_from];
299
300 cross_section_lower_upper_corner.second[d] =
301 boundary_points.second[index_to_write_from];
303
304 return BoundingBox<spacedim - 1, Number>(cross_section_lower_upper_corner);
305}
306
307
309template <int spacedim, typename Number>
312 const Point<spacedim, Number> &point) const
313{
314 auto unit = point;
315 const auto diag = boundary_points.second - boundary_points.first;
316 unit -= boundary_points.first;
317 for (unsigned int d = 0; d < spacedim; ++d)
318 unit[d] /= diag[d];
319 return unit;
320}
321
323
324template <int spacedim, typename Number>
327 const Point<spacedim, Number> &point) const
328{
329 auto real = boundary_points.first;
330 const auto diag = boundary_points.second - boundary_points.first;
331 for (unsigned int d = 0; d < spacedim; ++d)
332 real[d] += diag[d] * point[d];
333 return real;
334}
335
336
337
338template <int spacedim, typename Number>
339Number
341 const Point<spacedim, Number> &point,
342 const unsigned int direction) const
344 const Number p1 = lower_bound(direction);
345 const Number p2 = upper_bound(direction);
346
347 if (point[direction] > p2)
348 return point[direction] - p2;
349 else if (point[direction] < p1)
350 return p1 - point[direction];
351 else
352 return -std::min(point[direction] - p1, p2 - point[direction]);
353}
355
356
357template <int spacedim, typename Number>
358Number
360 const Point<spacedim, Number> &point) const
361{
362 // calculate vector of orthogonal signed distances
363 std::array<Number, spacedim> distances;
364 for (unsigned int d = 0; d < spacedim; ++d)
365 distances[d] = signed_distance(point, d);
366
367 // determine the number of positive signed distances
368 const unsigned int n_positive_signed_distances =
369 std::count_if(distances.begin(), distances.end(), [](const auto &a) {
370 return a > 0.0;
371 });
373 if (n_positive_signed_distances <= 1)
374 // point is inside of bounding box (0: all signed distances are
375 // negative; find the index with the smallest absolute value)
376 // or next to a face (1: all signed distances are negative
377 // but one; find this index)
378 return *std::max_element(distances.begin(), distances.end());
379 else
380 // point is next to a corner (2D/3D: all signed distances are
381 // positive) or a line (3D: all signed distances are positive
382 // but one) -> take the l2-norm of all positive signed distances
383 return std::sqrt(std::accumulate(distances.begin(),
384 distances.end(),
385 0.0,
386 [](const auto &a, const auto &b) {
387 return a + (b > 0 ? b * b : 0.0);
388 }));
389}
390
391
392
393template <int dim, typename Number>
396{
397 std::pair<Point<dim, Number>, Point<dim, Number>> lower_upper_corner;
398 for (unsigned int i = 0; i < dim; ++i)
399 {
400 lower_upper_corner.second[i] = 1;
401 }
402 return BoundingBox<dim, Number>(lower_upper_corner);
403}
404
405
406#include "bounding_box.inst"
NeighborType
BoundingBox< dim, Number > create_unit_bounding_box()
std::pair< Point< spacedim, Number >, Point< spacedim, Number > > boundary_points
BoundingBox< 1, Number > bounds(const unsigned int direction) const
bool has_overlap_with(const BoundingBox< spacedim, Number > &other_bbox, const double tolerance=std::numeric_limits< Number >::epsilon()) const
Point< spacedim, Number > center() const
Number lower_bound(const unsigned int direction) const
Number signed_distance(const Point< spacedim, Number > &point, const unsigned int direction) const
void merge_with(const BoundingBox< spacedim, Number > &other_bbox)
std::pair< Point< spacedim, Number >, Point< spacedim, Number > > & get_boundary_points()
bool point_inside(const Point< spacedim, Number > &p, const double tolerance=std::numeric_limits< Number >::epsilon()) const
double volume() const
Point< spacedim, Number > real_to_unit(const Point< spacedim, Number > &point) const
Number side_length(const unsigned int direction) const
BoundingBox< spacedim, Number > child(const unsigned int index) const
NeighborType get_neighbor_type(const BoundingBox< spacedim, Number > &other_bbox, const double tolerance=std::numeric_limits< Number >::epsilon()) const
Point< spacedim, Number > vertex(const unsigned int index) const
BoundingBox< spacedim - 1, Number > cross_section(const unsigned int direction) const
Number upper_bound(const unsigned int direction) const
Point< spacedim, Number > unit_to_real(const Point< spacedim, Number > &point) const
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > center
#define AssertIndexRange(index, range)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static Point< dim > unit_cell_vertex(const unsigned int vertex)