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 500 : RayTracingMeshOutput::validParams()
31 : {
32 500 : InputParameters params = FileOutput::validParams();
33 :
34 1000 : params.addRequiredParam<UserObjectName>("study", "The RayTracingStudy to get the segments from");
35 :
36 1000 : params.addParam<bool>("output_data", false, "Whether or not to also output the Ray's data");
37 1000 : params.addParam<std::vector<std::string>>("output_data_names",
38 : "The names of specific data to output");
39 1000 : params.addParam<bool>("output_data_nodal",
40 1000 : 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 1000 : params.addParam<bool>(
45 1000 : "output_aux_data", false, "Whether or not to also output the Ray's aux data");
46 1000 : 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 500 : "ray_id intersections");
51 1000 : params.addParam<MultiMooseEnum>("output_properties", props, "Which Ray properties to output");
52 :
53 500 : return params;
54 500 : }
55 :
56 256 : RayTracingMeshOutput::RayTracingMeshOutput(const InputParameters & params)
57 : : FileOutput(params),
58 : UserObjectInterface(this),
59 256 : _study(getUserObject<RayTracingStudy>("study")),
60 512 : _output_data(getParam<bool>("output_data")),
61 256 : _output_data_names(isParamValid("output_data_names")
62 300 : ? &getParam<std::vector<std::string>>("output_data_names")
63 : : nullptr),
64 512 : _output_data_nodal(getParam<bool>("output_data_nodal")),
65 512 : _output_aux_data(getParam<bool>("output_aux_data")),
66 512 : _output_aux_data_names(isParamValid("output_aux_data_names")
67 282 : ? &getParam<std::vector<std::string>>("output_aux_data_names")
68 : : nullptr),
69 256 : _ray_id_var(invalid_uint),
70 256 : _intersections_var(invalid_uint),
71 256 : _pid_var(invalid_uint),
72 256 : _processor_crossings_var(invalid_uint),
73 256 : _trajectory_changes_var(invalid_uint),
74 512 : _segmented_rays(false)
75 : {
76 256 : 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 254 : 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 252 : 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 250 : if (_output_data_nodal && !_study.segmentsOnCacheTraces())
93 2 : paramError("output_data_nodal", "Not supported when study segments_on_cache_traces = false");
94 248 : if (_output_data && _output_data_names)
95 2 : paramError("output_data", "Cannot be used in addition to 'output_data_names'; choose one");
96 246 : 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 244 : }
100 :
101 : std::string
102 629 : RayTracingMeshOutput::filename()
103 : {
104 : // Append the file extension on the base file name
105 629 : std::ostringstream output;
106 629 : output << _file_base << fileExtension();
107 :
108 : // Add the _000x extension to the file
109 629 : if (_file_num > 0)
110 175 : output << "-s" << std::setw(_padding) << std::setprecision(0) << std::setfill('0') << std::right
111 175 : << _file_num;
112 :
113 : // Return the filename
114 629 : return output.str();
115 629 : }
116 :
117 : void
118 391 : RayTracingMeshOutput::output()
119 : {
120 : // Do we even have any traces?
121 391 : auto num_segments = _study.getCachedTraces().size();
122 391 : _communicator.sum(num_segments);
123 391 : 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 389 : buildBoundingBoxes();
128 :
129 : // Build the _segment_mesh
130 389 : buildSegmentMesh();
131 :
132 : // Setup the system to store the Ray field data
133 389 : setupEquationSystem();
134 :
135 : // Fill the field data
136 385 : fillFields();
137 :
138 : // And output
139 385 : 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 385 : _sys = nullptr;
146 : _segment_mesh = nullptr;
147 385 : }
148 :
149 : void
150 389 : RayTracingMeshOutput::buildIDMap()
151 : {
152 778 : 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 3031 : for (const auto & trace_data : _study.getCachedTraces())
157 : {
158 2642 : const auto find = local_ray_needed_nodes.find(trace_data._ray_id);
159 2642 : if (find == local_ray_needed_nodes.end())
160 2601 : local_ray_needed_nodes[trace_data._ray_id] = neededNodes(trace_data);
161 : else
162 : {
163 41 : auto & value = find->second;
164 82 : 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 389 : if (local_ray_needed_nodes.size())
171 : {
172 362 : auto & root_entry = send_needed_nodes[0];
173 2963 : for (const auto & id_nodes_pair : local_ray_needed_nodes)
174 2601 : 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 362 : [&global_needed_nodes, &pid_received_ids](
185 : processor_id_type pid, const std::vector<std::pair<RayID, dof_id_type>> & pairs)
186 : {
187 362 : auto & pid_received_entry = pid_received_ids[pid];
188 2963 : for (const auto & id_max_pair : pairs)
189 : {
190 2601 : const RayID ray_id = id_max_pair.first;
191 2601 : const dof_id_type max = id_max_pair.second;
192 :
193 2601 : const auto find = global_needed_nodes.find(ray_id);
194 2601 : if (find == global_needed_nodes.end())
195 2329 : global_needed_nodes.emplace(ray_id, max);
196 : else
197 : {
198 272 : auto & current_max = find->second;
199 272 : current_max = std::max(current_max, max);
200 : }
201 :
202 2601 : pid_received_entry.insert(ray_id);
203 : }
204 389 : };
205 389 : 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 2718 : for (auto & pair : global_needed_nodes)
212 : {
213 2329 : const RayID ray_id = pair.first;
214 2329 : const dof_id_type max_nodes = pair.second;
215 :
216 2329 : global_ids.emplace(ray_id, std::make_pair(current_node_id, current_elem_id));
217 :
218 2329 : current_node_id += max_nodes;
219 2329 : if (max_nodes > 1)
220 1591 : current_elem_id += max_nodes - 1;
221 : else // stationary
222 738 : 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 389 : _max_node_id = current_node_id;
227 389 : 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 751 : for (auto & pid_ids_pair : pid_received_ids)
233 : {
234 362 : 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 2963 : for (const RayID ray_id : ray_ids)
239 : {
240 2601 : const auto & global_entry = global_ids.at(ray_id);
241 2601 : const dof_id_type node_id = global_entry.first;
242 2601 : const dof_id_type elem_id = global_entry.second;
243 2601 : 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 362 : [&](processor_id_type,
250 : const std::vector<std::tuple<RayID, dof_id_type, dof_id_type>> & tuples)
251 : {
252 2963 : for (const auto & tuple : tuples)
253 : {
254 2601 : const RayID ray_id = std::get<0>(tuple);
255 2601 : const dof_id_type node_id = std::get<1>(tuple);
256 2601 : const dof_id_type elem_id = std::get<2>(tuple);
257 2601 : _ray_starting_id_map.emplace(ray_id, std::make_pair(node_id, elem_id));
258 : }
259 389 : };
260 :
261 : _ray_starting_id_map.clear();
262 389 : Parallel::push_parallel_vector_data(comm(), send_ids, append_ids);
263 389 : }
264 :
265 : void
266 389 : RayTracingMeshOutput::buildSegmentMesh()
267 : {
268 778 : TIME_SECTION("buildSegmentMesh", 3, "Building RayTracing Mesh Output");
269 :
270 389 : _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 3031 : for (const auto & entry : _study.getCachedTraces())
277 : {
278 : const auto num_segments = entry.numSegments();
279 2642 : num_nodes += num_segments + 1;
280 2642 : if (num_segments >= 1)
281 : {
282 1880 : _segmented_rays = true;
283 1880 : num_elems += num_segments;
284 : }
285 : else if (entry.stationary())
286 738 : num_elems += 1;
287 : }
288 :
289 389 : comm().max(_segmented_rays);
290 :
291 : // Build the segment mesh
292 : mooseAssert(!_segment_mesh, "Not cleared");
293 389 : if (_communicator.size() == 1)
294 124 : _segment_mesh = std::make_unique<ReplicatedMesh>(_communicator, _mesh_ptr->dimension());
295 : else
296 : {
297 265 : _segment_mesh = std::make_unique<DistributedMesh>(_communicator, _mesh_ptr->dimension());
298 265 : _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 389 : _segment_mesh->reserve_nodes(num_nodes);
308 389 : _segment_mesh->reserve_elem(num_elems);
309 :
310 389 : 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 3031 : 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 2642 : if (trace_data.numSegments() == 0 && !trace_data.stationary())
320 24 : continue;
321 :
322 : dof_id_type node_id, elem_id;
323 2618 : 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 2618 : _segment_mesh->add_point(trace_data._point_data[0]._point, node_id, processor_id());
329 2618 : 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 738 : auto elem = _segment_mesh->add_elem(Elem::build_with_id(NODEELEM, elem_id));
336 : elem->processor_id(processor_id());
337 738 : elem->set_unique_id(_max_node_id + elem_id++);
338 738 : 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 8850 : for (std::size_t i = 1; i < trace_data._point_data.size(); ++i)
346 : {
347 6970 : 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 6970 : Node * node = _segment_mesh->add_point(point, node_id, processor_id());
352 6970 : 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 6970 : Elem * elem = _segment_mesh->add_elem(Elem::build_with_id(EDGE2, elem_id));
357 : elem->processor_id(processor_id());
358 6970 : elem->set_unique_id(_max_node_id + elem_id++);
359 6970 : elem->set_node(0, last_node);
360 6970 : elem->set_node(1, node);
361 :
362 : // Set neighbor links
363 6970 : 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 389 : 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 1583 : for (const auto & trace_data : _study.getCachedTraces())
388 : {
389 : const auto num_segments = trace_data.numSegments();
390 :
391 1318 : if (num_segments == 0)
392 270 : continue;
393 :
394 : dof_id_type start_node_id, start_elem_id;
395 1048 : startingIDs(trace_data, start_node_id, start_elem_id);
396 :
397 : // Another part of the trace for this Ray happened before this one
398 1148 : if ((_study.segmentsOnCacheTraces() ? trace_data._intersections
399 100 : : (unsigned long int)(trace_data._processor_crossings +
400 100 : trace_data._trajectory_changes)) != 0)
401 : {
402 : // The element before start_elem (may be nullptr, meaning it didn't trace on this proc)
403 299 : 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 299 : if (!previous_elem)
409 : {
410 299 : Elem * start_elem = _segment_mesh->elem_ptr(start_elem_id);
411 299 : start_elem->set_neighbor(0, const_cast<RemoteElem *>(remote_elem));
412 : start_elem->node_ptr(0)->invalidate_processor_id();
413 299 : 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 1048 : if (!trace_data._last)
422 : {
423 : // The element after end_elem (may be nullptr, meaning it didn't happen on this proc)
424 289 : 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 289 : if (!next_elem)
428 : {
429 289 : Elem * end_elem = _segment_mesh->elem_ptr(start_elem_id + num_segments - 1);
430 289 : end_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
431 : end_elem->node_ptr(1)->invalidate_processor_id();
432 289 : 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 853 : for (const Node * node : need_node_owners)
441 1471 : for (const auto & neighbor_pid_bbox_pair : _inflated_neighbor_bboxes)
442 : {
443 883 : const auto neighbor_pid = neighbor_pid_bbox_pair.first;
444 883 : const auto & neighbor_bbox = neighbor_pid_bbox_pair.second;
445 883 : if (neighbor_bbox.contains_point(*node))
446 598 : 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 183 : [this, &confirm_node_owners_sends](processor_id_type pid,
457 : const std::vector<dof_id_type> & incoming_node_ids)
458 : {
459 183 : auto & confirm_pid = confirm_node_owners_sends[pid];
460 781 : for (const auto & node_id : incoming_node_ids)
461 : {
462 598 : Node * node = _segment_mesh->query_node_ptr(node_id);
463 598 : if (node)
464 : {
465 : mooseAssert(!node->valid_processor_id(), "Should be invalid");
466 558 : node->processor_id() = std::min(pid, processor_id());
467 558 : confirm_pid.push_back(node_id);
468 : }
469 : }
470 448 : };
471 :
472 : // Ship the nodes that need owners for and pick an owner when we receive
473 265 : Parallel::push_parallel_vector_data(
474 265 : _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 853 : for (Node * node : need_node_owners)
480 588 : if (!node->valid_processor_id())
481 30 : node->processor_id() = processor_id();
482 :
483 : // We're done!
484 265 : _segment_mesh->prepare_for_use();
485 389 : }
486 :
487 : void
488 389 : RayTracingMeshOutput::setupEquationSystem()
489 : {
490 778 : TIME_SECTION("setupEquationSystem", 3, "Setting Up Ray Tracing MeshOutput Equation System");
491 :
492 778 : _es = std::make_unique<EquationSystems>(*_segment_mesh);
493 389 : _sys = &_es->add_system<libMesh::ExplicitSystem>("sys");
494 389 : _data_vars.clear();
495 389 : _aux_data_vars.clear();
496 :
497 : // Add variables for the basic properties if enabled
498 389 : _ray_id_var = invalid_uint;
499 389 : _intersections_var = invalid_uint;
500 389 : _pid_var = invalid_uint;
501 389 : _processor_crossings_var = invalid_uint;
502 389 : _trajectory_changes_var = invalid_uint;
503 1382 : for (auto & prop : _pars.get<MultiMooseEnum>("output_properties"))
504 993 : switch (prop)
505 : {
506 253 : case 0: // ray_id (stationary and segments)
507 253 : _ray_id_var = _sys->add_variable("ray_id", CONSTANT, MONOMIAL);
508 253 : break;
509 365 : case 1: // intersections (segments only)
510 365 : if (_segmented_rays)
511 365 : _intersections_var = _sys->add_variable("intersections", CONSTANT, MONOMIAL);
512 : break;
513 125 : case 2: // pid (stationary and segments)
514 125 : _pid_var = _sys->add_variable("pid", CONSTANT, MONOMIAL);
515 125 : break;
516 125 : case 3: // processor_crossings (segments only)
517 125 : if (_segmented_rays)
518 125 : _processor_crossings_var = _sys->add_variable("processor_crossings", CONSTANT, MONOMIAL);
519 : break;
520 125 : case 4: // trajectory_changes (segments only)
521 125 : if (_segmented_rays)
522 125 : _trajectory_changes_var = _sys->add_variable("trajectory_changes", CONSTANT, MONOMIAL);
523 : break;
524 0 : default:
525 0 : mooseError("Invalid property");
526 : }
527 :
528 776 : 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 776 : aux ? (_output_aux_data ? &_study.rayAuxDataNames() : _output_aux_data_names)
535 389 : : (_output_data ? &_study.rayDataNames() : _output_data_names);
536 :
537 : // Nothing to output
538 642 : if (!from_names)
539 : return vars;
540 :
541 162 : const auto output_data_nodal = aux ? false : _output_data_nodal;
542 :
543 344 : for (const auto & name : *from_names)
544 : {
545 : const auto data_index =
546 186 : aux ? _study.getRayAuxDataIndex(name, true) : _study.getRayDataIndex(name, true);
547 186 : 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 306 : const std::string var_prefix = aux ? "aux_" : "";
555 182 : if (output_data_nodal)
556 68 : _sys->add_variable(var_prefix + name, FIRST, LAGRANGE);
557 : else
558 296 : _sys->add_variable(var_prefix + name, CONSTANT, MONOMIAL);
559 182 : const auto var_num = _sys->variable_number(var_prefix + name);
560 182 : vars.emplace_back(data_index, var_num);
561 : }
562 :
563 : return vars;
564 0 : };
565 :
566 776 : _data_vars = get_data_vars(false);
567 772 : _aux_data_vars = get_data_vars(true);
568 :
569 : // All done
570 385 : _es->init();
571 385 : }
572 :
573 : void
574 385 : RayTracingMeshOutput::fillFields()
575 : {
576 770 : TIME_SECTION("fillFields", 3, "Filling RayTracing MeshOutput Fields");
577 :
578 : // Helper for setting the solution (if enabled)
579 : const auto set_solution =
580 40018 : [this](const DofObject * const dof, const unsigned int var, const Real value)
581 : {
582 : mooseAssert(dof, "Nullptr dof");
583 40018 : if (var != invalid_uint)
584 : {
585 19092 : const auto dof_number = dof->dof_number(_sys->number(), var, 0);
586 19092 : _sys->solution->set(dof_number, value);
587 : }
588 40403 : };
589 :
590 : // Helper for filling data and aux data (if enabled)
591 : const auto fill_data =
592 15526 : [&set_solution](const DofObject * const dof, const auto & vars, const auto & data)
593 : {
594 : mooseAssert(dof, "Nullptr dof");
595 19238 : for (const auto & [data_index, var_num] : vars)
596 3712 : set_solution(dof, var_num, data[data_index]);
597 15911 : };
598 :
599 3023 : for (const auto & trace_data : _study.getCachedTraces())
600 : {
601 2638 : 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 2638 : if (!trace_data.numSegments() && !trace_data.stationary())
606 24 : continue;
607 :
608 : dof_id_type node_id, elem_id;
609 2614 : startingIDs(trace_data, node_id, elem_id);
610 :
611 2614 : 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 2614 : if (_output_data_nodal && _study.hasRayData() && !trace_data.stationary())
616 118 : 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 10318 : 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 7704 : set_solution(elem, _ray_id_var, trace_data._ray_id);
625 7704 : set_solution(elem, _pid_var, processor_id());
626 : // Elemental properties that only apply to segments
627 : if (!trace_data.stationary())
628 : {
629 6966 : set_solution(elem, _intersections_var, intersection++);
630 6966 : set_solution(elem, _processor_crossings_var, trace_data._processor_crossings);
631 6966 : 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 7704 : if (_output_data_nodal && !trace_data.stationary())
636 400 : fill_data(elem->node_ptr(1), _data_vars, point_data._data);
637 : else
638 7304 : fill_data(elem, _data_vars, point_data._data);
639 : // Aux data fields
640 7704 : 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 385 : _sys->solution->close();
649 385 : _sys->update();
650 385 : }
651 :
652 : void
653 389 : RayTracingMeshOutput::buildBoundingBoxes()
654 : {
655 778 : TIME_SECTION("buildBoundingBoxes", 3, "Building Bounding Boxes for RayTracing Mesh Output");
656 :
657 : // Not used in the one proc case
658 389 : if (_communicator.size() == 1)
659 : return;
660 :
661 : // Local bounding box
662 265 : _bbox = MeshTools::create_local_bounding_box(_mesh_ptr->getMesh());
663 :
664 : // Gather the bounding boxes of all processors
665 265 : std::vector<std::pair<Point, Point>> bb_points = {static_cast<std::pair<Point, Point>>(_bbox)};
666 265 : _communicator.allgather(bb_points, true);
667 265 : _inflated_bboxes.resize(_communicator.size());
668 1010 : for (processor_id_type pid = 0; pid < _communicator.size(); ++pid)
669 : {
670 745 : BoundingBox pid_bbox = static_cast<BoundingBox>(bb_points[pid]);
671 745 : pid_bbox.scale(0.01);
672 : _inflated_bboxes[pid] = pid_bbox;
673 : }
674 :
675 : // Find intersecting (neighbor) bounding boxes
676 265 : _inflated_neighbor_bboxes.clear();
677 1010 : for (processor_id_type pid = 0; pid < _communicator.size(); ++pid)
678 745 : if (pid != processor_id())
679 : {
680 : // Insert if the searched processor's bbox intersects my bbox
681 480 : const auto & pid_bbox = _inflated_bboxes[pid];
682 480 : if (_bbox.intersects(pid_bbox))
683 470 : _inflated_neighbor_bboxes.emplace_back(pid, pid_bbox);
684 : }
685 389 : }
686 :
687 : void
688 6280 : RayTracingMeshOutput::startingIDs(const TraceData & trace_data,
689 : dof_id_type & start_node_id,
690 : dof_id_type & start_elem_id) const
691 : {
692 6280 : 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 4804 : : (_study.segmentsOnCacheTraces()
698 4804 : ? trace_data._intersections
699 300 : : (trace_data._processor_crossings + trace_data._trajectory_changes));
700 :
701 6280 : start_node_id = begin_node_id + offset;
702 6280 : start_elem_id = begin_elem_id + offset;
703 6280 : }
704 :
705 : dof_id_type
706 2642 : RayTracingMeshOutput::neededNodes(const TraceData & trace_data) const
707 : {
708 2642 : if (_study.segmentsOnCacheTraces())
709 2522 : return trace_data._intersections + trace_data._point_data.size();
710 120 : return trace_data._processor_crossings + trace_data._trajectory_changes +
711 120 : trace_data._point_data.size();
712 : }
|
