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/enum_to_string.h"
28 : #include "libmesh/face_c0polygon.h"
29 : #include "libmesh/int_range.h"
30 : #include "libmesh/mesh_base.h"
31 : #include "libmesh/parallel_elem.h"
32 : #include "libmesh/parallel_mesh.h"
33 : #include "libmesh/remote_elem.h"
34 :
35 : // Helper functions in anonymous namespace
36 :
37 : namespace
38 : {
39 : using namespace libMesh;
40 :
41 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
42 : static const unsigned int header_size = 12;
43 : #else
44 : static const unsigned int header_size = 11;
45 : #endif
46 :
47 : #ifndef NDEBUG
48 : // Currently this constant is only used for debugging.
49 : static const largest_id_type elem_magic_header = 987654321;
50 : #endif
51 : }
52 :
53 :
54 : namespace libMesh
55 : {
56 :
57 : namespace Parallel
58 : {
59 :
60 : template <>
61 : unsigned int
62 219293 : Packing<const Elem *>::packed_size (std::vector<largest_id_type>::const_iterator in)
63 : {
64 : #ifndef NDEBUG
65 108932 : const largest_id_type packed_header = *in++;
66 108932 : libmesh_assert_equal_to (packed_header, elem_magic_header);
67 : #endif
68 :
69 : // int 0: level
70 : const unsigned int level =
71 219293 : cast_int<unsigned int>(*in);
72 :
73 : // int 4: element type
74 219293 : const int typeint = cast_int<int>(*(in+4));
75 108932 : libmesh_assert_greater_equal (typeint, 0);
76 108932 : libmesh_assert_less (typeint, INVALID_ELEM);
77 : const ElemType type =
78 108932 : cast_int<ElemType>(typeint);
79 :
80 219293 : unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
81 219293 : unsigned int n_sides = Elem::type_to_n_sides_map[type];
82 219293 : unsigned int n_edges = Elem::type_to_n_edges_map[type];
83 108932 : unsigned int variable_topology_size = 0;
84 :
85 219293 : if (n_nodes == invalid_uint)
86 : {
87 8 : libmesh_error_msg_if(type != C0POLYGON && type != C0POLYHEDRON,
88 : "Unsupported variable-topology element "
89 : << Utility::enum_to_string(type));
90 :
91 4 : auto topology_in = in + header_size;
92 8 : n_nodes = cast_int<unsigned int>(*topology_in++);
93 8 : n_sides = cast_int<unsigned int>(*topology_in++);
94 8 : n_edges = cast_int<unsigned int>(*topology_in++);
95 4 : variable_topology_size = 3;
96 :
97 8 : if (type == C0POLYHEDRON)
98 : {
99 2 : topology_in += n_nodes;
100 36 : for (unsigned int s = 0; s != n_sides; ++s)
101 : {
102 : const unsigned int n_side_nodes =
103 32 : cast_int<unsigned int>(*topology_in++);
104 16 : topology_in += n_side_nodes;
105 32 : variable_topology_size += n_side_nodes + 1;
106 : }
107 : }
108 : }
109 :
110 219293 : const unsigned int pre_indexing_size =
111 219293 : header_size + variable_topology_size + n_nodes + n_sides*2;
112 :
113 : const unsigned int indexing_size =
114 219293 : DofObject::unpackable_indexing_size(in+pre_indexing_size);
115 :
116 : // We communicate if we are on the boundary or not
117 108932 : unsigned int total_packed_bc_data = 1;
118 219293 : largest_id_type on_boundary = *(in + pre_indexing_size + indexing_size);
119 :
120 219293 : if (on_boundary)
121 : {
122 : // Extracting if the children are allowed on the boundary
123 12299 : total_packed_bc_data++;
124 25123 : largest_id_type allow_children_on_boundary = *(in + pre_indexing_size + indexing_size + 1);
125 :
126 : // For now, children are only supported on sides, the nodes and shell faces are
127 : // treated using the top parents only
128 25123 : if (level == 0 || allow_children_on_boundary)
129 : {
130 140295 : for (unsigned int s = 0; s != n_sides; ++s)
131 : {
132 : const int n_bcs = cast_int<int>
133 171304 : (*(in + pre_indexing_size + indexing_size +
134 171304 : total_packed_bc_data++));
135 56132 : libmesh_assert_greater_equal (n_bcs, 0);
136 115172 : total_packed_bc_data += n_bcs;
137 : }
138 : }
139 : }
140 219293 : if (level == 0)
141 : {
142 1442783 : for (unsigned int e = 0; e != n_edges; ++e)
143 : {
144 : const int n_bcs = cast_int<int>
145 1857996 : (*(in + pre_indexing_size + indexing_size +
146 1857996 : total_packed_bc_data++));
147 617354 : libmesh_assert_greater_equal (n_bcs, 0);
148 1240642 : total_packed_bc_data += n_bcs;
149 : }
150 :
151 606423 : for (unsigned short sf=0; sf != 2; ++sf)
152 : {
153 : const int n_bcs = cast_int<int>
154 605442 : (*(in + pre_indexing_size + indexing_size +
155 605442 : total_packed_bc_data++));
156 201160 : libmesh_assert_greater_equal (n_bcs, 0);
157 404282 : total_packed_bc_data += n_bcs;
158 : }
159 : }
160 :
161 : return
162 : #ifndef NDEBUG
163 : 1 + // Account for magic header
164 : #endif
165 219293 : pre_indexing_size + indexing_size + total_packed_bc_data;
166 : }
167 :
168 :
169 :
170 : template <>
171 : unsigned int
172 53956 : Packing<const Elem *>::packed_size (std::vector<largest_id_type>::iterator in)
173 : {
174 53956 : return packed_size(std::vector<largest_id_type>::const_iterator(in));
175 : }
176 :
177 :
178 :
179 : template <>
180 : unsigned int
181 216346 : Packing<const Elem *>::packable_size (const Elem * const & elem,
182 : const MeshBase * mesh)
183 : {
184 107950 : unsigned int variable_topology_size = 0;
185 216346 : if (elem->type() == C0POLYGON || elem->type() == C0POLYHEDRON)
186 : {
187 : // Store the dynamic node, side, and edge counts.
188 4 : variable_topology_size = 3;
189 :
190 : // A polygon's node ordering fully specifies its topology. A
191 : // polyhedron additionally needs each side's local node indices.
192 8 : if (elem->type() == C0POLYHEDRON)
193 37 : for (auto s : elem->side_index_range())
194 32 : variable_topology_size +=
195 40 : 1 + cast_int<unsigned int>(elem->nodes_on_side(s).size());
196 : }
197 :
198 : // We always communicate if we are on a boundary or not
199 107950 : unsigned int total_packed_bcs = 1;
200 216346 : const unsigned int n_sides = elem->n_sides();
201 :
202 107950 : largest_id_type on_boundary = 0;
203 1084406 : for (auto s : elem->side_index_range())
204 837757 : if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s))
205 : {
206 11822 : on_boundary = 1;
207 11822 : break;
208 : }
209 :
210 216346 : if (on_boundary)
211 : {
212 : // In this case we need another entry to check if we allow children on the boundary.
213 : // We only allow children on sides, edges and sheel faces are treated normally using
214 : // their top parents.
215 11822 : total_packed_bcs++;
216 23691 : if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side())
217 : {
218 23691 : total_packed_bcs += n_sides;
219 130861 : for (unsigned int s = 0; s != n_sides; ++s)
220 107170 : total_packed_bcs +=
221 107170 : mesh->get_boundary_info().n_raw_boundary_ids(elem,s);
222 : }
223 : }
224 :
225 216346 : if (elem->level() == 0)
226 : {
227 200538 : const unsigned int n_edges = elem->n_edges();
228 200538 : total_packed_bcs += n_edges;
229 1425474 : for (unsigned int e = 0; e != n_edges; ++e)
230 1224936 : total_packed_bcs +=
231 1224936 : mesh->get_boundary_info().n_edge_boundary_ids(elem,e);
232 :
233 200538 : total_packed_bcs += 2; // shellfaces
234 601614 : for (unsigned short sf=0; sf != 2; ++sf)
235 401076 : total_packed_bcs +=
236 401076 : mesh->get_boundary_info().n_shellface_boundary_ids(elem,sf);
237 : }
238 :
239 : return
240 : #ifndef NDEBUG
241 : 1 + // add an int for the magic header when testing
242 : #endif
243 216346 : header_size + variable_topology_size + elem->n_nodes() + n_sides*2 +
244 216346 : elem->packed_indexing_size() + total_packed_bcs;
245 : }
246 :
247 :
248 :
249 : template <>
250 : unsigned int
251 30872 : Packing<const Elem *>::packable_size (const Elem * const & elem,
252 : const DistributedMesh * mesh)
253 : {
254 30872 : return packable_size(elem, static_cast<const MeshBase *>(mesh));
255 : }
256 :
257 :
258 :
259 : template <>
260 : unsigned int
261 0 : Packing<const Elem *>::packable_size (const Elem * const & elem,
262 : const ParallelMesh * mesh)
263 : {
264 0 : return packable_size(elem, static_cast<const MeshBase *>(mesh));
265 : }
266 :
267 :
268 :
269 : template <>
270 : void
271 162280 : Packing<const Elem *>::pack (const Elem * const & elem,
272 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
273 : const MeshBase * mesh)
274 : {
275 53956 : libmesh_assert(elem);
276 :
277 : #ifndef NDEBUG
278 53956 : *data_out++ = elem_magic_header;
279 : #endif
280 :
281 : #ifdef LIBMESH_ENABLE_AMR
282 162280 : *data_out++ = (static_cast<largest_id_type>(elem->level()));
283 162280 : *data_out++ = (static_cast<largest_id_type>(elem->p_level()));
284 :
285 : // Encode both the refinement flag and whether the element has
286 : // children together. This coding is unambiguous because our
287 : // refinement state encoding starts at 0 and ends at
288 : // INVALID_REFINEMENTSTATE
289 : largest_id_type refinement_info =
290 162280 : static_cast<largest_id_type>(elem->refinement_flag());
291 53956 : if (elem->has_children())
292 2790 : refinement_info +=
293 : static_cast<largest_id_type>(Elem::INVALID_REFINEMENTSTATE) + 1;
294 53956 : *data_out++ = (refinement_info);
295 :
296 162280 : *data_out++ = (static_cast<largest_id_type>(elem->p_refinement_flag()));
297 : #else
298 : *data_out++ = (0);
299 : *data_out++ = (0);
300 : *data_out++ = (0);
301 : *data_out++ = (0);
302 : #endif
303 162280 : *data_out++ = (static_cast<largest_id_type>(elem->type()));
304 162280 : *data_out++ = (elem->processor_id());
305 162280 : *data_out++ = (elem->subdomain_id());
306 162280 : *data_out++ = (elem->id());
307 :
308 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
309 162280 : if (elem->valid_unique_id())
310 162280 : *data_out++ = (static_cast<largest_id_type>(elem->unique_id()));
311 : else
312 : // OK to send invalid unique id, we must not own this DOF
313 0 : *data_out++ = (static_cast<largest_id_type>(DofObject::invalid_unique_id));
314 : #endif
315 :
316 : #ifdef LIBMESH_ENABLE_AMR
317 : // use parent_ID of invalid_id to indicate a level 0 element
318 162280 : if (elem->level() == 0)
319 : {
320 100008 : *data_out++ =(DofObject::invalid_id);
321 100008 : *data_out++ =(DofObject::invalid_id);
322 : }
323 : else
324 : {
325 15808 : *data_out++ =(elem->parent()->id());
326 15808 : *data_out++ =(elem->parent()->which_child_am_i(elem));
327 : }
328 : #else
329 : *data_out++ = (DofObject::invalid_id);
330 : *data_out++ = (DofObject::invalid_id);
331 : #endif
332 :
333 219383 : if ((elem->dim() < LIBMESH_DIM) &&
334 57103 : elem->interior_parent())
335 84 : *data_out++ =(elem->interior_parent()->id());
336 : else
337 107856 : *data_out++ =(DofObject::invalid_id);
338 :
339 : const bool has_variable_topology =
340 162280 : elem->type() == C0POLYGON || elem->type() == C0POLYHEDRON;
341 53956 : if (has_variable_topology)
342 : {
343 6 : *data_out++ = elem->n_nodes();
344 6 : *data_out++ = elem->n_sides();
345 6 : *data_out++ = elem->n_edges();
346 : }
347 :
348 1194439 : for (const Node & node : elem->node_ref_range())
349 978203 : *data_out++ = node.id();
350 :
351 162280 : if (elem->type() == C0POLYHEDRON)
352 28 : for (auto s : elem->side_index_range())
353 : {
354 : const std::vector<unsigned int> side_nodes =
355 32 : elem->nodes_on_side(s);
356 40 : *data_out++ = side_nodes.size();
357 132 : for (const auto node : side_nodes)
358 144 : *data_out++ = node;
359 : }
360 :
361 : // Add the id of and the side for any return link from each neighbor
362 887411 : for (auto neigh : elem->neighbor_ptr_range())
363 : {
364 671175 : if (neigh)
365 : {
366 124418 : *data_out++ = (neigh->id());
367 124418 : if (neigh == remote_elem)
368 5970 : *data_out++ = (DofObject::invalid_id);
369 : else
370 115455 : *data_out++ = neigh->which_neighbor_am_i(elem);
371 : }
372 : else
373 : {
374 364294 : *data_out++ = (DofObject::invalid_id);
375 364294 : *data_out++ = (DofObject::invalid_id);
376 : }
377 : }
378 :
379 : // Add any DofObject indices
380 162280 : elem->pack_indexing(data_out);
381 :
382 : // We check if this is a boundary cell. We use the raw
383 : // IDs because we also communicate the parents which
384 : // will bring their associated IDs
385 162280 : largest_id_type on_boundary = 0;
386 826722 : for (auto s : elem->side_index_range())
387 628185 : if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s))
388 : {
389 17699 : on_boundary = 1;
390 17699 : break;
391 : }
392 :
393 53956 : *data_out++ = on_boundary;
394 :
395 162280 : if (on_boundary)
396 : {
397 17699 : *data_out++ = mesh->get_boundary_info().is_children_on_boundary_side();
398 : // Again, only do this if we allow children to hold boundary sides, the edges and
399 : // shell faces are treated normally using their top parents
400 17699 : if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side())
401 : {
402 11788 : std::vector<boundary_id_type> bcs;
403 103711 : for (auto s : elem->side_index_range())
404 : {
405 80118 : mesh->get_boundary_info().raw_boundary_ids(elem, s, bcs);
406 :
407 133570 : *data_out++ =(bcs.size());
408 :
409 103230 : for (const auto & bid : bcs)
410 30840 : *data_out++ = bid;
411 : }
412 : }
413 : }
414 :
415 : // If this is a coarse element,
416 : // Add any element side boundary condition ids
417 162280 : if (elem->level() == 0)
418 : {
419 100008 : std::vector<boundary_id_type> bcs;
420 1118945 : for (auto e : elem->edge_index_range())
421 : {
422 918517 : mesh->get_boundary_info().edge_boundary_ids(elem, e, bcs);
423 :
424 1531049 : *data_out++ =(bcs.size());
425 :
426 918661 : for (const auto & bid : bcs)
427 192 : *data_out++ = bid;
428 : }
429 :
430 451272 : for (unsigned short sf=0; sf != 2; ++sf)
431 : {
432 300848 : mesh->get_boundary_info().shellface_boundary_ids(elem, sf, bcs);
433 :
434 501688 : *data_out++ =(bcs.size());
435 :
436 310832 : for (const auto & bid : bcs)
437 13312 : *data_out++ = bid;
438 : }
439 : }
440 162280 : }
441 :
442 :
443 :
444 : template <>
445 : void
446 23070 : Packing<const Elem *>::pack (const Elem * const & elem,
447 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
448 : const DistributedMesh * mesh)
449 : {
450 23070 : pack(elem, data_out, static_cast<const MeshBase*>(mesh));
451 23070 : }
452 :
453 :
454 :
455 : template <>
456 : void
457 0 : Packing<const Elem *>::pack (const Elem * const & elem,
458 : std::back_insert_iterator<std::vector<largest_id_type>> data_out,
459 : const ParallelMesh * mesh)
460 : {
461 0 : pack(elem, data_out, static_cast<const MeshBase*>(mesh));
462 0 : }
463 :
464 :
465 :
466 : // FIXME - this needs serious work to be 64-bit compatible
467 : template <>
468 : Elem *
469 165337 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
470 : MeshBase * mesh)
471 : {
472 : #ifndef NDEBUG
473 54976 : const std::vector<largest_id_type>::const_iterator original_in = in;
474 :
475 54976 : const largest_id_type incoming_header = *in++;
476 54976 : libmesh_assert_equal_to (incoming_header, elem_magic_header);
477 : #endif
478 :
479 : // int 0: level
480 : const unsigned int level =
481 165337 : cast_int<unsigned int>(*in++);
482 :
483 : #ifdef LIBMESH_ENABLE_AMR
484 : // int 1: p level
485 : const unsigned int p_level =
486 165337 : cast_int<unsigned int>(*in++);
487 :
488 : // int 2: refinement flag and encoded has_children
489 165337 : const int rflag = cast_int<int>(*in++);
490 : const int invalid_rflag =
491 54976 : cast_int<int>(Elem::INVALID_REFINEMENTSTATE);
492 54976 : libmesh_assert_greater_equal (rflag, 0);
493 :
494 54976 : libmesh_assert_less (rflag, invalid_rflag*2+1);
495 :
496 54976 : const bool has_children = (rflag > invalid_rflag);
497 :
498 165337 : const Elem::RefinementState refinement_flag = has_children ?
499 3192 : cast_int<Elem::RefinementState>(rflag - invalid_rflag - 1) :
500 53912 : cast_int<Elem::RefinementState>(rflag);
501 :
502 : // int 3: p refinement flag
503 165337 : const int pflag = cast_int<int>(*in++);
504 54976 : libmesh_assert_greater_equal (pflag, 0);
505 54976 : libmesh_assert_less (pflag, Elem::INVALID_REFINEMENTSTATE);
506 : const Elem::RefinementState p_refinement_flag =
507 54976 : cast_int<Elem::RefinementState>(pflag);
508 : #else
509 : in += 3;
510 : #endif // LIBMESH_ENABLE_AMR
511 :
512 : // int 4: element type
513 165337 : const int typeint = cast_int<int>(*in++);
514 54976 : libmesh_assert_greater_equal (typeint, 0);
515 54976 : libmesh_assert_less (typeint, INVALID_ELEM);
516 : const ElemType type =
517 54976 : cast_int<ElemType>(typeint);
518 :
519 165337 : unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
520 165337 : unsigned int n_sides = Elem::type_to_n_sides_map[type];
521 165337 : unsigned int n_edges = Elem::type_to_n_edges_map[type];
522 :
523 : // int 5: processor id
524 : const processor_id_type processor_id =
525 165337 : 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 165337 : cast_int<subdomain_id_type>(*in++);
532 :
533 : // int 7: dof object id
534 : const dof_id_type id =
535 165337 : 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 165337 : 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 165337 : ? static_cast<dof_id_type>(*in++)
553 13200 : : 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 165337 : ? static_cast<unsigned int>(*in++)
563 13200 : : 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 165337 : 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 165337 : if (n_nodes == invalid_uint)
576 : {
577 6 : libmesh_error_msg_if(type != C0POLYGON && type != C0POLYHEDRON,
578 : "Unsupported variable-topology element "
579 : << Utility::enum_to_string(type));
580 :
581 6 : n_nodes = cast_int<unsigned int>(*in++);
582 6 : n_sides = cast_int<unsigned int>(*in++);
583 6 : n_edges = cast_int<unsigned int>(*in++);
584 :
585 6 : if (type == C0POLYGON)
586 3 : libmesh_error_msg_if
587 : (n_nodes < 3 || n_sides != n_nodes || n_edges != n_nodes,
588 : "Invalid packed C0POLYGON topology with "
589 : << n_nodes << " nodes, " << n_sides << " sides, and "
590 : << n_edges << " edges");
591 : else
592 3 : libmesh_error_msg_if
593 : (n_nodes < 4 || n_sides < 4,
594 : "Invalid packed C0POLYHEDRON topology with "
595 : << n_nodes << " nodes and " << n_sides << " sides");
596 : }
597 :
598 54976 : const auto node_ids_in = in;
599 110361 : in += n_nodes;
600 :
601 109952 : std::vector<std::vector<unsigned int>> polyhedron_side_nodes;
602 165337 : if (type == C0POLYHEDRON)
603 : {
604 3 : polyhedron_side_nodes.resize(n_sides);
605 4 : std::vector<bool> node_seen(n_nodes, false);
606 1 : unsigned int next_new_node = 0;
607 27 : for (auto & side_nodes : polyhedron_side_nodes)
608 : {
609 : const unsigned int n_side_nodes =
610 24 : cast_int<unsigned int>(*in++);
611 24 : libmesh_error_msg_if(n_side_nodes < 3,
612 : "Cannot unpack a C0POLYHEDRON side with only "
613 : << n_side_nodes << " nodes");
614 24 : side_nodes.resize(n_side_nodes);
615 132 : for (auto & node : side_nodes)
616 : {
617 108 : node = cast_int<unsigned int>(*in++);
618 108 : libmesh_error_msg_if(node >= n_nodes,
619 : "Invalid local node " << node
620 : << " in packed C0POLYHEDRON side");
621 :
622 144 : if (!node_seen[node])
623 : {
624 36 : libmesh_error_msg_if
625 : (node != next_new_node,
626 : "Packed C0POLYHEDRON side connectivity first encounters "
627 : "local node " << node << " where local node "
628 : << next_new_node << " was expected");
629 12 : node_seen[node] = true;
630 36 : ++next_new_node;
631 : }
632 : }
633 : }
634 :
635 3 : libmesh_error_msg_if
636 : (next_new_node != n_nodes && next_new_node + 1 != n_nodes,
637 : "Packed C0POLYHEDRON sides reference " << next_new_node
638 : << " vertices, but the element has " << n_nodes
639 : << " nodes; at most one unreferenced midpoint node is permitted");
640 : }
641 :
642 165337 : Elem * elem = mesh->query_elem_ptr(id);
643 :
644 : // if we already have this element, make sure its
645 : // properties match, and update any missing neighbor
646 : // links, but then go on
647 165337 : if (elem)
648 : {
649 9690 : libmesh_assert_equal_to (elem->level(), level);
650 9690 : libmesh_assert_equal_to (elem->id(), id);
651 : //#ifdef LIBMESH_ENABLE_UNIQUE_ID
652 : // No check for unique id sanity
653 : //#endif
654 9690 : libmesh_assert_equal_to (elem->processor_id(), processor_id);
655 9690 : libmesh_assert_equal_to (elem->subdomain_id(), subdomain_id);
656 9690 : libmesh_assert_equal_to (elem->type(), type);
657 :
658 29061 : if (type == C0POLYGON || type == C0POLYHEDRON)
659 : {
660 0 : libmesh_error_msg_if
661 : (elem->type() != type ||
662 : elem->n_nodes() != n_nodes ||
663 : elem->n_sides() != n_sides ||
664 : elem->n_edges() != n_edges,
665 : "Existing " << Utility::enum_to_string(type)
666 : << " topology does not match its packed topology");
667 :
668 0 : for (unsigned int n = 0; n != n_nodes; ++n)
669 0 : libmesh_error_msg_if
670 : (elem->node_id(n) !=
671 : cast_int<dof_id_type>(*(node_ids_in + n)),
672 : "Existing " << Utility::enum_to_string(type)
673 : << " local node " << n
674 : << " does not match its packed node");
675 :
676 0 : if (type == C0POLYHEDRON)
677 0 : for (auto s : elem->side_index_range())
678 0 : libmesh_error_msg_if
679 : (elem->nodes_on_side(s) != polyhedron_side_nodes[s],
680 : "Existing C0POLYHEDRON side " << s
681 : << " does not match its packed topology");
682 : }
683 :
684 9690 : libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
685 9690 : libmesh_assert_equal_to (elem->n_sides(), n_sides);
686 9690 : libmesh_assert_equal_to (elem->n_edges(), n_edges);
687 :
688 : #ifndef NDEBUG
689 : // All our nodes should be correct
690 105990 : for (unsigned int i=0; i != n_nodes; ++i)
691 96300 : libmesh_assert(elem->node_id(i) ==
692 : cast_int<dof_id_type>(*(node_ids_in + i)));
693 : #endif
694 :
695 : #ifdef LIBMESH_ENABLE_AMR
696 9690 : libmesh_assert_equal_to (elem->refinement_flag(), refinement_flag);
697 9690 : libmesh_assert_equal_to (elem->has_children(), has_children);
698 :
699 : #ifdef DEBUG
700 9690 : if (elem->active())
701 : {
702 8712 : libmesh_assert_equal_to (elem->p_level(), p_level);
703 8712 : libmesh_assert_equal_to (elem->p_refinement_flag(), p_refinement_flag);
704 : }
705 : #endif
706 :
707 9690 : libmesh_assert (!level || elem->parent() != nullptr);
708 9690 : libmesh_assert (!level || elem->parent()->id() == parent_id);
709 9690 : libmesh_assert (!level || elem->parent()->child_ptr(which_child_am_i) == elem);
710 : #endif
711 : // Our interior_parent link should be "close to" correct - we
712 : // may have to update it, but we can check for some
713 : // inconsistencies.
714 : {
715 : // If the sending processor sees no interior_parent here, we'd
716 : // better agree.
717 29061 : if (interior_parent_id == DofObject::invalid_id)
718 : {
719 28977 : if (elem->dim() < LIBMESH_DIM)
720 2623 : libmesh_assert (!(elem->interior_parent()));
721 : }
722 :
723 : // If the sending processor has a remote_elem interior_parent,
724 : // then all we know is that we'd better have *some*
725 : // interior_parent
726 84 : else if (interior_parent_id == remote_elem->id())
727 : {
728 0 : libmesh_assert(elem->interior_parent());
729 : }
730 : else
731 : {
732 : Elem * ip =
733 84 : mesh->interior_mesh().query_elem_ptr(interior_parent_id);
734 :
735 : // The sending processor sees an interior parent here, so
736 : // if we don't have that interior element, then we'd
737 : // better have a remote_elem signifying that fact.
738 84 : if (!ip)
739 0 : libmesh_assert_equal_to (elem->interior_parent(), remote_elem);
740 : else
741 : {
742 : // The sending processor has an interior_parent here,
743 : // and we have that element, but that does *NOT* mean
744 : // we're already linking to it. Perhaps we initially
745 : // received elem from a processor on which the
746 : // interior_parent link was remote?
747 28 : libmesh_assert(elem->interior_parent() == ip ||
748 : elem->interior_parent() == remote_elem);
749 :
750 : // If the link was originally remote, update it
751 84 : if (elem->interior_parent() == remote_elem)
752 : {
753 0 : elem->set_interior_parent(ip);
754 : }
755 : }
756 : }
757 : }
758 :
759 : // Our neighbor links should be "close to" correct - we may have
760 : // to update a remote_elem link, and we can check for possible
761 : // inconsistencies along the way.
762 : //
763 : // Even for subactive elements, we'll try to keep neighbor links
764 : // in good shape now, if only so any future find_neighbors() is
765 : // idempotent.
766 173823 : for (auto n : elem->side_index_range())
767 : {
768 : const dof_id_type neighbor_id =
769 144762 : cast_int<dof_id_type>(*in++);
770 :
771 : const dof_id_type neighbor_side =
772 144762 : cast_int<dof_id_type>(*in++);
773 :
774 : // If the sending processor sees a domain boundary here,
775 : // we'd better agree ... unless all we see is a remote_elem?
776 : // In that case maybe we just couldn't keep up with a user's
777 : // delete_elem. Let's trust them.
778 144762 : if (neighbor_id == DofObject::invalid_id)
779 : {
780 14634 : const Elem * my_neigh = elem->neighbor_ptr(n);
781 14634 : if (my_neigh == remote_elem)
782 1 : elem->set_neighbor(n, nullptr);
783 : else
784 4883 : libmesh_assert (!my_neigh);
785 14634 : continue;
786 4866 : }
787 :
788 : // If the sending processor has a remote_elem neighbor here,
789 : // then all we know is that we'd better *not* have a domain
790 : // boundary ... except that maybe it's the *sending*
791 : // processor who missed a delete_elem we saw.
792 130128 : if (neighbor_id == remote_elem->id())
793 : {
794 : // At this level of the code we can't even assert in
795 : // cases where the neighbor should know what they're
796 : // talking about, so skip it.
797 :
798 : // libmesh_assert(elem->neighbor_ptr(n));
799 6104 : continue;
800 : }
801 :
802 120979 : Elem * neigh = mesh->query_elem_ptr(neighbor_id);
803 :
804 : // The sending processor sees a neighbor here, so if we
805 : // don't have that neighboring element, then we'd better
806 : // have a remote_elem signifying that fact.
807 120979 : if (!neigh)
808 : {
809 2442 : libmesh_assert_equal_to (elem->neighbor_ptr(n), remote_elem);
810 4884 : continue;
811 : }
812 :
813 : // The sending processor has a neighbor here, and we have
814 : // that element, but that does *NOT* mean we're already
815 : // linking to it. Perhaps we initially received both elem
816 : // and neigh from processors on which their mutual link was
817 : // remote?
818 37895 : libmesh_assert(elem->neighbor_ptr(n) == neigh ||
819 : elem->neighbor_ptr(n) == remote_elem);
820 :
821 : // If the link was originally remote, we should update it,
822 : // and make sure the appropriate parts of its family link
823 : // back to us.
824 151548 : if (elem->neighbor_ptr(n) == remote_elem)
825 : {
826 0 : elem->set_neighbor(n, neigh);
827 : }
828 : else
829 37895 : libmesh_assert(elem->subactive() ||
830 : neigh->level() < elem->level() ||
831 : neigh->neighbor_ptr(neighbor_side) == elem);
832 :
833 113653 : if (neighbor_side != libMesh::invalid_uint)
834 113653 : elem->make_links_to_me_local(n, neighbor_side);
835 : }
836 :
837 : // Our p level and refinement flags should be "close to" correct
838 : // if we're not an active element - we might have a p level
839 : // increased or decreased by changes in remote_elem children.
840 : //
841 : // But if we have remote_elem children, then we shouldn't be
842 : // doing a projection on this inactive element on this
843 : // processor, so we won't need correct p settings. Couldn't
844 : // hurt to update, though.
845 : #ifdef LIBMESH_ENABLE_AMR
846 38751 : if (elem->processor_id() != mesh->processor_id())
847 : {
848 : // Do this simultaneously; otherwise we can get a false
849 : // positive when a hack_p_level or set_p_refineemnt_flag
850 : // assertion sees inconsistency between an old flag and new
851 : // value or vice-versa
852 4106 : elem->hack_p_level_and_refinement_flag(p_level, p_refinement_flag);
853 : }
854 : #endif // LIBMESH_ENABLE_AMR
855 :
856 : // FIXME: We should add some debug mode tests to ensure that the
857 : // encoded indexing and boundary conditions are consistent.
858 : }
859 : else
860 : {
861 : // We don't already have the element, so we need to create it.
862 :
863 : // Find the parent if necessary
864 90990 : Elem * parent = nullptr;
865 : #ifdef LIBMESH_ENABLE_AMR
866 : // Find a child element's parent
867 136276 : if (level > 0)
868 : {
869 : // Note that we must be very careful to construct the send
870 : // connectivity so that parents are encountered before
871 : // children. If we get here and can't find the parent that
872 : // is a fatal error.
873 1056 : parent = mesh->elem_ptr(parent_id);
874 : }
875 : // Or assert that the sending processor sees no parent
876 : else
877 44934 : libmesh_assert_equal_to (parent_id, DofObject::invalid_id);
878 : #else
879 : // No non-level-0 elements without AMR
880 : libmesh_assert_equal_to (level, 0);
881 : #endif
882 :
883 136276 : if (type == C0POLYGON)
884 3 : elem = std::make_unique<C0Polygon>(n_nodes, parent).release();
885 136273 : else if (type == C0POLYHEDRON)
886 : {
887 5 : std::vector<std::shared_ptr<Polygon>> sides(n_sides);
888 27 : for (auto s : index_range(sides))
889 : {
890 16 : const auto & side_nodes = polyhedron_side_nodes[s];
891 : auto side = std::make_shared<C0Polygon>
892 40 : (cast_int<unsigned int>(side_nodes.size()));
893 132 : for (auto n : index_range(side_nodes))
894 : {
895 : const dof_id_type node_id =
896 : cast_int<dof_id_type>
897 180 : (*(node_ids_in + side_nodes[n]));
898 108 : side->set_node(n, mesh->node_ptr(node_id));
899 : }
900 16 : sides[s] = std::move(side);
901 : }
902 :
903 3 : std::unique_ptr<Node> generated_mid_node;
904 : auto polyhedron = std::make_unique<C0Polyhedron>
905 3 : (sides, generated_mid_node, parent);
906 :
907 3 : libmesh_error_msg_if
908 : (polyhedron->n_nodes() != n_nodes,
909 : "Packed C0POLYHEDRON has " << n_nodes
910 : << " nodes, but reconstructing its topology produced "
911 : << polyhedron->n_nodes() << " nodes");
912 :
913 3 : if (generated_mid_node)
914 : {
915 : const dof_id_type mid_node_id =
916 3 : cast_int<dof_id_type>(*(node_ids_in + n_nodes - 1));
917 4 : polyhedron->set_node(n_nodes - 1,
918 3 : mesh->node_ptr(mid_node_id));
919 : }
920 :
921 1 : elem = polyhedron.release();
922 1 : }
923 : else
924 136270 : elem = Elem::build(type,parent).release();
925 45286 : libmesh_assert (elem);
926 :
927 : #ifdef LIBMESH_ENABLE_AMR
928 136276 : if (level != 0)
929 : {
930 : // Since this is a newly created element, the parent must
931 : // have previously thought of this child as a remote element.
932 352 : libmesh_assert_equal_to (parent->child_ptr(which_child_am_i), remote_elem);
933 :
934 1056 : parent->add_child(elem, which_child_am_i);
935 : }
936 :
937 : // Assign the refinement flags and levels
938 136276 : elem->set_p_level(p_level);
939 45286 : elem->set_refinement_flag(refinement_flag);
940 45286 : elem->set_p_refinement_flag(p_refinement_flag);
941 45286 : libmesh_assert_equal_to (elem->level(), level);
942 :
943 : // If this element should have children, assign remote_elem to
944 : // all of them for now, for consistency. Later unpacked
945 : // elements may overwrite that.
946 136276 : if (has_children)
947 : {
948 258 : const unsigned int nc = elem->n_children();
949 1206 : for (unsigned int c=0; c != nc; ++c)
950 948 : elem->add_child(const_cast<RemoteElem *>(remote_elem), c);
951 : }
952 :
953 : #endif // LIBMESH_ENABLE_AMR
954 :
955 : // Assign the IDs
956 136276 : elem->subdomain_id() = subdomain_id;
957 136276 : elem->processor_id() = processor_id;
958 136276 : elem->set_id() = id;
959 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
960 45286 : elem->set_unique_id(unique_id);
961 : #endif
962 :
963 : // Assign the connectivity
964 45286 : libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
965 :
966 136276 : if (type == C0POLYHEDRON)
967 : {
968 3 : libmesh_error_msg_if(elem->n_sides() != n_sides,
969 : "Packed C0POLYHEDRON has " << n_sides
970 : << " sides, but reconstructing its topology produced "
971 : << elem->n_sides() << " sides");
972 3 : libmesh_error_msg_if(elem->n_edges() != n_edges,
973 : "Packed C0POLYHEDRON has " << n_edges
974 : << " edges, but reconstructing its topology produced "
975 : << elem->n_edges() << " edges");
976 :
977 42 : for (unsigned int n = 0; n != n_nodes; ++n)
978 65 : libmesh_error_msg_if
979 : (elem->node_id(n) !=
980 : cast_int<dof_id_type>(*(node_ids_in + n)),
981 : "Packed C0POLYHEDRON local node ordering was not preserved");
982 : }
983 : else
984 855117 : for (unsigned int n=0; n != n_nodes; n++)
985 : elem->set_node
986 957894 : (n, mesh->node_ptr
987 957894 : (cast_int<dof_id_type>(*(node_ids_in + n))));
988 :
989 : // Set interior_parent if found
990 : {
991 : // We may be unpacking an element that was a ghost element on the
992 : // sender, in which case the element's interior_parent may not be
993 : // known by the packed element. We'll have to set such
994 : // interior_parents to remote_elem ourselves and wait for a
995 : // later packed element to give us better information.
996 136276 : if (interior_parent_id == remote_elem->id())
997 : {
998 : elem->set_interior_parent
999 0 : (const_cast<RemoteElem *>(remote_elem));
1000 : }
1001 136276 : else if (interior_parent_id != DofObject::invalid_id)
1002 : {
1003 : // If we don't have the interior parent element, then it's
1004 : // a remote_elem until we get it.
1005 : Elem * ip =
1006 84 : mesh->interior_mesh().query_elem_ptr(interior_parent_id);
1007 84 : if (!ip )
1008 : elem->set_interior_parent
1009 0 : (const_cast<RemoteElem *>(remote_elem));
1010 : else
1011 84 : elem->set_interior_parent(ip);
1012 : }
1013 : }
1014 :
1015 676789 : for (auto n : elem->side_index_range())
1016 : {
1017 : const dof_id_type neighbor_id =
1018 540513 : cast_int<dof_id_type>(*in++);
1019 :
1020 : const dof_id_type neighbor_side =
1021 540513 : cast_int<dof_id_type>(*in++);
1022 :
1023 540513 : if (neighbor_id == DofObject::invalid_id)
1024 357422 : continue;
1025 :
1026 : // We may be unpacking an element that was a ghost element on the
1027 : // sender, in which case the element's neighbors may not all be
1028 : // known by the packed element. We'll have to set such
1029 : // neighbors to remote_elem ourselves and wait for a later
1030 : // packed element to give us better information.
1031 4544 : if (neighbor_id == remote_elem->id())
1032 : {
1033 96 : elem->set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
1034 96 : continue;
1035 : }
1036 :
1037 : // If we don't have the neighbor element, then it's a
1038 : // remote_elem until we get it.
1039 4448 : Elem * neigh = mesh->query_elem_ptr(neighbor_id);
1040 4448 : if (!neigh)
1041 : {
1042 2089 : elem->set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
1043 2089 : continue;
1044 : }
1045 :
1046 : // If we have the neighbor element, then link to it, and
1047 : // make sure any appropriate parts of its family link back
1048 : // to us.
1049 2359 : elem->set_neighbor(n, neigh);
1050 :
1051 2359 : if (neighbor_side != libMesh::invalid_uint)
1052 2359 : elem->make_links_to_me_local(n, neighbor_side);
1053 : }
1054 :
1055 136276 : elem->unpack_indexing(in);
1056 :
1057 136276 : mesh->add_elem(elem);
1058 : }
1059 :
1060 165337 : in += elem->packed_indexing_size();
1061 :
1062 : // We check if this is cell holds a boundary ID or not
1063 165337 : auto on_boundary = *in++;
1064 165337 : if (on_boundary)
1065 : {
1066 : // Only treat the sides with caution. This is because we might hold boundary IDs
1067 : // on the sides of the children. This is not supported for edges and shell faces, thus
1068 : // they are treated assuming that only top parents can hold the IDs.
1069 19229 : auto children_on_boundary = *in++;
1070 19229 : if (elem->level() == 0 || children_on_boundary)
1071 : {
1072 107735 : for (auto s : elem->side_index_range())
1073 : {
1074 : const boundary_id_type num_bcs =
1075 88506 : cast_int<boundary_id_type>(*in++);
1076 :
1077 114870 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1078 8778 : mesh->get_boundary_info().add_side
1079 35142 : (elem, s, cast_int<boundary_id_type>(*in++));
1080 : }
1081 : }
1082 : }
1083 :
1084 : // If this is a coarse element,
1085 : // add any element side or edge boundary condition ids
1086 165337 : if (level == 0)
1087 : {
1088 1086794 : for (auto e : elem->edge_index_range())
1089 : {
1090 : const boundary_id_type num_bcs =
1091 934657 : cast_int<boundary_id_type>(*in++);
1092 :
1093 934801 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1094 48 : mesh->get_boundary_info().add_edge
1095 192 : (elem, e, cast_int<boundary_id_type>(*in++));
1096 : }
1097 :
1098 456411 : for (unsigned short sf=0; sf != 2; ++sf)
1099 : {
1100 : const boundary_id_type num_bcs =
1101 304274 : cast_int<boundary_id_type>(*in++);
1102 :
1103 314258 : for (boundary_id_type bc_it=0; bc_it < num_bcs; bc_it++)
1104 3328 : mesh->get_boundary_info().add_shellface
1105 13312 : (elem, sf, cast_int<boundary_id_type>(*in++));
1106 : }
1107 : }
1108 :
1109 : // Return the new element
1110 220313 : return elem;
1111 55385 : }
1112 :
1113 :
1114 :
1115 : template <>
1116 : Elem *
1117 23070 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
1118 : DistributedMesh * mesh)
1119 : {
1120 23070 : return unpack(in, static_cast<MeshBase*>(mesh));
1121 : }
1122 :
1123 :
1124 :
1125 : template <>
1126 : Elem *
1127 0 : Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
1128 : ParallelMesh * mesh)
1129 : {
1130 0 : return unpack(in, static_cast<MeshBase*>(mesh));
1131 : }
1132 :
1133 : } // namespace Parallel
1134 :
1135 : } // namespace libMesh
|
