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 : // Moose includes
11 : #include "RayTracing.h"
12 : #include "LineSegment.h"
13 : #include "MooseError.h"
14 :
15 : #include "libmesh/plane.h"
16 : #include "libmesh/point.h"
17 : #include "libmesh/mesh.h"
18 : #include "libmesh/elem.h"
19 :
20 : using namespace libMesh;
21 :
22 : namespace Moose
23 : {
24 :
25 : /**
26 : * Figure out which (if any) side of an Elem is intersected by a line.
27 : *
28 : * @param elem The elem to search
29 : * @param not_side Sides to _not_ search (Use -1 if you want to search all sides)
30 : * @param intersection_point If an intersection is found this will be filled with the x,y,z position
31 : * of that intersection
32 : * @return The side that is intersected by the line. Will return -1 if it doesn't intersect any
33 : * side
34 : */
35 : int
36 498 : sideIntersectedByLine(const Elem * elem,
37 : std::vector<int> & not_side,
38 : const LineSegment & line_segment,
39 : Point & intersection_point)
40 : {
41 498 : unsigned int n_sides = elem->n_sides();
42 :
43 : // Whether or not they intersect
44 498 : bool intersect = false;
45 :
46 498 : unsigned int dim = elem->dim();
47 :
48 1661 : for (unsigned int i = 0; i < n_sides; i++)
49 : {
50 : // Don't search the "not_side"
51 : // Note: A linear search is fine here because this vector is going to be < n_sides
52 1575 : if (std::find(not_side.begin(), not_side.end(), static_cast<int>(i)) != not_side.end())
53 405 : continue;
54 :
55 : // Get a simplified side element
56 1170 : std::unique_ptr<const Elem> side_elem = elem->side_ptr(i);
57 :
58 1170 : if (dim == 3)
59 : {
60 : // Make a plane out of the first three nodes on the side
61 572 : Plane plane(side_elem->point(0), side_elem->point(1), side_elem->point(2));
62 :
63 : // See if they intersect
64 572 : intersect = line_segment.intersect(plane, intersection_point);
65 572 : }
66 598 : else if (dim == 2)
67 : {
68 : // Make a Line Segment out of the first two nodes on the side
69 345 : LineSegment side_segment(side_elem->point(0), side_elem->point(1));
70 :
71 : // See if they intersect
72 345 : intersect = line_segment.intersect(side_segment, intersection_point);
73 345 : }
74 : else // 1D
75 : {
76 : // See if the line segment contains the point
77 253 : intersect = line_segment.contains_point(side_elem->point(0));
78 :
79 : // If it does then save off that one point as the intersection point
80 253 : if (intersect)
81 209 : intersection_point = side_elem->point(0);
82 : }
83 :
84 1170 : if (intersect)
85 : {
86 412 : if (side_elem->contains_point(intersection_point))
87 : {
88 412 : const Elem * neighbor = elem->neighbor_ptr(i);
89 :
90 : // If this side is on a boundary, let's do another search and see if we can find a better
91 : // candidate
92 412 : if (!neighbor)
93 : {
94 22 : not_side.push_back(i); // Make sure we don't find this side again
95 :
96 22 : int better_side = sideIntersectedByLine(elem, not_side, line_segment, intersection_point);
97 :
98 22 : if (better_side != -1)
99 11 : return better_side;
100 : }
101 :
102 401 : return i;
103 : }
104 : }
105 1170 : }
106 :
107 : // Didn't find one
108 86 : return -1;
109 : }
110 :
111 : /**
112 : * Returns the side number for elem that neighbor is on
113 : *
114 : * Returns -1 if the neighbor can't be found to be a neighbor
115 : */
116 : int
117 390 : sideNeighborIsOn(const Elem * elem, const Elem * neighbor)
118 : {
119 390 : unsigned int n_sides = elem->n_sides();
120 :
121 471 : for (unsigned int i = 0; i < n_sides; i++)
122 : {
123 471 : if (elem->neighbor_ptr(i) == neighbor)
124 390 : return i;
125 : }
126 :
127 0 : return -1;
128 : }
129 :
130 : /**
131 : * Recursively find all elements intersected by a line segment
132 : *
133 : * Works by moving from one element to the next _through_ the side of the current element.
134 : * This means that (other than for the first element) there is always an incoming_side
135 : * that is the reason we ended up in this element in the first place. Search all the _other_
136 : * sides to see if there is a next element...
137 : *
138 : * @param line_segment the LineSegment to intersect
139 : * @param current_elem The current element that needs to be searched
140 : * @param incoming_side The side of the current element that was intersected by the LineSegment that
141 : * brought us here
142 : * @param intersected_elems The output
143 : * @param segments Line segments for the path across each element
144 : */
145 : void
146 476 : recursivelyFindElementsIntersectedByLine(const LineSegment & line_segment,
147 : const Elem * current_elem,
148 : int incoming_side,
149 : const Point & incoming_point,
150 : std::vector<Elem *> & intersected_elems,
151 : std::vector<LineSegment> & segments)
152 : {
153 476 : Point intersection_point;
154 :
155 476 : std::vector<int> not_side(1, incoming_side);
156 :
157 : // Find the side of this element that the LineSegment intersects... while ignoring the incoming
158 : // side (we don't want to move backward!)
159 : int intersected_side =
160 476 : sideIntersectedByLine(current_elem, not_side, line_segment, intersection_point);
161 :
162 476 : if (intersected_side != -1) // -1 means that we didn't find any side
163 : {
164 : // Get the neighbor on that side
165 401 : const Elem * neighbor = current_elem->neighbor_ptr(intersected_side);
166 :
167 401 : if (neighbor)
168 : {
169 : // Add it to the list
170 390 : intersected_elems.push_back(const_cast<Elem *>(neighbor));
171 :
172 : // Add the line segment across the element to the segments list
173 390 : segments.push_back(LineSegment(incoming_point, intersection_point));
174 :
175 : // Note: This is finding the side the current_elem is on for the neighbor. That's the
176 : // "incoming_side" for the neighbor
177 390 : int incoming_side = sideNeighborIsOn(neighbor, current_elem);
178 :
179 : // Recurse
180 390 : recursivelyFindElementsIntersectedByLine(
181 : line_segment, neighbor, incoming_side, intersection_point, intersected_elems, segments);
182 : }
183 : else // Add the final segment
184 11 : segments.push_back(LineSegment(incoming_point, line_segment.end()));
185 : }
186 : else // Add the final segment
187 75 : segments.push_back(LineSegment(incoming_point, line_segment.end()));
188 :
189 : // Finished... return out!
190 952 : return;
191 476 : }
192 :
193 : void
194 86 : elementsIntersectedByLine(const Point & p0,
195 : const Point & p1,
196 : const MeshBase & /*mesh*/,
197 : const PointLocatorBase & point_locator,
198 : std::vector<Elem *> & intersected_elems,
199 : std::vector<LineSegment> & segments)
200 : {
201 : // Make sure our list is clear
202 86 : intersected_elems.clear();
203 :
204 : // Find the starting element
205 86 : const Elem * first_elem = point_locator(p0);
206 :
207 : // Quick return if can't even locate the first element.
208 86 : if (!first_elem)
209 0 : return;
210 :
211 86 : intersected_elems.push_back(const_cast<Elem *>(first_elem));
212 :
213 : // Make a LineSegment object out of our two points for ease:
214 86 : LineSegment line_segment = LineSegment(p0, p1);
215 :
216 : // Find 'em!
217 86 : recursivelyFindElementsIntersectedByLine(
218 : line_segment, first_elem, -1, p0, intersected_elems, segments);
219 86 : }
220 : }
|
