www.mooseframework.org
RayTracing.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 // Moose includes
16 #include "RayTracing.h"
17 #include "LineSegment.h"
18 #include "MooseError.h"
19 
20 #include "libmesh/plane.h"
21 #include "libmesh/point.h"
22 #include "libmesh/mesh.h"
23 #include "libmesh/elem.h"
24 
25 namespace Moose
26 {
27 
38 int
39 sideIntersectedByLine(const Elem * elem,
40  std::vector<int> & not_side,
41  const LineSegment & line_segment,
42  Point & intersection_point)
43 {
44  unsigned int n_sides = elem->n_sides();
45 
46  // Whether or not they intersect
47  bool intersect = false;
48 
49  unsigned int dim = elem->dim();
50 
51  for (unsigned int i = 0; i < n_sides; i++)
52  {
53  // Don't search the "not_side"
54  // Note: A linear search is fine here because this vector is going to be < n_sides
55  if (std::find(not_side.begin(), not_side.end(), static_cast<int>(i)) != not_side.end())
56  continue;
57 
58  // Get a simplified side element
59  std::unique_ptr<Elem> side_elem = elem->side(i);
60 
61  if (dim == 3)
62  {
63  // Make a plane out of the first three nodes on the side
64  Plane plane(side_elem->point(0), side_elem->point(1), side_elem->point(2));
65 
66  // See if they intersect
67  intersect = line_segment.intersect(plane, intersection_point);
68  }
69  else if (dim == 2)
70  {
71  // Make a Line Segment out of the first two nodes on the side
72  LineSegment side_segment(side_elem->point(0), side_elem->point(1));
73 
74  // See if they intersect
75  intersect = line_segment.intersect(side_segment, intersection_point);
76  }
77  else // 1D
78  {
79  // See if the line segment contains the point
80  intersect = line_segment.contains_point(side_elem->point(0));
81 
82  // If it does then save off that one point as the intersection point
83  if (intersect)
84  intersection_point = side_elem->point(0);
85  }
86 
87  if (intersect)
88  {
89  if (side_elem->contains_point(intersection_point))
90  {
91  const Elem * neighbor = elem->neighbor_ptr(i);
92 
93  // If this side is on a boundary, let's do another search and see if we can find a better
94  // candidate
95  if (!neighbor)
96  {
97  not_side.push_back(i); // Make sure we don't find this side again
98 
99  int better_side = sideIntersectedByLine(elem, not_side, line_segment, intersection_point);
100 
101  if (better_side != -1)
102  return better_side;
103  }
104 
105  return i;
106  }
107  }
108  }
109 
110  // Didn't find one
111  return -1;
112 }
113 
119 int
120 sideNeighborIsOn(const Elem * elem, const Elem * neighbor)
121 {
122  unsigned int n_sides = elem->n_sides();
123 
124  for (unsigned int i = 0; i < n_sides; i++)
125  {
126  if (elem->neighbor_ptr(i) == neighbor)
127  return i;
128  }
129 
130  return -1;
131 }
132 
148 void
150  const Elem * current_elem,
151  int incoming_side,
152  const Point & incoming_point,
153  std::vector<Elem *> & intersected_elems,
154  std::vector<LineSegment> & segments)
155 {
156  Point intersection_point;
157 
158  std::vector<int> not_side(1, incoming_side);
159 
160  // Find the side of this element that the LineSegment intersects... while ignoring the incoming
161  // side (we don't want to move backward!)
162  int intersected_side =
163  sideIntersectedByLine(current_elem, not_side, line_segment, intersection_point);
164 
165  if (intersected_side != -1) // -1 means that we didn't find any side
166  {
167  // Get the neighbor on that side
168  const Elem * neighbor = current_elem->neighbor_ptr(intersected_side);
169 
170  if (neighbor)
171  {
172  // Add it to the list
173  intersected_elems.push_back(const_cast<Elem *>(neighbor));
174 
175  // Add the line segment across the element to the segments list
176  segments.push_back(LineSegment(incoming_point, intersection_point));
177 
178  // Note: This is finding the side the current_elem is on for the neighbor. That's the
179  // "incoming_side" for the neighbor
180  int incoming_side = sideNeighborIsOn(neighbor, current_elem);
181 
182  // Recurse
184  line_segment, neighbor, incoming_side, intersection_point, intersected_elems, segments);
185  }
186  else // Add the final segment
187  segments.push_back(LineSegment(incoming_point, line_segment.end()));
188  }
189  else // Add the final segment
190  segments.push_back(LineSegment(incoming_point, line_segment.end()));
191 
192  // Finished... return out!
193  return;
194 }
195 
196 void
197 elementsIntersectedByLine(const Point & p0,
198  const Point & p1,
199  const MeshBase & /*mesh*/,
200  const PointLocatorBase & point_locator,
201  std::vector<Elem *> & intersected_elems,
202  std::vector<LineSegment> & segments)
203 {
204  // Make sure our list is clear
205  intersected_elems.clear();
206 
207  // Find the starting element
208  const Elem * first_elem = point_locator(p0);
209 
210  // Quick return if can't even locate the first element.
211  if (!first_elem)
212  return;
213 
214  intersected_elems.push_back(const_cast<Elem *>(first_elem));
215 
216  // Make a LineSegment object out of our two points for ease:
217  LineSegment line_segment = LineSegment(p0, p1);
218 
219  // Find 'em!
221  line_segment, first_elem, -1, p0, intersected_elems, segments);
222 }
223 }
int sideNeighborIsOn(const Elem *elem, const Elem *neighbor)
Returns the side number for elem that neighbor is on.
Definition: RayTracing.C:120
The LineSegment class is used by the LineMaterialSamplerBase class and for some ray tracing stuff...
Definition: LineSegment.h:33
bool contains_point(const Point &p) const
Determines whether a point is in a line segment or not.
Definition: LineSegment.C:62
bool intersect(const Plane &pl, Point &intersect_p) const
Definition: LineSegment.C:69
void elementsIntersectedByLine(const Point &p0, const Point &p1, const MeshBase &mesh, const PointLocatorBase &point_locator, std::vector< Elem * > &intersected_elems, std::vector< LineSegment > &segments)
Find all of the elements intersected by a line.
Definition: RayTracing.C:197
void recursivelyFindElementsIntersectedByLine(const LineSegment &line_segment, const Elem *current_elem, int incoming_side, const Point &incoming_point, std::vector< Elem * > &intersected_elems, std::vector< LineSegment > &segments)
Recursively find all elements intersected by a line segment.
Definition: RayTracing.C:149
Definition: Moose.h:84
const Point & end() const
Ending of the line segment.
Definition: LineSegment.h:71
int sideIntersectedByLine(const Elem *elem, std::vector< int > &not_side, const LineSegment &line_segment, Point &intersection_point)
Figure out which (if any) side of an Elem is intersected by a line.
Definition: RayTracing.C:39