Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 :
9 : #include "RadialGreensConvolution.h"
10 : #include "ThreadedRadialGreensConvolutionLoop.h"
11 : #include "NonlinearSystemBase.h"
12 : #include "FEProblemBase.h"
13 :
14 : // Remove after idaholab/moose#26102
15 : #include "MagpieNanoflann.h"
16 :
17 : #include "libmesh/nanoflann.hpp"
18 : #include "libmesh/parallel_algebra.h"
19 : #include "libmesh/mesh_tools.h"
20 :
21 : #include <list>
22 : #include <iterator>
23 : #include <algorithm>
24 :
25 : registerMooseObject("MagpieApp", RadialGreensConvolution);
26 :
27 : // specialization for PointListAdaptor<RadialGreensConvolution::QPData>
28 : template <>
29 : const Point &
30 258769660 : PointListAdaptor<RadialGreensConvolution::QPData>::getPoint(
31 : const RadialGreensConvolution::QPData & item) const
32 : {
33 258769660 : return item._q_point;
34 : }
35 :
36 : InputParameters
37 54 : RadialGreensConvolution::validParams()
38 : {
39 54 : InputParameters params = ElementUserObject::validParams();
40 54 : params.addClassDescription("Perform a radial Green's function convolution");
41 108 : params.addCoupledVar("v", "Variable to gather");
42 108 : params.addRequiredParam<FunctionName>(
43 : "function",
44 : "Green's function (distance is substituted for x) without geometrical attenuation");
45 108 : params.addRequiredParam<Real>("r_cut", "Cut-off radius for the Green's function");
46 108 : params.addParam<bool>("normalize", false, "Normalize the Green's function integral to one");
47 :
48 : // we run this object once at the beginning of the timestep by default
49 108 : params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
50 :
51 : // make sure we always have geometric point neighbors ghosted
52 108 : params.addRelationshipManager("ElementPointNeighborLayers",
53 : Moose::RelationshipManagerType::GEOMETRIC |
54 : Moose::RelationshipManagerType::ALGEBRAIC);
55 :
56 54 : return params;
57 0 : }
58 :
59 30 : RadialGreensConvolution::RadialGreensConvolution(const InputParameters & parameters)
60 : : ElementUserObject(parameters),
61 30 : _v(coupledValue("v")),
62 30 : _v_var(coupled("v")),
63 30 : _function(getFunction("function")),
64 60 : _r_cut(getParam<Real>("r_cut")),
65 60 : _normalize(getParam<bool>("normalize")),
66 30 : _dim(_mesh.dimension()),
67 30 : _correction_integral(100),
68 30 : _dof_map(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).dofMap()),
69 30 : _update_communication_lists(false),
70 30 : _my_pid(processor_id()),
71 30 : _perf_meshchanged(registerTimedSection("meshChanged", 3)),
72 30 : _perf_updatelists(registerTimedSection("updateCommunicationLists", 3)),
73 90 : _perf_finalize(registerTimedSection("finalize", 2))
74 : {
75 : // collect mesh periodicity data
76 90 : for (unsigned int i = 0; i < _dim; ++i)
77 : {
78 60 : _periodic[i] = _mesh.isRegularOrthogonal() && _mesh.isTranslatedPeriodic(_v_var, i);
79 :
80 60 : _periodic_min[i] = _mesh.getMinInDimension(i);
81 60 : _periodic_max[i] = _mesh.getMaxInDimension(i);
82 60 : _periodic_vector[i](i) = _mesh.dimensionWidth(i);
83 :
84 : // we could allow this, but then we'd have to search over more than just the nearest periodic
85 : // neighbors
86 60 : if (_periodic[i] && 2.0 * _r_cut > _periodic_vector[i](i))
87 0 : paramError(
88 : "r_cut",
89 : "The cut-off radius cannot be larger than half the periodic size of the simulation cell");
90 : }
91 30 : }
92 :
93 : void
94 30 : RadialGreensConvolution::initialSetup()
95 : {
96 : // Get a pointer to the PeriodicBoundaries buried in libMesh
97 30 : _pbs = _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).dofMap().get_periodic_boundaries();
98 :
99 : // set up processor boundary node list
100 30 : meshChanged();
101 30 : }
102 :
103 : void
104 30 : RadialGreensConvolution::initialize()
105 : {
106 : _qp_data.clear();
107 :
108 : /* compute the dimensional correction
109 : *
110 : * |---__
111 : * | / .
112 : * | R/ | .
113 : * | / h| \ r_cut
114 : * |______|__|____
115 : * r
116 : *
117 : * Out of plane (2D) and out of line (1D) 3D contributions to the convolution
118 : * are not naturally captured by the 1D or 2D numerical integration of QP neighborhood.
119 : * At every Qp an integration in the h direction needs to be performed.
120 : * We can pre-tabulate this integral.
121 : */
122 30 : if (_dim < 3)
123 : {
124 25 : const Real dr = _r_cut / _correction_integral.size();
125 2525 : for (unsigned int i = 0; i < _correction_integral.size(); ++i)
126 : {
127 2500 : const Real r = i * dr;
128 2500 : const Real h = std::sqrt(_r_cut * _r_cut - r * r);
129 2500 : const unsigned int hsteps = std::ceil(h / dr);
130 2500 : const Real dh = h / hsteps;
131 :
132 2500 : _correction_integral[i] = 0.0;
133 :
134 : // for r = 0 we skip the first bin of the integral!
135 2500 : if (i == 0)
136 25 : _zero_dh = dh;
137 :
138 201205 : for (unsigned int j = i > 0 ? 0 : 1; j < hsteps; ++j)
139 : {
140 198705 : const Real h1 = j * dh;
141 198705 : const Real h2 = h1 + dh;
142 :
143 : // we evaluate the Green's function at the geometric middle of the height interval
144 : // and assume it to be constant.
145 198705 : const Real R2 = r * r + h1 * h2;
146 198705 : const Real G = _function.value(_t, Point(std::sqrt(R2), 0.0, 0.0));
147 :
148 : // We integrate the geometric attenuation analytically over the interval [h1, h2).
149 198705 : _correction_integral[i] += G * attenuationIntegral(h1, h2, r, _dim);
150 : }
151 : }
152 : }
153 30 : }
154 :
155 : Real
156 329625 : RadialGreensConvolution::attenuationIntegral(Real h1, Real h2, Real r, unsigned int dim) const
157 : {
158 329625 : if (dim == 2)
159 : {
160 : // in 2D we need to return the above and below plane contribution (hence 2.0 * ...)
161 289765 : if (r > 0.0)
162 : // integral 1/(h^2+r^2) dh = 1/r*atan(h/r)|
163 156985 : return 2.0 * 1.0 / (4.0 * libMesh::pi * r) * (std::atan(h2 / r) - std::atan(h1 / r));
164 : else
165 : // integral 1/h^2 dh = -1/h|
166 132780 : return 2.0 * 1.0 / (4.0 * libMesh::pi) * (-1.0 / h2 + 1.0 / h1);
167 : }
168 : else
169 : // dim = 1, we multiply the attenuation by 2pi*h (1/2 * 2pi/4pi = 0.25)
170 : // integral h/(h^2+r^2) dh = 1/2 * ln(h^2+r^2)|
171 39860 : return 0.25 * (std::log(h2 * h2 + r * r) - std::log(h1 * h1 + r * r));
172 : }
173 :
174 : void
175 34260 : RadialGreensConvolution::execute()
176 : {
177 34260 : auto id = _current_elem->id();
178 :
179 : // collect all QP data
180 177180 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
181 142920 : _qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _v[qp]);
182 :
183 : // make sure the result map entry for the current element is sized correctly
184 : auto i = _convolution.find(id);
185 34260 : if (i == _convolution.end())
186 68520 : _convolution.insert(std::make_pair(id, std::vector<Real>(_qrule->n_points())));
187 : else
188 0 : i->second.resize(_qrule->n_points());
189 34260 : }
190 :
191 : void
192 24 : RadialGreensConvolution::finalize()
193 : {
194 24 : TIME_SECTION(_perf_finalize);
195 :
196 : // the first chunk of data is always the local data - remember its size
197 24 : unsigned int local_size = _qp_data.size();
198 :
199 : // if normalization is requested we need to integrate the current variable field
200 24 : Real source_integral = 0.0;
201 24 : if (_normalize)
202 : {
203 142944 : for (const auto & qpd : _qp_data)
204 142920 : source_integral += qpd._volume * qpd._value;
205 : gatherSum(source_integral);
206 : }
207 :
208 : // communicate the qp data list if n_proc > 1
209 24 : if (_app.n_processors() > 1)
210 : {
211 : // !!!!!!!!!!!
212 : // !!CAREFUL!! Is it guaranteed that _qp_data is in the same order if the mesh has not changed?
213 : // According to @friedmud it is not guaranteed if threads are used
214 : // !!!!!!!!!!!
215 :
216 : // update after mesh changes and/or if a displaced problem exists
217 12 : if (_update_communication_lists || _fe_problem.getDisplacedProblem() ||
218 : libMesh::n_threads() > 1)
219 12 : updateCommunicationLists();
220 :
221 : // sparse send data (processor ID,)
222 : std::vector<std::size_t> non_zero_comm;
223 36 : for (auto i = beginIndex(_communication_lists); i < _communication_lists.size(); ++i)
224 24 : if (!_communication_lists[i].empty())
225 12 : non_zero_comm.push_back(i);
226 :
227 : // data structures for sparse point to point communication
228 12 : std::vector<std::vector<QPData>> send(non_zero_comm.size());
229 12 : std::vector<Parallel::Request> send_requests(non_zero_comm.size());
230 12 : Parallel::MessageTag send_tag = _communicator.get_unique_tag(4711);
231 : std::vector<QPData> receive;
232 :
233 12 : const auto item_type = StandardType<QPData>(&(_qp_data[0]));
234 :
235 : #if 0
236 : // output local qp locations
237 : // _console << name() << ' ' << receive.size() << '\n' << name() << std::flush;
238 : for (auto item : _qp_data)
239 : _console << name() << ' ' << _my_pid << ' '
240 : << item._q_point(0) << ' '
241 : << item._q_point(1) << ' '
242 : << item._q_point(2) << std::flush;
243 : #endif
244 :
245 : // fill buffer and send structures
246 24 : for (auto i = beginIndex(non_zero_comm); i < non_zero_comm.size(); ++i)
247 : {
248 12 : const auto pid = non_zero_comm[i];
249 : const auto & list = _communication_lists[pid];
250 :
251 : // fill send buffer for transfer to pid
252 12 : send[i].reserve(list.size());
253 15204 : for (const auto & item : list)
254 : {
255 15192 : send[i].push_back(_qp_data[item]);
256 : #if 0
257 : // output sent qp locations
258 : _console << name() << ' '
259 : << _qp_data[item]._q_point(0) << ' '
260 : << _qp_data[item]._q_point(1) << ' '
261 : << _qp_data[item]._q_point(2) << ' '
262 : << pid << std::flush;
263 : #endif
264 : }
265 :
266 : // issue non-blocking send
267 12 : _communicator.send(pid, send[i], send_requests[i], send_tag);
268 : }
269 :
270 : // receive messages - we assume that we receive as many messages as we send!
271 24 : for (auto i = beginIndex(non_zero_comm); i < non_zero_comm.size(); ++i)
272 : {
273 : // inspect incoming message
274 12 : Parallel::Status status(_communicator.probe(Parallel::any_source, send_tag));
275 : const auto source_pid = TIMPI::cast_int<processor_id_type>(status.source());
276 : const auto message_size = status.size(item_type);
277 :
278 : // resize receive buffer accordingly and receive data
279 12 : receive.resize(message_size);
280 12 : _communicator.receive(source_pid, receive, send_tag);
281 :
282 : #if 0
283 : // output received qp locations
284 : // _console << name() << ' ' << receive.size() << '\n' << name() << std::flush;
285 : for (auto item : receive)
286 : _console << name() << ' ' << source_pid << ' '
287 : << item._q_point(0) << ' '
288 : << item._q_point(1) << ' '
289 : << item._q_point(2) << std::flush;
290 : #endif
291 :
292 : // append communicated data
293 12 : _qp_data.insert(_qp_data.end(), receive.begin(), receive.end());
294 : }
295 :
296 : // wait until all send requests are at least buffered and we can destroy
297 : // the send buffers by going out of scope
298 12 : Parallel::wait(send_requests);
299 12 : }
300 :
301 : // build KD-Tree using data we just received
302 : const unsigned int max_leaf_size = 20; // slightly affects runtime
303 : auto point_list = PointListAdaptor<QPData>(_qp_data.begin(), _qp_data.end());
304 24 : _kd_tree = std::make_unique<KDTreeType>(
305 48 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
306 :
307 : mooseAssert(_kd_tree != nullptr, "KDTree was not properly initialized.");
308 24 : _kd_tree->buildIndex();
309 :
310 : // result map entry
311 : const auto end_it = _convolution.end();
312 : auto it = end_it;
313 :
314 : // build thread loop functor
315 24 : ThreadedRadialGreensConvolutionLoop rgcl(*this);
316 :
317 : // run threads
318 : auto local_range_begin = _qp_data.begin();
319 24 : auto local_range_end = local_range_begin;
320 : std::advance(local_range_end, local_size);
321 48 : Threads::parallel_reduce(QPDataRange(local_range_begin, local_range_end), rgcl);
322 :
323 24 : Real convolution_integral = rgcl.convolutionIntegral();
324 :
325 : // if normalization is requested we need to communicate the convolution result
326 : // and normalize the result entries
327 24 : if (_normalize)
328 : {
329 : gatherSum(convolution_integral);
330 :
331 : // we may not need to or may not be able to normalize
332 24 : if (convolution_integral == 0.0)
333 : {
334 0 : if (source_integral == 0.0)
335 : return;
336 0 : mooseError("Unable to normalize Green's function. Is it all zero?");
337 : }
338 :
339 : // normalize result entries
340 34284 : for (auto & re : _convolution)
341 177180 : for (auto & ri : re.second)
342 142920 : ri *= source_integral / convolution_integral;
343 : }
344 :
345 : // make it a differential result
346 : it = end_it;
347 142944 : for (unsigned int i = 0; i < local_size; ++i)
348 : {
349 142920 : const auto & local_qp = _qp_data[i];
350 :
351 : // Look up result map iterator only if we enter a new element. this saves a bunch
352 : // of map lookups because same element entries are consecutive in the _qp_data vector.
353 142920 : if (it == end_it || it->first != local_qp._elem_id)
354 34260 : it = _convolution.find(local_qp._elem_id);
355 :
356 142920 : it->second[local_qp._qp] -= local_qp._volume * local_qp._value;
357 : }
358 : }
359 :
360 : void
361 6 : RadialGreensConvolution::threadJoin(const UserObject & y)
362 : {
363 : const RadialGreensConvolution & uo = static_cast<const RadialGreensConvolution &>(y);
364 6 : _qp_data.insert(_qp_data.begin(), uo._qp_data.begin(), uo._qp_data.end());
365 6 : _convolution.insert(uo._convolution.begin(), uo._convolution.end());
366 6 : }
367 :
368 : void
369 822 : RadialGreensConvolution::insertNotLocalPointNeighbors(dof_id_type node)
370 : {
371 : mooseAssert(!_nodes_to_elem_map[node].empty(), "Node not found in _nodes_to_elem_map");
372 :
373 4570 : for (const auto * elem : _nodes_to_elem_map[node])
374 3748 : if (elem->processor_id() != _my_pid && elem->active())
375 1874 : _point_neighbors.insert(elem);
376 822 : }
377 :
378 : void
379 20280 : RadialGreensConvolution::insertNotLocalPeriodicPointNeighbors(dof_id_type node,
380 : const Node * reference)
381 : {
382 : mooseAssert(!_nodes_to_elem_map[node].empty(), "Node not found in _nodes_to_elem_map");
383 :
384 20280 : const Node * first = _mesh.nodePtr(node);
385 70080 : for (const auto * elem : _nodes_to_elem_map[node])
386 49800 : if (elem->processor_id() != _my_pid && elem->active())
387 2731 : _periodic_point_neighbors.emplace(elem, first, reference);
388 20280 : }
389 :
390 : void
391 12088 : RadialGreensConvolution::findNotLocalPeriodicPointNeighbors(const Node * first)
392 : {
393 12088 : auto iters = _periodic_node_map.equal_range(first->id());
394 : auto it = iters.first;
395 :
396 : // no periodic copies
397 12088 : if (it == iters.second)
398 10136 : return;
399 :
400 : // insert first periodic neighbor
401 12088 : insertNotLocalPeriodicPointNeighbors(it->second, first);
402 : ++it;
403 :
404 : // usual case, node was on the face of a periodic boundary (and not on its edge or corner)
405 12088 : if (it == iters.second)
406 : return;
407 :
408 : // number of periodic directions
409 : unsigned int periodic_dirs = 1;
410 : std::set<dof_id_type> nodes;
411 1952 : nodes.insert(iters.first->second); // probably not necessary
412 :
413 : // insert remaining periodic copies
414 : do
415 : {
416 2048 : nodes.insert(it->second);
417 : it++;
418 2048 : periodic_dirs++;
419 2048 : } while (it != iters.second);
420 :
421 : // now jump periodic_dirs-1 times from those nodes to their periodic copies
422 4000 : for (unsigned int i = 1; i < periodic_dirs; ++i)
423 : {
424 : std::set<dof_id_type> new_nodes;
425 6720 : for (auto node : nodes)
426 : {
427 : // periodic copies of the set members we already inserted
428 : auto new_iters = _periodic_node_map.equal_range(node);
429 :
430 : // insert the ids of the periodic copies of the node
431 14976 : for (it = new_iters.first; it != new_iters.second; ++it)
432 10304 : new_nodes.insert(it->second);
433 : }
434 2048 : nodes.insert(new_nodes.begin(), new_nodes.end());
435 : }
436 :
437 : // add all jumped to nodes
438 10144 : for (auto node : nodes)
439 8192 : insertNotLocalPeriodicPointNeighbors(node, first);
440 : }
441 :
442 : void
443 30 : RadialGreensConvolution::meshChanged()
444 : {
445 30 : TIME_SECTION(_perf_meshchanged);
446 :
447 : // get underlying libMesh mesh
448 30 : auto & mesh = _mesh.getMesh();
449 :
450 : // get a fresh point locator
451 60 : _point_locator = _mesh.getPointLocator();
452 :
453 : // rebuild periodic node map (this is the heaviest part by far)
454 30 : _mesh.buildPeriodicNodeMap(_periodic_node_map, _v_var, _pbs);
455 :
456 : // Build a new node to element map
457 : _nodes_to_elem_map.clear();
458 30 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), _nodes_to_elem_map);
459 :
460 : // clear point neighbor set
461 : _point_neighbors.clear();
462 :
463 : // my processor id
464 : const auto my_pid = processor_id();
465 :
466 : // iterate over active local elements
467 30 : const auto end = mesh.active_local_elements_end();
468 91420 : for (auto it = mesh.active_local_elements_begin(); it != end; ++it)
469 : // find a face that faces either a boundary (nullptr) or a different processor
470 232240 : for (unsigned int s = 0; s < (*it)->n_sides(); ++s)
471 : {
472 186560 : const auto * neighbor = (*it)->neighbor_ptr(s);
473 186560 : if (neighbor)
474 : {
475 182312 : if (neighbor->processor_id() != my_pid)
476 : {
477 : // add all face node touching elements directly
478 3650 : for (unsigned int n = 0; n < (*it)->n_nodes(); ++n)
479 1644 : if ((*it)->is_node_on_side(n, s))
480 822 : insertNotLocalPointNeighbors((*it)->node_id(n));
481 : }
482 : }
483 : else
484 : {
485 : // find periodic node neighbors and
486 52600 : for (unsigned int n = 0; n < (*it)->n_nodes(); ++n)
487 24176 : if ((*it)->is_node_on_side(n, s))
488 12088 : findNotLocalPeriodicPointNeighbors((*it)->node_ptr(n));
489 : }
490 : }
491 :
492 : // request communication list update
493 30 : _update_communication_lists = true;
494 30 : }
495 :
496 : void
497 12 : RadialGreensConvolution::updateCommunicationLists()
498 : {
499 12 : TIME_SECTION(_perf_updatelists);
500 :
501 : // clear communication lists
502 12 : _communication_lists.clear();
503 12 : _communication_lists.resize(n_processors());
504 :
505 : // build KD-Tree using local qpoint data
506 : const unsigned int max_leaf_size = 20; // slightly affects runtime
507 : auto point_list = PointListAdaptor<QPData>(_qp_data.begin(), _qp_data.end());
508 : auto kd_tree = std::make_unique<KDTreeType>(
509 12 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
510 : mooseAssert(kd_tree != nullptr, "KDTree was not properly initialized.");
511 12 : kd_tree->buildIndex();
512 :
513 : std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
514 : nanoflann::SearchParameters search_params;
515 :
516 : // iterate over non periodic point neighbor elements
517 374 : for (auto elem : _point_neighbors)
518 : {
519 : ret_matches.clear();
520 362 : Point centroid = elem->vertex_average();
521 362 : const Real r_cut2 = _r_cut + elem->hmax() / 2.0;
522 362 : kd_tree->radiusSearch(&(centroid(0)), r_cut2 * r_cut2, ret_matches, search_params);
523 89036 : for (auto & match : ret_matches)
524 88674 : _communication_lists[elem->processor_id()].insert(match.first);
525 : }
526 :
527 : // iterate over periodic point neighbor elements
528 1014 : for (auto tuple : _periodic_point_neighbors)
529 : {
530 : const auto * elem = std::get<0>(tuple);
531 : const auto * first = std::get<1>(tuple);
532 : const auto * second = std::get<2>(tuple);
533 :
534 1002 : Point centroid = elem->vertex_average() - (*first - *second);
535 :
536 1002 : const Real r_cut2 = _r_cut + elem->hmax() / 2.0;
537 :
538 : ret_matches.clear();
539 1002 : kd_tree->radiusSearch(&(centroid(0)), r_cut2 * r_cut2, ret_matches, search_params);
540 259656 : for (auto & match : ret_matches)
541 258654 : _communication_lists[elem->processor_id()].insert(match.first);
542 : }
543 :
544 : // request fulfilled
545 12 : _update_communication_lists = false;
546 24 : }
547 :
548 : namespace TIMPI
549 : {
550 :
551 36 : StandardType<RadialGreensConvolution::QPData>::StandardType(
552 36 : const RadialGreensConvolution::QPData * example)
553 : {
554 : // We need an example for MPI_Address to use
555 36 : static const RadialGreensConvolution::QPData p;
556 36 : if (!example)
557 : example = &p;
558 :
559 : #ifdef LIBMESH_HAVE_MPI
560 :
561 : // Get the sub-data-types, and make sure they live long enough
562 : // to construct the derived type
563 36 : StandardType<Point> d1(&example->_q_point);
564 36 : StandardType<dof_id_type> d2(&example->_elem_id);
565 36 : StandardType<short> d3(&example->_qp);
566 36 : StandardType<Real> d4(&example->_volume);
567 36 : StandardType<Real> d5(&example->_value);
568 :
569 : MPI_Datatype types[] = {
570 36 : (data_type)d1, (data_type)d2, (data_type)d3, (data_type)d4, (data_type)d5};
571 36 : int blocklengths[] = {1, 1, 1, 1, 1};
572 : MPI_Aint displs[5], start;
573 :
574 36 : libmesh_call_mpi(MPI_Get_address(const_cast<RadialGreensConvolution::QPData *>(example), &start));
575 36 : libmesh_call_mpi(MPI_Get_address(const_cast<Point *>(&example->_q_point), &displs[0]));
576 36 : libmesh_call_mpi(MPI_Get_address(const_cast<dof_id_type *>(&example->_elem_id), &displs[1]));
577 36 : libmesh_call_mpi(MPI_Get_address(const_cast<short *>(&example->_qp), &displs[2]));
578 36 : libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_volume), &displs[3]));
579 36 : libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_value), &displs[4]));
580 :
581 216 : for (std::size_t i = 0; i < 5; ++i)
582 180 : displs[i] -= start;
583 :
584 : // create a prototype structure
585 : MPI_Datatype tmptype;
586 36 : libmesh_call_mpi(MPI_Type_create_struct(5, blocklengths, displs, types, &tmptype));
587 36 : libmesh_call_mpi(MPI_Type_commit(&tmptype));
588 :
589 : // resize the structure type to account for padding, if any
590 36 : libmesh_call_mpi(
591 : MPI_Type_create_resized(tmptype, 0, sizeof(RadialGreensConvolution::QPData), &_datatype));
592 36 : libmesh_call_mpi(MPI_Type_free(&tmptype));
593 :
594 : this->commit();
595 :
596 : #endif // LIBMESH_HAVE_MPI
597 36 : }
598 : } // namespace TIMPI
599 :
600 0 : StandardType<RadialGreensConvolution::QPData>::StandardType(
601 0 : const StandardType<RadialGreensConvolution::QPData> & t)
602 0 : : TIMPI::DataType()
603 : {
604 : #ifdef LIBMESH_HAVE_MPI
605 0 : libmesh_call_mpi(MPI_Type_dup(t._datatype, &_datatype));
606 : #endif
607 0 : }
|
