Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2017 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
18
19#include <limits>
20#include <numeric>
21
23
24template <int spacedim, typename Number>
25bool
27 const double tolerance) const
28{
29 for (unsigned int i = 0; i < spacedim; ++i)
30 {
31 // Bottom left-top right convention: the point is outside if it's smaller
32 // than the first or bigger than the second boundary point The bounding
33 // box is defined as a closed set
34 if ((p[i] <
35 this->boundary_points.first[i] - tolerance * side_length(i)) ||
36 (p[i] > this->boundary_points.second[i] + tolerance * side_length(i)))
37 return false;
38 }
39 return true;
40}
41
42
43
44template <int spacedim, typename Number>
45void
47 const BoundingBox<spacedim, Number> &other_bbox)
48{
49 for (unsigned int i = 0; i < spacedim; ++i)
50 {
51 this->boundary_points.first[i] =
52 std::min(this->boundary_points.first[i],
53 other_bbox.boundary_points.first[i]);
54 this->boundary_points.second[i] =
55 std::max(this->boundary_points.second[i],
56 other_bbox.boundary_points.second[i]);
57 }
58}
59
60template <int spacedim, typename Number>
61bool
63 const BoundingBox<spacedim, Number> &other_bbox,
64 const double tolerance) const
65{
66 for (unsigned int i = 0; i < spacedim; ++i)
67 {
68 // testing if the boxes are close enough to intersect
69 if ((other_bbox.boundary_points.second[i] <
70 this->boundary_points.first[i] - tolerance * side_length(i)) ||
71 (other_bbox.boundary_points.first[i] >
72 this->boundary_points.second[i] + tolerance * side_length(i)))
73 return false;
74 }
75 return true;
76}
77
78template <int spacedim, typename Number>
81 const BoundingBox<spacedim, Number> &other_bbox,
82 const double tolerance) const
83{
84 if (!has_overlap_with(other_bbox, tolerance))
86
87 if (spacedim == 1)
88 {
89 // In dimension 1 if the two bounding box are neighbors
90 // we can merge them
92 }
93 else
94 {
95 const std::array<Point<spacedim, Number>, 2> bbox1 = {
96 {this->get_boundary_points().first,
97 this->get_boundary_points().second}};
98 const std::array<Point<spacedim, Number>, 2> bbox2 = {
99 {other_bbox.get_boundary_points().first,
100 other_bbox.get_boundary_points().second}};
101
102 // The boxes intersect: we need to understand now how they intersect.
103 // We begin by computing the intersection:
104 std::array<double, spacedim> intersect_bbox_min;
105 std::array<double, spacedim> intersect_bbox_max;
106 for (unsigned int d = 0; d < spacedim; ++d)
107 {
108 intersect_bbox_min[d] = std::max(bbox1[0][d], bbox2[0][d]);
109 intersect_bbox_max[d] = std::min(bbox1[1][d], bbox2[1][d]);
110 }
111
112 // Finding the intersection's dimension
113 int intersect_dim = spacedim;
114 for (unsigned int d = 0; d < spacedim; ++d)
115 if (std::abs(intersect_bbox_min[d] - intersect_bbox_max[d]) <=
116 tolerance * (std::abs(intersect_bbox_min[d]) +
117 std::abs(intersect_bbox_max[d])))
118 --intersect_dim;
119
120 if (intersect_dim == 0 || intersect_dim == spacedim - 2)
122
123 // Checking the two mergeable cases: first if the boxes are aligned so
124 // that they can be merged
125 unsigned int not_align_1 = 0, not_align_2 = 0;
126 bool same_direction = true;
127 for (unsigned int d = 0; d < spacedim; ++d)
128 {
129 if (std::abs(bbox2[0][d] - bbox1[0][d]) >
130 tolerance * (std::abs(bbox2[0][d]) + std::abs(bbox1[0][d])))
131 ++not_align_1;
132 if (std::abs(bbox1[1][d] - bbox2[1][d]) >
133 tolerance * (std::abs(bbox1[1][d]) + std::abs(bbox2[1][d])))
134 ++not_align_2;
135 if (not_align_1 != not_align_2)
136 {
137 same_direction = false;
138 break;
139 }
140 }
141
142 if (not_align_1 <= 1 && not_align_2 <= 1 && same_direction)
144
145 // Second: one box is contained/equal to the other
146 if ((this->point_inside(bbox2[0]) && this->point_inside(bbox2[1])) ||
147 (other_bbox.point_inside(bbox1[0], tolerance) &&
148 other_bbox.point_inside(bbox1[1], tolerance)))
150
151 // Degenerate and mergeable cases have been found, it remains:
153 }
154}
155
156
157
158template <int spacedim, typename Number>
159double
161{
162 double vol = 1.0;
163 for (unsigned int i = 0; i < spacedim; ++i)
164 vol *= (this->boundary_points.second[i] - this->boundary_points.first[i]);
165 return vol;
166}
167
168
169
170template <int spacedim, typename Number>
171Number
172BoundingBox<spacedim, Number>::lower_bound(const unsigned int direction) const
173{
174 AssertIndexRange(direction, spacedim);
175
176 return boundary_points.first[direction];
177}
178
179
180
181template <int spacedim, typename Number>
182Number
183BoundingBox<spacedim, Number>::upper_bound(const unsigned int direction) const
184{
185 AssertIndexRange(direction, spacedim);
186
187 return boundary_points.second[direction];
188}
189
190
191
192template <int spacedim, typename Number>
195{
197 for (unsigned int i = 0; i < spacedim; ++i)
198 point[i] = .5 * (boundary_points.first[i] + boundary_points.second[i]);
199
200 return point;
201}
202
203
204
205template <int spacedim, typename Number>
207BoundingBox<spacedim, Number>::bounds(const unsigned int direction) const
209 AssertIndexRange(direction, spacedim);
210
211 std::pair<Point<1, Number>, Point<1, Number>> lower_upper_bounds;
212 lower_upper_bounds.first[0] = lower_bound(direction);
213 lower_upper_bounds.second[0] = upper_bound(direction);
214
215 return BoundingBox<1, Number>(lower_upper_bounds);
217
218
219
220template <int spacedim, typename Number>
221Number
222BoundingBox<spacedim, Number>::side_length(const unsigned int direction) const
223{
224 AssertIndexRange(direction, spacedim);
225
226 return boundary_points.second[direction] - boundary_points.first[direction];
227}
228
229
230
231template <int spacedim, typename Number>
233BoundingBox<spacedim, Number>::vertex(const unsigned int index) const
234{
236
237 const Point<spacedim> unit_cell_vertex =
239
241 for (unsigned int i = 0; i < spacedim; ++i)
242 point[i] = boundary_points.first[i] + side_length(i) * unit_cell_vertex[i];
243
244 return point;
245}
246
247
248
249template <int spacedim, typename Number>
251BoundingBox<spacedim, Number>::child(const unsigned int index) const
252{
254
255 // Vertex closest to child.
256 const Point<spacedim, Number> parent_vertex = vertex(index);
257 const Point<spacedim, Number> parent_center = center();
258
259 const Point<spacedim> upper_corner_unit_cell =
262
263 const Point<spacedim> lower_corner_unit_cell =
265
266 std::pair<Point<spacedim, Number>, Point<spacedim, Number>>
267 child_lower_upper_corner;
268 for (unsigned int i = 0; i < spacedim; ++i)
269 {
270 const double child_side_length = side_length(i) / 2;
271
272 const double child_center = (parent_center[i] + parent_vertex[i]) / 2;
273
274 child_lower_upper_corner.first[i] =
275 child_center + child_side_length * (lower_corner_unit_cell[i] - .5);
276 child_lower_upper_corner.second[i] =
277 child_center + child_side_length * (upper_corner_unit_cell[i] - .5);
278 }
280 return BoundingBox<spacedim, Number>(child_lower_upper_corner);
281}
282
283
284
285template <int spacedim, typename Number>
286BoundingBox<spacedim - 1, Number>
287BoundingBox<spacedim, Number>::cross_section(const unsigned int direction) const
288{
289 AssertIndexRange(direction, spacedim);
290
291 std::pair<Point<spacedim - 1, Number>, Point<spacedim - 1, Number>>
292 cross_section_lower_upper_corner;
293 for (unsigned int d = 0; d < spacedim - 1; ++d)
294 {
295 const int index_to_write_from =
296 internal::coordinate_to_one_dim_higher<spacedim - 1>(direction, d);
298 cross_section_lower_upper_corner.first[d] =
299 boundary_points.first[index_to_write_from];
300
301 cross_section_lower_upper_corner.second[d] =
302 boundary_points.second[index_to_write_from];
304
305 return BoundingBox<spacedim - 1, Number>(cross_section_lower_upper_corner);
306}
307
308
310template <int spacedim, typename Number>
313 const Point<spacedim, Number> &point) const
314{
315 auto unit = point;
316 const auto diag = boundary_points.second - boundary_points.first;
317 unit -= boundary_points.first;
318 for (unsigned int d = 0; d < spacedim; ++d)
319 unit[d] /= diag[d];
320 return unit;
321}
322
324
325template <int spacedim, typename Number>
328 const Point<spacedim, Number> &point) const
329{
330 auto real = boundary_points.first;
331 const auto diag = boundary_points.second - boundary_points.first;
332 for (unsigned int d = 0; d < spacedim; ++d)
333 real[d] += diag[d] * point[d];
334 return real;
335}
336
337
338
339template <int spacedim, typename Number>
340Number
342 const Point<spacedim, Number> &point,
343 const unsigned int direction) const
345 const Number p1 = lower_bound(direction);
346 const Number p2 = upper_bound(direction);
347
348 if (point[direction] > p2)
349 return point[direction] - p2;
350 else if (point[direction] < p1)
351 return p1 - point[direction];
352 else
353 return -std::min(point[direction] - p1, p2 - point[direction]);
354}
356
357
358template <int spacedim, typename Number>
359Number
361 const Point<spacedim, Number> &point) const
362{
363 // calculate vector of orthogonal signed distances
364 std::array<Number, spacedim> distances;
365 for (unsigned int d = 0; d < spacedim; ++d)
366 distances[d] = signed_distance(point, d);
367
368 // determine the number of positive signed distances
369 const unsigned int n_positive_signed_distances =
370 std::count_if(distances.begin(), distances.end(), [](const auto &a) {
371 return a > 0.0;
372 });
374 if (n_positive_signed_distances <= 1)
375 // point is inside of bounding box (0: all signed distances are
376 // negative; find the index with the smallest absolute value)
377 // or next to a face (1: all signed distances are negative
378 // but one; find this index)
379 return *std::max_element(distances.begin(), distances.end());
380 else
381 // point is next to a corner (2D/3D: all signed distances are
382 // positive) or a line (3D: all signed distances are positive
383 // but one) -> take the l2-norm of all positive signed distances
384 return std::sqrt(std::accumulate(distances.begin(),
385 distances.end(),
386 0.0,
387 [](const auto &a, const auto &b) {
388 return a + (b > 0 ? b * b : 0.0);
389 }));
390}
391
392
393
394template <int dim, typename Number>
397{
398 std::pair<Point<dim, Number>, Point<dim, Number>> lower_upper_corner;
399 for (unsigned int i = 0; i < dim; ++i)
400 {
401 lower_upper_corner.second[i] = 1;
402 }
403 return BoundingBox<dim, Number>(lower_upper_corner);
404}
405
406
407#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:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)