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 : // TODO: remove ignore warnings once std::tuple<> is instantiated as a StandardType
11 : // Using push_parallel_vector_data with std::tuple<> leads to a -Wextra error
12 : // https://github.com/libMesh/TIMPI/issues/52
13 : #include "libmesh/ignore_warnings.h"
14 : #include "RayTracingMeshOutput.h"
15 : #include "libmesh/restore_warnings.h"
16 :
17 : // Local Includes
18 : #include "RayTracingStudy.h"
19 : #include "TraceData.h"
20 :
21 : // libMesh includes
22 : #include "libmesh/mesh_tools.h"
23 : #include "libmesh/parallel_algebra.h"
24 : #include "libmesh/parallel_sync.h"
25 : #include "libmesh/remote_elem.h"
26 :
27 : using namespace libMesh;
28 :
29 : InputParameters
30 874 : RayTracingMeshOutput::validParams()
31 : {
32 874 : InputParameters params = FileOutput::validParams();
33 :
34 1748 : params.addRequiredParam<UserObjectName>("study", "The RayTracingStudy to get the segments from");
35 :
36 1748 : params.addParam<bool>("output_data", false, "Whether or not to also output the Ray's data");
37 1748 : params.addParam<std::vector<std::string>>("output_data_names",
38 : "The names of specific data to output");
39 1748 : params.addParam<bool>("output_data_nodal",
40 1748 : false,
41 : "Whether or not to output the Ray's data in a nodal sense, in which the "
42 : "data is interpolated linearly across segments");
43 :
44 1748 : params.addParam<bool>(
45 1748 : "output_aux_data", false, "Whether or not to also output the Ray's aux data");
46 1748 : params.addParam<std::vector<std::string>>("output_aux_data_names",
47 : "The names of specific aux data to output");
48 :
49 : MultiMooseEnum props("ray_id intersections pid processor_crossings trajectory_changes",
50 874 : "ray_id intersections");
51 1748 : params.addParam<MultiMooseEnum>("output_properties", props, "Which Ray properties to output");
52 :
53 874 : return params;
54 874 : }
55 :
56 443 : RayTracingMeshOutput::RayTracingMeshOutput(const InputParameters & params)
57 : : FileOutput(params),
58 : UserObjectInterface(this),
59 443 : _study(getUserObject<RayTracingStudy>("study")),
60 886 : _output_data(getParam<bool>("output_data")),
61 443 : _output_data_names(isParamValid("output_data_names")
62 527 : ? &getParam<std::vector<std::string>>("output_data_names")
63 : : nullptr),
64 886 : _output_data_nodal(getParam<bool>("output_data_nodal")),
65 886 : _output_aux_data(getParam<bool>("output_aux_data")),
66 886 : _output_aux_data_names(isParamValid("output_aux_data_names")
67 489 : ? &getParam<std::vector<std::string>>("output_aux_data_names")
68 : : nullptr),
69 443 : _ray_id_var(invalid_uint),
70 443 : _intersections_var(invalid_uint),
71 443 : _pid_var(invalid_uint),
72 443 : _processor_crossings_var(invalid_uint),
73 443 : _trajectory_changes_var(invalid_uint),
74 886 : _segmented_rays(false)
75 : {
76 443 : if ((_output_data || _output_data_names) && !_study.dataOnCacheTraces())
77 2 : mooseError("In order to output Ray data in output '",
78 : name(),
79 : "', the RayTracingStudy '",
80 : _study.name(),
81 : "' must set data_on_cache_traces = true");
82 441 : if ((_output_aux_data || _output_aux_data_names) && !_study.auxDataOnCacheTraces())
83 2 : mooseError("In order to output Ray aux data in output '",
84 : name(),
85 : "', the RayTracingStudy '",
86 : _study.name(),
87 : "' must set aux_data_on_cache_traces = true");
88 439 : if (_output_data_nodal && !_output_data && !_output_data_names)
89 2 : paramError("output_data_nodal",
90 : "Cannot be used unless there is data to output; in addition set either "
91 : "'output_data' or 'output_data_names");
92 437 : if (_output_data_nodal && !_study.segmentsOnCacheTraces())
93 2 : paramError("output_data_nodal", "Not supported when study segments_on_cache_traces = false");
94 435 : if (_output_data && _output_data_names)
95 2 : paramError("output_data", "Cannot be used in addition to 'output_data_names'; choose one");
96 433 : if (_output_aux_data && _output_aux_data_names)
97 2 : paramError("output_aux_data",
98 : "Cannot be used in addition to 'output_aux_data_names'; choose one");
99 431 : }
100 :
101 : std::string
102 1025 : RayTracingMeshOutput::filename()
103 : {
104 : // Append the .e extension on the base file name
105 1025 : std::ostringstream output;
106 1025 : output << _file_base << ".e";
107 :
108 : // Add the _000x extension to the file
109 1025 : if (_file_num > 0)
110 283 : output << "-s" << std::setw(_padding) << std::setprecision(0) << std::setfill('0') << std::right
111 283 : << _file_num;
112 :
113 : // Return the filename
114 1025 : return output.str();
115 1025 : }
116 :
117 : void
118 600 : RayTracingMeshOutput::output()
119 : {
120 : // Do we even have any traces?
121 600 : auto num_segments = _study.getCachedTraces().size();
122 600 : _communicator.sum(num_segments);
123 600 : if (!num_segments)
124 2 : mooseError("No cached trace segments were found in the study '", _study.name(), "'.");
125 :
126 : // Build the _inflated_neighbor_bboxes
127 598 : buildBoundingBoxes();
128 :
129 : // Build the _segment_mesh
130 598 : buildSegmentMesh();
131 :
132 : // Setup the system to store the Ray field data
133 598 : setupEquationSystem();
134 :
135 : // Fill the field data
136 594 : fillFields();
137 :
138 : // And output
139 594 : outputMesh();
140 :
141 : // Done with these
142 : // We don't necessarily need to create a new mesh every time, but it's easier than
143 : // checking if the Rays have changed from last time we built a mesh
144 : _es = nullptr;
145 594 : _sys = nullptr;
146 : _segment_mesh = nullptr;
147 594 : }
148 :
149 : void
150 598 : RayTracingMeshOutput::buildIDMap()
151 : {
152 1196 : TIME_SECTION("buildIDMap", 3, "Building RayTracing ID Map");
153 :
154 : // Build the maximum number of nodes required to represent each one of my local Rays
155 : std::map<RayID, dof_id_type> local_ray_needed_nodes;
156 4569 : for (const auto & trace_data : _study.getCachedTraces())
157 : {
158 3971 : const auto find = local_ray_needed_nodes.find(trace_data._ray_id);
159 3971 : if (find == local_ray_needed_nodes.end())
160 3861 : local_ray_needed_nodes[trace_data._ray_id] = neededNodes(trace_data);
161 : else
162 : {
163 110 : auto & value = find->second;
164 220 : value = std::max(value, neededNodes(trace_data));
165 : }
166 : }
167 :
168 : // Fill all of my local Ray maxima to be sent to processor 0
169 : std::map<processor_id_type, std::vector<std::pair<RayID, dof_id_type>>> send_needed_nodes;
170 598 : if (local_ray_needed_nodes.size())
171 : {
172 563 : auto & root_entry = send_needed_nodes[0];
173 4424 : for (const auto & id_nodes_pair : local_ray_needed_nodes)
174 3861 : root_entry.emplace_back(id_nodes_pair);
175 : }
176 :
177 : // The global map of ray -> required nodes needed to be filled on processor 0
178 : std::map<RayID, dof_id_type> global_needed_nodes;
179 : // Keep track of what Ray IDs we received from what processors so we know who to send back to
180 : std::unordered_map<processor_id_type, std::set<RayID>> pid_received_ids;
181 : // Take the required nodes for each Ray and determine on processor 0 the global maximum
182 : // nodes needed to represent each Ray
183 : const auto append_global_max =
184 563 : [&global_needed_nodes, &pid_received_ids](
185 : processor_id_type pid, const std::vector<std::pair<RayID, dof_id_type>> & pairs)
186 : {
187 563 : auto & pid_received_entry = pid_received_ids[pid];
188 4424 : for (const auto & id_max_pair : pairs)
189 : {
190 3861 : const RayID ray_id = id_max_pair.first;
191 3861 : const dof_id_type max = id_max_pair.second;
192 :
193 3861 : const auto find = global_needed_nodes.find(ray_id);
194 3861 : if (find == global_needed_nodes.end())
195 3388 : global_needed_nodes.emplace(ray_id, max);
196 : else
197 : {
198 473 : auto & current_max = find->second;
199 473 : current_max = std::max(current_max, max);
200 : }
201 :
202 3861 : pid_received_entry.insert(ray_id);
203 : }
204 598 : };
205 598 : Parallel::push_parallel_vector_data(comm(), send_needed_nodes, append_global_max);
206 :
207 : // Decide on the starting representative starting node and elem ID for each Ray
208 : std::map<RayID, std::pair<dof_id_type, dof_id_type>> global_ids;
209 : dof_id_type current_node_id = 0;
210 : dof_id_type current_elem_id = 0;
211 3986 : for (auto & pair : global_needed_nodes)
212 : {
213 3388 : const RayID ray_id = pair.first;
214 3388 : const dof_id_type max_nodes = pair.second;
215 :
216 3388 : global_ids.emplace(ray_id, std::make_pair(current_node_id, current_elem_id));
217 :
218 3388 : current_node_id += max_nodes;
219 3388 : if (max_nodes > 1)
220 2281 : current_elem_id += max_nodes - 1;
221 : else // stationary
222 1107 : current_elem_id += 1;
223 : }
224 :
225 : // Share the max node IDs so that we have a starting point for unique IDs for elems
226 598 : _max_node_id = current_node_id;
227 598 : comm().max(_max_node_id);
228 :
229 : // Fill the starting ID information to each processor that needs it from processor 0
230 : std::unordered_map<processor_id_type, std::vector<std::tuple<RayID, dof_id_type, dof_id_type>>>
231 : send_ids;
232 1161 : for (auto & pid_ids_pair : pid_received_ids)
233 : {
234 563 : const processor_id_type pid = pid_ids_pair.first;
235 : const std::set<RayID> & ray_ids = pid_ids_pair.second;
236 :
237 : auto & pid_send = send_ids[pid];
238 4424 : for (const RayID ray_id : ray_ids)
239 : {
240 3861 : const auto & global_entry = global_ids.at(ray_id);
241 3861 : const dof_id_type node_id = global_entry.first;
242 3861 : const dof_id_type elem_id = global_entry.second;
243 3861 : pid_send.emplace_back(ray_id, node_id, elem_id);
244 : }
245 : }
246 :
247 : // Take the starting ID information from processor 0 and store it locally
248 : const auto append_ids =
249 563 : [&](processor_id_type,
250 : const std::vector<std::tuple<RayID, dof_id_type, dof_id_type>> & tuples)
251 : {
252 4424 : for (const auto & tuple : tuples)
253 : {
254 3861 : const RayID ray_id = std::get<0>(tuple);
255 3861 : const dof_id_type node_id = std::get<1>(tuple);
256 3861 : const dof_id_type elem_id = std::get<2>(tuple);
257 3861 : _ray_starting_id_map.emplace(ray_id, std::make_pair(node_id, elem_id));
258 : }
259 598 : };
260 :
261 : _ray_starting_id_map.clear();
262 598 : Parallel::push_parallel_vector_data(comm(), send_ids, append_ids);
263 598 : }
264 :
265 : void
266 598 : RayTracingMeshOutput::buildSegmentMesh()
267 : {
268 1196 : TIME_SECTION("buildSegmentMesh", 3, "Building RayTracing Mesh Output");
269 :
270 598 : _segmented_rays = false;
271 :
272 : // Tally nodes and elems for the local mesh ahead of time so we can reserve
273 : // Each segment requires an element, and we need one more node than elems
274 : dof_id_type num_nodes = 0;
275 : dof_id_type num_elems = 0;
276 4569 : for (const auto & entry : _study.getCachedTraces())
277 : {
278 : const auto num_segments = entry.numSegments();
279 3971 : num_nodes += num_segments + 1;
280 3971 : if (num_segments >= 1)
281 : {
282 2829 : _segmented_rays = true;
283 2829 : num_elems += num_segments;
284 : }
285 : else if (entry.stationary())
286 1107 : num_elems += 1;
287 : }
288 :
289 598 : comm().max(_segmented_rays);
290 :
291 : // Build the segment mesh
292 : mooseAssert(!_segment_mesh, "Not cleared");
293 598 : if (_communicator.size() == 1)
294 124 : _segment_mesh = std::make_unique<ReplicatedMesh>(_communicator, _mesh_ptr->dimension());
295 : else
296 : {
297 474 : _segment_mesh = std::make_unique<DistributedMesh>(_communicator, _mesh_ptr->dimension());
298 474 : _segment_mesh->set_distributed();
299 : }
300 : // We set neighbor links
301 : _segment_mesh->allow_find_neighbors(false);
302 : // Don't renumber so that we can remain consistent between processor counts
303 : _segment_mesh->allow_renumbering(false);
304 : // As we're building segments for just this partiton, we're partitioning on our own
305 : _segment_mesh->skip_partitioning(true);
306 : // Reserve nodes and elems
307 598 : _segment_mesh->reserve_nodes(num_nodes);
308 598 : _segment_mesh->reserve_elem(num_elems);
309 :
310 598 : buildIDMap();
311 :
312 : // Decide on a starting ID for each processor
313 : // Build a mesh segment for each local Ray segment
314 : // Each one of these objects represents a single Ray's path through this processor
315 4569 : for (const auto & trace_data : _study.getCachedTraces())
316 : {
317 : // If we have no segments, this is a trace that skimmed the corner of a processor boundary and
318 : // didn't contribute anything
319 3971 : if (trace_data.numSegments() == 0 && !trace_data.stationary())
320 35 : continue;
321 :
322 : dof_id_type node_id, elem_id;
323 3936 : startingIDs(trace_data, node_id, elem_id);
324 :
325 : // Add the start point
326 : mooseAssert(!_segment_mesh->query_node_ptr(node_id), "Node already exists");
327 : Node * last_node =
328 3936 : _segment_mesh->add_point(trace_data._point_data[0]._point, node_id, processor_id());
329 3936 : last_node->set_unique_id(node_id++);
330 :
331 : // Stationary, add a NodeElem
332 : if (trace_data.stationary())
333 : {
334 : mooseAssert(!_segment_mesh->query_elem_ptr(elem_id), "Elem already exists");
335 1107 : auto elem = _segment_mesh->add_elem(Elem::build_with_id(NODEELEM, elem_id));
336 : elem->processor_id(processor_id());
337 1107 : elem->set_unique_id(_max_node_id + elem_id++);
338 1107 : elem->set_node(0, last_node);
339 : }
340 : // Not stationary; add a point and element for each segment
341 : else
342 : {
343 : Elem * last_elem = nullptr;
344 :
345 12976 : for (std::size_t i = 1; i < trace_data._point_data.size(); ++i)
346 : {
347 10147 : const auto & point = trace_data._point_data[i]._point;
348 :
349 : // Add next point on the trace
350 : mooseAssert(!_segment_mesh->query_node_ptr(node_id), "Node already exists");
351 10147 : Node * node = _segment_mesh->add_point(point, node_id, processor_id());
352 10147 : node->set_unique_id(node_id++);
353 :
354 : // Build a segment from this point to the last
355 : mooseAssert(!_segment_mesh->query_elem_ptr(elem_id), "Elem already exists");
356 10147 : Elem * elem = _segment_mesh->add_elem(Elem::build_with_id(EDGE2, elem_id));
357 : elem->processor_id(processor_id());
358 10147 : elem->set_unique_id(_max_node_id + elem_id++);
359 10147 : elem->set_node(0, last_node);
360 10147 : elem->set_node(1, node);
361 :
362 : // Set neighbor links
363 10147 : if (last_elem)
364 : {
365 : elem->set_neighbor(0, last_elem);
366 : last_elem->set_neighbor(1, elem);
367 : }
368 :
369 : last_elem = elem;
370 : last_node = node;
371 : }
372 : }
373 : }
374 :
375 : // If the mesh is replicated, everything that follows is unnecessary. Prepare and be done.
376 598 : if (_segment_mesh->is_replicated())
377 : {
378 124 : _segment_mesh->prepare_for_use();
379 : return;
380 : }
381 :
382 : // Find the Nodes on processor boundaries that we need to decide on owners for. Also set up
383 : // remote_elem links for neighbors we clearly don't have. We don't build the neighbor maps at
384 : // all, so all we have are nullptr and remote_elem. We can only do all of this after the mesh is
385 : // built because the mesh isn't necessarily built in the order Rays are traced
386 : std::vector<Node *> need_node_owners;
387 3121 : for (const auto & trace_data : _study.getCachedTraces())
388 : {
389 : const auto num_segments = trace_data.numSegments();
390 :
391 2647 : if (num_segments == 0)
392 650 : continue;
393 :
394 : dof_id_type start_node_id, start_elem_id;
395 1997 : startingIDs(trace_data, start_node_id, start_elem_id);
396 :
397 : // Another part of the trace for this Ray happened before this one
398 2117 : if ((_study.segmentsOnCacheTraces() ? trace_data._intersections
399 120 : : (unsigned long int)(trace_data._processor_crossings +
400 120 : trace_data._trajectory_changes)) != 0)
401 : {
402 : // The element before start_elem (may be nullptr, meaning it didn't trace on this proc)
403 560 : Elem * previous_elem = _segment_mesh->query_elem_ptr(start_elem_id - 1);
404 :
405 : // We don't have the previous element segment, so it exists on another processor
406 : // Set the remote_elem and mark that we need to find out who owns the first node on this
407 : // trace
408 560 : if (!previous_elem)
409 : {
410 560 : Elem * start_elem = _segment_mesh->elem_ptr(start_elem_id);
411 560 : start_elem->set_neighbor(0, const_cast<RemoteElem *>(remote_elem));
412 : start_elem->node_ptr(0)->invalidate_processor_id();
413 560 : need_node_owners.push_back(start_elem->node_ptr(0));
414 : }
415 : }
416 :
417 : // If you have multiple sets of segments on a single processor, it is possible
418 : // to have part of the trace on another processor within the segments we know about.
419 : // In this case, we have to do some magic to the neighbor links so that the mesh
420 : // doesn't fail in assertions.
421 1997 : if (!trace_data._last)
422 : {
423 : // The element after end_elem (may be nullptr, meaning it didn't happen on this proc)
424 548 : Elem * next_elem = _segment_mesh->query_elem_ptr(start_elem_id + num_segments);
425 :
426 : // We have the next element from another trace, set neighbors as we own both
427 548 : if (!next_elem)
428 : {
429 548 : Elem * end_elem = _segment_mesh->elem_ptr(start_elem_id + num_segments - 1);
430 548 : end_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
431 : end_elem->node_ptr(1)->invalidate_processor_id();
432 548 : need_node_owners.push_back(end_elem->node_ptr(1));
433 : }
434 : }
435 : }
436 :
437 : // Sort through the neighboring bounding boxes and prepare requests to each processor that /may/
438 : // also have one of the nodes that we need to decide on an owner for
439 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> need_node_owners_sends;
440 1582 : for (const Node * node : need_node_owners)
441 2570 : for (const auto & neighbor_pid_bbox_pair : _inflated_neighbor_bboxes)
442 : {
443 1462 : const auto neighbor_pid = neighbor_pid_bbox_pair.first;
444 1462 : const auto & neighbor_bbox = neighbor_pid_bbox_pair.second;
445 1462 : if (neighbor_bbox.contains_point(*node))
446 1120 : need_node_owners_sends[neighbor_pid].push_back(node->id());
447 : }
448 :
449 : // Functor that takes in a set of incoming node IDs from a processor and decides on an owner.
450 : // For every node that we need an owner for, this should be called for said node on both
451 : // processors. Therefore, just pick the minimum processor ID as the owner. This should never
452 : // happen more than once for a single node because we're building 1D elems and we can only ever
453 : // have two procs that touch a node.
454 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> confirm_node_owners_sends;
455 : auto decide_node_owners_functor =
456 294 : [this, &confirm_node_owners_sends](processor_id_type pid,
457 : const std::vector<dof_id_type> & incoming_node_ids)
458 : {
459 294 : auto & confirm_pid = confirm_node_owners_sends[pid];
460 1414 : for (const auto & node_id : incoming_node_ids)
461 : {
462 1120 : Node * node = _segment_mesh->query_node_ptr(node_id);
463 1120 : if (node)
464 : {
465 : mooseAssert(!node->valid_processor_id(), "Should be invalid");
466 1072 : node->processor_id() = std::min(pid, processor_id());
467 1072 : confirm_pid.push_back(node_id);
468 : }
469 : }
470 768 : };
471 :
472 : // Ship the nodes that need owners for and pick an owner when we receive
473 474 : Parallel::push_parallel_vector_data(
474 474 : _communicator, need_node_owners_sends, decide_node_owners_functor);
475 :
476 : // At this point, any nodes that we actually need ghosts for should have their processor_ids
477 : // satisfied. Anything that is left isn't actually ghosted and is ours. Therefore, set
478 : // everything left that is invalid to be owned by this proc.
479 1582 : for (Node * node : need_node_owners)
480 1108 : if (!node->valid_processor_id())
481 36 : node->processor_id() = processor_id();
482 :
483 : // We're done!
484 474 : _segment_mesh->prepare_for_use();
485 474 : }
486 :
487 : void
488 598 : RayTracingMeshOutput::setupEquationSystem()
489 : {
490 1196 : TIME_SECTION("setupEquationSystem", 3, "Setting Up Ray Tracing MeshOutput Equation System");
491 :
492 1196 : _es = std::make_unique<EquationSystems>(*_segment_mesh);
493 598 : _sys = &_es->add_system<libMesh::ExplicitSystem>("sys");
494 598 : _data_vars.clear();
495 598 : _aux_data_vars.clear();
496 :
497 : // Add variables for the basic properties if enabled
498 598 : _ray_id_var = invalid_uint;
499 598 : _intersections_var = invalid_uint;
500 598 : _pid_var = invalid_uint;
501 598 : _processor_crossings_var = invalid_uint;
502 598 : _trajectory_changes_var = invalid_uint;
503 2008 : for (auto & prop : _pars.get<MultiMooseEnum>("output_properties"))
504 1410 : switch (prop)
505 : {
506 404 : case 0: // ray_id (stationary and segments)
507 404 : _ray_id_var = _sys->add_variable("ray_id", CONSTANT, MONOMIAL);
508 404 : break;
509 556 : case 1: // intersections (segments only)
510 556 : if (_segmented_rays)
511 556 : _intersections_var = _sys->add_variable("intersections", CONSTANT, MONOMIAL);
512 : break;
513 150 : case 2: // pid (stationary and segments)
514 150 : _pid_var = _sys->add_variable("pid", CONSTANT, MONOMIAL);
515 150 : break;
516 150 : case 3: // processor_crossings (segments only)
517 150 : if (_segmented_rays)
518 150 : _processor_crossings_var = _sys->add_variable("processor_crossings", CONSTANT, MONOMIAL);
519 : break;
520 150 : case 4: // trajectory_changes (segments only)
521 150 : if (_segmented_rays)
522 150 : _trajectory_changes_var = _sys->add_variable("trajectory_changes", CONSTANT, MONOMIAL);
523 : break;
524 0 : default:
525 0 : mooseError("Invalid property");
526 : }
527 :
528 1194 : const auto get_data_vars = [this](const bool aux)
529 : {
530 : // The data index -> variable result
531 : std::vector<std::pair<RayDataIndex, unsigned int>> vars;
532 :
533 : const auto from_names =
534 1194 : aux ? (_output_aux_data ? &_study.rayAuxDataNames() : _output_aux_data_names)
535 598 : : (_output_data ? &_study.rayDataNames() : _output_data_names);
536 :
537 : // Nothing to output
538 976 : if (!from_names)
539 : return vars;
540 :
541 264 : const auto output_data_nodal = aux ? false : _output_data_nodal;
542 :
543 566 : for (const auto & name : *from_names)
544 : {
545 : const auto data_index =
546 306 : aux ? _study.getRayAuxDataIndex(name, true) : _study.getRayDataIndex(name, true);
547 306 : if (data_index == Ray::INVALID_RAY_DATA_INDEX)
548 : {
549 6 : const std::string names_param = aux ? "output_aux_data_names" : "output_data_names";
550 6 : const std::string data_prefix = aux ? "aux " : "";
551 4 : paramError(names_param, "The ray ", data_prefix, "data '", name, "' is not registered");
552 : }
553 :
554 508 : const std::string var_prefix = aux ? "aux_" : "";
555 302 : if (output_data_nodal)
556 108 : _sys->add_variable(var_prefix + name, FIRST, LAGRANGE);
557 : else
558 496 : _sys->add_variable(var_prefix + name, CONSTANT, MONOMIAL);
559 302 : const auto var_num = _sys->variable_number(var_prefix + name);
560 302 : vars.emplace_back(data_index, var_num);
561 : }
562 :
563 : return vars;
564 0 : };
565 :
566 1194 : _data_vars = get_data_vars(false);
567 1190 : _aux_data_vars = get_data_vars(true);
568 :
569 : // All done
570 594 : _es->init();
571 594 : }
572 :
573 : void
574 594 : RayTracingMeshOutput::fillFields()
575 : {
576 1188 : TIME_SECTION("fillFields", 3, "Filling RayTracing MeshOutput Fields");
577 :
578 : // Helper for setting the solution (if enabled)
579 : const auto set_solution =
580 58413 : [this](const DofObject * const dof, const unsigned int var, const Real value)
581 : {
582 : mooseAssert(dof, "Nullptr dof");
583 58413 : if (var != invalid_uint)
584 : {
585 27342 : const auto dof_number = dof->dof_number(_sys->number(), var, 0);
586 27342 : _sys->solution->set(dof_number, value);
587 : }
588 59007 : };
589 :
590 : // Helper for filling data and aux data (if enabled)
591 : const auto fill_data =
592 22683 : [&set_solution](const DofObject * const dof, const auto & vars, const auto & data)
593 : {
594 : mooseAssert(dof, "Nullptr dof");
595 28167 : for (const auto & [data_index, var_num] : vars)
596 5484 : set_solution(dof, var_num, data[data_index]);
597 23277 : };
598 :
599 4561 : for (const auto & trace_data : _study.getCachedTraces())
600 : {
601 3967 : auto intersection = trace_data._intersections;
602 :
603 : // No segments and not stationary; means this ray bounced off
604 : // this processor and never actually moved on this processor
605 3967 : if (!trace_data.numSegments() && !trace_data.stationary())
606 35 : continue;
607 :
608 : dof_id_type node_id, elem_id;
609 3932 : startingIDs(trace_data, node_id, elem_id);
610 :
611 3932 : const Elem * elem = _segment_mesh->elem_ptr(elem_id);
612 :
613 : // Fill first node's nodal data if we need it; the loop that follows will handle
614 : // the rest of the nodes
615 3932 : if (_output_data_nodal && _study.hasRayData() && !trace_data.stationary())
616 183 : fill_data(elem->node_ptr(0), _data_vars, trace_data._point_data[0]._data);
617 :
618 : const std::size_t start = trace_data.stationary() ? 0 : 1;
619 15182 : for (const auto i : make_range(start, trace_data._point_data.size()))
620 : {
621 : mooseAssert(elem, "Nullptr elem");
622 :
623 : // Elemental ID and pid
624 11250 : set_solution(elem, _ray_id_var, trace_data._ray_id);
625 11250 : set_solution(elem, _pid_var, processor_id());
626 : // Elemental properties that only apply to segments
627 : if (!trace_data.stationary())
628 : {
629 10143 : set_solution(elem, _intersections_var, intersection++);
630 10143 : set_solution(elem, _processor_crossings_var, trace_data._processor_crossings);
631 10143 : set_solution(elem, _trajectory_changes_var, trace_data._trajectory_changes);
632 : }
633 : const auto & point_data = trace_data._point_data[i];
634 : // Data fields
635 11250 : if (_output_data_nodal && !trace_data.stationary())
636 570 : fill_data(elem->node_ptr(1), _data_vars, point_data._data);
637 : else
638 10680 : fill_data(elem, _data_vars, point_data._data);
639 : // Aux data fields
640 11250 : fill_data(elem, _aux_data_vars, point_data._aux_data);
641 :
642 : // Advance to the next element
643 : if (!trace_data.stationary())
644 : elem = elem->neighbor_ptr(1);
645 : }
646 : }
647 :
648 594 : _sys->solution->close();
649 594 : _sys->update();
650 594 : }
651 :
652 : void
653 598 : RayTracingMeshOutput::buildBoundingBoxes()
654 : {
655 1196 : TIME_SECTION("buildBoundingBoxes", 3, "Building Bounding Boxes for RayTracing Mesh Output");
656 :
657 : // Not used in the one proc case
658 598 : if (_communicator.size() == 1)
659 : return;
660 :
661 : // Local bounding box
662 474 : _bbox = MeshTools::create_local_bounding_box(_mesh_ptr->getMesh());
663 :
664 : // Gather the bounding boxes of all processors
665 474 : std::vector<std::pair<Point, Point>> bb_points = {static_cast<std::pair<Point, Point>>(_bbox)};
666 474 : _communicator.allgather(bb_points, true);
667 474 : _inflated_bboxes.resize(_communicator.size());
668 1680 : for (processor_id_type pid = 0; pid < _communicator.size(); ++pid)
669 : {
670 1206 : BoundingBox pid_bbox = static_cast<BoundingBox>(bb_points[pid]);
671 1206 : pid_bbox.scale(0.01);
672 : _inflated_bboxes[pid] = pid_bbox;
673 : }
674 :
675 : // Find intersecting (neighbor) bounding boxes
676 474 : _inflated_neighbor_bboxes.clear();
677 1680 : for (processor_id_type pid = 0; pid < _communicator.size(); ++pid)
678 1206 : if (pid != processor_id())
679 : {
680 : // Insert if the searched processor's bbox intersects my bbox
681 732 : const auto & pid_bbox = _inflated_bboxes[pid];
682 732 : if (_bbox.intersects(pid_bbox))
683 720 : _inflated_neighbor_bboxes.emplace_back(pid, pid_bbox);
684 : }
685 474 : }
686 :
687 : void
688 9865 : RayTracingMeshOutput::startingIDs(const TraceData & trace_data,
689 : dof_id_type & start_node_id,
690 : dof_id_type & start_elem_id) const
691 : {
692 9865 : const auto [begin_node_id, begin_elem_id] = _ray_starting_id_map.at(trace_data._ray_id);
693 :
694 : const auto offset =
695 : trace_data.stationary()
696 : ? 1
697 7651 : : (_study.segmentsOnCacheTraces()
698 7651 : ? trace_data._intersections
699 360 : : (trace_data._processor_crossings + trace_data._trajectory_changes));
700 :
701 9865 : start_node_id = begin_node_id + offset;
702 9865 : start_elem_id = begin_elem_id + offset;
703 9865 : }
704 :
705 : dof_id_type
706 3971 : RayTracingMeshOutput::neededNodes(const TraceData & trace_data) const
707 : {
708 3971 : if (_study.segmentsOnCacheTraces())
709 3827 : return trace_data._intersections + trace_data._point_data.size();
710 144 : return trace_data._processor_crossings + trace_data._trajectory_changes +
711 144 : trace_data._point_data.size();
712 : }
|
