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 "ViewFactorRayStudy.h"
11 :
12 : // Local includes
13 : #include "ViewFactorRayBC.h"
14 : #include "GeometryUtils.h"
15 :
16 : // libMesh includes
17 : #include "libmesh/parallel_algebra.h"
18 : #include "libmesh/parallel_sync.h"
19 : #include "libmesh/enum_quadrature_type.h"
20 : #include "libmesh/fe_base.h"
21 : #include "libmesh/quadrature.h"
22 :
23 : // Ray tracing includes
24 : #include "ReflectRayBC.h"
25 : #include "RayTracingPackingUtils.h"
26 :
27 : using namespace libMesh;
28 :
29 : registerMooseObject("HeatTransferApp", ViewFactorRayStudy);
30 :
31 : InputParameters
32 330 : ViewFactorRayStudy::validParams()
33 : {
34 330 : auto params = RayTracingStudy::validParams();
35 :
36 660 : params.addRequiredParam<std::vector<BoundaryName>>(
37 : "boundary", "The list of boundaries where view factors are desired");
38 :
39 : MooseEnum qorders("CONSTANT FIRST SECOND THIRD FOURTH FIFTH SIXTH SEVENTH EIGHTH NINTH TENTH "
40 : "ELEVENTH TWELFTH THIRTEENTH FOURTEENTH FIFTEENTH SIXTEENTH SEVENTEENTH "
41 : "EIGHTTEENTH NINTEENTH TWENTIETH",
42 660 : "CONSTANT");
43 660 : params.addParam<MooseEnum>("face_order", qorders, "The face quadrature rule order");
44 :
45 660 : MooseEnum qtypes("GAUSS GRID", "GRID");
46 660 : params.addParam<MooseEnum>("face_type", qtypes, "The face quadrature type");
47 :
48 660 : MooseEnum convention("positive=0 negative=1", "positive");
49 660 : params.addParam<MooseEnum>(
50 : "internal_convention",
51 : convention,
52 : "The convention for spawning rays from internal sidesets; denotes the sign of the dot "
53 : "product between a ray and the internal sideset side normal");
54 :
55 660 : params.addParam<unsigned int>(
56 : "polar_quad_order",
57 660 : 16,
58 : "Order of the polar quadrature [polar angle is between ray and normal]. Must be even.");
59 660 : params.addParam<unsigned int>(
60 : "azimuthal_quad_order",
61 660 : 8,
62 : "Order of the azimuthal quadrature per quadrant [azimuthal angle is measured in "
63 : "a plane perpendicular to the normal].");
64 :
65 : // Shouldn't ever need RayKernels for view factors
66 330 : params.set<bool>("ray_kernel_coverage_check") = false;
67 330 : params.suppressParameter<bool>("ray_kernel_coverage_check");
68 :
69 : // So that the study executes before the RayTracingViewFactor
70 330 : params.set<bool>("force_preaux") = true;
71 330 : params.suppressParameter<bool>("force_preaux");
72 :
73 : // Need to use internal sidesets
74 330 : params.set<bool>("use_internal_sidesets") = true;
75 330 : params.suppressParameter<bool>("use_internal_sidesets");
76 :
77 : // Don't verify Rays in opt mode by default - it's expensive
78 330 : params.set<bool>("verify_rays") = false;
79 :
80 : // No need to use Ray registration
81 330 : params.set<bool>("_use_ray_registration") = false;
82 : // Do not need to bank Rays on completion
83 330 : params.set<bool>("_bank_rays_on_completion") = false;
84 :
85 330 : params.addClassDescription(
86 : "This ray study is used to compute view factors in cavities with obstruction. It sends out "
87 : "rays from surfaces bounding the radiation cavity into a set of directions determined by an "
88 : "angular quadrature. The rays are tracked and view factors are computed by determining the "
89 : "surface where the ray dies.");
90 330 : return params;
91 330 : }
92 :
93 165 : ViewFactorRayStudy::ViewFactorRayStudy(const InputParameters & parameters)
94 : : RayTracingStudy(parameters),
95 165 : _bnd_ids_vec(_mesh.getBoundaryIDs(getParam<std::vector<BoundaryName>>("boundary"))),
96 165 : _bnd_ids(_bnd_ids_vec.begin(), _bnd_ids_vec.end()),
97 330 : _internal_convention(getParam<MooseEnum>("internal_convention")),
98 165 : _ray_index_start_bnd_id(registerRayAuxData("start_bnd_id")),
99 165 : _ray_index_start_total_weight(registerRayAuxData("start_total_weight")),
100 165 : _fe_face(FEBase::build(_mesh.dimension(), FEType(CONSTANT, MONOMIAL))),
101 660 : _q_face(QBase::build(Moose::stringToEnum<QuadratureType>(getParam<MooseEnum>("face_type")),
102 165 : _mesh.dimension() - 1,
103 165 : Moose::stringToEnum<Order>(getParam<MooseEnum>("face_order")))),
104 165 : _is_3d(_mesh.dimension() == 3),
105 330 : _threaded_vf_info(libMesh::n_threads())
106 : {
107 165 : _fe_face->attach_quadrature_rule(_q_face.get());
108 : _fe_face->get_xyz();
109 :
110 : // create angular quadrature
111 165 : if (!_is_3d)
112 : {
113 : // In 2D, we integrate over angle theta instead of mu = cos(theta)
114 : // The integral over theta is approximated using a Gauss Legendre
115 : // quadrature. The integral we need to approximate is given by:
116 : //
117 : // int_{-pi/2}^{pi/2} cos(theta) d theta
118 : //
119 : // We get abscissae x and weight w for range of integration
120 : // from 0 to 1 and then rescale it to the integration range
121 : //
122 : std::vector<Real> x;
123 : std::vector<Real> w;
124 63 : RayTracingAngularQuadrature::gaussLegendre(
125 126 : 2 * getParam<unsigned int>("polar_quad_order"), x, w);
126 :
127 63 : _2d_aq_angles.resize(x.size());
128 63 : _2d_aq_weights.resize(x.size());
129 4511 : for (unsigned int j = 0; j < x.size(); ++j)
130 : {
131 4448 : _2d_aq_angles[j] = (2 * x[j] - 1) * M_PI / 2;
132 4448 : _2d_aq_weights[j] = w[j] * M_PI;
133 : }
134 63 : _num_dir = _2d_aq_angles.size();
135 63 : }
136 : else
137 : {
138 204 : _3d_aq = std::make_unique<RayTracingAngularQuadrature>(
139 204 : _mesh.dimension(),
140 : getParam<unsigned int>("polar_quad_order"),
141 102 : 4 * getParam<unsigned int>("azimuthal_quad_order"),
142 204 : /* mu_min = */ 0,
143 204 : /* mu_max = */ 1);
144 :
145 102 : _num_dir = _3d_aq->numDirections();
146 : }
147 165 : }
148 :
149 : void
150 165 : ViewFactorRayStudy::initialSetup()
151 : {
152 165 : RayTracingStudy::initialSetup();
153 :
154 : // We optimized away RayKernels, so don't allow them
155 165 : if (hasRayKernels(/* tid = */ 0))
156 0 : mooseError("Not compatible with RayKernels.");
157 :
158 : // RayBC coverage checks (at least one ViewFactorRayBC and optionally a ReflectRayBC
159 : // on ONLY external boundaries).
160 : std::vector<RayBoundaryConditionBase *> ray_bcs;
161 165 : RayTracingStudy::getRayBCs(ray_bcs, 0);
162 : unsigned int vf_bc_count = 0;
163 360 : for (RayBoundaryConditionBase * rbc : ray_bcs)
164 : {
165 195 : auto view_factor_bc = dynamic_cast<ViewFactorRayBC *>(rbc);
166 195 : if (view_factor_bc)
167 : {
168 165 : ++vf_bc_count;
169 :
170 165 : if (!view_factor_bc->hasBoundary(_bnd_ids))
171 0 : mooseError("The boundary restriction of ",
172 : rbc->type(),
173 : " '",
174 : rbc->name(),
175 : "' does not match 'boundary'");
176 : }
177 : else
178 : {
179 30 : auto reflect_bc = dynamic_cast<ReflectRayBC *>(rbc);
180 30 : if (reflect_bc)
181 : {
182 30 : if (reflect_bc->hasBoundary(_bnd_ids))
183 0 : mooseError("The boundaries applied in ReflectRayBC '",
184 : rbc->name(),
185 : "' cannot include any of the boundaries in ",
186 : type());
187 :
188 71 : for (const BoundaryID internal_bnd_id : getInternalSidesets())
189 41 : if (reflect_bc->hasBoundary(internal_bnd_id))
190 0 : mooseError("The ReflectRayBC '",
191 : rbc->name(),
192 : "' is defined on an internal boundary (",
193 : internal_bnd_id,
194 : ").\n\n",
195 : "This is not allowed for view factor computation.");
196 : }
197 : else
198 0 : mooseError("Does not support the ",
199 : rbc->type(),
200 : " ray boundary condition.\nSupported RayBCs: ReflectRayBC and ViewFactorRayBC.");
201 : }
202 195 : if (vf_bc_count != 1)
203 0 : mooseError("Requires one and only one ViewFactorRayBC.");
204 : }
205 165 : }
206 :
207 : void
208 122 : ViewFactorRayStudy::preExecuteStudy()
209 : {
210 : // Clear and zero the view factor maps we're about to accumulate into for each thread
211 266 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
212 : {
213 144 : _threaded_vf_info[tid].clear();
214 1107 : for (const BoundaryID from_id : _bnd_ids)
215 7680 : for (const BoundaryID to_id : _bnd_ids)
216 6717 : _threaded_vf_info[tid][from_id][to_id] = 0;
217 : }
218 :
219 122 : generateStartElems();
220 122 : }
221 :
222 : void
223 122 : ViewFactorRayStudy::postExecuteStudy()
224 : {
225 : // Finalize the cumulative _vf_info;
226 : _vf_info.clear();
227 942 : for (const BoundaryID from_id : _bnd_ids)
228 6568 : for (const BoundaryID to_id : _bnd_ids)
229 : {
230 5748 : Real & entry = _vf_info[from_id][to_id];
231 :
232 : // Zero before summing
233 5748 : entry = 0;
234 :
235 : // Sum over threads
236 12465 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
237 13434 : entry += _threaded_vf_info[tid][from_id][to_id];
238 :
239 : // Sum over processors
240 5748 : _communicator.sum(entry);
241 : }
242 122 : }
243 :
244 : void
245 122 : ViewFactorRayStudy::generateRays()
246 : {
247 244 : TIME_SECTION("generateRays", 3, "ViewFactorRayStudy Generating Rays");
248 :
249 : // Determine number of Rays and points to allocate space before generation and for output
250 : std::size_t num_local_rays = 0;
251 : std::size_t num_local_start_points = 0;
252 3170 : for (const auto & start_elem : _start_elems)
253 : {
254 3048 : num_local_start_points += start_elem._points.size();
255 3048 : num_local_rays += start_elem._points.size() * _num_dir;
256 : }
257 :
258 : // Print out totals while we're here
259 122 : std::size_t num_total_points = num_local_start_points;
260 122 : std::size_t num_total_rays = num_local_rays;
261 122 : _communicator.sum(num_total_points);
262 122 : _communicator.sum(num_total_rays);
263 122 : _console << "ViewFactorRayStudy generated " << num_total_points
264 122 : << " points with an angular quadrature of " << _num_dir
265 122 : << " directions per point requiring " << num_total_rays << " rays" << std::endl;
266 :
267 : // Reserve space in the buffer ahead of time before we fill it
268 122 : reserveRayBuffer(num_local_rays);
269 :
270 : Point direction;
271 122 : unsigned int num_rays_skipped = 0;
272 :
273 : // loop through all starting points and spawn rays from each for each point and angle
274 3170 : for (const auto & start_elem : _start_elems)
275 : {
276 : // Get normal for the element we're starting on
277 : auto inward_normal =
278 3048 : getSideNormal(start_elem._start_elem, start_elem._incoming_side, /* tid = */ 0);
279 : // We actually want the normal of the original element (remember that we may swap starting
280 : // elements to the element on the other face per requirements of the ray tracer)
281 3048 : if (start_elem._start_elem != start_elem._elem)
282 : inward_normal *= -1;
283 : // Lastly, if the boundary is external and the internal convention is positive, we must
284 : // switch the normal because our AQ uses the inward normal
285 3048 : if (_internal_convention == 0 &&
286 3048 : !start_elem._start_elem->neighbor_ptr(start_elem._incoming_side))
287 : inward_normal *= -1;
288 :
289 : // Rotation for the quadrature to align with the normal; in 3D we can do all of this
290 : // once up front using the 3D aq object. For 2D, we will do it within the direction loop
291 3048 : if (_is_3d)
292 2664 : _3d_aq->rotate(inward_normal);
293 :
294 : // Loop through all points and then all directions
295 59208 : for (std::size_t start_i = 0; start_i < start_elem._points.size(); ++start_i)
296 21684960 : for (std::size_t l = 0; l < _num_dir; ++l)
297 : {
298 : // Get direction of the ray; in 3D we already rotated, in 2D we rotate here
299 21628800 : if (_is_3d)
300 21485952 : direction = _3d_aq->getDirection(l);
301 : else
302 : {
303 142848 : const Real sin_theta = std::sin(_2d_aq_angles[l]);
304 142848 : const Real cos_theta = std::cos(_2d_aq_angles[l]);
305 142848 : direction(0) = cos_theta * inward_normal(0) - sin_theta * inward_normal(1);
306 142848 : direction(1) = sin_theta * inward_normal(0) + cos_theta * inward_normal(1);
307 142848 : direction(2) = 0;
308 : }
309 :
310 : // Angular weight function differs in 2D/3D
311 : // 2D: the quadrature abscissae are the angles between direction & normal.
312 : // The integrand is the cosine of that angle
313 : // 3D: the quadrature abscissae are the azimuthal angle phi and the cosine of the angle
314 : // between normal and direction (= mu). The integrand is mu in that case.
315 21628800 : const auto awf = _is_3d ? inward_normal * direction * _3d_aq->getTotalWeight(l)
316 142848 : : std::cos(_2d_aq_angles[l]) * _2d_aq_weights[l];
317 21628800 : const auto start_weight = start_elem._weights[start_i] * awf;
318 :
319 : // Skip the ray if it exists the domain through the non-planar side it is starting from.
320 : // We do not expect there are any neighbor elements to track it on if it exits the
321 : // non-planar side.
322 : bool intersection_found = false;
323 21485952 : if (_is_3d && start_elem._start_elem &&
324 43114752 : !start_elem._start_elem->neighbor_ptr(start_elem._incoming_side) &&
325 18586368 : sideIsNonPlanar(start_elem._start_elem, start_elem._incoming_side))
326 : {
327 : // Find edge on side that is 'in front' of the future ray
328 : Point intersection_point(std::numeric_limits<Real>::max(), -1, -1);
329 129600 : const auto side_elem = start_elem._start_elem->side_ptr(start_elem._incoming_side);
330 : Point proj_dir;
331 321246 : for (const auto edge_i : side_elem->side_index_range())
332 : {
333 317898 : const auto edge_1 = side_elem->side_ptr(edge_i);
334 : // Project direction onto (start_point, node 1, node 2)
335 : const auto d1 = *edge_1->node_ptr(0) - start_elem._points[start_i];
336 : const auto d2 = *edge_1->node_ptr(1) - start_elem._points[start_i];
337 317898 : const auto d1_unit = d1.unit();
338 317898 : const auto d2_unit = d2.unit();
339 : // If the starting point is aligned with the edge, it wont cross it
340 317898 : if (MooseUtils::absoluteFuzzyEqual(std::abs(d1_unit * d2_unit), 1))
341 0 : continue;
342 317898 : const auto normal = (d1_unit.cross(d2_unit)).unit();
343 :
344 : // One of the nodes must be in front of the start point following the direction
345 317898 : if (d1 * direction < 0 && d2 * direction < 0)
346 65610 : continue;
347 :
348 252288 : proj_dir = (direction - (direction * normal) * normal).unit();
349 :
350 : // Only the side of interest will have the projected direction in between d1 and d2
351 252288 : if ((proj_dir * d2_unit > d1_unit * d2_unit) &&
352 : (proj_dir * d1_unit > d1_unit * d2_unit))
353 : {
354 126252 : const auto dist = geom_utils::distanceFromLine(
355 : start_elem._points[start_i], *edge_1->node_ptr(0), *edge_1->node_ptr(1));
356 : // Ortho-normalize the base on the plane
357 : intersection_point = start_elem._points[start_i] + dist * proj_dir;
358 : intersection_found = true;
359 : break;
360 : }
361 317898 : }
362 :
363 : // Skip the ray if it goes out of the element
364 129600 : const auto grazing_dir = (intersection_point - start_elem._points[start_i]).unit();
365 129600 : if (intersection_found && inward_normal * direction < inward_normal * grazing_dir)
366 : {
367 2196 : num_rays_skipped++;
368 : continue;
369 : }
370 129600 : }
371 :
372 : // Acquire a Ray and fill with the starting information
373 21626604 : std::shared_ptr<Ray> ray = acquireRay();
374 21626604 : ray->setStart(
375 21626604 : start_elem._points[start_i], start_elem._start_elem, start_elem._incoming_side);
376 21626604 : ray->setStartingDirection(direction);
377 21626604 : ray->auxData(_ray_index_start_bnd_id) = start_elem._bnd_id;
378 21626604 : ray->auxData(_ray_index_start_total_weight) = start_weight;
379 :
380 : // Move the Ray into the buffer to be traced
381 21626604 : moveRayToBuffer(ray);
382 : }
383 : }
384 122 : if (num_rays_skipped)
385 3 : mooseInfo(num_rays_skipped,
386 : " rays were skipped as they exited the mesh at their starting point through "
387 : "non-planar sides.");
388 122 : }
389 :
390 : void
391 22971762 : ViewFactorRayStudy::addToViewFactorInfo(Real value,
392 : const BoundaryID from_id,
393 : const BoundaryID to_id,
394 : const THREAD_ID tid)
395 : {
396 : mooseAssert(currentlyPropagating(), "Can only be called during Ray tracing");
397 : mooseAssert(_threaded_vf_info[tid].count(from_id),
398 : "Threaded view factor info does not have from boundary");
399 : mooseAssert(_threaded_vf_info[tid][from_id].count(to_id),
400 : "Threaded view factor info does not have from -> to boundary");
401 :
402 22971762 : _threaded_vf_info[tid][from_id][to_id] += value;
403 22971762 : }
404 :
405 : Real
406 5748 : ViewFactorRayStudy::viewFactorInfo(const BoundaryID from_id, const BoundaryID to_id) const
407 : {
408 : auto it = _vf_info.find(from_id);
409 5748 : if (it == _vf_info.end())
410 0 : mooseError("From boundary id ", from_id, " not in view factor map.");
411 :
412 : auto itt = it->second.find(to_id);
413 5748 : if (itt == it->second.end())
414 0 : mooseError("From boundary id ", from_id, " to boundary_id ", to_id, " not in view factor map.");
415 5748 : return itt->second;
416 : }
417 :
418 : void
419 122 : ViewFactorRayStudy::generateStartElems()
420 : {
421 : const auto & points = _fe_face->get_xyz();
422 : const auto & weights = _fe_face->get_JxW();
423 :
424 : // Clear before filling
425 122 : _start_elems.clear();
426 :
427 : // Starting elements we have that are on the wrong side of an internal boundary
428 : std::unordered_map<processor_id_type, std::vector<StartElem>> send_start_map;
429 :
430 : // Get all possible points on the user defined boundaries on this proc
431 8719 : for (const BndElement * belem : *_mesh.getBoundaryElementRange())
432 : {
433 8597 : const Elem * elem = belem->_elem;
434 8597 : const auto side = belem->_side;
435 8597 : const auto bnd_id = belem->_bnd_id;
436 :
437 : // Skip if we don't own you
438 8597 : if (elem->processor_id() != _pid)
439 5549 : continue;
440 :
441 : // Skip if the boundary id isn't one we're looking for
442 3201 : if (!_bnd_ids.count(bnd_id))
443 3201 : continue;
444 :
445 : // Sanity check on QGRID not working on some types
446 3048 : if (_q_face->type() == libMesh::QGRID && elem->type() == TET4)
447 0 : mooseError(
448 : "Cannot use GRID quadrature type with tetrahedral elements in ViewFactorRayStudy '",
449 0 : _name,
450 : "'");
451 :
452 : // The elem/side that we will actually start the trace from
453 : // (this may change on internal sidesets)
454 3048 : const Elem * start_elem = elem;
455 3048 : auto start_side = side;
456 :
457 : // Reinit this face for points
458 3048 : _fe_face->reinit(elem, side);
459 :
460 : // See if this boundary is internal
461 : const Elem * neighbor = elem->neighbor_ptr(side);
462 3048 : if (neighbor)
463 : {
464 : if (!neighbor->active())
465 0 : mooseError(type(), " does not work with adaptivity");
466 :
467 : // With the positive convention, the Rays that we want to spawn from internal boundaries
468 : // have positive dot products with the outward normal on the side. The ray-tracer requires
469 : // that we provide an element incoming side that is actually incoming (the dot product with
470 : // the direction and the normal is negative). Therefore, switch the physical trace to start
471 : // from the other element and the corresponding side
472 663 : if (_internal_convention == 0)
473 : {
474 663 : start_elem = neighbor;
475 663 : start_side = neighbor->which_neighbor_am_i(elem);
476 : }
477 : }
478 :
479 : // If we own the true starting elem, add to our start info. Otherwise, package the
480 : // start info to be sent to the processor that will actually start this trace
481 3048 : const auto start_pid = start_elem->processor_id();
482 3048 : auto & add_to = _pid ? _start_elems : send_start_map[start_pid];
483 3048 : add_to.emplace_back(elem, start_elem, start_side, bnd_id, points, weights);
484 : }
485 :
486 : // If the internal convention is positive, we may have points that we switched to another
487 : // element for the actual trace, so communicate those to the processors that will
488 : // actually be starting them
489 122 : if (_internal_convention == 0)
490 : {
491 : // Functor that takes in StartElems and appends them to our local list
492 91 : auto append_start_elems = [this](processor_id_type, const std::vector<StartElem> & start_elems)
493 : {
494 91 : _start_elems.reserve(_start_elems.size() + start_elems.size());
495 2463 : for (const StartElem & start_elem : start_elems)
496 2372 : _start_elems.emplace_back(start_elem);
497 213 : };
498 :
499 : // Communicate and act on data
500 122 : Parallel::push_parallel_packed_range(_communicator, send_start_map, this, append_start_elems);
501 : }
502 122 : }
503 :
504 : namespace libMesh
505 : {
506 : namespace Parallel
507 : {
508 :
509 : unsigned int
510 34 : Packing<ViewFactorRayStudy::StartElem>::packing_size(const std::size_t num_points)
511 : {
512 : // Number of points, elem_id, start_elem_id, incoming_side, bnd_id
513 : unsigned int total_size = 5;
514 : // Points
515 34 : total_size += num_points * 3;
516 : // Weights
517 34 : total_size += num_points;
518 :
519 34 : return total_size;
520 : }
521 :
522 : unsigned int
523 17 : Packing<ViewFactorRayStudy::StartElem>::packed_size(typename std::vector<Real>::const_iterator in)
524 : {
525 17 : const std::size_t num_points = *in++;
526 17 : return packing_size(num_points);
527 : }
528 :
529 : unsigned int
530 17 : Packing<ViewFactorRayStudy::StartElem>::packable_size(
531 : const ViewFactorRayStudy::StartElem & start_elem, const void *)
532 : {
533 : mooseAssert(start_elem._points.size() == start_elem._weights.size(), "Size mismatch");
534 17 : return packing_size(start_elem._points.size());
535 : }
536 :
537 : template <>
538 : ViewFactorRayStudy::StartElem
539 17 : Packing<ViewFactorRayStudy::StartElem>::unpack(std::vector<Real>::const_iterator in,
540 : ViewFactorRayStudy * study)
541 : {
542 : // StartElem to fill into
543 : ViewFactorRayStudy::StartElem start_elem;
544 :
545 : // Number of points
546 17 : const std::size_t num_points = static_cast<std::size_t>(*in++);
547 :
548 : // Elem id
549 17 : RayTracingPackingUtils::unpack(start_elem._elem, *in++, &study->meshBase());
550 :
551 : // Start elem id
552 17 : RayTracingPackingUtils::unpack(start_elem._start_elem, *in++, &study->meshBase());
553 :
554 : // Incoming side
555 17 : start_elem._incoming_side = static_cast<unsigned short>(*in++);
556 :
557 : // Boundary ID
558 17 : start_elem._bnd_id = static_cast<BoundaryID>(*in++);
559 :
560 : // Points
561 17 : start_elem._points.resize(num_points);
562 362 : for (std::size_t i = 0; i < num_points; ++i)
563 : {
564 345 : start_elem._points[i](0) = *in++;
565 345 : start_elem._points[i](1) = *in++;
566 345 : start_elem._points[i](2) = *in++;
567 : }
568 :
569 : // Weights
570 17 : start_elem._weights.resize(num_points);
571 362 : for (std::size_t i = 0; i < num_points; ++i)
572 345 : start_elem._weights[i] = *in++;
573 :
574 17 : return start_elem;
575 : }
576 :
577 : template <>
578 : void
579 17 : Packing<ViewFactorRayStudy::StartElem>::pack(const ViewFactorRayStudy::StartElem & start_elem,
580 : std::back_insert_iterator<std::vector<Real>> data_out,
581 : const ViewFactorRayStudy * study)
582 : {
583 : // Number of points
584 17 : data_out = static_cast<buffer_type>(start_elem._points.size());
585 :
586 : // Elem id
587 34 : data_out = RayTracingPackingUtils::pack<buffer_type>(start_elem._elem, &study->meshBase());
588 :
589 : // Start elem id
590 34 : data_out = RayTracingPackingUtils::pack<buffer_type>(start_elem._start_elem, &study->meshBase());
591 :
592 : // Incoming side
593 17 : data_out = static_cast<buffer_type>(start_elem._incoming_side);
594 :
595 : // Boundary id
596 17 : data_out = static_cast<buffer_type>(start_elem._bnd_id);
597 :
598 : // Points
599 362 : for (const auto & point : start_elem._points)
600 : {
601 : data_out = point(0);
602 : data_out = point(1);
603 : data_out = point(2);
604 : }
605 :
606 : // Weights
607 : std::copy(start_elem._weights.begin(), start_elem._weights.end(), data_out);
608 17 : }
609 :
610 : } // namespace Parallel
611 :
612 : } // namespace libMesh
|
