Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // C++ includes
21 :
22 : // Local includes
23 : #include "libmesh/boundary_info.h"
24 : #include "libmesh/cell_c0polyhedron.h"
25 : #include "libmesh/distributed_mesh.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/face_c0polygon.h"
28 : #include "libmesh/int_range.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/parallel_elem.h"
31 : #include "libmesh/parallel_mesh.h"
32 : #include "libmesh/remote_elem.h"
33 :
34 : // Helper functions in anonymous namespace
35 :
36 : namespace
37 : {
38 : using namespace libMesh;
39 :
40 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
41 : static const unsigned int header_size = 12;
42 : #else
43 : static const unsigned int header_size = 11;
44 : #endif
45 :
46 : #ifndef NDEBUG
47 : // Currently this constant is only used for debugging.
48 : static const largest_id_type elem_magic_header = 987654321;
49 : #endif
50 : }
51 :
52 :
53 : namespace libMesh
54 : {
55 :
56 : namespace Parallel
57 : {
58 :
59 : template <>
60 : unsigned int
61 80124442 : Packing<const Elem *>::packed_size (std::vector<largest_id_type>::const_iterator in)
62 : {
63 : #ifndef NDEBUG
64 108932 : const largest_id_type packed_header = *in++;
65 108932 : libmesh_assert_equal_to (packed_header, elem_magic_header);
66 : #endif
67 :
68 : // int 0: level
69 : const unsigned int level =
70 80124442 : cast_int<unsigned int>(*in);
71 :
72 : // int 4: element type
73 80124442 : const int typeint = cast_int<int>(*(in+4));
74 108932 : libmesh_assert_greater_equal (typeint, 0);
75 108932 : libmesh_assert_less (typeint, INVALID_ELEM);
76 : const ElemType type =
77 108932 : cast_int<ElemType>(typeint);
78 :
79 80124442 : unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
80 80124442 : unsigned int n_sides = Elem::type_to_n_sides_map[type];
81 80124442 : unsigned int n_edges = Elem::type_to_n_edges_map[type];
82 108932 : unsigned int variable_topology_size = 0;
83 : // No Elem exists yet, so use the static-count sentinel corresponding
84 : // to Elem::runtime_topology() in type_to_n_nodes_map.
85 108932 : const bool has_runtime_topology = (n_nodes == invalid_uint);
86 :
87 80124442 : if (has_runtime_topology)
88 : {
89 4 : auto topology_in = in + header_size;
90 244 : n_nodes = cast_int<unsigned int>(*topology_in++);
91 244 : n_sides = cast_int<unsigned int>(*topology_in++);
92 244 : n_edges = cast_int<unsigned int>(*topology_in++);
93 4 : variable_topology_size = 3;
94 :
95 244 : if (Elem::type_to_dim_map[type] == 3)
96 : {
97 2 : topology_in += n_nodes;
98 1387 : for (unsigned int s = 0; s != n_sides; ++s)
99 : {
100 : const unsigned int n_side_nodes =
101 1210 : cast_int<unsigned int>(*topology_in++);
102 16 : topology_in += n_side_nodes;
103 1210 : variable_topology_size += n_side_nodes + 1;
104 : }
105 : }
106 : }
107 :
108 80124442 : const unsigned int pre_indexing_size =
109 80124442 : header_size + variable_topology_size + n_nodes + n_sides*2;
110 :
111 : const unsigned int indexing_size =
112 80124442 : DofObject::unpackable_indexing_size(in+pre_indexing_size);
113 :
114 : // We communicate if we are on the boundary or not
115 108932 : unsigned int total_packed_bc_data = 1;
116 80124442 : largest_id_type on_boundary = *(in + pre_indexing_size + indexing_size);
117 :
118 80124442 : if (on_boundary)
119 : {
120 : // Extracting if the children are allowed on the boundary
121 12299 : total_packed_bc_data++;
122 12568486 : largest_id_type allow_children_on_boundary = *(in + pre_indexing_size + indexing_size + 1);
123 :
124 : // For now, children are only supported on sides, the nodes and shell faces are
125 : // treated using the top parents only
126 12568486 : if (level == 0 || allow_children_on_boundary)
127 : {
128 66925652 : for (unsigned int s = 0; s != n_sides; ++s)
129 : {
130 : const int n_bcs = cast_int<int>
131 54413298 : (*(in + pre_indexing_size + indexing_size +
132 54413298 : total_packed_bc_data++));
133 56132 : libmesh_assert_greater_equal (n_bcs, 0);
134 54357166 : total_packed_bc_data += n_bcs;
135 : }
136 : }
137 : }
138 80124442 : if (level == 0)
139 : {
140 210843056 : for (unsigned int e = 0; e != n_edges; ++e)
141 : {
142 : const int n_bcs = cast_int<int>
143 183202101 : (*(in + pre_indexing_size + indexing_size +
144 183202101 : total_packed_bc_data++));
145 617354 : libmesh_assert_greater_equal (n_bcs, 0);
146 182584747 : total_packed_bc_data += n_bcs;
147 : }
148 :
149 84774927 : for (unsigned short sf=0; sf != 2; ++sf)
150 : {
151 : const int n_bcs = cast_int<int>
152 56717778 : (*(in + pre_indexing_size + indexing_size +
153 56717778 : total_packed_bc_data++));
154 201160 : libmesh_assert_greater_equal (n_bcs, 0);
155 56516618 : total_packed_bc_data += n_bcs;
156 : }
157 : }
158 :
159 : return
160 : #ifndef NDEBUG
161 : 1 + // Account for magic header
162 : #endif
163 80124442 : pre_indexing_size + indexing_size + total_packed_bc_data;
164 : }
165 :
166 :
167 :
168 : template <>
169 : unsigned int
170 53956 : Packing<const Elem *>::packed_size (std::vector<largest_id_type>::iterator in)
171 : {
172 53956 : return packed_size(std::vector<largest_id_type>::const_iterator(in));
173 : }
174 :
175 :
176 :
177 : template <>
178 : unsigned int
179 24357242 : Packing<const Elem *>::packable_size (const Elem * const & elem,
180 : const MeshBase * mesh)
181 : {
182 107950 : unsigned int variable_topology_size = 0;
183 24357242 : if (elem->runtime_topology())
184 : {
185 : // Store the dynamic node, side, and edge counts.
186 4 : variable_topology_size = 3;
187 :
188 : // A polygon's node ordering fully specifies its topology. A
189 : // polyhedron additionally needs each side's local node indices.
190 99 : if (elem->dim() == 3)
191 618 : for (auto s : elem->side_index_range())
192 536 : variable_topology_size +=
193 544 : 1 + cast_int<unsigned int>(elem->nodes_on_side(s).size());
194 : }
195 :
196 : // We always communicate if we are on a boundary or not
197 107950 : unsigned int total_packed_bcs = 1;
198 24357242 : const unsigned int n_sides = elem->n_sides();
199 :
200 107950 : largest_id_type on_boundary = 0;
201 108539259 : for (auto s : elem->side_index_range())
202 89463841 : if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s))
203 : {
204 11822 : on_boundary = 1;
205 11822 : break;
206 : }
207 :
208 24357242 : if (on_boundary)
209 : {
210 : // In this case we need another entry to check if we allow children on the boundary.
211 : // We only allow children on sides, edges and sheel faces are treated normally using
212 : // their top parents.
213 11822 : total_packed_bcs++;
214 5335818 : if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side())
215 : {
216 5335818 : total_packed_bcs += n_sides;
217 27994226 : for (unsigned int s = 0; s != n_sides; ++s)
218 22658408 : total_packed_bcs +=
219 22658408 : mesh->get_boundary_info().n_raw_boundary_ids(elem,s);
220 : }
221 : }
222 :
223 24357242 : if (elem->level() == 0)
224 : {
225 12872209 : const unsigned int n_edges = elem->n_edges();
226 12872209 : total_packed_bcs += n_edges;
227 97357661 : for (unsigned int e = 0; e != n_edges; ++e)
228 84485452 : total_packed_bcs +=
229 84485452 : mesh->get_boundary_info().n_edge_boundary_ids(elem,e);
230 :
231 12872209 : total_packed_bcs += 2; // shellfaces
232 38616627 : for (unsigned short sf=0; sf != 2; ++sf)
233 25744418 : total_packed_bcs +=
234 25744418 : mesh->get_boundary_info().n_shellface_boundary_ids(elem,sf);
235 : }
236 :
237 : return
238 : #ifndef NDEBUG
239 : 1 + // add an int for the magic header when testing
240 : #endif
241 24357242 : header_size + variable_topology_size + elem->n_nodes() + n_sides*2 +
242 24357242 : elem->packed_indexing_size() + total_packed_bcs;
243 : }
244 :
245 :
246 :
247 : template <>
248 : unsigned int
249 16345716 : Packing<const Elem *>::packable_size (const Elem * const & elem,
250 : const DistributedMesh * mesh)
251 : {
252 16345716 : return packable_size(elem, static_cast<const MeshBase *>(mesh));
253 : }
254 :
255 :
256 :
257 : template <>
258 : unsigned int
259 0 : Packing<const Elem *>::packable_size (const Elem * const & elem,
260 : const ParallelMesh * mesh)
261 : {
262 0 : return packable_size(elem, static_cast<const MeshBase *>(mesh));
263 : }
264 :
265 :
266 :
267 : template <>
268 : void
269 24302506 : Packing<const Elem *>::pack (const Elem * const & elem,
270 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
271 : const MeshBase * mesh)
272 : {
273 53956 : libmesh_assert(elem);
274 :
275 : #ifndef NDEBUG
276 53956 : *data_out++ = elem_magic_header;
277 : #endif
278 :
279 : #ifdef LIBMESH_ENABLE_AMR
280 24302506 : *data_out++ = (static_cast<largest_id_type>(elem->level()));
281 24302506 : *data_out++ = (static_cast<largest_id_type>(elem->p_level()));
282 :
283 : // Encode both the refinement flag and whether the element has
284 : // children together. This coding is unambiguous because our
285 : // refinement state encoding starts at 0 and ends at
286 : // INVALID_REFINEMENTSTATE
287 : largest_id_type refinement_info =
288 24302506 : static_cast<largest_id_type>(elem->refinement_flag());
289 53956 : if (elem->has_children())
290 3164915 : refinement_info +=
291 : static_cast<largest_id_type>(Elem::INVALID_REFINEMENTSTATE) + 1;
292 53956 : *data_out++ = (refinement_info);
293 :
294 24302506 : *data_out++ = (static_cast<largest_id_type>(elem->p_refinement_flag()));
295 : #else
296 : *data_out++ = (0);
297 : *data_out++ = (0);
298 : *data_out++ = (0);
299 : *data_out++ = (0);
300 : #endif
301 24302506 : *data_out++ = (static_cast<largest_id_type>(elem->type()));
302 24302506 : *data_out++ = (elem->processor_id());
303 24302506 : *data_out++ = (elem->subdomain_id());
304 24302506 : *data_out++ = (elem->id());
305 :
306 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
307 24302506 : if (elem->valid_unique_id())
308 24302506 : *data_out++ = (static_cast<largest_id_type>(elem->unique_id()));
309 : else
310 : // OK to send invalid unique id, we must not own this DOF
311 0 : *data_out++ = (static_cast<largest_id_type>(DofObject::invalid_unique_id));
312 : #endif
313 :
314 : #ifdef LIBMESH_ENABLE_AMR
315 : // use parent_ID of invalid_id to indicate a level 0 element
316 24302506 : if (elem->level() == 0)
317 : {
318 100008 : *data_out++ =(DofObject::invalid_id);
319 100008 : *data_out++ =(DofObject::invalid_id);
320 : }
321 : else
322 : {
323 11485033 : *data_out++ =(elem->parent()->id());
324 11485033 : *data_out++ =(elem->parent()->which_child_am_i(elem));
325 : }
326 : #else
327 : *data_out++ = (DofObject::invalid_id);
328 : *data_out++ = (DofObject::invalid_id);
329 : #endif
330 :
331 31410921 : if ((elem->dim() < LIBMESH_DIM) &&
332 7108415 : elem->interior_parent())
333 4588 : *data_out++ =(elem->interior_parent()->id());
334 : else
335 107856 : *data_out++ =(DofObject::invalid_id);
336 :
337 24302506 : const bool has_variable_topology = elem->runtime_topology();
338 24302506 : if (has_variable_topology)
339 : {
340 97 : *data_out++ = elem->n_nodes();
341 97 : *data_out++ = elem->n_sides();
342 97 : *data_out++ = elem->n_edges();
343 : }
344 :
345 260421249 : for (const Node & node : elem->node_ref_range())
346 236064787 : *data_out++ = node.id();
347 :
348 24302506 : if (has_variable_topology && elem->dim() == 3) // AKA, is c0polyhedron
349 609 : for (auto s : elem->side_index_range())
350 : {
351 : const std::vector<unsigned int> side_nodes =
352 536 : elem->nodes_on_side(s);
353 1048 : *data_out++ = side_nodes.size();
354 2736 : for (const auto node : side_nodes)
355 4344 : *data_out++ = node;
356 : }
357 :
358 : // Add the id of and the side for any return link from each neighbor
359 127896497 : for (auto neigh : elem->neighbor_ptr_range())
360 : {
361 103540035 : if (neigh)
362 : {
363 93366572 : *data_out++ = (neigh->id());
364 93366572 : if (neigh == remote_elem)
365 5970 : *data_out++ = (DofObject::invalid_id);
366 : else
367 85547068 : *data_out++ = neigh->which_neighbor_am_i(elem);
368 : }
369 : else
370 : {
371 364294 : *data_out++ = (DofObject::invalid_id);
372 364294 : *data_out++ = (DofObject::invalid_id);
373 : }
374 : }
375 :
376 : // Add any DofObject indices
377 24302506 : elem->pack_indexing(data_out);
378 :
379 : // We check if this is a boundary cell. We use the raw
380 : // IDs because we also communicate the parents which
381 : // will bring their associated IDs
382 24302506 : largest_id_type on_boundary = 0;
383 108279826 : for (auto s : elem->side_index_range())
384 89252632 : if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s))
385 : {
386 5329268 : on_boundary = 1;
387 5329268 : break;
388 : }
389 :
390 53956 : *data_out++ = on_boundary;
391 :
392 24302506 : if (on_boundary)
393 : {
394 5329268 : *data_out++ = mesh->get_boundary_info().is_children_on_boundary_side();
395 : // Again, only do this if we allow children to hold boundary sides, the edges and
396 : // shell faces are treated normally using their top parents
397 5329268 : if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side())
398 : {
399 11788 : std::vector<boundary_id_type> bcs;
400 27964300 : for (auto s : elem->side_index_range())
401 : {
402 22629138 : mesh->get_boundary_info().raw_boundary_ids(elem, s, bcs);
403 :
404 45231610 : *data_out++ =(bcs.size());
405 :
406 28541672 : for (const auto & bid : bcs)
407 11809684 : *data_out++ = bid;
408 : }
409 : }
410 : }
411 :
412 : // If this is a coarse element,
413 : // Add any element side boundary condition ids
414 24302506 : if (elem->level() == 0)
415 : {
416 100008 : std::vector<boundary_id_type> bcs;
417 97047796 : for (auto e : elem->edge_index_range())
418 : {
419 84176367 : mesh->get_boundary_info().edge_boundary_ids(elem, e, bcs);
420 :
421 168046749 : *data_out++ =(bcs.size());
422 :
423 84177023 : for (const auto & bid : bcs)
424 1216 : *data_out++ = bid;
425 : }
426 :
427 38464275 : for (unsigned short sf=0; sf != 2; ++sf)
428 : {
429 25642850 : mesh->get_boundary_info().shellface_boundary_ids(elem, sf, bcs);
430 :
431 51185692 : *data_out++ =(bcs.size());
432 :
433 25710610 : for (const auto & bid : bcs)
434 128864 : *data_out++ = bid;
435 : }
436 : }
437 24302506 : }
438 :
439 :
440 :
441 : template <>
442 : void
443 16337244 : Packing<const Elem *>::pack (const Elem * const & elem,
444 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
445 : const DistributedMesh * mesh)
446 : {
447 16337244 : pack(elem, data_out, static_cast<const MeshBase*>(mesh));
448 16337244 : }
449 :
450 :
451 :
452 : template <>
453 : void
454 0 : Packing<const Elem *>::pack (const Elem * const & elem,
455 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
456 : const ParallelMesh * mesh)
457 : {
458 0 : pack(elem, data_out, static_cast<const MeshBase*>(mesh));
459 0 : }
460 :
461 :
462 :
463 : // FIXME - this needs serious work to be 64-bit compatible
464 : template <>
465 : Elem *
466 80070486 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
467 : MeshBase * mesh)
468 : {
469 : #ifndef NDEBUG
470 54976 : const std::vector<largest_id_type>::const_iterator original_in = in;
471 :
472 54976 : const largest_id_type incoming_header = *in++;
473 54976 : libmesh_assert_equal_to (incoming_header, elem_magic_header);
474 : #endif
475 :
476 : // int 0: level
477 : const unsigned int level =
478 80070486 : cast_int<unsigned int>(*in++);
479 :
480 : #ifdef LIBMESH_ENABLE_AMR
481 : // int 1: p level
482 : const unsigned int p_level =
483 80070486 : cast_int<unsigned int>(*in++);
484 :
485 : // int 2: refinement flag and encoded has_children
486 80070486 : const int rflag = cast_int<int>(*in++);
487 : const int invalid_rflag =
488 54976 : cast_int<int>(Elem::INVALID_REFINEMENTSTATE);
489 54976 : libmesh_assert_greater_equal (rflag, 0);
490 :
491 54976 : libmesh_assert_less (rflag, invalid_rflag*2+1);
492 :
493 54976 : const bool has_children = (rflag > invalid_rflag);
494 :
495 80070486 : const Elem::RefinementState refinement_flag = has_children ?
496 15270495 : cast_int<Elem::RefinementState>(rflag - invalid_rflag - 1) :
497 53912 : cast_int<Elem::RefinementState>(rflag);
498 :
499 : // int 3: p refinement flag
500 80070486 : const int pflag = cast_int<int>(*in++);
501 54976 : libmesh_assert_greater_equal (pflag, 0);
502 54976 : libmesh_assert_less (pflag, Elem::INVALID_REFINEMENTSTATE);
503 : const Elem::RefinementState p_refinement_flag =
504 54976 : cast_int<Elem::RefinementState>(pflag);
505 : #else
506 : in += 3;
507 : #endif // LIBMESH_ENABLE_AMR
508 :
509 : // int 4: element type
510 80070486 : const int typeint = cast_int<int>(*in++);
511 54976 : libmesh_assert_greater_equal (typeint, 0);
512 54976 : libmesh_assert_less (typeint, INVALID_ELEM);
513 : const ElemType type =
514 54976 : cast_int<ElemType>(typeint);
515 :
516 80070486 : unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
517 80070486 : unsigned int n_sides = Elem::type_to_n_sides_map[type];
518 54976 : unsigned int n_edges = Elem::type_to_n_edges_map[type];
519 : // No Elem exists yet, so use the static-count sentinel corresponding
520 : // to Elem::runtime_topology().
521 54976 : const bool has_runtime_topology = (n_nodes == invalid_uint);
522 :
523 : // int 5: processor id
524 : const processor_id_type processor_id =
525 80070486 : cast_int<processor_id_type>(*in++);
526 54976 : libmesh_assert (processor_id < mesh->n_processors() ||
527 : processor_id == DofObject::invalid_processor_id);
528 :
529 : // int 6: subdomain id
530 : const subdomain_id_type subdomain_id =
531 80070486 : cast_int<subdomain_id_type>(*in++);
532 :
533 : // int 7: dof object id
534 : const dof_id_type id =
535 80070486 : cast_int<dof_id_type>(*in++);
536 54976 : libmesh_assert_not_equal_to (id, DofObject::invalid_id);
537 :
538 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
539 : // int 8: dof object unique id
540 : const unique_id_type unique_id =
541 80070486 : cast_int<unique_id_type>(*in++);
542 : #endif
543 :
544 : #ifdef LIBMESH_ENABLE_AMR
545 : // int 9: parent dof object id.
546 : // Note: If level==0, then (*in) == invalid_id. In
547 : // this case, the equality check in cast_int<unsigned>(*in) will
548 : // never succeed. Therefore, we should only attempt the more
549 : // rigorous cast verification in cases where level != 0.
550 : const dof_id_type parent_id =
551 54976 : (level == 0)
552 80070486 : ? static_cast<dof_id_type>(*in++)
553 51862181 : : cast_int<dof_id_type>(*in++);
554 54976 : libmesh_assert (level == 0 || parent_id != DofObject::invalid_id);
555 54976 : libmesh_assert (level != 0 || parent_id == DofObject::invalid_id);
556 :
557 : // int 10: local child id
558 : // Note: If level==0, then which_child_am_i is not valid, so don't
559 : // do the more rigorous cast verification.
560 : const unsigned int which_child_am_i =
561 54976 : (level == 0)
562 80070486 : ? static_cast<unsigned int>(*in++)
563 51862181 : : cast_int<unsigned int>(*in++);
564 : #else
565 : in += 2;
566 : #endif // LIBMESH_ENABLE_AMR
567 :
568 : const dof_id_type interior_parent_id =
569 80070486 : static_cast<dof_id_type>(*in++);
570 :
571 : // Make sure we don't miscount above when adding the "magic" header
572 : // plus the real data header
573 54976 : libmesh_assert_equal_to (in - original_in, header_size + 1);
574 :
575 80070486 : if (has_runtime_topology)
576 : {
577 242 : n_nodes = cast_int<unsigned int>(*in++);
578 242 : n_sides = cast_int<unsigned int>(*in++);
579 2 : n_edges = cast_int<unsigned int>(*in++);
580 :
581 2 : if (Elem::type_to_dim_map[type] == 2)
582 : {
583 1 : libmesh_assert_less (2, n_sides);
584 1 : libmesh_assert_equal_to (n_nodes % n_sides, 0);
585 1 : libmesh_assert_equal_to (n_edges, n_sides);
586 : }
587 1 : else if (Elem::type_to_dim_map[type] == 3)
588 : {
589 1 : libmesh_assert_less (3, n_nodes);
590 1 : libmesh_assert_less (3, n_sides);
591 : }
592 : }
593 54976 : libmesh_ignore(n_edges); // unused outside dbg/devel
594 :
595 54976 : const auto node_ids_in = in;
596 80015510 : in += n_nodes;
597 :
598 109952 : std::vector<std::vector<unsigned int>> polyhedron_side_nodes;
599 80070486 : if (has_runtime_topology && Elem::type_to_dim_map[type] == 3)
600 : {
601 176 : polyhedron_side_nodes.resize(n_sides);
602 : #ifndef NDEBUG
603 2 : std::vector<bool> node_seen(n_nodes, false);
604 1 : unsigned int next_new_node = 0;
605 : #endif
606 1378 : for (auto & side_nodes : polyhedron_side_nodes)
607 : {
608 : const unsigned int n_side_nodes =
609 1202 : cast_int<unsigned int>(*in++);
610 8 : libmesh_assert_less (2, n_side_nodes);
611 1202 : side_nodes.resize(n_side_nodes);
612 6302 : for (auto & node : side_nodes)
613 : {
614 5136 : node = cast_int<unsigned int>(*in++);
615 36 : libmesh_assert_less (node, n_nodes);
616 :
617 : #ifndef NDEBUG
618 36 : if (type == C0POLYHEDRON && !node_seen[node])
619 : {
620 12 : libmesh_assert_equal_to (node, next_new_node);
621 12 : node_seen[node] = true;
622 12 : ++next_new_node;
623 : }
624 : #endif
625 : }
626 : }
627 :
628 : #ifndef NDEBUG
629 1 : if (type == C0POLYHEDRON)
630 1 : libmesh_assert (next_new_node == n_nodes ||
631 : next_new_node + 1 == n_nodes);
632 : #endif
633 : }
634 :
635 80070486 : Elem * elem = mesh->query_elem_ptr(id);
636 :
637 : // if we already have this element, make sure its
638 : // properties match, and update any missing neighbor
639 : // links, but then go on
640 80070486 : if (elem)
641 : {
642 9690 : libmesh_assert_equal_to (elem->level(), level);
643 9690 : libmesh_assert_equal_to (elem->id(), id);
644 : //#ifdef LIBMESH_ENABLE_UNIQUE_ID
645 : // No check for unique id sanity
646 : //#endif
647 9690 : libmesh_assert_equal_to (elem->processor_id(), processor_id);
648 9690 : libmesh_assert_equal_to (elem->subdomain_id(), subdomain_id);
649 9690 : libmesh_assert_equal_to (elem->type(), type);
650 9690 : libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
651 9690 : libmesh_assert_equal_to (elem->n_sides(), n_sides);
652 9690 : libmesh_assert_equal_to (elem->n_edges(), n_edges);
653 :
654 : #ifndef NDEBUG
655 9690 : if (elem->runtime_topology() && elem->dim() == 3)
656 0 : for (auto s : elem->side_index_range())
657 0 : libmesh_assert(elem->nodes_on_side(s) ==
658 : polyhedron_side_nodes[s]);
659 : #endif
660 :
661 : #ifndef NDEBUG
662 : // All our nodes should be correct
663 105990 : for (unsigned int i=0; i != n_nodes; ++i)
664 96300 : libmesh_assert_equal_to
665 : (elem->node_id(i),
666 : cast_int<dof_id_type>(*(node_ids_in + i)));
667 : #endif
668 :
669 : #ifdef LIBMESH_ENABLE_AMR
670 9690 : libmesh_assert_equal_to (elem->refinement_flag(), refinement_flag);
671 9690 : libmesh_assert_equal_to (elem->has_children(), has_children);
672 :
673 : #ifdef DEBUG
674 9690 : if (elem->active())
675 : {
676 8712 : libmesh_assert_equal_to (elem->p_level(), p_level);
677 8712 : libmesh_assert_equal_to (elem->p_refinement_flag(), p_refinement_flag);
678 : }
679 : #endif
680 :
681 9690 : libmesh_assert (!level || elem->parent() != nullptr);
682 9690 : libmesh_assert (!level || elem->parent()->id() == parent_id);
683 9690 : libmesh_assert (!level || elem->parent()->child_ptr(which_child_am_i) == elem);
684 : #endif
685 : // Our interior_parent link should be "close to" correct - we
686 : // may have to update it, but we can check for some
687 : // inconsistencies.
688 : {
689 : // If the sending processor sees no interior_parent here, we'd
690 : // better agree.
691 49954485 : if (interior_parent_id == DofObject::invalid_id)
692 : {
693 49936788 : if (elem->dim() < LIBMESH_DIM)
694 2623 : libmesh_assert (!(elem->interior_parent()));
695 : }
696 :
697 : // If the sending processor has a remote_elem interior_parent,
698 : // then all we know is that we'd better have *some*
699 : // interior_parent
700 17697 : else if (interior_parent_id == remote_elem->id())
701 : {
702 0 : libmesh_assert(elem->interior_parent());
703 : }
704 : else
705 : {
706 : Elem * ip =
707 17697 : mesh->interior_mesh().query_elem_ptr(interior_parent_id);
708 :
709 : // The sending processor sees an interior parent here, so
710 : // if we don't have that interior element, then we'd
711 : // better have a remote_elem signifying that fact.
712 17697 : if (!ip)
713 0 : libmesh_assert_equal_to (elem->interior_parent(), remote_elem);
714 : else
715 : {
716 : // The sending processor has an interior_parent here,
717 : // and we have that element, but that does *NOT* mean
718 : // we're already linking to it. Perhaps we initially
719 : // received elem from a processor on which the
720 : // interior_parent link was remote?
721 28 : libmesh_assert(elem->interior_parent() == ip ||
722 : elem->interior_parent() == remote_elem);
723 :
724 : // If the link was originally remote, update it
725 17138 : if (elem->interior_parent() == remote_elem)
726 : {
727 0 : elem->set_interior_parent(ip);
728 : }
729 : }
730 : }
731 : }
732 :
733 : // Our neighbor links should be "close to" correct - we may have
734 : // to update a remote_elem link, and we can check for possible
735 : // inconsistencies along the way.
736 : //
737 : // Even for subactive elements, we'll try to keep neighbor links
738 : // in good shape now, if only so any future find_neighbors() is
739 : // idempotent.
740 260843130 : for (auto n : elem->side_index_range())
741 : {
742 : const dof_id_type neighbor_id =
743 210888645 : cast_int<dof_id_type>(*in++);
744 :
745 : const dof_id_type neighbor_side =
746 210888645 : cast_int<dof_id_type>(*in++);
747 :
748 : // If the sending processor sees a domain boundary here,
749 : // we'd better agree ... unless all we see is a remote_elem?
750 : // In that case maybe we just couldn't keep up with a user's
751 : // delete_elem. Let's trust them.
752 210888645 : if (neighbor_id == DofObject::invalid_id)
753 : {
754 20013611 : const Elem * my_neigh = elem->neighbor_ptr(n);
755 20013611 : if (my_neigh == remote_elem)
756 1 : elem->set_neighbor(n, nullptr);
757 : else
758 4883 : libmesh_assert (!my_neigh);
759 20013611 : continue;
760 20003843 : }
761 :
762 : // If the sending processor has a remote_elem neighbor here,
763 : // then all we know is that we'd better *not* have a domain
764 : // boundary ... except that maybe it's the *sending*
765 : // processor who missed a delete_elem we saw.
766 190875034 : if (neighbor_id == remote_elem->id())
767 : {
768 : // At this level of the code we can't even assert in
769 : // cases where the neighbor should know what they're
770 : // talking about, so skip it.
771 :
772 : // libmesh_assert(elem->neighbor_ptr(n));
773 18715306 : continue;
774 : }
775 :
776 172156683 : Elem * neigh = mesh->query_elem_ptr(neighbor_id);
777 :
778 : // The sending processor sees a neighbor here, so if we
779 : // don't have that neighboring element, then we'd better
780 : // have a remote_elem signifying that fact.
781 172156683 : if (!neigh)
782 : {
783 2442 : libmesh_assert_equal_to (elem->neighbor_ptr(n), remote_elem);
784 5785889 : continue;
785 : }
786 :
787 : // The sending processor has a neighbor here, and we have
788 : // that element, but that does *NOT* mean we're already
789 : // linking to it. Perhaps we initially received both elem
790 : // and neigh from processors on which their mutual link was
791 : // remote?
792 37895 : libmesh_assert(elem->neighbor_ptr(n) == neigh ||
793 : elem->neighbor_ptr(n) == remote_elem);
794 :
795 : // If the link was originally remote, we should update it,
796 : // and make sure the appropriate parts of its family link
797 : // back to us.
798 166406247 : if (elem->neighbor_ptr(n) == remote_elem)
799 : {
800 0 : elem->set_neighbor(n, neigh);
801 : }
802 : else
803 37895 : libmesh_assert(elem->subactive() ||
804 : neigh->level() < elem->level() ||
805 : neigh->neighbor_ptr(neighbor_side) == elem);
806 :
807 166368352 : if (neighbor_side != libMesh::invalid_uint)
808 166361849 : elem->make_links_to_me_local(n, neighbor_side);
809 : }
810 :
811 : // Our p level and refinement flags should be "close to" correct
812 : // if we're not an active element - we might have a p level
813 : // increased or decreased by changes in remote_elem children.
814 : //
815 : // But if we have remote_elem children, then we shouldn't be
816 : // doing a projection on this inactive element on this
817 : // processor, so we won't need correct p settings. Couldn't
818 : // hurt to update, though.
819 : #ifdef LIBMESH_ENABLE_AMR
820 49964175 : if (elem->processor_id() != mesh->processor_id())
821 : {
822 : // Do this simultaneously; otherwise we can get a false
823 : // positive when a hack_p_level or set_p_refineemnt_flag
824 : // assertion sees inconsistency between an old flag and new
825 : // value or vice-versa
826 4106 : elem->hack_p_level_and_refinement_flag(p_level, p_refinement_flag);
827 : }
828 : #endif // LIBMESH_ENABLE_AMR
829 :
830 : // FIXME: We should add some debug mode tests to ensure that the
831 : // encoded indexing and boundary conditions are consistent.
832 : }
833 : else
834 : {
835 : // We don't already have the element, so we need to create it.
836 :
837 : // Find the parent if necessary
838 30070715 : Elem * parent = nullptr;
839 : #ifdef LIBMESH_ENABLE_AMR
840 : // Find a child element's parent
841 30116001 : if (level > 0)
842 : {
843 : // Note that we must be very careful to construct the send
844 : // connectivity so that parents are encountered before
845 : // children. If we get here and can't find the parent that
846 : // is a fatal error.
847 25281983 : parent = mesh->elem_ptr(parent_id);
848 : }
849 : // Or assert that the sending processor sees no parent
850 : else
851 44934 : libmesh_assert_equal_to (parent_id, DofObject::invalid_id);
852 : #else
853 : // No non-level-0 elements without AMR
854 : libmesh_assert_equal_to (level, 0);
855 : #endif
856 :
857 30116001 : if (type == C0POLYGON)
858 59 : elem = std::make_unique<C0Polygon>(n_nodes, parent).release();
859 30115942 : else if (type == C0POLYHEDRON)
860 : {
861 108 : std::vector<std::shared_ptr<Polygon>> sides(n_sides);
862 860 : for (auto s : index_range(sides))
863 : {
864 16 : const auto & side_nodes = polyhedron_side_nodes[s];
865 : auto side = std::make_shared<C0Polygon>
866 770 : (cast_int<unsigned int>(side_nodes.size()));
867 4006 : for (auto n : index_range(side_nodes))
868 : {
869 : const dof_id_type node_id =
870 : cast_int<dof_id_type>
871 3324 : (*(node_ids_in + side_nodes[n]));
872 3252 : side->set_node(n, mesh->node_ptr(node_id));
873 : }
874 16 : sides[s] = std::move(side);
875 : }
876 :
877 106 : std::unique_ptr<Node> generated_mid_node;
878 : auto polyhedron = std::make_unique<C0Polyhedron>
879 106 : (sides, generated_mid_node, parent);
880 :
881 1 : libmesh_assert_equal_to (polyhedron->n_nodes(), n_nodes);
882 :
883 106 : if (generated_mid_node)
884 : {
885 : const dof_id_type mid_node_id =
886 59 : cast_int<dof_id_type>(*(node_ids_in + n_nodes - 1));
887 60 : polyhedron->set_node(n_nodes - 1,
888 59 : mesh->node_ptr(mid_node_id));
889 : }
890 :
891 1 : elem = polyhedron.release();
892 104 : }
893 : else
894 30115836 : elem = Elem::build(type,parent).release();
895 45286 : libmesh_assert (elem);
896 :
897 : #ifdef LIBMESH_ENABLE_AMR
898 30116001 : if (level != 0)
899 : {
900 : // Since this is a newly created element, the parent must
901 : // have previously thought of this child as a remote element.
902 352 : libmesh_assert_equal_to (parent->child_ptr(which_child_am_i), remote_elem);
903 :
904 25281983 : parent->add_child(elem, which_child_am_i);
905 : }
906 :
907 : // Assign the refinement flags and levels
908 30116001 : elem->set_p_level(p_level);
909 45286 : elem->set_refinement_flag(refinement_flag);
910 45286 : elem->set_p_refinement_flag(p_refinement_flag);
911 45286 : libmesh_assert_equal_to (elem->level(), level);
912 :
913 : // If this element should have children, assign remote_elem to
914 : // all of them for now, for consistency. Later unpacked
915 : // elements may overwrite that.
916 30116001 : if (has_children)
917 : {
918 6020952 : const unsigned int nc = elem->n_children();
919 30373338 : for (unsigned int c=0; c != nc; ++c)
920 24352386 : elem->add_child(const_cast<RemoteElem *>(remote_elem), c);
921 : }
922 :
923 : #endif // LIBMESH_ENABLE_AMR
924 :
925 : // Assign the IDs
926 30116001 : elem->subdomain_id() = subdomain_id;
927 30116001 : elem->processor_id() = processor_id;
928 30116001 : elem->set_id() = id;
929 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
930 45286 : elem->set_unique_id(unique_id);
931 : #endif
932 :
933 : // Assign the connectivity
934 45286 : libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
935 45286 : libmesh_assert_equal_to (elem->n_sides(), n_sides);
936 45286 : libmesh_assert_equal_to (elem->n_edges(), n_edges);
937 :
938 30116001 : if (!elem->runtime_topology() || elem->dim() != 3)
939 170692281 : for (unsigned int n=0; n != n_nodes; n++)
940 : elem->set_node
941 140815436 : (n, mesh->node_ptr
942 140815436 : (cast_int<dof_id_type>(*(node_ids_in + n))));
943 :
944 : #ifndef NDEBUG
945 284349 : for (unsigned int n = 0; n != n_nodes; ++n)
946 239063 : libmesh_assert_equal_to
947 : (elem->node_id(n),
948 : cast_int<dof_id_type>(*(node_ids_in + n)));
949 : #endif
950 :
951 : // Set interior_parent if found
952 : {
953 : // We may be unpacking an element that was a ghost element on the
954 : // sender, in which case the element's interior_parent may not be
955 : // known by the packed element. We'll have to set such
956 : // interior_parents to remote_elem ourselves and wait for a
957 : // later packed element to give us better information.
958 30116001 : if (interior_parent_id == remote_elem->id())
959 : {
960 : elem->set_interior_parent
961 0 : (const_cast<RemoteElem *>(remote_elem));
962 : }
963 30116001 : else if (interior_parent_id != DofObject::invalid_id)
964 : {
965 : // If we don't have the interior parent element, then it's
966 : // a remote_elem until we get it.
967 : Elem * ip =
968 11209 : mesh->interior_mesh().query_elem_ptr(interior_parent_id);
969 11209 : if (!ip )
970 : elem->set_interior_parent
971 1065 : (const_cast<RemoteElem *>(remote_elem));
972 : else
973 10144 : elem->set_interior_parent(ip);
974 : }
975 : }
976 :
977 151256434 : for (auto n : elem->side_index_range())
978 : {
979 : const dof_id_type neighbor_id =
980 121140433 : cast_int<dof_id_type>(*in++);
981 :
982 : const dof_id_type neighbor_side =
983 121140433 : cast_int<dof_id_type>(*in++);
984 :
985 121140433 : if (neighbor_id == DofObject::invalid_id)
986 14133694 : continue;
987 :
988 : // We may be unpacking an element that was a ghost element on the
989 : // sender, in which case the element's neighbors may not all be
990 : // known by the packed element. We'll have to set such
991 : // neighbors to remote_elem ourselves and wait for a later
992 : // packed element to give us better information.
993 106828192 : if (neighbor_id == remote_elem->id())
994 : {
995 3899205 : elem->set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
996 3899205 : continue;
997 : }
998 :
999 : // If we don't have the neighbor element, then it's a
1000 : // remote_elem until we get it.
1001 102928987 : Elem * neigh = mesh->query_elem_ptr(neighbor_id);
1002 102928987 : if (!neigh)
1003 : {
1004 48650132 : elem->set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
1005 48650132 : continue;
1006 : }
1007 :
1008 : // If we have the neighbor element, then link to it, and
1009 : // make sure any appropriate parts of its family link back
1010 : // to us.
1011 54278855 : elem->set_neighbor(n, neigh);
1012 :
1013 54278855 : if (neighbor_side != libMesh::invalid_uint)
1014 54277703 : elem->make_links_to_me_local(n, neighbor_side);
1015 : }
1016 :
1017 30116001 : elem->unpack_indexing(in);
1018 :
1019 30116001 : mesh->add_elem(elem);
1020 : }
1021 :
1022 80070486 : in += elem->packed_indexing_size();
1023 :
1024 : // We check if this is cell holds a boundary ID or not
1025 80070486 : auto on_boundary = *in++;
1026 80070486 : if (on_boundary)
1027 : {
1028 : // Only treat the sides with caution. This is because we might hold boundary IDs
1029 : // on the sides of the children. This is not supported for edges and shell faces, thus
1030 : // they are treated assuming that only top parents can hold the IDs.
1031 12562592 : auto children_on_boundary = *in++;
1032 12562592 : if (elem->level() == 0 || children_on_boundary)
1033 : {
1034 66893092 : for (auto s : elem->side_index_range())
1035 : {
1036 : const boundary_id_type num_bcs =
1037 54330500 : cast_int<boundary_id_type>(*in++);
1038 :
1039 69470893 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1040 8778 : mesh->get_boundary_info().add_side
1041 15149171 : (elem, s, cast_int<boundary_id_type>(*in++));
1042 : }
1043 : }
1044 : }
1045 :
1046 : // If this is a coarse element,
1047 : // add any element side or edge boundary condition ids
1048 80070486 : if (level == 0)
1049 : {
1050 210487067 : for (auto e : elem->edge_index_range())
1051 : {
1052 : const boundary_id_type num_bcs =
1053 182278762 : cast_int<boundary_id_type>(*in++);
1054 :
1055 182281946 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1056 48 : mesh->get_boundary_info().add_edge
1057 3232 : (elem, e, cast_int<boundary_id_type>(*in++));
1058 : }
1059 :
1060 84624915 : for (unsigned short sf=0; sf != 2; ++sf)
1061 : {
1062 : const boundary_id_type num_bcs =
1063 56416610 : cast_int<boundary_id_type>(*in++);
1064 :
1065 56677518 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1066 3328 : mesh->get_boundary_info().add_shellface
1067 264236 : (elem, sf, cast_int<boundary_id_type>(*in++));
1068 : }
1069 : }
1070 :
1071 : // Return the new element
1072 80125462 : return elem;
1073 79960534 : }
1074 :
1075 :
1076 :
1077 : template <>
1078 : Elem *
1079 16337244 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
1080 : DistributedMesh * mesh)
1081 : {
1082 16337244 : return unpack(in, static_cast<MeshBase*>(mesh));
1083 : }
1084 :
1085 :
1086 :
1087 : template <>
1088 : Elem *
1089 0 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
1090 : ParallelMesh * mesh)
1091 : {
1092 0 : return unpack(in, static_cast<MeshBase*>(mesh));
1093 : }
1094 :
1095 : } // namespace Parallel
1096 :
1097 : } // namespace libMesh
|
