1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 // Local includes
21 #include "libmesh/boundary_info.h"
22 #include "libmesh/ghosting_functor.h"
23 #include "libmesh/unstructured_mesh.h"
24 #include "libmesh/libmesh_logging.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/mesh_tools.h" // For n_levels
27 #include "libmesh/parallel.h"
28 #include "libmesh/remote_elem.h"
29 #include "libmesh/namebased_io.h"
30 #include "libmesh/partitioner.h"
31 #include "libmesh/enum_order.h"
32 #include "libmesh/mesh_communication.h"
33
34 // C++ includes
35 #include <fstream>
36 #include <sstream>
37 #include <iomanip>
38 #include <unordered_map>
39
40 // C includes
41 #include <sys/types.h> // for pid_t
42 #include <unistd.h> // for getpid(), unlink()
43
44 namespace libMesh
45 {
46
47
48 // ------------------------------------------------------------
49 // UnstructuredMesh class member functions
UnstructuredMesh(const Parallel::Communicator & comm_in,unsigned char d)50 UnstructuredMesh::UnstructuredMesh (const Parallel::Communicator & comm_in,
51 unsigned char d) :
52 MeshBase (comm_in,d)
53 {
54 libmesh_assert (libMesh::initialized());
55 }
56
57
58
copy_nodes_and_elements(const UnstructuredMesh & other_mesh,const bool skip_find_neighbors,dof_id_type element_id_offset,dof_id_type node_id_offset,unique_id_type unique_id_offset)59 void UnstructuredMesh::copy_nodes_and_elements(const UnstructuredMesh & other_mesh,
60 const bool skip_find_neighbors,
61 dof_id_type element_id_offset,
62 dof_id_type node_id_offset,
63 unique_id_type
64 #ifdef LIBMESH_ENABLE_UNIQUE_ID
65 unique_id_offset
66 #endif
67 )
68 {
69 LOG_SCOPE("copy_nodes_and_elements()", "UnstructuredMesh");
70
71 std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
72 extra_int_maps = this->merge_extra_integer_names(other_mesh);
73
74 const unsigned int n_old_node_ints = extra_int_maps.second.size(),
75 n_new_node_ints = _node_integer_names.size(),
76 n_old_elem_ints = extra_int_maps.first.size(),
77 n_new_elem_ints = _elem_integer_names.size();
78
79 // If we are partitioned into fewer parts than the incoming mesh,
80 // then we need to "wrap" the other Mesh's processor ids to fit
81 // within our range. This can happen, for example, while stitching
82 // ReplicatedMeshes with small numbers of elements in parallel...
83 bool wrap_proc_ids = (_n_parts < other_mesh._n_parts);
84
85 // We're assuming our subclass data needs no copy
86 libmesh_assert_equal_to (_is_prepared, other_mesh._is_prepared);
87
88 // We're assuming the other mesh has proper element number ordering,
89 // so that we add parents before their children, and that the other
90 // mesh is consistently partitioned.
91 #ifdef DEBUG
92 MeshTools::libmesh_assert_valid_amr_elem_ids(other_mesh);
93 MeshTools::libmesh_assert_valid_procids<Node>(other_mesh);
94 #endif
95
96 //Copy in Nodes
97 {
98 //Preallocate Memory if necessary
99 this->reserve_nodes(other_mesh.n_nodes());
100
101 for (const auto & oldn : other_mesh.node_ptr_range())
102 {
103 processor_id_type added_pid = cast_int<processor_id_type>
104 (wrap_proc_ids ? oldn->processor_id() % _n_parts : oldn->processor_id());
105
106 // Add new nodes in old node Point locations
107 Node * newn =
108 this->add_point(*oldn,
109 oldn->id() + node_id_offset,
110 added_pid);
111
112 newn->add_extra_integers(n_new_node_ints);
113 for (unsigned int i = 0; i != n_old_node_ints; ++i)
114 newn->set_extra_integer(extra_int_maps.second[i],
115 oldn->get_extra_integer(i));
116
117 #ifdef LIBMESH_ENABLE_UNIQUE_ID
118 newn->set_unique_id(oldn->unique_id() + unique_id_offset);
119 #endif
120 }
121 }
122
123 //Copy in Elements
124 {
125 //Preallocate Memory if necessary
126 this->reserve_elem(other_mesh.n_elem());
127
128 // Declare a map linking old and new elements, needed to copy the neighbor lists
129 typedef std::unordered_map<const Elem *, Elem *> map_type;
130 map_type old_elems_to_new_elems, ip_map;
131
132 // Loop over the elements
133 for (const auto & old : other_mesh.element_ptr_range())
134 {
135 // Build a new element
136 Elem * newparent = old->parent() ?
137 this->elem_ptr(old->parent()->id() + element_id_offset) :
138 nullptr;
139 auto el = Elem::build(old->type(), newparent);
140
141 el->subdomain_id() = old->subdomain_id();
142
143 // Hold off on trying to set the interior parent because we may actually
144 // add lower dimensional elements before their interior parents
145 if (el->dim() < other_mesh.mesh_dimension() && old->interior_parent())
146 ip_map[old] = el.get();
147
148 #ifdef LIBMESH_ENABLE_AMR
149 if (old->has_children())
150 for (unsigned int c = 0, nc = old->n_children(); c != nc; ++c)
151 if (old->child_ptr(c) == remote_elem)
152 el->add_child(const_cast<RemoteElem *>(remote_elem), c);
153
154 //Create the parent's child pointers if necessary
155 if (newparent)
156 {
157 unsigned int oldc = old->parent()->which_child_am_i(old);
158 newparent->add_child(el.get(), oldc);
159 }
160
161 // Copy the refinement flags
162 el->set_refinement_flag(old->refinement_flag());
163
164 // Use hack_p_level since we may not have sibling elements
165 // added yet
166 el->hack_p_level(old->p_level());
167
168 el->set_p_refinement_flag(old->p_refinement_flag());
169 #endif // #ifdef LIBMESH_ENABLE_AMR
170
171 //Assign all the nodes
172 for (auto i : el->node_index_range())
173 el->set_node(i) =
174 this->node_ptr(old->node_id(i) + node_id_offset);
175
176 // And start it off with the same processor id (mod _n_parts).
177 el->processor_id() = cast_int<processor_id_type>
178 (wrap_proc_ids ? old->processor_id() % _n_parts : old->processor_id());
179
180 // Give it the same element and unique ids
181 el->set_id(old->id() + element_id_offset);
182
183 el->add_extra_integers(n_new_elem_ints);
184 for (unsigned int i = 0; i != n_old_elem_ints; ++i)
185 el->set_extra_integer(extra_int_maps.first[i],
186 old->get_extra_integer(i));
187
188 #ifdef LIBMESH_ENABLE_UNIQUE_ID
189 el->set_unique_id(old->unique_id() + unique_id_offset);
190 #endif
191
192 el->set_mapping_type(old->mapping_type());
193 el->set_mapping_data(old->mapping_data());
194
195 //Hold onto it
196 if (!skip_find_neighbors)
197 {
198 for (auto s : old->side_index_range())
199 if (old->neighbor_ptr(s) == remote_elem)
200 el->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
201 this->add_elem(std::move(el));
202 }
203 else
204 {
205 Elem * new_el = this->add_elem(std::move(el));
206 old_elems_to_new_elems[old] = new_el;
207 }
208 }
209
210 for (auto & elem_pair : ip_map)
211 elem_pair.second->set_interior_parent(
212 this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset));
213
214 // Loop (again) over the elements to fill in the neighbors
215 if (skip_find_neighbors)
216 {
217 old_elems_to_new_elems[remote_elem] = const_cast<RemoteElem*>(remote_elem);
218
219 for (const auto & old_elem : other_mesh.element_ptr_range())
220 {
221 Elem * new_elem = old_elems_to_new_elems[old_elem];
222 for (auto s : old_elem->side_index_range())
223 {
224 const Elem * old_neighbor = old_elem->neighbor_ptr(s);
225 Elem * new_neighbor = old_elems_to_new_elems[old_neighbor];
226 new_elem->set_neighbor(s, new_neighbor);
227 }
228 }
229 }
230 }
231
232 //Finally prepare the new Mesh for use. Keep the same numbering and
233 //partitioning for now.
234 this->allow_renumbering(false);
235 this->allow_remote_element_removal(false);
236 this->allow_find_neighbors(!skip_find_neighbors);
237
238 // We should generally be able to skip *all* partitioning here
239 // because we're only adding one already-consistent mesh to another.
240 this->skip_partitioning(true);
241
242 this->prepare_for_use();
243
244 //But in the long term, use the same renumbering and partitioning
245 //policies as our source mesh.
246 this->allow_find_neighbors(other_mesh.allow_find_neighbors());
247 this->allow_renumbering(other_mesh.allow_renumbering());
248 this->allow_remote_element_removal(other_mesh.allow_remote_element_removal());
249 this->skip_partitioning(other_mesh._skip_all_partitioning);
250 this->skip_noncritical_partitioning(other_mesh._skip_noncritical_partitioning);
251 }
252
253
254
~UnstructuredMesh()255 UnstructuredMesh::~UnstructuredMesh ()
256 {
257 // this->clear (); // Nothing to clear at this level
258
259 libmesh_exceptionless_assert (!libMesh::closed());
260 }
261
262
263
264
265
find_neighbors(const bool reset_remote_elements,const bool reset_current_list)266 void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
267 const bool reset_current_list)
268 {
269 // We might actually want to run this on an empty mesh
270 // (e.g. the boundary mesh for a nonexistent bcid!)
271 // libmesh_assert_not_equal_to (this->n_nodes(), 0);
272 // libmesh_assert_not_equal_to (this->n_elem(), 0);
273
274 // This function must be run on all processors at once
275 parallel_object_only();
276
277 LOG_SCOPE("find_neighbors()", "Mesh");
278
279 //TODO:[BSK] This should be removed later?!
280 if (reset_current_list)
281 for (const auto & e : this->element_ptr_range())
282 for (auto s : e->side_index_range())
283 if (e->neighbor_ptr(s) != remote_elem || reset_remote_elements)
284 e->set_neighbor(s, nullptr);
285
286 // Find neighboring elements by first finding elements
287 // with identical side keys and then check to see if they
288 // are neighbors
289 {
290 // data structures -- Use the hash_multimap if available
291 typedef unsigned int key_type;
292 typedef std::pair<Elem *, unsigned char> val_type;
293 typedef std::unordered_multimap<key_type, val_type> map_type;
294
295 // A map from side keys to corresponding elements & side numbers
296 map_type side_to_elem_map;
297
298 // Pull objects out of the loop to reduce heap operations
299 std::unique_ptr<Elem> my_side, their_side;
300
301 for (const auto & element : this->element_ptr_range())
302 {
303 for (auto ms : element->side_index_range())
304 {
305 next_side:
306 // If we haven't yet found a neighbor on this side, try.
307 // Even if we think our neighbor is remote, that
308 // information may be out of date.
309 if (element->neighbor_ptr(ms) == nullptr ||
310 element->neighbor_ptr(ms) == remote_elem)
311 {
312 // Get the key for the side of this element
313 const unsigned int key = element->key(ms);
314
315 // Look for elements that have an identical side key
316 auto bounds = side_to_elem_map.equal_range(key);
317
318 // May be multiple keys, check all the possible
319 // elements which _might_ be neighbors.
320 if (bounds.first != bounds.second)
321 {
322 // Get the side for this element
323 element->side_ptr(my_side, ms);
324
325 // Look at all the entries with an equivalent key
326 while (bounds.first != bounds.second)
327 {
328 // Get the potential element
329 Elem * neighbor = bounds.first->second.first;
330
331 // Get the side for the neighboring element
332 const unsigned int ns = bounds.first->second.second;
333 neighbor->side_ptr(their_side, ns);
334 //libmesh_assert(my_side.get());
335 //libmesh_assert(their_side.get());
336
337 // If found a match with my side
338 //
339 // We need special tests here for 1D:
340 // since parents and children have an equal
341 // side (i.e. a node), we need to check
342 // ns != ms, and we also check level() to
343 // avoid setting our neighbor pointer to
344 // any of our neighbor's descendants
345 if ((*my_side == *their_side) &&
346 (element->level() == neighbor->level()) &&
347 ((element->dim() != 1) || (ns != ms)))
348 {
349 // So share a side. Is this a mixed pair
350 // of subactive and active/ancestor
351 // elements?
352 // If not, then we're neighbors.
353 // If so, then the subactive's neighbor is
354
355 if (element->subactive() ==
356 neighbor->subactive())
357 {
358 // an element is only subactive if it has
359 // been coarsened but not deleted
360 element->set_neighbor (ms,neighbor);
361 neighbor->set_neighbor(ns,element);
362 }
363 else if (element->subactive())
364 {
365 element->set_neighbor(ms,neighbor);
366 }
367 else if (neighbor->subactive())
368 {
369 neighbor->set_neighbor(ns,element);
370 }
371 side_to_elem_map.erase (bounds.first);
372
373 // get out of this nested crap
374 goto next_side;
375 }
376
377 ++bounds.first;
378 }
379 }
380
381 // didn't find a match...
382 // Build the map entry for this element
383 side_to_elem_map.emplace
384 (key, std::make_pair(element, cast_int<unsigned char>(ms)));
385 }
386 }
387 }
388 }
389
390 #ifdef LIBMESH_ENABLE_AMR
391
392 /**
393 * Here we look at all of the child elements which
394 * don't already have valid neighbors.
395 *
396 * If a child element has a nullptr neighbor it is
397 * either because it is on the boundary or because
398 * its neighbor is at a different level. In the
399 * latter case we must get the neighbor from the
400 * parent.
401 *
402 * If a child element has a remote_elem neighbor
403 * on a boundary it shares with its parent, that
404 * info may have become out-dated through coarsening
405 * of the neighbor's parent. In this case, if the
406 * parent's neighbor is active then the child should
407 * share it.
408 *
409 * Furthermore, that neighbor better be active,
410 * otherwise we missed a child somewhere.
411 *
412 *
413 * We also need to look through children ordered by increasing
414 * refinement level in order to add new interior_parent() links in
415 * boundary elements which have just been generated by refinement,
416 * and fix links in boundary elements whose previous
417 * interior_parent() has just been coarsened away.
418 */
419 const unsigned int n_levels = MeshTools::n_levels(*this);
420 for (unsigned int level = 1; level < n_levels; ++level)
421 {
422 for (auto & current_elem : as_range(level_elements_begin(level),
423 level_elements_end(level)))
424 {
425 libmesh_assert(current_elem);
426 Elem * parent = current_elem->parent();
427 libmesh_assert(parent);
428 const unsigned int my_child_num = parent->which_child_am_i(current_elem);
429
430 for (auto s : current_elem->side_index_range())
431 {
432 if (current_elem->neighbor_ptr(s) == nullptr ||
433 (current_elem->neighbor_ptr(s) == remote_elem &&
434 parent->is_child_on_side(my_child_num, s)))
435 {
436 Elem * neigh = parent->neighbor_ptr(s);
437
438 // If neigh was refined and had non-subactive children
439 // made remote earlier, then our current elem should
440 // actually have one of those remote children as a
441 // neighbor
442 if (neigh &&
443 (neigh->ancestor() ||
444 // If neigh has subactive children which should have
445 // matched as neighbors of the current element but
446 // did not, then those likewise must be remote
447 // children.
448 (current_elem->subactive() && neigh->has_children() &&
449 (neigh->level()+1) == current_elem->level())))
450 {
451 #ifdef DEBUG
452 // Let's make sure that "had children made remote"
453 // situation is actually the case
454 libmesh_assert(neigh->has_children());
455 bool neigh_has_remote_children = false;
456 for (auto & child : neigh->child_ref_range())
457 if (&child == remote_elem)
458 neigh_has_remote_children = true;
459 libmesh_assert(neigh_has_remote_children);
460
461 // And let's double-check that we don't have
462 // a remote_elem neighboring an active local element
463 if (current_elem->active())
464 libmesh_assert_not_equal_to (current_elem->processor_id(),
465 this->processor_id());
466 #endif // DEBUG
467 neigh = const_cast<RemoteElem *>(remote_elem);
468 }
469 // If neigh and current_elem are more than one level
470 // apart, figuring out whether we have a remote
471 // neighbor here becomes much harder.
472 else if (neigh && (current_elem->subactive() &&
473 neigh->has_children()))
474 {
475 // Find the deepest descendant of neigh which
476 // we could consider for a neighbor. If we run
477 // out of neigh children, then that's our
478 // neighbor. If we find a potential neighbor
479 // with remote_children and we don't find any
480 // potential neighbors among its non-remote
481 // children, then our neighbor must be remote.
482 while (neigh != remote_elem &&
483 neigh->has_children())
484 {
485 bool found_neigh = false;
486 for (unsigned int c = 0, nc = neigh->n_children();
487 !found_neigh && c != nc; ++c)
488 {
489 Elem * child = neigh->child_ptr(c);
490 if (child == remote_elem)
491 continue;
492 for (auto ncn : child->neighbor_ptr_range())
493 {
494 if (ncn != remote_elem &&
495 ncn->is_ancestor_of(current_elem))
496 {
497 neigh = ncn;
498 found_neigh = true;
499 break;
500 }
501 }
502 }
503 if (!found_neigh)
504 neigh = const_cast<RemoteElem *>(remote_elem);
505 }
506 }
507 current_elem->set_neighbor(s, neigh);
508 #ifdef DEBUG
509 if (neigh != nullptr && neigh != remote_elem)
510 // We ignore subactive elements here because
511 // we don't care about neighbors of subactive element.
512 if ((!neigh->active()) && (!current_elem->subactive()))
513 {
514 libMesh::err << "On processor " << this->processor_id()
515 << std::endl;
516 libMesh::err << "Bad element ID = " << current_elem->id()
517 << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
518 libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
519 << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
520 libMesh::err << "Bad element size = " << current_elem->hmin()
521 << ", Bad neighbor size = " << neigh->hmin() << std::endl;
522 libMesh::err << "Bad element center = " << current_elem->centroid()
523 << ", Bad neighbor center = " << neigh->centroid() << std::endl;
524 libMesh::err << "ERROR: "
525 << (current_elem->active()?"Active":"Ancestor")
526 << " Element at level "
527 << current_elem->level() << std::endl;
528 libMesh::err << "with "
529 << (parent->active()?"active":
530 (parent->subactive()?"subactive":"ancestor"))
531 << " parent share "
532 << (neigh->subactive()?"subactive":"ancestor")
533 << " neighbor at level " << neigh->level()
534 << std::endl;
535 NameBasedIO(*this).write ("bad_mesh.gmv");
536 libmesh_error_msg("Problematic mesh written to bad_mesh.gmv.");
537 }
538 #endif // DEBUG
539 }
540 }
541
542 // We can skip to the next element if we're full-dimension
543 // and therefore don't have any interior parents
544 if (current_elem->dim() >= LIBMESH_DIM)
545 continue;
546
547 // We have no interior parents unless we can find one later
548 current_elem->set_interior_parent(nullptr);
549
550 Elem * pip = parent->interior_parent();
551
552 if (!pip)
553 continue;
554
555 // If there's no interior_parent children, whether due to a
556 // remote element or a non-conformity, then there's no
557 // children to search.
558 if (pip == remote_elem || pip->active())
559 {
560 current_elem->set_interior_parent(pip);
561 continue;
562 }
563
564 // For node comparisons we'll need a sensible tolerance
565 Real node_tolerance = current_elem->hmin() * TOLERANCE;
566
567 // Otherwise our interior_parent should be a child of our
568 // parent's interior_parent.
569 for (auto & child : pip->child_ref_range())
570 {
571 // If we have a remote_elem, that might be our
572 // interior_parent. We'll set it provisionally now and
573 // keep trying to find something better.
574 if (&child == remote_elem)
575 {
576 current_elem->set_interior_parent
577 (const_cast<RemoteElem *>(remote_elem));
578 continue;
579 }
580
581 bool child_contains_our_nodes = true;
582 for (auto & n : current_elem->node_ref_range())
583 {
584 bool child_contains_this_node = false;
585 for (auto & cn : child.node_ref_range())
586 if (cn.absolute_fuzzy_equals
587 (n, node_tolerance))
588 {
589 child_contains_this_node = true;
590 break;
591 }
592 if (!child_contains_this_node)
593 {
594 child_contains_our_nodes = false;
595 break;
596 }
597 }
598 if (child_contains_our_nodes)
599 {
600 current_elem->set_interior_parent(&child);
601 break;
602 }
603 }
604
605 // We should have found *some* interior_parent at this
606 // point, whether semilocal or remote.
607 libmesh_assert(current_elem->interior_parent());
608 }
609 }
610
611 #endif // AMR
612
613
614 #ifdef DEBUG
615 MeshTools::libmesh_assert_valid_neighbors(*this,
616 !reset_remote_elements);
617 MeshTools::libmesh_assert_valid_amr_interior_parents(*this);
618 #endif
619 }
620
621
622
read(const std::string & name,void *,bool skip_renumber_nodes_and_elements,bool skip_find_neighbors)623 void UnstructuredMesh::read (const std::string & name,
624 void *,
625 bool skip_renumber_nodes_and_elements,
626 bool skip_find_neighbors)
627 {
628 // Set the skip_renumber_nodes_and_elements flag on all processors
629 // if necessary.
630 // This ensures that renumber_nodes_and_elements is *not* called
631 // during prepare_for_use() for certain types of mesh files.
632 // This is required in cases where there is an associated solution
633 // file which expects a certain ordering of the nodes.
634 if (name.rfind(".gmv") + 4 == name.size())
635 this->allow_renumbering(false);
636
637 NameBasedIO(*this).read(name);
638
639 if (skip_renumber_nodes_and_elements)
640 {
641 // Use MeshBase::allow_renumbering() yourself instead.
642 libmesh_deprecated();
643 this->allow_renumbering(false);
644 }
645
646 // Done reading the mesh. Now prepare it for use.
647 const bool old_allow_find_neighbors = this->allow_find_neighbors();
648 this->allow_find_neighbors(!skip_find_neighbors);
649 this->prepare_for_use();
650 this->allow_find_neighbors(old_allow_find_neighbors);
651 }
652
653
654
write(const std::string & name)655 void UnstructuredMesh::write (const std::string & name)
656 {
657 LOG_SCOPE("write()", "Mesh");
658
659 NameBasedIO(*this).write(name);
660 }
661
662
663
write(const std::string & name,const std::vector<Number> & v,const std::vector<std::string> & vn)664 void UnstructuredMesh::write (const std::string & name,
665 const std::vector<Number> & v,
666 const std::vector<std::string> & vn)
667 {
668 LOG_SCOPE("write()", "Mesh");
669
670 NameBasedIO(*this).write_nodal_data(name, v, vn);
671 }
672
673
674
675
676
create_pid_mesh(UnstructuredMesh & pid_mesh,const processor_id_type pid)677 void UnstructuredMesh::create_pid_mesh(UnstructuredMesh & pid_mesh,
678 const processor_id_type pid) const
679 {
680
681 // Issue a warning if the number the number of processors
682 // currently available is less that that requested for
683 // partitioning. This is not necessarily an error since
684 // you may run on one processor and still partition the
685 // mesh into several partitions.
686 #ifdef DEBUG
687 if (this->n_processors() < pid)
688 {
689 libMesh::out << "WARNING: You are creating a "
690 << "mesh for a processor id (="
691 << pid
692 << ") greater than "
693 << "the number of processors available for "
694 << "the calculation. (="
695 << this->n_processors()
696 << ")."
697 << std::endl;
698 }
699 #endif
700
701 this->create_submesh (pid_mesh,
702 this->active_pid_elements_begin(pid),
703 this->active_pid_elements_end(pid));
704 }
705
706
707
708
709
710
711
create_submesh(UnstructuredMesh & new_mesh,const const_element_iterator & it,const const_element_iterator & it_end)712 void UnstructuredMesh::create_submesh (UnstructuredMesh & new_mesh,
713 const const_element_iterator & it,
714 const const_element_iterator & it_end) const
715 {
716 // Just in case the subdomain_mesh already has some information
717 // in it, get rid of it.
718 new_mesh.clear();
719
720 // If we're not serial, our submesh isn't either.
721 // There are no remote elements to delete on an empty mesh, but
722 // calling the method to do so marks the mesh as parallel.
723 if (!this->is_serial())
724 new_mesh.delete_remote_elements();
725
726 // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
727 // This may happen if the user accidentally passes the original mesh into
728 // this function! We will check this by making sure we did not just
729 // clear ourself.
730 libmesh_assert_not_equal_to (this->n_nodes(), 0);
731 libmesh_assert_not_equal_to (this->n_elem(), 0);
732
733 // Container to catch boundary IDs handed back by BoundaryInfo
734 std::vector<boundary_id_type> bc_ids;
735
736 // Put any extra integers on the new mesh too
737 new_mesh.merge_extra_integer_names(*this);
738 const unsigned int n_node_ints = _node_integer_names.size(),
739 n_elem_ints = _elem_integer_names.size();
740
741 for (const auto & old_elem : as_range(it, it_end))
742 {
743 // Add an equivalent element type to the new_mesh.
744 // Copy ids for this element.
745 auto uelem = Elem::build(old_elem->type());
746 uelem->set_id() = old_elem->id();
747 #ifdef LIBMESH_ENABLE_UNIQUE_ID
748 uelem->set_unique_id(old_elem->unique_id());
749 #endif
750 uelem->subdomain_id() = old_elem->subdomain_id();
751 uelem->processor_id() = old_elem->processor_id();
752
753 uelem->add_extra_integers(n_elem_ints);
754 for (unsigned int i = 0; i != n_elem_ints; ++i)
755 uelem->set_extra_integer(i, old_elem->get_extra_integer(i));
756
757 uelem->set_mapping_type(old_elem->mapping_type());
758 uelem->set_mapping_data(old_elem->mapping_data());
759
760 Elem * new_elem = new_mesh.add_elem(std::move(uelem));
761
762 libmesh_assert(new_elem);
763
764 // Loop over the nodes on this element.
765 for (auto n : old_elem->node_index_range())
766 {
767 const dof_id_type this_node_id = old_elem->node_id(n);
768
769 // Add this node to the new mesh if it's not there already
770 if (!new_mesh.query_node_ptr(this_node_id))
771 {
772 Node * newn =
773 new_mesh.add_point (old_elem->point(n),
774 this_node_id,
775 old_elem->node_ptr(n)->processor_id());
776
777 newn->add_extra_integers(n_node_ints);
778 for (unsigned int i = 0; i != n_node_ints; ++i)
779 newn->set_extra_integer(i, old_elem->node_ptr(n)->get_extra_integer(i));
780
781 #ifdef LIBMESH_ENABLE_UNIQUE_ID
782 newn->set_unique_id(old_elem->node_ptr(n)->unique_id());
783 #endif
784 }
785
786 // Define this element's connectivity on the new mesh
787 new_elem->set_node(n) = new_mesh.node_ptr(this_node_id);
788 }
789
790 // Maybe add boundary conditions for this element
791 for (auto s : old_elem->side_index_range())
792 {
793 this->get_boundary_info().boundary_ids(old_elem, s, bc_ids);
794 new_mesh.get_boundary_info().add_side (new_elem, s, bc_ids);
795 }
796 } // end loop over elements
797
798 // Prepare the new_mesh for use
799 new_mesh.prepare_for_use();
800 }
801
802
803
804 #ifdef LIBMESH_ENABLE_AMR
contract()805 bool UnstructuredMesh::contract ()
806 {
807 LOG_SCOPE ("contract()", "Mesh");
808
809 // Flag indicating if this call actually changes the mesh
810 bool mesh_changed = false;
811
812 #ifdef DEBUG
813 for (const auto & elem : this->element_ptr_range())
814 libmesh_assert(elem->active() || elem->subactive() || elem->ancestor());
815 #endif
816
817 // Loop over the elements.
818 for (auto & elem : this->element_ptr_range())
819 {
820 // Delete all the subactive ones
821 if (elem->subactive())
822 {
823 // No level-0 element should be subactive.
824 // Note that we CAN'T test elem->level(), as that
825 // touches elem->parent()->dim(), and elem->parent()
826 // might have already been deleted!
827 libmesh_assert(elem->parent());
828
829 // Delete the element
830 // This just sets a pointer to nullptr, and doesn't
831 // invalidate any iterators
832 this->delete_elem(elem);
833
834 // the mesh has certainly changed
835 mesh_changed = true;
836 }
837 else
838 {
839 // Compress all the active ones
840 if (elem->active())
841 elem->contract();
842 else
843 libmesh_assert (elem->ancestor());
844 }
845 }
846
847 // Strip any newly-created nullptr voids out of the element array
848 this->renumber_nodes_and_elements();
849
850 // FIXME: Need to understand why deleting subactive children
851 // invalidates the point locator. For now we will clear it explicitly
852 this->clear_point_locator();
853
854 // Allow our GhostingFunctor objects to reinit if necessary.
855 for (auto & gf : as_range(this->ghosting_functors_begin(),
856 this->ghosting_functors_end()))
857 {
858 libmesh_assert(gf);
859 gf->mesh_reinit();
860 }
861
862 return mesh_changed;
863 }
864 #endif // #ifdef LIBMESH_ENABLE_AMR
865
866
867
all_first_order()868 void UnstructuredMesh::all_first_order ()
869 {
870 /*
871 * when the mesh is not prepared,
872 * at least renumber the nodes and
873 * elements, so that the node ids
874 * are correct
875 */
876 if (!this->_is_prepared)
877 this->renumber_nodes_and_elements ();
878
879 START_LOG("all_first_order()", "Mesh");
880
881 /**
882 * Prepare to identify (and then delete) a bunch of no-longer-used nodes.
883 */
884 std::vector<bool> node_touched_by_me(this->max_node_id(), false);
885
886 // Loop over the high-ordered elements.
887 // First make sure they _are_ indeed high-order, and then replace
888 // them with an equivalent first-order element.
889 for (auto & so_elem : element_ptr_range())
890 {
891 libmesh_assert(so_elem);
892
893 /*
894 * build the first-order equivalent, add to
895 * the new_elements list.
896 */
897 auto lo_elem = Elem::build
898 (Elem::first_order_equivalent_type
899 (so_elem->type()), so_elem->parent());
900
901 const unsigned short n_sides = so_elem->n_sides();
902
903 for (unsigned short s=0; s != n_sides; ++s)
904 if (so_elem->neighbor_ptr(s) == remote_elem)
905 lo_elem->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
906
907 #ifdef LIBMESH_ENABLE_AMR
908 /*
909 * Reset the parent links of any child elements
910 */
911 if (so_elem->has_children())
912 for (unsigned int c = 0, nc = so_elem->n_children(); c != nc; ++c)
913 {
914 Elem * child = so_elem->child_ptr(c);
915 if (child != remote_elem)
916 child->set_parent(lo_elem.get());
917 lo_elem->add_child(child, c);
918 }
919
920 /*
921 * Reset the child link of any parent element
922 */
923 if (so_elem->parent())
924 {
925 unsigned int c =
926 so_elem->parent()->which_child_am_i(so_elem);
927 lo_elem->parent()->replace_child(lo_elem.get(), c);
928 }
929
930 /*
931 * Copy as much data to the new element as makes sense
932 */
933 lo_elem->set_p_level(so_elem->p_level());
934 lo_elem->set_refinement_flag(so_elem->refinement_flag());
935 lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
936 #endif
937
938 libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
939
940 /*
941 * By definition the vertices of the linear and
942 * second order element are identically numbered.
943 * transfer these.
944 */
945 for (unsigned int v=0, snv=so_elem->n_vertices(); v < snv; v++)
946 {
947 lo_elem->set_node(v) = so_elem->node_ptr(v);
948 node_touched_by_me[lo_elem->node_id(v)] = true;
949 }
950
951 /*
952 * find_neighbors relies on remote_elem neighbor links being
953 * properly maintained.
954 */
955 for (unsigned short s=0; s != n_sides; s++)
956 {
957 if (so_elem->neighbor_ptr(s) == remote_elem)
958 lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
959 }
960
961 /**
962 * If the second order element had any boundary conditions they
963 * should be transferred to the first-order element. The old
964 * boundary conditions will be removed from the BoundaryInfo
965 * data structure by insert_elem.
966 */
967 this->get_boundary_info().copy_boundary_ids
968 (this->get_boundary_info(), so_elem, lo_elem.get());
969
970 /*
971 * The new first-order element is ready.
972 * Inserting it into the mesh will replace and delete
973 * the second-order element.
974 */
975 lo_elem->set_id(so_elem->id());
976 #ifdef LIBMESH_ENABLE_UNIQUE_ID
977 lo_elem->set_unique_id(so_elem->unique_id());
978 #endif
979 lo_elem->processor_id() = so_elem->processor_id();
980 lo_elem->subdomain_id() = so_elem->subdomain_id();
981
982 const unsigned int nei = so_elem->n_extra_integers();
983 lo_elem->add_extra_integers(nei);
984 for (unsigned int i=0; i != nei; ++i)
985 lo_elem->set_extra_integer(i, so_elem->get_extra_integer(i));
986
987 // This is probably moot but shouldn't hurt
988 lo_elem->set_mapping_type(so_elem->mapping_type());
989 lo_elem->set_mapping_data(so_elem->mapping_data());
990
991 this->insert_elem(std::move(lo_elem));
992 }
993
994 // Deleting nodes does not invalidate iterators, so this is safe.
995 for (const auto & node : this->node_ptr_range())
996 if (!node_touched_by_me[node->id()])
997 this->delete_node(node);
998
999 // If crazy people applied boundary info to non-vertices and then
1000 // deleted those non-vertices, we should make sure their boundary id
1001 // caches are correct.
1002 this->get_boundary_info().regenerate_id_sets();
1003
1004 STOP_LOG("all_first_order()", "Mesh");
1005
1006 // On hanging nodes that used to also be second order nodes, we
1007 // might now have an invalid nodal processor_id()
1008 Partitioner::set_node_processor_ids(*this);
1009
1010 // delete or renumber nodes if desired
1011 this->prepare_for_use();
1012 }
1013
1014
1015
all_second_order(const bool full_ordered)1016 void UnstructuredMesh::all_second_order (const bool full_ordered)
1017 {
1018 // This function must be run on all processors at once
1019 parallel_object_only();
1020
1021 /*
1022 * when the mesh is not prepared,
1023 * at least renumber the nodes and
1024 * elements, so that the node ids
1025 * are correct
1026 */
1027 if (!this->_is_prepared)
1028 this->renumber_nodes_and_elements ();
1029
1030 /*
1031 * If the mesh is empty
1032 * then we have nothing to do
1033 */
1034 if (!this->n_elem())
1035 return;
1036
1037 /*
1038 * If the mesh is already second order
1039 * then we have nothing to do.
1040 * We have to test for this in a round-about way to avoid
1041 * a bug on distributed parallel meshes with more processors
1042 * than elements.
1043 */
1044 bool already_second_order = false;
1045 if (this->elements_begin() != this->elements_end() &&
1046 (*(this->elements_begin()))->default_order() != FIRST)
1047 already_second_order = true;
1048 this->comm().max(already_second_order);
1049 if (already_second_order)
1050 return;
1051
1052 START_LOG("all_second_order()", "Mesh");
1053
1054 /*
1055 * this map helps in identifying second order
1056 * nodes. Namely, a second-order node:
1057 * - edge node
1058 * - face node
1059 * - bubble node
1060 * is uniquely defined through a set of adjacent
1061 * vertices. This set of adjacent vertices is
1062 * used to identify already added higher-order
1063 * nodes. We are safe to use node id's since we
1064 * make sure that these are correctly numbered.
1065 */
1066 std::map<std::vector<dof_id_type>, Node *> adj_vertices_to_so_nodes;
1067
1068 /*
1069 * The maximum number of new second order nodes we might be adding,
1070 * for use when picking unique unique_id values later
1071 */
1072 unsigned int max_new_nodes_per_elem;
1073
1074 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1075 unique_id_type max_unique_id = this->parallel_max_unique_id();
1076 #endif
1077
1078 /*
1079 * for speed-up of the \p add_point() method, we
1080 * can reserve memory. Guess the number of additional
1081 * nodes for different dimensions
1082 */
1083 switch (this->mesh_dimension())
1084 {
1085 case 1:
1086 /*
1087 * in 1D, there can only be order-increase from Edge2
1088 * to Edge3. Something like 1/2 of n_nodes() have
1089 * to be added
1090 */
1091 max_new_nodes_per_elem = 3 - 2;
1092 this->reserve_nodes(static_cast<unsigned int>
1093 (1.5*static_cast<double>(this->n_nodes())));
1094 break;
1095
1096 case 2:
1097 /*
1098 * in 2D, either refine from Tri3 to Tri6 (double the nodes)
1099 * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
1100 */
1101 max_new_nodes_per_elem = 9 - 4;
1102 this->reserve_nodes(static_cast<unsigned int>
1103 (2*static_cast<double>(this->n_nodes())));
1104 break;
1105
1106
1107 case 3:
1108 /*
1109 * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
1110 * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
1111 * quite some nodes, and since we do not want to overburden the memory by
1112 * a too conservative guess, use the lower bound
1113 */
1114 max_new_nodes_per_elem = 27 - 8;
1115 this->reserve_nodes(static_cast<unsigned int>
1116 (2.5*static_cast<double>(this->n_nodes())));
1117 break;
1118
1119 default:
1120 // Hm?
1121 libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1122 }
1123
1124
1125
1126 /*
1127 * form a vector that will hold the node id's of
1128 * the vertices that are adjacent to the son-th
1129 * second-order node. Pull this outside of the
1130 * loop so that silly compilers don't repeatedly
1131 * create and destroy the vector.
1132 */
1133 std::vector<dof_id_type> adjacent_vertices_ids;
1134
1135 /**
1136 * On distributed meshes we currently only support unpartitioned
1137 * meshes (where we'll add every node in sync) or
1138 * completely-partitioned meshes (where we'll sync nodes later);
1139 * let's keep track to make sure we're not in any in-between state.
1140 */
1141 dof_id_type n_unpartitioned_elem = 0,
1142 n_partitioned_elem = 0;
1143 const processor_id_type my_pid = this->processor_id();
1144
1145 /**
1146 * Loop over the low-ordered elements in the _elements vector.
1147 * First make sure they _are_ indeed low-order, and then replace
1148 * them with an equivalent second-order element. Don't
1149 * forget to delete the low-order element, or else it will leak!
1150 */
1151 for (auto & lo_elem : element_ptr_range())
1152 {
1153 // make sure it is linear order
1154 libmesh_error_msg_if(lo_elem->default_order() != FIRST,
1155 "ERROR: This is not a linear element: type=" << lo_elem->type());
1156
1157 // this does _not_ work for refined elements
1158 libmesh_assert_equal_to (lo_elem->level (), 0);
1159
1160 const processor_id_type lo_pid = lo_elem->processor_id();
1161
1162 if (lo_pid == DofObject::invalid_processor_id)
1163 ++n_unpartitioned_elem;
1164 else
1165 ++n_partitioned_elem;
1166
1167 /*
1168 * build the second-order equivalent, add to
1169 * the new_elements list. Note that this here
1170 * is the only point where \p full_ordered
1171 * is necessary. The remaining code works well
1172 * for either type of second-order equivalent, e.g.
1173 * Hex20 or Hex27, as equivalents for Hex8
1174 */
1175 auto so_elem =
1176 Elem::build (Elem::second_order_equivalent_type(lo_elem->type(),
1177 full_ordered));
1178
1179 libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
1180
1181
1182 /*
1183 * By definition the vertices of the linear and
1184 * second order element are identically numbered.
1185 * transfer these.
1186 */
1187 for (unsigned int v=0, lnv=lo_elem->n_vertices(); v < lnv; v++)
1188 so_elem->set_node(v) = lo_elem->node_ptr(v);
1189
1190 /*
1191 * Now handle the additional mid-side nodes. This
1192 * is simply handled through a map that remembers
1193 * the already-added nodes. This map maps the global
1194 * ids of the vertices (that uniquely define this
1195 * higher-order node) to the new node.
1196 * Notation: son = second-order node
1197 */
1198 const unsigned int son_begin = so_elem->n_vertices();
1199 const unsigned int son_end = so_elem->n_nodes();
1200
1201 for (unsigned int son=son_begin; son<son_end; son++)
1202 {
1203 const unsigned int n_adjacent_vertices =
1204 so_elem->n_second_order_adjacent_vertices(son);
1205
1206 adjacent_vertices_ids.resize(n_adjacent_vertices);
1207
1208 for (unsigned int v=0; v<n_adjacent_vertices; v++)
1209 adjacent_vertices_ids[v] =
1210 so_elem->node_id( so_elem->second_order_adjacent_vertex(son,v) );
1211
1212 /*
1213 * \p adjacent_vertices_ids is now in order of the current
1214 * side. sort it, so that comparisons with the
1215 * \p adjacent_vertices_ids created through other elements'
1216 * sides can match
1217 */
1218 std::sort(adjacent_vertices_ids.begin(),
1219 adjacent_vertices_ids.end());
1220
1221
1222 // does this set of vertices already have a mid-node added?
1223 auto pos = adj_vertices_to_so_nodes.equal_range (adjacent_vertices_ids);
1224
1225 // no, not added yet
1226 if (pos.first == pos.second)
1227 {
1228 /*
1229 * for this set of vertices, there is no
1230 * second_order node yet. Add it.
1231 *
1232 * compute the location of the new node as
1233 * the average over the adjacent vertices.
1234 */
1235 Point new_location = this->point(adjacent_vertices_ids[0]);
1236 for (unsigned int v=1; v<n_adjacent_vertices; v++)
1237 new_location += this->point(adjacent_vertices_ids[v]);
1238
1239 new_location /= static_cast<Real>(n_adjacent_vertices);
1240
1241 /* Add the new point to the mesh.
1242 *
1243 * If we are on a serialized mesh, then we're doing this
1244 * all in sync, and the node processor_id will be
1245 * consistent between processors.
1246 *
1247 * If we are on a distributed mesh, we can fix
1248 * inconsistent processor ids later, but only if every
1249 * processor gives new nodes a *locally* consistent
1250 * processor id, so we'll give the new node the
1251 * processor id of an adjacent element for now and then
1252 * we'll update that later if appropriate.
1253 */
1254 Node * so_node = this->add_point
1255 (new_location, DofObject::invalid_id, lo_pid);
1256
1257 /* Come up with a unique unique_id for a potentially new
1258 * node. On a distributed mesh we don't yet know what
1259 * processor_id will definitely own it, so we can't let
1260 * the pid determine the unique_id. But we're not
1261 * adding unpartitioned nodes in sync, so we can't let
1262 * the mesh autodetermine a unique_id for a new
1263 * unpartitioned node either. So we have to pick unique
1264 * unique_id values manually.
1265 *
1266 * We don't have to pick the *same* unique_id value as
1267 * will be picked on other processors, though; we'll
1268 * sync up each node later. We just need to make sure
1269 * we don't duplicate any unique_id that might be chosen
1270 * by the same process elsewhere.
1271 */
1272 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1273 unique_id_type new_unique_id = max_unique_id +
1274 max_new_nodes_per_elem * lo_elem->id() +
1275 son - son_begin;
1276
1277 so_node->set_unique_id(new_unique_id);
1278 #endif
1279 libmesh_ignore(max_new_nodes_per_elem);
1280
1281 /*
1282 * insert the new node with its defining vertex
1283 * set into the map, and relocate pos to this
1284 * new entry, so that the so_elem can use
1285 * \p pos for inserting the node
1286 */
1287 adj_vertices_to_so_nodes.emplace_hint
1288 (pos.first, adjacent_vertices_ids, so_node);
1289
1290 so_elem->set_node(son) = so_node;
1291 }
1292 // yes, already added.
1293 else
1294 {
1295 Node * so_node = pos.first->second;
1296 libmesh_assert(so_node);
1297
1298 so_elem->set_node(son) = so_node;
1299
1300 // We need to ensure that the processor who should own a
1301 // node *knows* they own the node. And because
1302 // Node::choose_processor_id() may depend on Node id,
1303 // which may not yet be authoritative, we still have to
1304 // use a dumb-but-id-independent partitioning heuristic.
1305 processor_id_type chosen_pid =
1306 std::min (so_node->processor_id(), lo_pid);
1307
1308 // Plus, if we just discovered that we own this node,
1309 // then on a distributed mesh we need to make sure to
1310 // give it a valid id, not just a placeholder id!
1311 if (!this->is_replicated() &&
1312 so_node->processor_id() != my_pid &&
1313 chosen_pid == my_pid)
1314 this->own_node(*so_node);
1315
1316 so_node->processor_id() = chosen_pid;
1317 }
1318 }
1319
1320 /*
1321 * find_neighbors relies on remote_elem neighbor links being
1322 * properly maintained.
1323 */
1324 for (auto s : lo_elem->side_index_range())
1325 {
1326 if (lo_elem->neighbor_ptr(s) == remote_elem)
1327 so_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1328 }
1329
1330 /**
1331 * If the linear element had any boundary conditions they
1332 * should be transferred to the second-order element. The old
1333 * boundary conditions will be removed from the BoundaryInfo
1334 * data structure by insert_elem.
1335 *
1336 * Also, prepare_for_use() will reconstruct most of our neighbor
1337 * links, but if we have any remote_elem links in a distributed
1338 * mesh, they need to be preserved. We do that in the same loop
1339 * here.
1340 */
1341 this->get_boundary_info().copy_boundary_ids
1342 (this->get_boundary_info(), lo_elem, so_elem.get());
1343
1344 /*
1345 * The new second-order element is ready.
1346 * Inserting it into the mesh will replace and delete
1347 * the first-order element.
1348 */
1349 so_elem->set_id(lo_elem->id());
1350 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1351 so_elem->set_unique_id(lo_elem->unique_id());
1352 #endif
1353 so_elem->processor_id() = lo_pid;
1354 so_elem->subdomain_id() = lo_elem->subdomain_id();
1355
1356 const unsigned int nei = so_elem->n_extra_integers();
1357 so_elem->add_extra_integers(nei);
1358 for (unsigned int i=0; i != nei; ++i)
1359 so_elem->set_extra_integer(i, lo_elem->get_extra_integer(i));
1360
1361 // This might not help anything but shouldn't hurt.
1362 so_elem->set_mapping_type(lo_elem->mapping_type());
1363 so_elem->set_mapping_data(lo_elem->mapping_data());
1364
1365 this->insert_elem(std::move(so_elem));
1366 }
1367
1368 // we can clear the map
1369 adj_vertices_to_so_nodes.clear();
1370
1371 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1372 const unique_id_type new_max_unique_id = max_unique_id +
1373 max_new_nodes_per_elem * this->n_elem();
1374 this->set_next_unique_id(new_max_unique_id);
1375 #endif
1376
1377
1378
1379 STOP_LOG("all_second_order()", "Mesh");
1380
1381 // On a DistributedMesh our ghost node processor ids may be bad,
1382 // the ids of nodes touching remote elements may be inconsistent,
1383 // unique_ids of newly added non-local nodes remain unset, and our
1384 // partitioning of new nodes may not be well balanced.
1385 //
1386 // make_nodes_parallel_consistent() will fix all this.
1387 if (!this->is_replicated())
1388 {
1389 dof_id_type max_unpartitioned_elem = n_unpartitioned_elem;
1390 this->comm().max(max_unpartitioned_elem);
1391 if (max_unpartitioned_elem)
1392 {
1393 // We'd better be effectively serialized here. In theory we
1394 // could support more complicated cases but in practice we
1395 // only support "completely partitioned" and/or "serialized"
1396 if (!this->comm().verify(n_unpartitioned_elem) ||
1397 !this->comm().verify(n_partitioned_elem) ||
1398 !this->is_serial())
1399 libmesh_not_implemented();
1400 }
1401 else
1402 {
1403 MeshCommunication().make_nodes_parallel_consistent (*this);
1404 }
1405 }
1406
1407 // renumber nodes, elements etc
1408 this->prepare_for_use();
1409 }
1410
1411 } // namespace libMesh
1412