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 "ClaimRays.h"
11 :
12 : // Local includes
13 : #include "RayTracingStudy.h"
14 :
15 : // libMesh includes
16 : #include "libmesh/elem.h"
17 : #include "libmesh/parallel_algebra.h"
18 : #include "libmesh/parallel_sync.h"
19 : #include "libmesh/enum_point_locator_type.h"
20 : #include "libmesh/mesh_tools.h"
21 :
22 : using namespace libMesh;
23 :
24 3424 : ClaimRays::ClaimRays(RayTracingStudy & study,
25 : const std::vector<std::shared_ptr<Ray>> & rays,
26 : std::vector<std::shared_ptr<Ray>> & local_rays,
27 3424 : const bool do_exchange)
28 : : ParallelObject(study.comm()),
29 : MeshChangedInterface(study.parameters()),
30 3424 : _mesh(study.mesh()),
31 3424 : _pid(comm().rank()),
32 3424 : _do_exchange(do_exchange),
33 3424 : _study(study),
34 3424 : _parallel_study(*_study.parallelStudy()),
35 3424 : _rays(rays),
36 3424 : _local_rays(local_rays),
37 3424 : _needs_init(true)
38 : {
39 3424 : }
40 :
41 : void
42 841 : ClaimRays::claim()
43 : {
44 841 : if (_needs_init)
45 : {
46 841 : init();
47 841 : _needs_init = false;
48 : }
49 :
50 841 : preClaim();
51 :
52 : // Clear these as we're about to fill
53 841 : _local_rays.clear();
54 :
55 : // Grab the point locator
56 1682 : _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, _mesh.getMesh());
57 841 : _point_locator->enable_out_of_mesh_mode();
58 :
59 : // Exchange: filter Rays into processors that _may_ claim them
60 : std::unordered_map<processor_id_type, std::vector<std::shared_ptr<Ray>>> rays_to_send;
61 841 : if (_do_exchange)
62 576 : for (processor_id_type pid = 0; pid < comm().size(); ++pid)
63 368 : if (_pid != pid)
64 : {
65 : const BoundingBox & pid_bbox = inflatedBoundingBox(pid);
66 26080 : for (auto & ray : _rays)
67 25920 : if (pid_bbox.contains_point(ray->currentPoint()))
68 14468 : rays_to_send[pid].push_back(ray);
69 : }
70 :
71 : // Functor for possibly claiming a vector of Rays
72 : auto claim_functor =
73 991 : [&](processor_id_type /* pid */, const std::vector<std::shared_ptr<Ray>> & rays)
74 : {
75 60254 : for (auto & ray : rays)
76 59263 : possiblyClaim(ray);
77 1832 : };
78 :
79 : // Send the relevant Rays to everyone and then attempt to claim the ones that we receive
80 841 : if (_do_exchange)
81 208 : Parallel::push_parallel_packed_range(comm(), rays_to_send, &_parallel_study, claim_functor);
82 :
83 : // Attempt to claim the locally generated rays in _rays
84 841 : claim_functor(_pid, _rays);
85 :
86 : // Verify the claiming if the study so desires
87 841 : if (_study.verifyRays())
88 841 : verifyClaiming();
89 :
90 841 : postClaim();
91 841 : }
92 :
93 : void
94 59263 : ClaimRays::possiblyClaim(const std::shared_ptr<Ray> & ray)
95 : {
96 59263 : prePossiblyClaimRay(ray);
97 :
98 : const auto elem =
99 59263 : claimPoint(ray->currentPoint(), getID(ray), (*_point_locator)(ray->currentPoint()));
100 59263 : if (elem)
101 : {
102 43441 : _local_rays.push_back(ray);
103 43441 : postClaimRay(_local_rays.back(), elem);
104 : }
105 59263 : }
106 :
107 : const Elem *
108 59263 : ClaimRays::claimPoint(const Point & point, const RayID id, const Elem * elem)
109 : {
110 59263 : if (elem)
111 : {
112 : // Looking for smallest (even ID Ray) or largest (odd ID Ray) elem id
113 45257 : const bool smallest = id % 2 == 0;
114 :
115 : // Start with the element we found, as it is a valid candidate
116 : const Elem * extremum_elem = elem;
117 :
118 : // All point neighbors for this element
119 : mooseAssert(_elem_point_neighbors.count(elem->id()), "Not in point neighbor map");
120 45257 : const auto & neighbors = _elem_point_neighbors.at(elem->id());
121 :
122 : // Find element that matches the extremum criteria
123 367961 : for (const auto & neighbor : neighbors)
124 : {
125 : mooseAssert(neighbor->active(), "Inactive neighbor");
126 :
127 322704 : if ((smallest && neighbor->id() < extremum_elem->id()) || // satisfies
128 123934 : (!smallest && neighbor->id() > extremum_elem->id())) // ...one of the id checks
129 130675 : if (neighbor->contains_point(point)) // and also contains the point
130 18095 : extremum_elem = neighbor;
131 : }
132 :
133 : // Claim the object if we own the extremum elem
134 45257 : if (extremum_elem->processor_id() == _pid)
135 : {
136 : mooseAssert(extremum_elem->active(), "Inactive element");
137 : return extremum_elem;
138 : }
139 : }
140 :
141 : return nullptr;
142 : }
143 :
144 : void
145 43441 : ClaimRays::postClaimRay(std::shared_ptr<Ray> & ray, const Elem * elem)
146 : {
147 : mooseAssert(_mesh.queryElemPtr(elem->id()) == elem, "Mesh doesn't contain elem");
148 : mooseAssert(elem->active(), "Inactive element");
149 :
150 : // If the incoming side is set and is not incoming, or if it is not set at all, see
151 : // if we can find an incoming side that is valid.
152 : auto starting_incoming_side = RayTracingCommon::invalid_side;
153 43441 : if (!(!ray->invalidCurrentIncomingSide() &&
154 0 : _study.elemSide(*elem, ray->currentIncomingSide()).contains_point(ray->currentPoint()) &&
155 0 : _study.sideIsIncoming(elem, ray->currentIncomingSide(), ray->direction(), /* tid = */ 0)))
156 98729 : for (const auto s : elem->side_index_range())
157 144988 : if (_study.elemSide(*elem, s).contains_point(ray->currentPoint()) &&
158 46614 : _study.sideIsIncoming(elem, s, ray->direction(), /* tid = */ 0))
159 : {
160 : starting_incoming_side = s;
161 : break;
162 : }
163 :
164 43441 : ray->setStart(ray->currentPoint(), elem, starting_incoming_side);
165 43441 : }
166 :
167 : void
168 841 : ClaimRays::init()
169 : {
170 841 : buildBoundingBoxes();
171 841 : buildPointNeighbors();
172 841 : }
173 :
174 : void
175 208 : ClaimRays::meshChanged()
176 : {
177 208 : _needs_init = true;
178 208 : }
179 :
180 : void
181 841 : ClaimRays::buildBoundingBoxes()
182 : {
183 : // Local bounding box
184 841 : const auto bbox = MeshTools::create_local_bounding_box(_mesh.getMesh());
185 :
186 : // Gather the bounding boxes of all processors
187 841 : std::vector<std::pair<Point, Point>> bb_points = {static_cast<std::pair<Point, Point>>(bbox)};
188 841 : comm().allgather(bb_points, true);
189 :
190 : // Inflate the local bboxes by a bit and store
191 841 : _inflated_bboxes.resize(comm().size());
192 2362 : for (processor_id_type pid = 0; pid < comm().size(); ++pid)
193 : {
194 1521 : BoundingBox pid_bbox = static_cast<BoundingBox>(bb_points[pid]);
195 1521 : pid_bbox.scale(0.01);
196 : _inflated_bboxes[pid] = pid_bbox;
197 : }
198 841 : }
199 :
200 : void
201 841 : ClaimRays::buildPointNeighbors()
202 : {
203 : _elem_point_neighbors.clear();
204 841 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
205 :
206 132194 : for (const auto & elem : _mesh.getMesh().active_element_ptr_range())
207 : {
208 65256 : auto & fill = _elem_point_neighbors[elem->id()];
209 368530 : for (unsigned int v = 0; v < elem->n_vertices(); ++v)
210 : {
211 303274 : const auto & node = elem->node_ptr(v);
212 4429300 : for (const auto & neighbor_id : node_to_elem_map.at(node->id()))
213 : {
214 4126026 : if (neighbor_id == elem->id())
215 303274 : continue;
216 :
217 3822752 : const auto & neighbor = _mesh.elemPtr(neighbor_id);
218 3822752 : if (std::count(fill.begin(), fill.end(), neighbor) == 0)
219 2857192 : fill.emplace_back(neighbor);
220 : }
221 : }
222 841 : }
223 841 : }
224 :
225 : void
226 841 : ClaimRays::verifyClaiming()
227 : {
228 : // NOTE for all of the following: we use char here in place of bool.
229 : // This is because bool is not instantiated as a StandardType in
230 : // TIMPI due to the fun of std::vector<bool>
231 :
232 : // Map from Ray ID -> whether or not it was generated (false) or
233 : // claimed/possibly also generated (true)
234 : std::map<RayID, char> local_map;
235 : auto add_to_local_map =
236 1682 : [this, &local_map](const std::vector<std::shared_ptr<Ray>> & rays, const bool claimed_rays)
237 : {
238 89918 : for (const auto & ray : rays)
239 : {
240 88236 : const auto id = getID(ray);
241 :
242 : // Try to insert into the map
243 88236 : auto emplace_pair = local_map.emplace(id, claimed_rays);
244 :
245 : // If it already exists but has not been claimed yet, set it to being claimed
246 88236 : if (!emplace_pair.second && claimed_rays)
247 : {
248 : mooseAssert(!emplace_pair.first->second,
249 : "Ray was claimed more than once on a single processor");
250 33524 : emplace_pair.first->second = true;
251 : }
252 : }
253 2523 : };
254 :
255 : // Build the local_map
256 841 : add_to_local_map(_rays, false);
257 841 : add_to_local_map(_local_rays, true);
258 :
259 : // Build the structure to send the local generation/claiming information to rank 0
260 : std::map<processor_id_type, std::vector<std::pair<RayID, char>>> send_info;
261 841 : if (local_map.size())
262 : send_info.emplace(std::piecewise_construct,
263 841 : std::forward_as_tuple(0),
264 841 : std::forward_as_tuple(local_map.begin(), local_map.end()));
265 :
266 : // The mapping (filled on rank 0) from Ray ID -> (processor id, claiming status)
267 : std::map<RayID, std::vector<std::pair<processor_id_type, char>>> global_map;
268 :
269 : // Functor for receiving the generation/claiming information
270 841 : auto receive_functor = [&global_map](processor_id_type pid,
271 54712 : const std::vector<std::pair<RayID, char>> & id_claimed_pairs)
272 : {
273 55553 : for (const auto & id_claimed_pair : id_claimed_pairs)
274 54712 : global_map[id_claimed_pair.first].emplace_back(pid, id_claimed_pair.second);
275 841 : };
276 :
277 : // Send claiming information to rank 0
278 841 : Parallel::push_parallel_vector_data(comm(), send_info, receive_functor);
279 :
280 : // Rank 0 will make sure everything looks good
281 841 : if (_pid == 0)
282 43962 : for (const auto & id_pairs_pair : global_map)
283 : {
284 43441 : const RayID id = id_pairs_pair.first;
285 : const std::vector<std::pair<processor_id_type, char>> & pid_claimed_pairs =
286 : id_pairs_pair.second;
287 :
288 : std::vector<processor_id_type> claimed_pids;
289 98153 : for (const auto & pid_claimed_pair : pid_claimed_pairs)
290 54712 : if (pid_claimed_pair.second)
291 43441 : claimed_pids.push_back(pid_claimed_pair.first);
292 :
293 43441 : if (claimed_pids.size() == 0)
294 0 : _study.mooseError("Failed to claim the Ray with ID ", id);
295 : mooseAssert(claimed_pids.size() == 1, "Ray was claimed on multiple processors");
296 : }
297 841 : }
|
