Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "BoundingBoxIntersectionHelper.h"
11 :
12 : // libMesh includes
13 : #include "libmesh/bounding_box.h"
14 : #include "libmesh/face_quad4.h"
15 : #include "libmesh/cell_hex8.h"
16 :
17 : #include "TraceRayTools.h"
18 : #include "DebugRay.h"
19 :
20 : using namespace libMesh;
21 :
22 879 : BoundingBoxIntersectionHelper::BoundingBoxIntersectionHelper(const BoundingBox & bbox,
23 879 : const unsigned int dim)
24 879 : : _comm(), _mesh(std::make_unique<Mesh>(_comm, dim))
25 : {
26 : // Add the nodes that represent our bounding box
27 879 : _mesh->add_point(Point(bbox.min()(0), bbox.min()(1), bbox.min()(2)), 0, 0);
28 879 : _mesh->add_point(Point(bbox.max()(0), bbox.min()(1), bbox.min()(2)), 1, 0);
29 879 : if (dim > 1)
30 : {
31 855 : _mesh->add_point(Point(bbox.max()(0), bbox.max()(1), bbox.min()(2)), 2, 0);
32 855 : _mesh->add_point(Point(bbox.min()(0), bbox.max()(1), bbox.min()(2)), 3, 0);
33 :
34 855 : if (dim == 3)
35 : {
36 592 : _mesh->add_point(Point(bbox.min()(0), bbox.min()(1), bbox.max()(2)), 4, 0);
37 592 : _mesh->add_point(Point(bbox.max()(0), bbox.min()(1), bbox.max()(2)), 5, 0);
38 592 : _mesh->add_point(Point(bbox.max()(0), bbox.max()(1), bbox.max()(2)), 6, 0);
39 592 : _mesh->add_point(Point(bbox.min()(0), bbox.max()(1), bbox.max()(2)), 7, 0);
40 : }
41 : }
42 :
43 : // Build the element and set the nodes accordingly
44 : Elem * elem;
45 879 : if (dim == 1)
46 24 : elem = Elem::build(EDGE2).release();
47 855 : else if (dim == 2)
48 263 : elem = Elem::build(QUAD4).release();
49 : else
50 592 : elem = Elem::build(HEX8).release();
51 879 : elem = _mesh->add_elem(elem);
52 879 : elem->subdomain_id() = 0;
53 : elem->processor_id(0);
54 879 : elem->set_id() = 0;
55 6715 : for (unsigned int n = 0; n < _mesh->n_nodes(); ++n)
56 5836 : elem->set_node(n, _mesh->node_ptr(n));
57 :
58 : // All done with the mesh
59 : _mesh->skip_partitioning(true);
60 879 : _mesh->prepare_for_use();
61 :
62 : // This costs a little bit so why not cache it now
63 879 : _hmax = elem->hmax();
64 879 : }
65 :
66 : Point
67 834180 : BoundingBoxIntersectionHelper::intersection(const Point & point, const Point & direction) const
68 : {
69 : mooseAssert(std::abs(1.0 - direction.norm()) < TOLERANCE, "Unnormalized direction");
70 :
71 834180 : const Elem * elem = _mesh->elem_ptr(0);
72 :
73 : // We're going to check all possible intersections and keep the one that's furthest away
74 834180 : Point best_intersection_point = RayTracingCommon::invalid_point;
75 : Real best_intersection_distance = -std::numeric_limits<Real>::max();
76 :
77 : // For use in the intersection routines
78 : Point intersection_point;
79 : ElemExtrema intersected_extrema;
80 : Real intersection_distance;
81 : bool intersected;
82 :
83 : // Check intersection with sides
84 5814060 : for (const auto s : elem->side_index_range())
85 : {
86 4979880 : intersection_point = RayTracingCommon::invalid_point;
87 : intersected_extrema.invalidate();
88 :
89 4979880 : if (elem->dim() == 3) // 3D: Intersect HEX8
90 4929912 : intersected = TraceRayTools::sideIntersectedByLine<Hex8>(elem,
91 : point,
92 : direction,
93 : s,
94 : _hmax,
95 : intersection_point,
96 : intersection_distance,
97 : intersected_extrema,
98 4929912 : _hmax
99 : #ifdef DEBUG_RAY_INTERSECTIONS
100 : ,
101 : false
102 : #endif
103 : );
104 49968 : else if (elem->dim() == 2) // 2D: Intersect QUAD4
105 : {
106 49824 : intersected = TraceRayTools::sideIntersectedByLine<Quad4>(elem,
107 : point,
108 : direction,
109 : s,
110 : _hmax,
111 : intersection_point,
112 : intersection_distance,
113 : intersected_extrema,
114 49824 : _hmax
115 : #ifdef DEBUG_RAY_INTERSECTIONS
116 : ,
117 : false
118 : #endif
119 : );
120 : }
121 : else // 1D: see if side point is between point and end far away
122 : {
123 : // Dummy end point outside of the element
124 : const Point dummy_end = point + 1.01 * direction * _hmax;
125 :
126 144 : intersected = TraceRayTools::isWithinSegment(point, dummy_end, elem->point(s));
127 144 : if (intersected)
128 : {
129 144 : intersection_distance = (elem->point(s) - point).norm();
130 144 : intersection_point = elem->point(s);
131 : }
132 : }
133 :
134 4979880 : if (intersected && intersection_distance > best_intersection_distance)
135 : {
136 : best_intersection_distance = intersection_distance;
137 840619 : best_intersection_point = intersection_point;
138 : }
139 : }
140 :
141 834180 : return best_intersection_point;
142 : }
|
