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 : // Local includes
21 : #include "libmesh/boundary_info.h"
22 : #include "libmesh/ghosting_functor.h"
23 : #include "libmesh/ghost_point_neighbors.h"
24 : #include "libmesh/unstructured_mesh.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/elem_range.h"
28 : #include "libmesh/face_polygon.h"
29 : #include "libmesh/mesh_tools.h" // For n_levels
30 : #include "libmesh/parallel.h"
31 : #include "libmesh/remote_elem.h"
32 : #include "libmesh/namebased_io.h"
33 : #include "libmesh/partitioner.h"
34 : #include "libmesh/enum_order.h"
35 : #include "libmesh/mesh_communication.h"
36 : #include "libmesh/enum_to_string.h"
37 : #include "libmesh/mesh_serializer.h"
38 : #include "libmesh/utility.h"
39 :
40 : #ifdef LIBMESH_HAVE_NANOFLANN
41 : #include "libmesh/nanoflann.hpp"
42 : #endif
43 :
44 : // C++ includes
45 : #include <algorithm> // std::all_of
46 : #include <atomic>
47 : #include <fstream>
48 : #include <iomanip>
49 : #include <map>
50 : #include <sstream>
51 : #include <unordered_map>
52 :
53 : // for disjoint neighbors
54 : #include "libmesh/periodic_boundaries.h"
55 : #include "libmesh/periodic_boundary.h"
56 :
57 : namespace {
58 :
59 : using namespace libMesh;
60 :
61 : // Helper functions for all_second_order, all_complete_order
62 :
63 : std::map<std::vector<dof_id_type>, Node *>::iterator
64 35075642 : map_hi_order_node(unsigned int hon,
65 : const Elem & hi_elem,
66 : std::map<std::vector<dof_id_type>, Node *> & adj_vertices_to_ho_nodes)
67 : {
68 : /*
69 : * form a vector that will hold the node id's of
70 : * the vertices that are adjacent to the nth
71 : * higher-order node.
72 : */
73 : const unsigned int n_adjacent_vertices =
74 35075642 : hi_elem.n_second_order_adjacent_vertices(hon);
75 :
76 35075642 : std::vector<dof_id_type> adjacent_vertices_ids(n_adjacent_vertices);
77 :
78 117884639 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
79 85012075 : adjacent_vertices_ids[v] =
80 82808997 : hi_elem.node_id( hi_elem.second_order_adjacent_vertex(hon,v) );
81 :
82 : /*
83 : * \p adjacent_vertices_ids is now in order of the current
84 : * side. sort it, so that comparisons with the
85 : * \p adjacent_vertices_ids created through other elements'
86 : * sides can match
87 : */
88 35075642 : std::sort(adjacent_vertices_ids.begin(),
89 : adjacent_vertices_ids.end());
90 :
91 : // Does this set of vertices already have a mid-node added? If not
92 : // we'll want to add it.
93 36924666 : return adj_vertices_to_ho_nodes.try_emplace(adjacent_vertices_ids, nullptr).first;
94 : }
95 :
96 3739233 : void transfer_elem(Elem & lo_elem,
97 : std::unique_ptr<Elem> hi_elem,
98 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
99 : unique_id_type max_unique_id,
100 : unique_id_type max_new_nodes_per_elem,
101 : #endif
102 : UnstructuredMesh & mesh,
103 : std::map<std::vector<dof_id_type>, Node *> & adj_vertices_to_ho_nodes,
104 : std::unordered_map<Elem *, std::vector<Elem *>> & exterior_children_of)
105 : {
106 97314 : libmesh_assert_equal_to (lo_elem.n_vertices(), hi_elem->n_vertices());
107 :
108 194628 : const processor_id_type my_pid = mesh.processor_id();
109 3739233 : const processor_id_type lo_pid = lo_elem.processor_id();
110 :
111 : /*
112 : * Now handle the additional higher-order nodes. This
113 : * is simply handled through a map that remembers
114 : * the already-added nodes. This map maps the global
115 : * ids of the vertices (that uniquely define this
116 : * higher-order node) to the new node.
117 : * Notation: hon = high-order node
118 : */
119 3739233 : const unsigned int hon_begin = lo_elem.n_nodes();
120 3739233 : const unsigned int hon_end = hi_elem->n_nodes();
121 :
122 97314 : libmesh_assert_less (hon_begin, hon_end);
123 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
124 97314 : libmesh_assert_less_equal
125 : (hon_end-hon_begin, max_new_nodes_per_elem);
126 : #endif
127 :
128 38529312 : for (unsigned int hon=hon_begin; hon<hon_end; hon++)
129 : {
130 34790079 : auto pos = map_hi_order_node(hon, *hi_elem, adj_vertices_to_ho_nodes);
131 :
132 : // no, not added yet
133 34790079 : if (!pos->second)
134 : {
135 328828 : const auto & adjacent_vertices_ids = pos->first;
136 :
137 : /*
138 : * for this set of vertices, there is no
139 : * second_order node yet. Add it.
140 : *
141 : * compute the location of the new node as
142 : * the average over the adjacent vertices.
143 : */
144 328828 : Point new_location = 0;
145 44475304 : for (dof_id_type vertex_id : adjacent_vertices_ids)
146 32004330 : new_location += mesh.point(vertex_id);
147 :
148 12799802 : new_location /= static_cast<Real>(adjacent_vertices_ids.size());
149 :
150 : /* Add the new point to the mesh.
151 : *
152 : * If we are on a serialized mesh, then we're doing this
153 : * all in sync, and the node processor_id will be
154 : * consistent between processors.
155 : *
156 : * If we are on a distributed mesh, we can fix
157 : * inconsistent processor ids later, but only if every
158 : * processor gives new nodes a *locally* consistent
159 : * processor id, so we'll give the new node the
160 : * processor id of an adjacent element for now and then
161 : * we'll update that later if appropriate.
162 : */
163 : Node * hi_node = mesh.add_point
164 12470974 : (new_location, DofObject::invalid_id, lo_pid);
165 :
166 : /* Come up with a unique unique_id for a potentially new
167 : * node. On a distributed mesh we don't yet know what
168 : * processor_id will definitely own it, so we can't let
169 : * the pid determine the unique_id. But we're not
170 : * adding unpartitioned nodes in sync, so we can't let
171 : * the mesh autodetermine a unique_id for a new
172 : * unpartitioned node either. So we have to pick unique
173 : * unique_id values manually.
174 : *
175 : * We don't have to pick the *same* unique_id value as
176 : * will be picked on other processors, though; we'll
177 : * sync up each node later. We just need to make sure
178 : * we don't duplicate any unique_id that might be chosen
179 : * by the same process elsewhere.
180 : */
181 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
182 12470974 : unique_id_type new_unique_id = max_unique_id +
183 12799802 : max_new_nodes_per_elem * lo_elem.id() +
184 12470974 : hon - hon_begin;
185 :
186 328828 : hi_node->set_unique_id(new_unique_id);
187 : #endif
188 :
189 : /*
190 : * insert the new node with its defining vertex
191 : * set into the map, and relocate pos to this
192 : * new entry, so that the hi_elem can use
193 : * \p pos for inserting the node
194 : */
195 12470974 : pos->second = hi_node;
196 :
197 12470974 : hi_elem->set_node(hon, hi_node);
198 : }
199 : // yes, already added.
200 : else
201 : {
202 581922 : Node * hi_node = pos->second;
203 581922 : libmesh_assert(hi_node);
204 581922 : libmesh_assert_equal_to(mesh.node_ptr(hi_node->id()), hi_node);
205 :
206 22319105 : hi_elem->set_node(hon, hi_node);
207 :
208 : // We need to ensure that the processor who should own a
209 : // node *knows* they own the node. And because
210 : // Node::choose_processor_id() may depend on Node id,
211 : // which may not yet be authoritative, we still have to
212 : // use a dumb-but-id-independent partitioning heuristic.
213 : processor_id_type chosen_pid =
214 22319105 : std::min (hi_node->processor_id(), lo_pid);
215 :
216 : // Plus, if we just discovered that we own this node,
217 : // then on a distributed mesh we need to make sure to
218 : // give it a valid id, not just a placeholder id!
219 22319105 : if (!mesh.is_replicated() &&
220 22319105 : hi_node->processor_id() != my_pid &&
221 : chosen_pid == my_pid)
222 4130 : mesh.own_node(*hi_node);
223 :
224 22319105 : hi_node->processor_id() = chosen_pid;
225 : }
226 : }
227 :
228 : /*
229 : * find_neighbors relies on remote_elem neighbor links being
230 : * properly maintained. Our own code here relies on ordinary
231 : * neighbor links being properly maintained, so let's just keep
232 : * everything up to date.
233 : */
234 19142826 : for (auto s : lo_elem.side_index_range())
235 : {
236 15403593 : Elem * neigh = lo_elem.neighbor_ptr(s);
237 15403593 : if (!neigh)
238 14499984 : continue;
239 :
240 527991 : if (neigh != remote_elem)
241 : {
242 : // We don't support AMR even outside our own range yet.
243 30346 : libmesh_assert_equal_to (neigh->level(), 0);
244 :
245 494094 : const unsigned int ns = neigh->which_neighbor_am_i(&lo_elem);
246 30346 : libmesh_assert_not_equal_to(ns, libMesh::invalid_uint);
247 :
248 60692 : neigh->set_neighbor(ns, hi_elem.get());
249 : }
250 :
251 61012 : hi_elem->set_neighbor(s, neigh);
252 : }
253 :
254 : /**
255 : * If the old element has an interior_parent(), transfer it to the
256 : * new element ... and if the interior_parent itself might be
257 : * getting upgraded, make sure we later consider the new element to
258 : * be its exterior child, not the old element.
259 : */
260 3739233 : Elem * interior_p = lo_elem.interior_parent();
261 3739233 : if (interior_p)
262 0 : hi_elem->set_interior_parent(interior_p);
263 :
264 3739233 : if (auto parent_exterior_it = exterior_children_of.find(interior_p);
265 97314 : parent_exterior_it != exterior_children_of.end())
266 : {
267 0 : auto & exteriors = parent_exterior_it->second;
268 0 : for (std::size_t i : index_range(exteriors))
269 0 : if (exteriors[i] == &lo_elem)
270 : {
271 0 : exteriors[i] = hi_elem.get();
272 0 : break;
273 : }
274 : }
275 :
276 : /**
277 : * If we had interior_parent() links to the old element, transfer
278 : * them to the new element.
279 : */
280 3836547 : if (auto exterior_it = exterior_children_of.find(&lo_elem);
281 97314 : exterior_it != exterior_children_of.end())
282 : {
283 3739233 : for (Elem * exterior_elem : exterior_it->second)
284 : {
285 0 : libmesh_assert(exterior_elem->interior_parent() == &lo_elem);
286 0 : exterior_elem->set_interior_parent(hi_elem.get());
287 : }
288 : }
289 :
290 : /**
291 : * If the old element had any boundary conditions they
292 : * should be transferred to the second-order element. The old
293 : * boundary conditions will be removed from the BoundaryInfo
294 : * data structure by insert_elem.
295 : *
296 : * Also, prepare_for_use() will reconstruct most of our neighbor
297 : * links, but if we have any remote_elem links in a distributed
298 : * mesh, they need to be preserved. We do that in the same loop
299 : * here.
300 : */
301 97314 : mesh.get_boundary_info().copy_boundary_ids
302 3739233 : (mesh.get_boundary_info(), &lo_elem, hi_elem.get());
303 :
304 : /*
305 : * The new second-order element is ready.
306 : * Inserting it into the mesh will replace and delete
307 : * the first-order element.
308 : */
309 194628 : hi_elem->set_id(lo_elem.id());
310 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
311 194628 : hi_elem->set_unique_id(lo_elem.unique_id());
312 : #endif
313 :
314 3739233 : const unsigned int nei = lo_elem.n_extra_integers();
315 3739233 : hi_elem->add_extra_integers(nei);
316 3739552 : for (unsigned int i=0; i != nei; ++i)
317 319 : hi_elem->set_extra_integer(i, lo_elem.get_extra_integer(i));
318 :
319 3641919 : hi_elem->inherit_data_from(lo_elem);
320 :
321 3836547 : mesh.insert_elem(std::move(hi_elem));
322 3739233 : }
323 :
324 :
325 : template <typename ElemTypeConverter>
326 : void
327 44767 : all_increased_order_range (UnstructuredMesh & mesh,
328 : const SimpleRange<MeshBase::element_iterator> & range,
329 : const unsigned int max_new_nodes_per_elem,
330 : const ElemTypeConverter & elem_type_converter)
331 : {
332 : // This function must be run on all processors at once
333 1276 : timpi_parallel_only(mesh.comm());
334 :
335 : /*
336 : * The maximum number of new higher-order nodes we might be adding,
337 : * for use when picking unique unique_id values later. This variable
338 : * is not used unless unique ids are enabled, so libmesh_ignore() it
339 : * to avoid warnings in that case.
340 : */
341 1276 : libmesh_ignore(max_new_nodes_per_elem);
342 :
343 : /*
344 : * The mesh should at least be consistent enough for us to add new
345 : * nodes consistently.
346 : */
347 44767 : mesh.update_parallel_id_counts();
348 :
349 : /*
350 : * If the mesh is empty then we have nothing to do
351 : */
352 44767 : if (!mesh.n_elem())
353 3984 : return;
354 :
355 : // If every element in the range _on every proc_ is already of the
356 : // requested higher order then we have nothing to do. However, if
357 : // any proc has some lower-order elements in the range, then _all_
358 : // processors need to continue this function because it is
359 : // parallel_only().
360 : //
361 : // Note: std::all_of() returns true for an empty range, which can
362 : // happen for example in the DistributedMesh case when there are
363 : // more processors than elements. In the case of an empty range we
364 : // therefore set already_second_order to true on that proc.
365 109956 : auto is_higher_order = [&elem_type_converter](const Elem * elem) {
366 42966 : ElemType old_type = elem->type();
367 2194 : ElemType new_type = elem_type_converter(old_type);
368 42966 : return old_type == new_type;
369 : };
370 :
371 44767 : bool already_higher_order =
372 88258 : std::all_of(range.begin(), range.end(), is_higher_order);
373 :
374 : // Check with other processors and possibly return early
375 44767 : mesh.comm().min(already_higher_order);
376 44767 : if (already_higher_order)
377 116 : return;
378 :
379 : /*
380 : * this map helps in identifying higher order
381 : * nodes. Namely, a higher-order node:
382 : * - edge node
383 : * - face node
384 : * - bubble node
385 : * is uniquely defined through a set of adjacent
386 : * vertices. This set of adjacent vertices is
387 : * used to identify already added higher-order
388 : * nodes. We are safe to use node id's since we
389 : * make sure that these are correctly numbered.
390 : *
391 : * We lazily use an ordered map here to avoid having to implement a
392 : * good hash for vector<dof_id_type>
393 : */
394 2320 : std::map<std::vector<dof_id_type>, Node *> adj_vertices_to_ho_nodes;
395 :
396 : /*
397 : * This map helps us reset any interior_parent() values from the
398 : * lower order element to its higher order replacement. Unlike with
399 : * neighbor pointers, we don't have backlinks here, so we have to
400 : * iterate over the mesh to track forward links.
401 : */
402 2320 : std::unordered_map<Elem *, std::vector<Elem *>> exterior_children_of;
403 :
404 : /*
405 : * max_new_nodes_per_elem is the maximum number of new higher order
406 : * nodes we might be adding, for use when picking unique unique_id
407 : * values later. This variable is not used unless unique ids are
408 : * enabled.
409 : */
410 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
411 40783 : unique_id_type max_unique_id = mesh.parallel_max_unique_id();
412 : #endif
413 :
414 : /**
415 : * On distributed meshes we currently only support unpartitioned
416 : * meshes (where we'll add every node in sync) or
417 : * completely-partitioned meshes (where we'll sync nodes later);
418 : * let's keep track to make sure we're not in any in-between state.
419 : */
420 1160 : dof_id_type n_unpartitioned_elem = 0;
421 :
422 : /**
423 : * Loop over the elements in the given range. If any are
424 : * already at higher than first-order, track their higher-order
425 : * nodes in case we need them for neighboring elements later.
426 : *
427 : * In this way we can use this method to "fix up" a mesh which has
428 : * otherwise inconsistent neighbor pairs of lower and higher order
429 : * geometric elements.
430 : *
431 : * If any elements are not at the desired order yet, we need to
432 : * check their neighbors and even their edge neighbors for higher
433 : * order; we may need to share elements with a neighbor not in the
434 : * range.
435 : */
436 7476057 : auto track_if_necessary = [&adj_vertices_to_ho_nodes,
437 : &exterior_children_of,
438 198492 : &elem_type_converter](Elem * elem) {
439 3776586 : if (elem && elem != remote_elem)
440 : {
441 3776586 : if (elem->default_order() != FIRST)
442 305187 : for (unsigned int hon : make_range(elem->n_vertices(), elem->n_nodes()))
443 : {
444 285563 : auto pos = map_hi_order_node(hon, *elem, adj_vertices_to_ho_nodes);
445 299325 : pos->second = elem->node_ptr(hon);
446 : }
447 :
448 3776586 : const ElemType old_type = elem->type();
449 134918 : const ElemType new_type = elem_type_converter(old_type);
450 3776586 : if (old_type != new_type)
451 3756962 : exterior_children_of.emplace(elem, std::vector<Elem *>());
452 : }
453 : };
454 :
455 : // If we're in the common case then just track everything; otherwise
456 : // find point neighbors to track
457 123285 : if (range.begin() == mesh.elements_begin() &&
458 150266 : range.end() == mesh.elements_end())
459 : {
460 7311037 : for (auto & elem : range)
461 3734225 : track_if_necessary(elem);
462 : }
463 : else
464 : {
465 240 : GhostingFunctor::map_type point_neighbor_elements;
466 :
467 240 : GhostPointNeighbors point_neighbor_finder(mesh);
468 11880 : point_neighbor_finder(range.begin(), range.end(),
469 : mesh.n_processors(),
470 : point_neighbor_elements);
471 :
472 46401 : for (auto & [elem, coupling_map] : point_neighbor_elements)
473 : {
474 2168 : libmesh_ignore(coupling_map);
475 42361 : track_if_necessary(const_cast<Elem *>(elem));
476 : }
477 : }
478 :
479 : /**
480 : * Loop over all mesh elements to look for interior_parent links we
481 : * need to upgrade later.
482 : */
483 7570059 : for (auto & elem : mesh.element_ptr_range())
484 3928323 : if (auto exterior_map_it = exterior_children_of.find(elem->interior_parent());
485 101364 : exterior_map_it != exterior_children_of.end())
486 0 : exterior_map_it->second.push_back(elem);
487 :
488 : /**
489 : * Loop over the low-ordered elements in the _elements vector.
490 : * First make sure they _are_ indeed low-order, and then replace
491 : * them with an equivalent second-order element. Don't
492 : * forget to delete the low-order element, or else it will leak!
493 : */
494 7325559 : for (auto & lo_elem : range)
495 : {
496 : // Now we can skip the elements in the range that are already
497 : // higher-order.
498 3739724 : const ElemType old_type = lo_elem->type();
499 132546 : const ElemType new_type = elem_type_converter(old_type);
500 :
501 3739724 : if (old_type == new_type)
502 491 : continue;
503 :
504 : // this does _not_ work for refined elements
505 97314 : libmesh_assert_equal_to (lo_elem->level(), 0);
506 :
507 3739233 : if (lo_elem->processor_id() == DofObject::invalid_processor_id)
508 3636227 : ++n_unpartitioned_elem;
509 :
510 : /*
511 : * Build the higher-order equivalent; add to
512 : * the new_elements list.
513 : */
514 3739233 : auto ho_elem = Elem::build (new_type);
515 :
516 97314 : libmesh_assert_equal_to (lo_elem->n_vertices(), ho_elem->n_vertices());
517 :
518 : /*
519 : * By definition the initial nodes of the lower and higher order
520 : * element are identically numbered. Transfer these.
521 : */
522 19302574 : for (unsigned int v=0, lnn=lo_elem->n_nodes(); v < lnn; v++)
523 15978605 : ho_elem->set_node(v, lo_elem->node_ptr(v));
524 :
525 3933861 : transfer_elem(*lo_elem, std::move(ho_elem),
526 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
527 : max_unique_id, max_new_nodes_per_elem,
528 : #endif
529 : mesh, adj_vertices_to_ho_nodes,
530 : exterior_children_of);
531 : } // end for (auto & lo_elem : range)
532 :
533 : // we can clear the map at this point.
534 1160 : adj_vertices_to_ho_nodes.clear();
535 :
536 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
537 40783 : const unique_id_type new_max_unique_id = max_unique_id +
538 40783 : max_new_nodes_per_elem * mesh.n_elem();
539 40783 : mesh.set_next_unique_id(new_max_unique_id);
540 : #endif
541 :
542 : // On a DistributedMesh our ghost node processor ids may be bad,
543 : // the ids of nodes touching remote elements may be inconsistent,
544 : // unique_ids of newly added non-local nodes remain unset, and our
545 : // partitioning of new nodes may not be well balanced.
546 : //
547 : // make_nodes_parallel_consistent() will fix all this.
548 40783 : if (!mesh.is_replicated())
549 : {
550 35298 : dof_id_type max_unpartitioned_elem = n_unpartitioned_elem;
551 35298 : mesh.comm().max(max_unpartitioned_elem);
552 35298 : if (max_unpartitioned_elem)
553 : {
554 : // We'd better be effectively serialized here. In theory we
555 : // could support more complicated cases but for now we
556 : // only support "completely partitioned" and/or "serialized"
557 27917 : if (mesh.is_serial())
558 22 : libmesh_assert(mesh.comm().verify(n_unpartitioned_elem));
559 : else
560 0 : libmesh_not_implemented();
561 : }
562 : else
563 : {
564 7381 : MeshCommunication().make_nodes_parallel_consistent (mesh);
565 : }
566 : }
567 :
568 : // renumber nodes, repartition nodes, etc. We may no longer need a
569 : // find_neighbors() here since we're keeping neighbor links intact
570 : // ourselves, *except* that if we're not already prepared we may
571 : // have user code that was expecting this call to prepare neighbors.
572 2320 : const bool old_find_neighbors = mesh.allow_find_neighbors();
573 40783 : if (mesh.is_prepared())
574 264 : mesh.allow_find_neighbors(false);
575 40783 : mesh.prepare_for_use();
576 1160 : mesh.allow_find_neighbors(old_find_neighbors);
577 : }
578 :
579 :
580 : } // anonymous namespace
581 :
582 :
583 : namespace libMesh
584 : {
585 :
586 : // This class adapts a vector of Nodes (represented by a pair of a Point and a dof_id_type)
587 : // for use in a nanoflann KD-Tree
588 :
589 0 : class VectorOfNodesAdaptor
590 : {
591 : private:
592 : const std::vector<std::pair<Point, dof_id_type>> _nodes;
593 :
594 : public:
595 0 : VectorOfNodesAdaptor(const std::vector<std::pair<Point, dof_id_type>> & nodes) :
596 0 : _nodes(nodes)
597 0 : {}
598 :
599 : /**
600 : * Must return the number of data points
601 : */
602 0 : inline size_t kdtree_get_point_count() const { return _nodes.size(); }
603 :
604 : /**
605 : * \returns The dim'th component of the idx'th point in the class:
606 : * Since this is inlined and the "dim" argument is typically an immediate value, the
607 : * "if's" are actually solved at compile time.
608 : */
609 0 : inline Real kdtree_get_pt(const size_t idx, int dim) const
610 : {
611 0 : libmesh_assert_less (idx, _nodes.size());
612 0 : libmesh_assert_less (dim, 3);
613 :
614 0 : const Point & p(_nodes[idx].first);
615 :
616 0 : if (dim==0) return p(0);
617 0 : if (dim==1) return p(1);
618 0 : return p(2);
619 : }
620 :
621 : /*
622 : * Optional bounding-box computation
623 : */
624 : template <class BBOX>
625 0 : bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
626 : };
627 :
628 :
629 : // ------------------------------------------------------------
630 : // UnstructuredMesh class member functions
631 361822 : UnstructuredMesh::UnstructuredMesh (const Parallel::Communicator & comm_in,
632 361822 : unsigned char d) :
633 361822 : MeshBase (comm_in,d)
634 : {
635 10768 : libmesh_assert (libMesh::initialized());
636 361822 : }
637 :
638 :
639 :
640 37953 : UnstructuredMesh::UnstructuredMesh (const MeshBase & other_mesh) :
641 37953 : MeshBase (other_mesh)
642 : {
643 14158 : libmesh_assert (libMesh::initialized());
644 37953 : }
645 :
646 :
647 :
648 39728 : void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh,
649 : const bool skip_find_neighbors,
650 : dof_id_type element_id_offset,
651 : dof_id_type node_id_offset,
652 : unique_id_type
653 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
654 : unique_id_offset
655 : #endif
656 : ,
657 : std::unordered_map<subdomain_id_type, subdomain_id_type> *
658 : id_remapping,
659 : const bool skip_preparation)
660 : {
661 28416 : LOG_SCOPE("copy_nodes_and_elements()", "UnstructuredMesh");
662 :
663 : // If we're asked to skip all preparation, we should be skipping
664 : // find_neighbors specifically.
665 14208 : libmesh_assert(!skip_preparation || skip_find_neighbors);
666 :
667 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
668 54780 : extra_int_maps = this->merge_extra_integer_names(other_mesh);
669 :
670 39728 : const unsigned int n_old_node_ints = extra_int_maps.second.size(),
671 39728 : n_new_node_ints = _node_integer_names.size(),
672 39728 : n_old_elem_ints = extra_int_maps.first.size(),
673 39728 : n_new_elem_ints = _elem_integer_names.size();
674 :
675 : // If we are partitioned into fewer parts than the incoming mesh has
676 : // processors to handle, then we need to "wrap" the other Mesh's
677 : // processor ids to fit within our range. This can happen, for
678 : // example, while stitching meshes with small numbers of elements in
679 : // parallel...
680 15052 : bool wrap_proc_ids = (this->n_processors() <
681 15052 : other_mesh.n_partitions());
682 :
683 : // We're assuming the other mesh has proper element number ordering,
684 : // so that we add parents before their children, and that the other
685 : // mesh is consistently partitioned. We're not assuming that node
686 : // proc ids are topologically consistent, so we don't just
687 : // libmesh_assert_valid_procids.
688 : #ifdef DEBUG
689 14208 : MeshTools::libmesh_assert_valid_amr_elem_ids(other_mesh);
690 14208 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(other_mesh);
691 : #endif
692 :
693 : //Copy in Nodes
694 : {
695 : //Preallocate Memory if necessary
696 39728 : this->reserve_nodes(other_mesh.n_nodes());
697 :
698 10214302 : for (const auto & oldn : other_mesh.node_ptr_range())
699 : {
700 : processor_id_type added_pid = cast_int<processor_id_type>
701 6517955 : (wrap_proc_ids ? oldn->processor_id() % this->n_processors() : oldn->processor_id());
702 :
703 : // Add new nodes in old node Point locations
704 : Node * newn =
705 14627777 : this->add_point(*oldn,
706 6517955 : oldn->id() + node_id_offset,
707 2886856 : added_pid);
708 :
709 6517955 : newn->add_extra_integers(n_new_node_ints);
710 6786395 : for (unsigned int i = 0; i != n_old_node_ints; ++i)
711 280388 : newn->set_extra_integer(extra_int_maps.second[i],
712 268440 : oldn->get_extra_integer(i));
713 :
714 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
715 6517955 : newn->set_unique_id(oldn->unique_id() + unique_id_offset);
716 : #endif
717 24676 : }
718 : }
719 :
720 : //Copy in Elements
721 : {
722 : //Preallocate Memory if necessary
723 39728 : this->reserve_elem(other_mesh.n_elem());
724 :
725 : // Declare a map linking old and new elements, needed to copy the neighbor lists
726 : typedef std::unordered_map<const Elem *, Elem *> map_type;
727 28416 : map_type old_elems_to_new_elems, ip_map;
728 :
729 : // Loop over the elements
730 20942282 : for (const auto & old : other_mesh.element_ptr_range())
731 : {
732 : // Build a new element
733 11743819 : Elem * newparent = old->parent() ?
734 1054761 : this->elem_ptr(old->parent()->id() + element_id_offset) :
735 1673074 : nullptr;
736 13104383 : auto el = old->disconnected_clone();
737 1985584 : el->set_parent(newparent);
738 :
739 11431309 : subdomain_id_type sbd_id = old->subdomain_id();
740 11431309 : if (id_remapping)
741 : {
742 538 : auto remapping_it = id_remapping->find(sbd_id);
743 19099 : if (remapping_it != id_remapping->end())
744 568 : sbd_id = remapping_it->second;
745 : }
746 11431309 : el->subdomain_id() = sbd_id;
747 :
748 : // Hold off on trying to set the interior parent because we may actually
749 : // add lower dimensional elements before their interior parents
750 11431309 : if (old->interior_parent())
751 3008 : ip_map[old] = el.get();
752 :
753 : #ifdef LIBMESH_ENABLE_AMR
754 11431309 : if (old->has_children())
755 1388529 : for (unsigned int c = 0, nc = old->n_children(); c != nc; ++c)
756 1147510 : if (old->child_ptr(c) == remote_elem)
757 68325 : el->add_child(const_cast<RemoteElem *>(remote_elem), c);
758 :
759 : //Create the parent's child pointers if necessary
760 11431309 : if (newparent)
761 : {
762 1079165 : unsigned int oldc = old->parent()->which_child_am_i(old);
763 1054761 : newparent->add_child(el.get(), oldc);
764 : }
765 :
766 : // Copy the refinement flags
767 11431309 : el->set_refinement_flag(old->refinement_flag());
768 :
769 : // Use hack_p_level since we may not have sibling elements
770 : // added yet
771 1985584 : el->hack_p_level(old->p_level());
772 :
773 1985584 : el->set_p_refinement_flag(old->p_refinement_flag());
774 : #endif // #ifdef LIBMESH_ENABLE_AMR
775 :
776 : //Assign all the nodes
777 64150951 : for (auto i : el->node_index_range())
778 63940344 : el->set_node(i,
779 52719642 : this->node_ptr(old->node_id(i) + node_id_offset));
780 :
781 : // And start it off with the same processor id (mod _n_parts).
782 11431309 : el->processor_id() = cast_int<processor_id_type>
783 11431309 : (wrap_proc_ids ? old->processor_id() % this->n_processors() : old->processor_id());
784 :
785 : // Give it the same element and unique ids
786 11431309 : el->set_id(old->id() + element_id_offset);
787 :
788 11431309 : el->add_extra_integers(n_new_elem_ints);
789 11453810 : for (unsigned int i = 0; i != n_old_elem_ints; ++i)
790 28721 : el->set_extra_integer(extra_int_maps.first[i],
791 22501 : old->get_extra_integer(i));
792 :
793 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
794 11431309 : el->set_unique_id(old->unique_id() + unique_id_offset);
795 : #endif
796 :
797 : //Hold onto it
798 11431309 : if (!skip_find_neighbors)
799 : {
800 220630 : for (auto s : old->side_index_range())
801 182136 : if (old->neighbor_ptr(s) == remote_elem)
802 256 : el->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
803 46942 : this->add_elem(std::move(el));
804 : }
805 : else
806 : {
807 11698285 : Elem * new_el = this->add_elem(std::move(el));
808 11387183 : old_elems_to_new_elems[old] = new_el;
809 : }
810 9470401 : }
811 :
812 : // If the other_mesh had some interior parents, we may need to
813 : // copy those pointers (if they're to elements in a third mesh),
814 : // or create new equivalent pointers (if they're to elements we
815 : // just copied), or scream and die (if the other mesh had interior
816 : // parents from a third mesh but we already have interior parents
817 : // that aren't to that same third mesh.
818 39728 : if (!ip_map.empty())
819 : {
820 298 : std::atomic<bool> existing_interior_parents{false};
821 :
822 : Threads::parallel_for
823 298 : (this->element_stored_range(),
824 610 : [&existing_interior_parents](const ElemRange & range)
825 : {
826 10582 : for (Elem * elem : range)
827 10284 : if (elem->interior_parent())
828 : {
829 0 : existing_interior_parents = true;
830 0 : break;
831 : }
832 298 : });
833 :
834 : MeshBase * other_interior_mesh =
835 142 : const_cast<MeshBase *>(&other_mesh.interior_mesh());
836 :
837 : // If we don't already have interior parents, then we can just
838 : // use whatever interior_mesh we need for the incoming
839 : // elements.
840 298 : if (!existing_interior_parents)
841 : {
842 298 : if (other_interior_mesh == &other_mesh)
843 62 : this->set_interior_mesh(*this);
844 : else
845 80 : this->set_interior_mesh(*other_interior_mesh);
846 : }
847 :
848 298 : if (other_interior_mesh == &other_mesh &&
849 218 : _interior_mesh == this)
850 1538 : for (auto & elem_pair : ip_map)
851 1320 : elem_pair.second->set_interior_parent(
852 1382 : this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset));
853 80 : else if (other_interior_mesh == _interior_mesh)
854 1768 : for (auto & elem_pair : ip_map)
855 : {
856 1688 : Elem * ip = const_cast<Elem *>(elem_pair.first->interior_parent());
857 1688 : libmesh_assert(ip == remote_elem ||
858 : ip == other_interior_mesh->elem_ptr(ip->id()));
859 1688 : elem_pair.second->set_interior_parent(ip);
860 : }
861 : else
862 0 : libmesh_error_msg("Cannot copy boundary elements between meshes with different interior meshes");
863 : }
864 :
865 : // Loop (again) over the elements to fill in the neighbors
866 39728 : if (skip_find_neighbors)
867 : {
868 39444 : old_elems_to_new_elems[remote_elem] = const_cast<RemoteElem*>(remote_elem);
869 :
870 20856286 : for (const auto & old_elem : other_mesh.element_ptr_range())
871 : {
872 11387183 : Elem * new_elem = old_elems_to_new_elems[old_elem];
873 57160056 : for (auto s : old_elem->side_index_range())
874 : {
875 46697387 : const Elem * old_neighbor = old_elem->neighbor_ptr(s);
876 45461771 : Elem * new_neighbor = old_elems_to_new_elems[old_neighbor];
877 7957326 : new_elem->set_neighbor(s, new_neighbor);
878 : }
879 24408 : }
880 : }
881 : }
882 :
883 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
884 : // We set the unique ids of nodes after adding them to the mesh such that our value of
885 : // _next_unique_id may be wrong. So we amend that here
886 39728 : this->set_next_unique_id(other_mesh.parallel_max_unique_id() + unique_id_offset + 1);
887 : #endif
888 :
889 : // Finally, partially prepare the new Mesh for use, if that isn't
890 : // being skipped.
891 : // Even the default behavior here is for backwards compatibility,
892 : // and we don't want to prepare everything.
893 :
894 39728 : if (!skip_preparation)
895 : {
896 : // Keep the same numbering and partitioning and distribution
897 : // status for now, but save our original policies to restore
898 : // later.
899 100 : const bool allowed_renumbering = this->allow_renumbering();
900 100 : const bool allowed_find_neighbors = this->allow_find_neighbors();
901 100 : const bool allowed_elem_removal = this->allow_remote_element_removal();
902 100 : const bool allowed_detect_detect_interior_parents = this->allow_detect_interior_parents();
903 50 : this->allow_renumbering(false);
904 50 : this->allow_remote_element_removal(false);
905 50 : this->allow_find_neighbors(!skip_find_neighbors);
906 100 : this->allow_detect_interior_parents(other_mesh.allow_detect_interior_parents());
907 :
908 : // We should generally be able to skip *all* partitioning here
909 : // because we're only adding one already-consistent mesh to
910 : // another.
911 100 : const bool skipped_partitioning = this->skip_partitioning();
912 50 : this->skip_partitioning(true);
913 :
914 50 : const Preparation old_preparation = this->preparation();
915 1775 : this->prepare_for_use();
916 :
917 : //But in the long term, don't change our policies.
918 50 : this->allow_find_neighbors(allowed_find_neighbors);
919 50 : this->allow_renumbering(allowed_renumbering);
920 50 : this->allow_remote_element_removal(allowed_elem_removal);
921 50 : this->skip_partitioning(skipped_partitioning);
922 50 : this->allow_detect_interior_parents(allowed_detect_detect_interior_parents);
923 :
924 : // That prepare_for_use() call marked us as prepared, but we
925 : // specifically avoided some important preparation, so we might not
926 : // actually be prepared now.
927 1775 : if (skip_find_neighbors)
928 42 : this->unset_has_neighbor_ptrs();
929 :
930 50 : const Preparation other_preparation = other_mesh.preparation();
931 1775 : if (!old_preparation.is_partitioned ||
932 1775 : !other_preparation.is_partitioned)
933 0 : this->unset_is_partitioned();
934 1775 : if (!old_preparation.has_removed_orphaned_nodes ||
935 1775 : !other_preparation.has_removed_orphaned_nodes)
936 0 : this->unset_has_removed_orphaned_nodes();
937 1775 : if (!old_preparation.has_removed_remote_elements ||
938 1775 : !other_preparation.has_removed_remote_elements)
939 0 : this->unset_has_removed_remote_elements();
940 : }
941 :
942 : // In general we've just invalidated just about everything, and we'd
943 : // like to unset_is_prepared(), but specific use cases might know a
944 : // priori that they're still partitioned well, or that they've
945 : // copied in a disjoint mesh component and don't need new neighbor
946 : // pointers, or that they're not adding anything that would change
947 : // cached subdomain/element/boundary sets, etc., so we'll rely on
948 : // users of the "advanced" skip_preparation option to also set what
949 : // preparation they still need.
950 :
951 : // else
952 : // this->unset_is_prepared();
953 39728 : }
954 :
955 :
956 :
957 399775 : UnstructuredMesh::~UnstructuredMesh ()
958 : {
959 : // this->clear (); // Nothing to clear at this level
960 :
961 24926 : libmesh_exceptionless_assert (!libMesh::closed());
962 399775 : }
963 :
964 :
965 :
966 : namespace {
967 : /**
968 : * \returns \p true if element sides \p a and \p b are the same elem,
969 : * for the purpose of linking them as neighbors in find_neighbors().
970 : *
971 : * This is normally just Elem::operator==, which compares (sorted) node
972 : * ids. But operator== first requires the two sides to have the same
973 : * element type, so it never matches a Polyhedron's polygonal side
974 : * against the TRI3/QUAD4 side of an adjacent standard element (tet,
975 : * hex, ...) even when they are geometrically the same face. At such a
976 : * mixed interface -- exactly one side is a polygon -- we fall back to
977 : * comparing the vertex node ids, which is what low_order_key() already
978 : * keys on. Purely standard/standard and polygon/polygon pairs are left
979 : * entirely to operator==.
980 : */
981 131785006 : bool sides_are_the_same_elem(const Elem & a, const Elem & b)
982 : {
983 131785006 : if (a == b)
984 5113695 : return true;
985 :
986 : // The check above would have been sufficient if they were both
987 : // polygons or both not. If just one is a polygon, we need to check the nodes
988 146 : if (!a.runtime_topology() && !b.runtime_topology())
989 0 : return false;
990 :
991 146 : const unsigned int nv = a.n_vertices();
992 146 : if (nv != b.n_vertices())
993 0 : return false;
994 :
995 162 : std::vector<dof_id_type> a_ids(nv), b_ids(nv);
996 584 : for (unsigned int v = 0; v != nv; ++v)
997 : {
998 450 : a_ids[v] = a.node_id(v);
999 462 : b_ids[v] = b.node_id(v);
1000 : }
1001 146 : std::sort(a_ids.begin(), a_ids.end());
1002 146 : std::sort(b_ids.begin(), b_ids.end());
1003 146 : return a_ids == b_ids;
1004 : }
1005 : }
1006 :
1007 :
1008 :
1009 656699 : void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
1010 : const bool reset_current_list,
1011 : const bool assert_valid,
1012 : const bool check_non_remote)
1013 : {
1014 : // We might actually want to run this on an empty mesh
1015 : // (e.g. the boundary mesh for a nonexistent bcid!)
1016 : // libmesh_assert_not_equal_to (this->n_nodes(), 0);
1017 : // libmesh_assert_not_equal_to (this->n_elem(), 0);
1018 :
1019 : // This function must be run on all processors at once
1020 27992 : parallel_object_only();
1021 :
1022 27992 : LOG_SCOPE("find_neighbors()", "Mesh");
1023 :
1024 : //TODO:[BSK] This should be removed later?!
1025 656699 : if (reset_current_list)
1026 : Threads::parallel_for
1027 585947 : (this->element_stored_range(),
1028 1681361 : [reset_remote_elements](const ElemRange & range)
1029 : {
1030 71780318 : for (Elem * e : range)
1031 356149024 : for (auto s : e->side_index_range())
1032 290668733 : if (e->neighbor_ptr(s) != remote_elem || reset_remote_elements)
1033 11202686 : e->set_neighbor(s, nullptr);
1034 586453 : });
1035 :
1036 : // Find neighboring elements by first finding elements
1037 : // with identical side keys and then check to see if they
1038 : // are neighbors
1039 : {
1040 : // data structures -- Use the hash_multimap if available
1041 : typedef dof_id_type key_type;
1042 : typedef std::pair<Elem *, unsigned char> val_type;
1043 : typedef std::unordered_multimap<key_type, val_type> map_type;
1044 :
1045 : // A map from side keys to corresponding elements & side numbers
1046 55984 : map_type side_to_elem_map;
1047 :
1048 : // Pull objects out of the loop to reduce heap operations
1049 656699 : std::unique_ptr<Elem> my_side, their_side;
1050 :
1051 151419676 : for (const auto & element : this->element_ptr_range())
1052 : {
1053 389263631 : for (auto ms : element->side_index_range())
1054 : {
1055 442220579 : next_side:
1056 : // If we haven't yet found a neighbor on this side, try.
1057 : // Even if we think our neighbor is remote, that
1058 : // information may be out of date.
1059 : //
1060 : // If we're only checking remote neighbors, after a
1061 : // redistribution, then we'll skip the non-remote ones
1062 455769644 : if ((element->neighbor_ptr(ms) == nullptr && check_non_remote) ||
1063 159372739 : element->neighbor_ptr(ms) == remote_elem)
1064 : {
1065 : // Get the key for the side of this element. Use the
1066 : // low_order_key so we can find neighbors in
1067 : // mixed-order meshes if necessary.
1068 286726829 : const dof_id_type key = element->low_order_key(ms);
1069 :
1070 : // Look for elements that have an identical side key
1071 11222174 : auto bounds = side_to_elem_map.equal_range(key);
1072 :
1073 : // May be multiple keys, check all the possible
1074 : // elements which _might_ be neighbors.
1075 286726829 : if (bounds.first != bounds.second)
1076 : {
1077 : // Get the side for this element
1078 131666523 : element->side_ptr(my_side, ms);
1079 :
1080 : // Look at all the entries with an equivalent key
1081 131874550 : while (bounds.first != bounds.second)
1082 : {
1083 : // Get the potential element
1084 131785006 : Elem * neighbor = bounds.first->second.first;
1085 :
1086 : // Get the side for the neighboring element
1087 131785006 : const unsigned int ns = bounds.first->second.second;
1088 131785006 : neighbor->side_ptr(their_side, ns);
1089 : //libmesh_assert(my_side.get());
1090 : //libmesh_assert(their_side.get());
1091 :
1092 : // If found a match with my side
1093 : //
1094 : // In 1D, since parents and children have an
1095 : // equal side (i.e. a node) we need to check
1096 : // for matching level() to avoid setting our
1097 : // neighbor pointer to any of our neighbor's
1098 : // descendants.
1099 263570012 : if (sides_are_the_same_elem(*my_side, *their_side) &&
1100 131785006 : (element->level() == neighbor->level()))
1101 : {
1102 : // So share a side. Is this a mixed pair
1103 : // of subactive and active/ancestor
1104 : // elements?
1105 : // If not, then we're neighbors.
1106 : // If so, then the subactive's neighbor is
1107 :
1108 134179701 : if (element->subactive() ==
1109 131576979 : neighbor->subactive())
1110 : {
1111 : // an element is only subactive if it has
1112 : // been coarsened but not deleted
1113 7688869 : element->set_neighbor (ms,neighbor);
1114 131396042 : neighbor->set_neighbor(ns,element);
1115 : }
1116 180937 : else if (element->subactive())
1117 : {
1118 7644 : element->set_neighbor(ms,neighbor);
1119 : }
1120 99221 : else if (neighbor->subactive())
1121 : {
1122 14052 : neighbor->set_neighbor(ns,element);
1123 : }
1124 5107843 : side_to_elem_map.erase (bounds.first);
1125 :
1126 : // get out of this nested crap
1127 131576979 : goto next_side;
1128 : }
1129 :
1130 5860 : ++bounds.first;
1131 : }
1132 : }
1133 :
1134 : // didn't find a match...
1135 : // Build the map entry for this element
1136 : side_to_elem_map.emplace
1137 155149850 : (key, std::make_pair(element, cast_int<unsigned char>(ms)));
1138 : }
1139 : }
1140 614355 : }
1141 614355 : }
1142 :
1143 : #ifdef LIBMESH_ENABLE_PERIODIC
1144 : // Get the disjoint neighbor boundary pairs object (from periodic BCs)
1145 656699 : auto * db = this->get_disjoint_neighbor_boundary_pairs();
1146 :
1147 656699 : if (db)
1148 : {
1149 : // Obtain a point locator
1150 1085 : std::unique_ptr<PointLocatorBase> point_locator = this->sub_point_locator();
1151 :
1152 7202 : for (const auto & element : this->element_ptr_range())
1153 : {
1154 13322 : for (auto ms : element->side_index_range())
1155 : {
1156 : // Skip if this side already has a valid neighbor (including remote neighbors)
1157 10912 : if (element->neighbor_ptr(ms) != nullptr &&
1158 2047 : element->neighbor_ptr(ms) != remote_elem)
1159 2003 : continue;
1160 :
1161 32779 : for (const auto & [id, boundary_ptr] : *db)
1162 : {
1163 24210 : if (!this->get_boundary_info().has_boundary_id(element, ms, id))
1164 22747 : continue;
1165 :
1166 : unsigned int neigh_side;
1167 : const Elem * neigh =
1168 1491 : db->neighbor(id, *point_locator, element, ms, &neigh_side);
1169 :
1170 1463 : if (neigh && neigh != remote_elem && neigh != element)
1171 : {
1172 1463 : auto neigh_changeable = this->elem_ptr(neigh->id());
1173 1463 : element->set_neighbor(ms, neigh_changeable);
1174 1463 : neigh_changeable->set_neighbor(neigh_side, element);
1175 : }
1176 : }
1177 : }
1178 985 : }
1179 985 : }
1180 : #endif // LIBMESH_ENABLE_PERIODIC
1181 :
1182 : #ifdef LIBMESH_ENABLE_AMR
1183 :
1184 : /**
1185 : * Here we look at all of the child elements which
1186 : * don't already have valid neighbors.
1187 : *
1188 : * If a child element has a nullptr neighbor it is
1189 : * either because it is on the boundary or because
1190 : * its neighbor is at a different level. In the
1191 : * latter case we must get the neighbor from the
1192 : * parent.
1193 : *
1194 : * If a child element has a remote_elem neighbor
1195 : * on a boundary it shares with its parent, that
1196 : * info may have become out-dated through coarsening
1197 : * of the neighbor's parent. In this case, if the
1198 : * parent's neighbor is active then the child should
1199 : * share it.
1200 : *
1201 : * Furthermore, that neighbor better be active,
1202 : * otherwise we missed a child somewhere.
1203 : *
1204 : *
1205 : * We also need to look through children ordered by increasing
1206 : * refinement level in order to add new interior_parent() links in
1207 : * boundary elements which have just been generated by refinement,
1208 : * and fix links in boundary elements whose previous
1209 : * interior_parent() has just been coarsened away.
1210 : */
1211 656699 : const unsigned int n_levels = MeshTools::n_levels(*this);
1212 936292 : for (unsigned int level = 1; level < n_levels; ++level)
1213 : {
1214 2195590 : for (auto & current_elem : as_range(level_elements_begin(level),
1215 101926628 : level_elements_end(level)))
1216 : {
1217 1645384 : libmesh_assert(current_elem);
1218 50828777 : Elem * parent = current_elem->parent();
1219 1645384 : libmesh_assert(parent);
1220 50828777 : const unsigned int my_child_num = parent->which_child_am_i(current_elem);
1221 :
1222 252944373 : for (auto s : current_elem->side_index_range())
1223 : {
1224 210537240 : if (current_elem->neighbor_ptr(s) == nullptr ||
1225 193638569 : (current_elem->neighbor_ptr(s) == remote_elem &&
1226 1859896 : parent->is_child_on_side(my_child_num, s)))
1227 : {
1228 698600 : Elem * neigh = parent->neighbor_ptr(s);
1229 :
1230 : // If neigh was refined and had non-subactive children
1231 : // made remote earlier, then our current elem should
1232 : // actually have one of those remote children as a
1233 : // neighbor
1234 15079721 : if (neigh &&
1235 4476844 : (neigh->ancestor() ||
1236 : // If neigh has subactive children which should have
1237 : // matched as neighbors of the current element but
1238 : // did not, then those likewise must be remote
1239 : // children.
1240 4289835 : (current_elem->subactive() && neigh->has_children() &&
1241 189 : (neigh->level()+1) == current_elem->level())))
1242 : {
1243 : #ifdef DEBUG
1244 : // Let's make sure that "had children made remote"
1245 : // situation is actually the case
1246 0 : libmesh_assert(neigh->has_children());
1247 0 : bool neigh_has_remote_children = false;
1248 0 : for (auto & child : neigh->child_ref_range())
1249 0 : if (&child == remote_elem)
1250 0 : neigh_has_remote_children = true;
1251 0 : libmesh_assert(neigh_has_remote_children);
1252 :
1253 : // And let's double-check that we don't have
1254 : // a remote_elem neighboring an active local element
1255 0 : if (current_elem->active())
1256 0 : libmesh_assert_not_equal_to (current_elem->processor_id(),
1257 : this->processor_id());
1258 : #endif // DEBUG
1259 78693 : neigh = const_cast<RemoteElem *>(remote_elem);
1260 : }
1261 : // If neigh and current_elem are more than one level
1262 : // apart, figuring out whether we have a remote
1263 : // neighbor here becomes much harder.
1264 10648678 : else if (neigh && (current_elem->subactive() &&
1265 15800 : neigh->has_children()))
1266 : {
1267 : // Find the deepest descendant of neigh which
1268 : // we could consider for a neighbor. If we run
1269 : // out of neigh children, then that's our
1270 : // neighbor. If we find a potential neighbor
1271 : // with remote_children and we don't find any
1272 : // potential neighbors among its non-remote
1273 : // children, then our neighbor must be remote.
1274 0 : while (neigh != remote_elem &&
1275 0 : neigh->has_children())
1276 : {
1277 0 : bool found_neigh = false;
1278 0 : for (unsigned int c = 0, nc = neigh->n_children();
1279 0 : !found_neigh && c != nc; ++c)
1280 : {
1281 0 : Elem * child = neigh->child_ptr(c);
1282 0 : if (child == remote_elem)
1283 0 : continue;
1284 0 : for (auto ncn : child->neighbor_ptr_range())
1285 : {
1286 0 : if (ncn != remote_elem &&
1287 0 : ncn->is_ancestor_of(current_elem))
1288 : {
1289 0 : neigh = ncn;
1290 0 : found_neigh = true;
1291 0 : break;
1292 : }
1293 : }
1294 : }
1295 0 : if (!found_neigh)
1296 0 : neigh = const_cast<RemoteElem *>(remote_elem);
1297 : }
1298 : }
1299 10711571 : current_elem->set_neighbor(s, neigh);
1300 : #ifdef DEBUG
1301 464104 : if (neigh != nullptr && neigh != remote_elem)
1302 : // We ignore subactive elements here because
1303 : // we don't care about neighbors of subactive element.
1304 214304 : if ((!neigh->active()) && (!current_elem->subactive()))
1305 : {
1306 0 : libMesh::err << "On processor " << this->processor_id()
1307 0 : << std::endl;
1308 0 : libMesh::err << "Bad element ID = " << current_elem->id()
1309 0 : << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
1310 0 : libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
1311 0 : << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
1312 0 : libMesh::err << "Bad element size = " << current_elem->hmin()
1313 0 : << ", Bad neighbor size = " << neigh->hmin() << std::endl;
1314 0 : libMesh::err << "Bad element center = " << current_elem->vertex_average()
1315 0 : << ", Bad neighbor center = " << neigh->vertex_average() << std::endl;
1316 0 : libMesh::err << "ERROR: "
1317 0 : << (current_elem->active()?"Active":"Ancestor")
1318 0 : << " Element at level "
1319 0 : << current_elem->level() << std::endl;
1320 0 : libMesh::err << "with "
1321 0 : << (parent->active()?"active":
1322 0 : (parent->subactive()?"subactive":"ancestor"))
1323 0 : << " parent share "
1324 0 : << (neigh->subactive()?"subactive":"ancestor")
1325 0 : << " neighbor at level " << neigh->level()
1326 0 : << std::endl;
1327 0 : NameBasedIO(*this).write ("bad_mesh.gmv");
1328 0 : libmesh_error_msg("Problematic mesh written to bad_mesh.gmv.");
1329 : }
1330 : #endif // DEBUG
1331 : }
1332 : }
1333 :
1334 : // We can skip to the next element if we're full-dimension
1335 : // and therefore don't have any interior parents
1336 50828777 : if (current_elem->dim() >= LIBMESH_DIM)
1337 8453873 : continue;
1338 :
1339 : // We have no interior parents unless we can find one later
1340 43408960 : current_elem->set_interior_parent(nullptr);
1341 :
1342 43408960 : Elem * pip = parent->interior_parent();
1343 :
1344 43408960 : if (!pip)
1345 42758003 : continue;
1346 :
1347 : // If there's no interior_parent children, whether due to a
1348 : // remote element or a non-conformity, then there's no
1349 : // children to search.
1350 25419 : if (pip == remote_elem || pip->active())
1351 : {
1352 2052 : current_elem->set_interior_parent(pip);
1353 2052 : continue;
1354 : }
1355 :
1356 : // For node comparisons we'll need a sensible tolerance
1357 23367 : Real node_tolerance = current_elem->hmin() * TOLERANCE;
1358 :
1359 : // Otherwise our interior_parent should be a child of our
1360 : // parent's interior_parent.
1361 77739 : for (auto & child : pip->child_ref_range())
1362 : {
1363 : // If we have a remote_elem, that might be our
1364 : // interior_parent. We'll set it provisionally now and
1365 : // keep trying to find something better.
1366 76601 : if (&child == remote_elem)
1367 : {
1368 : current_elem->set_interior_parent
1369 4788 : (const_cast<RemoteElem *>(remote_elem));
1370 4788 : continue;
1371 : }
1372 :
1373 5860 : bool child_contains_our_nodes = true;
1374 145312 : for (auto & n : current_elem->node_ref_range())
1375 : {
1376 10176 : bool child_contains_this_node = false;
1377 795528 : for (auto & cn : child.node_ref_range())
1378 745944 : if (cn.absolute_fuzzy_equals
1379 714744 : (n, node_tolerance))
1380 : {
1381 6268 : child_contains_this_node = true;
1382 6268 : break;
1383 : }
1384 120075 : if (!child_contains_this_node)
1385 : {
1386 3908 : child_contains_our_nodes = false;
1387 3908 : break;
1388 : }
1389 : }
1390 71813 : if (child_contains_our_nodes)
1391 : {
1392 22229 : current_elem->set_interior_parent(&child);
1393 1952 : break;
1394 : }
1395 : }
1396 :
1397 : // We should have found *some* interior_parent at this
1398 : // point, whether semilocal or remote.
1399 1952 : libmesh_assert(current_elem->interior_parent());
1400 266057 : }
1401 : }
1402 : #endif // AMR
1403 :
1404 : #ifdef DEBUG
1405 27992 : if (assert_valid)
1406 : {
1407 27856 : MeshTools::libmesh_assert_valid_neighbors(*this,
1408 27856 : !reset_remote_elements);
1409 27856 : MeshTools::libmesh_assert_valid_amr_interior_parents(*this);
1410 : }
1411 : #else
1412 : libmesh_ignore(assert_valid);
1413 : #endif
1414 :
1415 656699 : this->_preparation.has_neighbor_ptrs = true;
1416 656699 : }
1417 :
1418 :
1419 :
1420 5776 : void UnstructuredMesh::read (const std::string & name,
1421 : void *,
1422 : bool skip_renumber_nodes_and_elements,
1423 : bool skip_find_neighbors,
1424 : bool skip_detect_interior_parents)
1425 : {
1426 : // Set the skip_renumber_nodes_and_elements flag on all processors
1427 : // if necessary.
1428 : // This ensures that renumber_nodes_and_elements is *not* called
1429 : // during prepare_for_use() for certain types of mesh files.
1430 : // This is required in cases where there is an associated solution
1431 : // file which expects a certain ordering of the nodes.
1432 5776 : if (Utility::ends_with(name, ".gmv"))
1433 0 : this->allow_renumbering(false);
1434 :
1435 5776 : NameBasedIO(*this).read(name);
1436 :
1437 5776 : if (skip_renumber_nodes_and_elements)
1438 : {
1439 : // Use MeshBase::allow_renumbering() yourself instead.
1440 : libmesh_deprecated();
1441 0 : this->allow_renumbering(false);
1442 : }
1443 :
1444 : // Done reading the mesh. Now prepare it for use.
1445 356 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
1446 356 : const bool old_allow_detect_interior_parents = this->allow_detect_interior_parents();
1447 :
1448 178 : this->allow_find_neighbors(!skip_find_neighbors);
1449 178 : this->allow_detect_interior_parents(!skip_detect_interior_parents);
1450 :
1451 5776 : this->prepare_for_use();
1452 :
1453 178 : this->allow_find_neighbors(old_allow_find_neighbors);
1454 178 : this->allow_detect_interior_parents(old_allow_detect_interior_parents);
1455 5776 : }
1456 :
1457 :
1458 :
1459 2938 : void UnstructuredMesh::write (const std::string & name) const
1460 : {
1461 598 : LOG_SCOPE("write()", "Mesh");
1462 :
1463 3536 : NameBasedIO(*this).write(name);
1464 2938 : }
1465 :
1466 :
1467 :
1468 0 : void UnstructuredMesh::write (const std::string & name,
1469 : const std::vector<Number> & v,
1470 : const std::vector<std::string> & vn) const
1471 : {
1472 0 : LOG_SCOPE("write()", "Mesh");
1473 :
1474 0 : NameBasedIO(*this).write_nodal_data(name, v, vn);
1475 0 : }
1476 :
1477 :
1478 :
1479 :
1480 :
1481 0 : void UnstructuredMesh::create_pid_mesh(UnstructuredMesh & pid_mesh,
1482 : const processor_id_type pid) const
1483 : {
1484 :
1485 : // Issue a warning if the number the number of processors
1486 : // currently available is less that that requested for
1487 : // partitioning. This is not necessarily an error since
1488 : // you may run on one processor and still partition the
1489 : // mesh into several partitions.
1490 : #ifdef DEBUG
1491 0 : if (this->n_processors() < pid)
1492 : {
1493 0 : libMesh::out << "WARNING: You are creating a "
1494 0 : << "mesh for a processor id (="
1495 0 : << pid
1496 0 : << ") greater than "
1497 0 : << "the number of processors available for "
1498 0 : << "the calculation. (="
1499 0 : << this->n_processors()
1500 0 : << ")."
1501 0 : << std::endl;
1502 : }
1503 : #endif
1504 :
1505 0 : this->create_submesh (pid_mesh,
1506 0 : this->active_pid_elements_begin(pid),
1507 0 : this->active_pid_elements_end(pid));
1508 0 : }
1509 :
1510 :
1511 :
1512 :
1513 :
1514 :
1515 :
1516 0 : void UnstructuredMesh::create_submesh (UnstructuredMesh & new_mesh,
1517 : const const_element_iterator & it,
1518 : const const_element_iterator & it_end) const
1519 : {
1520 : // Just in case the subdomain_mesh already has some information
1521 : // in it, get rid of it.
1522 0 : new_mesh.clear();
1523 :
1524 : // If we're not serial, our submesh isn't either.
1525 : // There are no remote elements to delete on an empty mesh, but
1526 : // calling the method to do so marks the mesh as parallel.
1527 0 : if (!this->is_serial())
1528 0 : new_mesh.delete_remote_elements();
1529 :
1530 : // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
1531 : // This may happen if the user accidentally passes the original mesh into
1532 : // this function! We will check this by making sure we did not just
1533 : // clear ourself.
1534 0 : libmesh_assert_not_equal_to (this->n_nodes(), 0);
1535 0 : libmesh_assert_not_equal_to (this->n_elem(), 0);
1536 :
1537 : // Container to catch boundary IDs handed back by BoundaryInfo
1538 0 : std::vector<boundary_id_type> bc_ids;
1539 :
1540 : // Put any extra integers on the new mesh too
1541 0 : new_mesh.merge_extra_integer_names(*this);
1542 0 : const unsigned int n_node_ints = _node_integer_names.size();
1543 :
1544 0 : for (const auto & old_elem : as_range(it, it_end))
1545 : {
1546 : // Add an equivalent element type to the new_mesh.
1547 : // disconnected_clone() copies ids, extra element integers, etc.
1548 0 : auto uelem = old_elem->disconnected_clone();
1549 0 : Elem * new_elem = new_mesh.add_elem(std::move(uelem));
1550 0 : libmesh_assert(new_elem);
1551 :
1552 : // Loop over the nodes on this element.
1553 0 : for (auto n : old_elem->node_index_range())
1554 : {
1555 0 : const dof_id_type this_node_id = old_elem->node_id(n);
1556 :
1557 : // Add this node to the new mesh if it's not there already
1558 0 : if (!new_mesh.query_node_ptr(this_node_id))
1559 : {
1560 : Node * newn =
1561 0 : new_mesh.add_point (old_elem->point(n),
1562 : this_node_id,
1563 0 : old_elem->node_ptr(n)->processor_id());
1564 :
1565 0 : newn->add_extra_integers(n_node_ints);
1566 0 : for (unsigned int i = 0; i != n_node_ints; ++i)
1567 0 : newn->set_extra_integer(i, old_elem->node_ptr(n)->get_extra_integer(i));
1568 :
1569 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1570 0 : newn->set_unique_id(old_elem->node_ptr(n)->unique_id());
1571 : #endif
1572 : }
1573 :
1574 : // Define this element's connectivity on the new mesh
1575 0 : new_elem->set_node(n, new_mesh.node_ptr(this_node_id));
1576 : }
1577 :
1578 : // Maybe add boundary conditions for this element
1579 0 : for (auto s : old_elem->side_index_range())
1580 : {
1581 0 : this->get_boundary_info().boundary_ids(old_elem, s, bc_ids);
1582 0 : new_mesh.get_boundary_info().add_side (new_elem, s, bc_ids);
1583 : }
1584 0 : } // end loop over elements
1585 :
1586 : // Prepare the new_mesh for use
1587 0 : new_mesh.prepare_for_use();
1588 0 : }
1589 :
1590 :
1591 :
1592 : #ifdef LIBMESH_ENABLE_AMR
1593 25897 : bool UnstructuredMesh::contract ()
1594 : {
1595 860 : LOG_SCOPE ("contract()", "Mesh");
1596 :
1597 : // Flag indicating if this call actually changes the mesh
1598 860 : bool mesh_changed = false;
1599 :
1600 : #ifdef DEBUG
1601 545354 : for (const auto & elem : this->element_ptr_range())
1602 544494 : libmesh_assert(elem->active() || elem->subactive() || elem->ancestor());
1603 : #endif
1604 :
1605 : // Loop over the elements.
1606 20051476 : for (auto & elem : this->element_ptr_range())
1607 : {
1608 : // Delete all the subactive ones
1609 10544765 : if (elem->subactive())
1610 : {
1611 : // No level-0 element should be subactive.
1612 : // Note that we CAN'T test elem->level(), as that
1613 : // touches elem->parent()->dim(), and elem->parent()
1614 : // might have already been deleted!
1615 71808 : libmesh_assert(elem->parent());
1616 :
1617 : // Delete the element
1618 : // This just sets a pointer to nullptr, and doesn't
1619 : // invalidate any iterators
1620 1499078 : this->delete_elem(elem);
1621 :
1622 : // the mesh has certainly changed
1623 71808 : mesh_changed = true;
1624 : }
1625 : else
1626 : {
1627 : // Compress all the active ones
1628 472686 : if (elem->active())
1629 6776920 : elem->contract();
1630 : else
1631 113362 : libmesh_assert (elem->ancestor());
1632 : }
1633 24177 : }
1634 :
1635 : // Strip any newly-created nullptr voids out of the element array
1636 25897 : this->renumber_nodes_and_elements();
1637 :
1638 : // FIXME: Need to understand why deleting subactive children
1639 : // invalidates the point locator. For now we will clear it explicitly
1640 25897 : this->clear_point_locator();
1641 :
1642 : // Allow our GhostingFunctor objects to reinit if necessary.
1643 27863 : for (auto & gf : as_range(this->ghosting_functors_begin(),
1644 116549 : this->ghosting_functors_end()))
1645 : {
1646 2826 : libmesh_assert(gf);
1647 86966 : gf->mesh_reinit();
1648 : }
1649 :
1650 26757 : return mesh_changed;
1651 : }
1652 : #endif // #ifdef LIBMESH_ENABLE_AMR
1653 :
1654 :
1655 :
1656 11401 : void UnstructuredMesh::all_first_order ()
1657 : {
1658 792 : LOG_SCOPE("all_first_order()", "Mesh");
1659 :
1660 : /**
1661 : * Prepare to identify (and then delete) a bunch of no-longer-used nodes.
1662 : */
1663 11797 : std::vector<bool> node_touched_by_me(this->max_node_id(), false);
1664 :
1665 : // Loop over the high-ordered elements.
1666 : // First make sure they _are_ indeed high-order, and then replace
1667 : // them with an equivalent first-order element.
1668 484366 : for (auto & so_elem : element_ptr_range())
1669 : {
1670 25880 : libmesh_assert(so_elem);
1671 :
1672 : /*
1673 : * build the first-order equivalent, add to
1674 : * the new_elements list.
1675 : */
1676 : auto lo_elem = Elem::build
1677 : (Elem::first_order_equivalent_type
1678 282740 : (so_elem->type()), so_elem->parent());
1679 :
1680 256860 : const unsigned short n_sides = so_elem->n_sides();
1681 :
1682 1225522 : for (unsigned short s=0; s != n_sides; ++s)
1683 1066562 : if (so_elem->neighbor_ptr(s) == remote_elem)
1684 0 : lo_elem->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1685 :
1686 : #ifdef LIBMESH_ENABLE_AMR
1687 : /*
1688 : * Reset the parent links of any child elements
1689 : */
1690 256860 : if (so_elem->has_children())
1691 376697 : for (unsigned int c = 0, nc = so_elem->n_children(); c != nc; ++c)
1692 : {
1693 295596 : Elem * child = so_elem->child_ptr(c);
1694 295596 : if (child != remote_elem)
1695 48208 : child->set_parent(lo_elem.get());
1696 295596 : lo_elem->add_child(child, c);
1697 : }
1698 :
1699 : /*
1700 : * Reset the child link of any parent element
1701 : */
1702 282740 : if (so_elem->parent())
1703 : {
1704 : unsigned int c =
1705 231063 : so_elem->parent()->which_child_am_i(so_elem);
1706 255167 : lo_elem->parent()->replace_child(lo_elem.get(), c);
1707 : }
1708 :
1709 : /*
1710 : * Copy as much data to the new element as makes sense
1711 : */
1712 282740 : lo_elem->set_p_level(so_elem->p_level());
1713 256860 : lo_elem->set_refinement_flag(so_elem->refinement_flag());
1714 51760 : lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
1715 : #endif
1716 :
1717 25880 : libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
1718 :
1719 : /*
1720 : * By definition the vertices of the linear and
1721 : * second order element are identically numbered.
1722 : * transfer these.
1723 : */
1724 1234462 : for (unsigned int v=0, snv=so_elem->n_vertices(); v < snv; v++)
1725 : {
1726 1075982 : lo_elem->set_node(v, so_elem->node_ptr(v));
1727 295140 : node_touched_by_me[lo_elem->node_id(v)] = true;
1728 : }
1729 :
1730 : /*
1731 : * find_neighbors relies on remote_elem neighbor links being
1732 : * properly maintained.
1733 : */
1734 1225522 : for (unsigned short s=0; s != n_sides; s++)
1735 : {
1736 1066562 : if (so_elem->neighbor_ptr(s) == remote_elem)
1737 0 : lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1738 : }
1739 :
1740 : /**
1741 : * If the second order element had any boundary conditions they
1742 : * should be transferred to the first-order element. The old
1743 : * boundary conditions will be removed from the BoundaryInfo
1744 : * data structure by insert_elem.
1745 : */
1746 25880 : this->get_boundary_info().copy_boundary_ids
1747 256860 : (this->get_boundary_info(), so_elem, lo_elem.get());
1748 :
1749 : /*
1750 : * The new first-order element is ready.
1751 : * Inserting it into the mesh will replace and delete
1752 : * the second-order element.
1753 : */
1754 256860 : lo_elem->set_id(so_elem->id());
1755 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1756 51760 : lo_elem->set_unique_id(so_elem->unique_id());
1757 : #endif
1758 :
1759 256860 : const unsigned int nei = so_elem->n_extra_integers();
1760 256860 : lo_elem->add_extra_integers(nei);
1761 259687 : for (unsigned int i=0; i != nei; ++i)
1762 2827 : lo_elem->set_extra_integer(i, so_elem->get_extra_integer(i));
1763 :
1764 256860 : lo_elem->inherit_data_from(*so_elem);
1765 :
1766 308620 : this->insert_elem(std::move(lo_elem));
1767 215709 : }
1768 :
1769 : // Deleting nodes does not invalidate iterators, so this is safe.
1770 1701128 : for (const auto & node : this->node_ptr_range())
1771 1014397 : if (!node_touched_by_me[node->id()])
1772 631294 : this->delete_node(node);
1773 :
1774 : // If crazy people applied boundary info to non-vertices and then
1775 : // deleted those non-vertices, we should make sure their boundary id
1776 : // caches are correct.
1777 11401 : this->get_boundary_info().regenerate_id_sets();
1778 :
1779 : // On hanging nodes that used to also be second order nodes, we
1780 : // might now have an invalid nodal processor_id()
1781 11401 : Partitioner::set_node_processor_ids(*this);
1782 :
1783 : // delete or renumber nodes if desired
1784 11401 : this->prepare_for_use();
1785 11401 : }
1786 :
1787 :
1788 :
1789 : void
1790 22980 : UnstructuredMesh::all_second_order_range (const SimpleRange<element_iterator> & range,
1791 : const bool full_ordered)
1792 : {
1793 666 : LOG_SCOPE("all_second_order_range()", "Mesh");
1794 :
1795 : /*
1796 : * The maximum number of new second order nodes we might be adding,
1797 : * for use when picking unique unique_id values later. This variable
1798 : * is not used unless unique ids are enabled.
1799 : */
1800 : unsigned int max_new_nodes_per_elem;
1801 :
1802 : /*
1803 : * For speed-up of the \p add_point() method, we
1804 : * can reserve memory. Guess the number of additional
1805 : * nodes based on the element spatial dimensions and the
1806 : * total number of nodes in the mesh as an upper bound.
1807 : */
1808 22980 : switch (this->mesh_dimension())
1809 : {
1810 923 : case 1:
1811 : /*
1812 : * in 1D, there can only be order-increase from Edge2
1813 : * to Edge3. Something like 1/2 of n_nodes() have
1814 : * to be added
1815 : */
1816 26 : max_new_nodes_per_elem = 3 - 2;
1817 1846 : this->reserve_nodes(static_cast<unsigned int>
1818 923 : (1.5*static_cast<double>(this->n_nodes())));
1819 897 : break;
1820 :
1821 3728 : case 2:
1822 : /*
1823 : * in 2D, either refine from Tri3 to Tri6 (double the nodes)
1824 : * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
1825 : */
1826 120 : max_new_nodes_per_elem = 9 - 4;
1827 7456 : this->reserve_nodes(static_cast<unsigned int>
1828 3728 : (2*static_cast<double>(this->n_nodes())));
1829 3608 : break;
1830 :
1831 :
1832 18329 : case 3:
1833 : /*
1834 : * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
1835 : * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
1836 : * quite some nodes, and since we do not want to overburden the memory by
1837 : * a too conservative guess, use the lower bound
1838 : */
1839 520 : max_new_nodes_per_elem = 27 - 8;
1840 36658 : this->reserve_nodes(static_cast<unsigned int>
1841 18329 : (2.5*static_cast<double>(this->n_nodes())));
1842 17809 : break;
1843 :
1844 0 : default:
1845 : // Hm?
1846 0 : libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1847 : }
1848 :
1849 : // All the real work is done in the helper function
1850 22980 : all_increased_order_range(*this, range, max_new_nodes_per_elem,
1851 93320 : [full_ordered](ElemType t) {
1852 2219083 : return Elem::second_order_equivalent_type(t, full_ordered);
1853 : });
1854 22980 : }
1855 :
1856 :
1857 :
1858 21787 : void UnstructuredMesh::all_complete_order_range(const SimpleRange<element_iterator> & range)
1859 : {
1860 610 : LOG_SCOPE("all_complete_order()", "Mesh");
1861 :
1862 : /*
1863 : * The maximum number of new higher-order nodes we might be adding,
1864 : * for use when picking unique unique_id values later. This variable
1865 : * is not used unless unique ids are enabled.
1866 : */
1867 : unsigned int max_new_nodes_per_elem;
1868 :
1869 : /*
1870 : * for speed-up of the \p add_point() method, we
1871 : * can reserve memory. Guess the number of additional
1872 : * nodes based on the element spatial dimensions and the
1873 : * total number of nodes in the mesh as an upper bound.
1874 : */
1875 21787 : switch (this->mesh_dimension())
1876 : {
1877 0 : case 1:
1878 : /*
1879 : * in 1D, there can only be order-increase from Edge2
1880 : * to Edge3. Something like 1/2 of n_nodes() have
1881 : * to be added
1882 : */
1883 0 : max_new_nodes_per_elem = 3 - 2;
1884 0 : this->reserve_nodes(static_cast<unsigned int>
1885 0 : (1.5*static_cast<double>(this->n_nodes())));
1886 0 : break;
1887 :
1888 1775 : case 2:
1889 : /*
1890 : * in 2D, we typically refine from Tri3 or Tri6 to Tri7 (2.3333
1891 : * or 1.1667 times the nodes) but might refine from Quad4 to
1892 : * Quad9 (2.25 times the nodes)
1893 : */
1894 50 : max_new_nodes_per_elem = 9 - 4;
1895 3550 : this->reserve_nodes(static_cast<unsigned int>
1896 1775 : (2*static_cast<double>(this->n_nodes())));
1897 1725 : break;
1898 :
1899 :
1900 20012 : case 3:
1901 : /*
1902 : * in 3D, we typically refine from Tet10 to Tet14 (factor = 1.4)
1903 : * but may go Hex8 to Hex27 or Tet4 to Tet14 (something > 3).
1904 : * Since in 3D there _are_ already quite some nodes, and since
1905 : * we do not want to overburden the memory by a too-conservative
1906 : * guess, use a moderate bound
1907 : */
1908 560 : max_new_nodes_per_elem = 27 - 8;
1909 40024 : this->reserve_nodes(static_cast<unsigned int>
1910 20012 : (2.5*static_cast<double>(this->n_nodes())));
1911 19452 : break;
1912 :
1913 0 : default:
1914 : // Hm?
1915 0 : libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1916 : }
1917 :
1918 : // All the real work is done in the helper function
1919 21787 : all_increased_order_range(*this, range, max_new_nodes_per_elem,
1920 146881 : [](ElemType t) {
1921 5340193 : return Elem::complete_order_equivalent_type(t);
1922 : });
1923 21787 : }
1924 :
1925 :
1926 : std::size_t
1927 1633 : UnstructuredMesh::stitch_meshes (const MeshBase & other_mesh,
1928 : boundary_id_type this_mesh_boundary_id,
1929 : boundary_id_type other_mesh_boundary_id,
1930 : Real tol,
1931 : bool clear_stitched_boundary_ids,
1932 : bool verbose,
1933 : bool use_binary_search,
1934 : bool enforce_all_nodes_match_on_boundaries,
1935 : bool merge_boundary_nodes_all_or_nothing,
1936 : bool remap_subdomain_ids,
1937 : bool prepare_after_stitching)
1938 : {
1939 92 : LOG_SCOPE("stitch_meshes()", "UnstructuredMesh");
1940 1633 : return stitching_helper(&other_mesh,
1941 : this_mesh_boundary_id,
1942 : other_mesh_boundary_id,
1943 : tol,
1944 : clear_stitched_boundary_ids,
1945 : verbose,
1946 : use_binary_search,
1947 : enforce_all_nodes_match_on_boundaries,
1948 : true,
1949 : merge_boundary_nodes_all_or_nothing,
1950 : remap_subdomain_ids,
1951 1533 : prepare_after_stitching);
1952 : }
1953 :
1954 :
1955 : std::size_t
1956 213 : UnstructuredMesh::stitch_surfaces (boundary_id_type boundary_id_1,
1957 : boundary_id_type boundary_id_2,
1958 : Real tol,
1959 : bool clear_stitched_boundary_ids,
1960 : bool verbose,
1961 : bool use_binary_search,
1962 : bool enforce_all_nodes_match_on_boundaries,
1963 : bool merge_boundary_nodes_all_or_nothing,
1964 : bool prepare_after_stitching)
1965 :
1966 : {
1967 213 : return stitching_helper(nullptr,
1968 : boundary_id_1,
1969 : boundary_id_2,
1970 : tol,
1971 : clear_stitched_boundary_ids,
1972 : verbose,
1973 : use_binary_search,
1974 : enforce_all_nodes_match_on_boundaries,
1975 : /* skip_find_neighbors = */ true,
1976 : merge_boundary_nodes_all_or_nothing,
1977 : /* remap_subdomain_ids = */ false,
1978 213 : prepare_after_stitching);
1979 : }
1980 :
1981 :
1982 : std::size_t
1983 1846 : UnstructuredMesh::stitching_helper (const MeshBase * other_mesh,
1984 : boundary_id_type this_mesh_boundary_id,
1985 : boundary_id_type other_mesh_boundary_id,
1986 : Real tol,
1987 : bool clear_stitched_boundary_ids,
1988 : bool verbose,
1989 : bool use_binary_search,
1990 : bool enforce_all_nodes_match_on_boundaries,
1991 : bool skip_find_neighbors,
1992 : bool merge_boundary_nodes_all_or_nothing,
1993 : bool remap_subdomain_ids,
1994 : bool prepare_after_stitching)
1995 : {
1996 : #ifdef DEBUG
1997 : // We rely on neighbor links here
1998 52 : MeshTools::libmesh_assert_valid_neighbors(*this);
1999 : #endif
2000 :
2001 52 : bool is_valid_disjoint_pair_to_stitch = false;
2002 :
2003 : #ifdef LIBMESH_ENABLE_PERIODIC
2004 1846 : auto * this_db = this->get_disjoint_neighbor_boundary_pairs();
2005 1846 : auto * other_db = (other_mesh ? other_mesh->get_disjoint_neighbor_boundary_pairs() : nullptr);
2006 : const bool have_disc_bdys =
2007 1846 : (this_db && !this_db->empty()) || (other_db && !other_db->empty());
2008 :
2009 52 : if (have_disc_bdys)
2010 : {
2011 10 : const boundary_id_type a = this_mesh_boundary_id;
2012 10 : const boundary_id_type b = other_mesh_boundary_id;
2013 :
2014 40 : auto get_pb = [](const PeriodicBoundaries * db, boundary_id_type id)
2015 : {
2016 539 : return db ? db->boundary(id) : nullptr;
2017 : };
2018 :
2019 : // this mesh
2020 355 : const auto * pb_this_a = get_pb(this_db, a);
2021 355 : const auto * pb_this_b = get_pb(this_db, b);
2022 10 : const bool in_this =
2023 355 : (pb_this_a && pb_this_a->pairedboundary == b) ||
2024 0 : (pb_this_b && pb_this_b->pairedboundary == a);
2025 :
2026 : // other mesh
2027 345 : const auto * pb_other_b = get_pb(other_db, b);
2028 10 : const auto * pb_other_a = get_pb(other_db, a);
2029 10 : const bool in_other =
2030 355 : (pb_other_b && pb_other_b->pairedboundary == a) ||
2031 0 : (pb_other_a && pb_other_a->pairedboundary == b);
2032 :
2033 : // Conflict conditions:
2034 : // Case 1: On "this" mesh, a or b exist but are not paired,
2035 : // while the other mesh pairs them.
2036 355 : if (!in_this && (pb_this_a || pb_this_b) && in_other)
2037 146 : libmesh_error_msg("Disjoint neighbor boundary pairing mismatch: on 'this' mesh, "
2038 : "boundary (" << a << " or " << b
2039 : << ") exists but is not paired; on 'other' mesh the pair is present.");
2040 :
2041 : // Case 2: On "other" mesh, a or b exist but are not paired,
2042 : // while this mesh pairs them.
2043 284 : if (!in_other && (pb_other_a || pb_other_b) && in_this)
2044 0 : libmesh_error_msg("Disjoint neighbor boundary pairing mismatch: on 'other' mesh, "
2045 : "boundary (" << a << " or " << b
2046 : << ") exists but is not paired; on 'this' mesh the pair is present.");
2047 :
2048 : // Legal conditions: either side has a correct pairing
2049 284 : if (in_this || in_other)
2050 6 : is_valid_disjoint_pair_to_stitch = true;
2051 : }
2052 : #endif // LIBMESH_ENABLE_PERIODIC
2053 :
2054 : // We can't even afford any unset neighbor links here.
2055 1775 : if (!this->is_prepared())
2056 71 : this->find_neighbors();
2057 :
2058 : // FIXME: make distributed mesh support efficient.
2059 : // Yes, we currently suck.
2060 1875 : MeshSerializer serialize(*this);
2061 :
2062 : // *Badly*.
2063 1725 : std::unique_ptr<MeshSerializer> serialize_other;
2064 1775 : if (other_mesh)
2065 : serialize_other = std::make_unique<MeshSerializer>
2066 3036 : (*const_cast<MeshBase *>(other_mesh));
2067 :
2068 102 : std::map<dof_id_type, dof_id_type> node_to_node_map, other_to_this_node_map; // The second is the inverse map of the first
2069 100 : std::map<dof_id_type, std::vector<dof_id_type>> node_to_elems_map;
2070 :
2071 : typedef dof_id_type key_type;
2072 : typedef std::pair<const Elem *, unsigned char> val_type;
2073 : typedef std::pair<key_type, val_type> key_val_pair;
2074 : typedef std::unordered_multimap<key_type, val_type> map_type;
2075 : // Mapping between all side keys in this mesh and elements+side numbers relevant to the boundary in this mesh as well.
2076 100 : map_type side_to_elem_map;
2077 :
2078 : // If there is only one mesh (i.e. other_mesh == nullptr), then loop over this mesh twice
2079 1775 : if (!other_mesh)
2080 : {
2081 6 : other_mesh = this;
2082 : }
2083 :
2084 1775 : if ((this_mesh_boundary_id != BoundaryInfo::invalid_id) &&
2085 50 : (other_mesh_boundary_id != BoundaryInfo::invalid_id))
2086 : {
2087 100 : LOG_SCOPE("stitch_meshes node merging", "UnstructuredMesh");
2088 :
2089 : // While finding nodes on the boundary, also find the minimum edge length
2090 : // of all faces on both boundaries. This will later be used in relative
2091 : // distance checks when stitching nodes.
2092 1775 : Real h_min = std::numeric_limits<Real>::max();
2093 50 : bool h_min_updated = false;
2094 :
2095 : // Loop below fills in these sets for the two meshes.
2096 100 : std::set<dof_id_type> this_boundary_node_ids, other_boundary_node_ids;
2097 :
2098 : // Pull objects out of the loop to reduce heap operations
2099 1775 : std::unique_ptr<const Elem> side;
2100 :
2101 : {
2102 : // Make temporary fixed-size arrays for loop
2103 1775 : boundary_id_type id_array[2] = {this_mesh_boundary_id, other_mesh_boundary_id};
2104 1775 : std::set<dof_id_type> * set_array[2] = {&this_boundary_node_ids, &other_boundary_node_ids};
2105 1775 : const MeshBase * mesh_array[2] = {this, other_mesh};
2106 :
2107 5325 : for (unsigned i=0; i<2; ++i)
2108 : {
2109 : // First we deal with node boundary IDs. We only enter
2110 : // this loop if we have at least one nodeset. Note that we
2111 : // do not attempt to make an h_min determination here.
2112 : // The h_min determination is done while looping over the
2113 : // Elems and checking their sides and edges for boundary
2114 : // information, below.
2115 3650 : if (mesh_array[i]->get_boundary_info().n_nodeset_conds() > 0)
2116 : {
2117 : // build_node_list() returns a vector of (node-id, bc-id) tuples
2118 324270 : for (const auto & t : mesh_array[i]->get_boundary_info().build_node_list())
2119 : {
2120 321204 : boundary_id_type node_bc_id = std::get<1>(t);
2121 321204 : if (node_bc_id == id_array[i])
2122 : {
2123 61415 : dof_id_type this_node_id = std::get<0>(t);
2124 61415 : set_array[i]->insert( this_node_id );
2125 : }
2126 : }
2127 : }
2128 :
2129 : // Container to catch boundary IDs passed back from BoundaryInfo.
2130 200 : std::vector<boundary_id_type> bc_ids;
2131 :
2132 : // Pointers to boundary NodeElems encountered while looping over the entire Mesh
2133 : // and checking side and edge boundary ids. The Nodes associated with NodeElems
2134 : // may be in a boundary nodeset, but not connected to any other Elems. In this
2135 : // case, we also consider the "minimum node separation distance" amongst all
2136 : // NodeElems when determining the relevant h_min value for this mesh.
2137 200 : std::vector<const Elem *> boundary_node_elems;
2138 :
2139 86902 : for (auto & el : mesh_array[i]->element_ptr_range())
2140 : {
2141 : // Now check whether elem has a face on the specified boundary
2142 267337 : for (auto side_id : el->side_index_range())
2143 : {
2144 : bool should_stitch_this_side =
2145 231450 : (el->neighbor_ptr(side_id) == nullptr) ||
2146 1420 : (is_valid_disjoint_pair_to_stitch &&
2147 1460 : mesh_array[i]->get_boundary_info().has_boundary_id(el, side_id, id_array[i]));
2148 :
2149 6340 : if (should_stitch_this_side)
2150 : {
2151 : // Get *all* boundary IDs on this side, not just the first one!
2152 99718 : mesh_array[i]->get_boundary_info().boundary_ids (el, side_id, bc_ids);
2153 :
2154 99718 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2155 : {
2156 16472 : el->build_side_ptr(side, side_id);
2157 112535 : for (auto & n : side->node_ref_range())
2158 96063 : set_array[i]->insert(n.id());
2159 :
2160 16472 : h_min = std::min(h_min, side->hmin());
2161 464 : h_min_updated = true;
2162 :
2163 : // This side is on the boundary, add its information to side_to_elem
2164 16472 : if (skip_find_neighbors && (i==0))
2165 : {
2166 8236 : key_type key = el->low_order_key(side_id);
2167 232 : val_type val;
2168 8236 : val.first = el;
2169 232 : val.second = cast_int<unsigned char>(side_id);
2170 :
2171 8236 : key_val_pair kvp;
2172 8236 : kvp.first = key;
2173 232 : kvp.second = val;
2174 232 : side_to_elem_map.insert (kvp);
2175 : }
2176 : }
2177 :
2178 : // Also, check the edges on this side. We don't have to worry about
2179 : // updating neighbor info in this case since elements don't store
2180 : // neighbor info on edges.
2181 1211294 : for (auto edge_id : el->edge_index_range())
2182 : {
2183 1111576 : if (el->is_edge_on_side(edge_id, side_id))
2184 : {
2185 : // Get *all* boundary IDs on this edge, not just the first one!
2186 368490 : mesh_array[i]->get_boundary_info().edge_boundary_ids (el, edge_id, bc_ids);
2187 :
2188 368490 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2189 : {
2190 0 : std::unique_ptr<const Elem> edge (el->build_edge_ptr(edge_id));
2191 0 : for (auto & n : edge->node_ref_range())
2192 0 : set_array[i]->insert( n.id() );
2193 :
2194 0 : h_min = std::min(h_min, edge->hmin());
2195 0 : h_min_updated = true;
2196 0 : }
2197 : }
2198 : } // end for (edge_id)
2199 : } // end if (should_stitch_this_side)
2200 : } // end for (side_id)
2201 :
2202 : // Alternatively, is this a boundary NodeElem? If so,
2203 : // add it to a list of NodeElems that will later be
2204 : // used to set h_min based on the minimum node
2205 : // separation distance between all pairs of boundary
2206 : // NodeElems.
2207 41109 : if (el->type() == NODEELEM)
2208 : {
2209 2700 : mesh_array[i]->get_boundary_info().boundary_ids(el->node_ptr(0), bc_ids);
2210 2628 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2211 : {
2212 2556 : boundary_node_elems.push_back(el);
2213 :
2214 : // Debugging:
2215 : // libMesh::out << "Elem " << el->id() << " is a NodeElem on boundary " << id_array[i] << std::endl;
2216 : }
2217 : } // end if (el->type() == NODEELEM)
2218 3350 : } // end for (el)
2219 :
2220 : // Compute the minimum node separation distance amongst
2221 : // all boundary NodeElem pairs.
2222 : {
2223 200 : const auto N = boundary_node_elems.size();
2224 6106 : for (auto node_elem_i : make_range(N))
2225 47286 : for (auto node_elem_j : make_range(node_elem_i+1, N))
2226 : {
2227 : Real node_sep =
2228 47250 : (boundary_node_elems[node_elem_i]->point(0) - boundary_node_elems[node_elem_j]->point(0)).norm();
2229 :
2230 : // We only want to consider non-coincident
2231 : // boundary NodeElem pairs when determining the
2232 : // minimum node separation distance.
2233 44730 : if (node_sep > 0.)
2234 : {
2235 44730 : h_min = std::min(h_min, node_sep);
2236 1260 : h_min_updated = true;
2237 : }
2238 : } // end for (node_elem_j)
2239 : } // end minimum NodeElem separation scope
2240 : } // end for (i)
2241 : } // end scope
2242 :
2243 1775 : if (verbose)
2244 : {
2245 14 : libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
2246 14 : << "This mesh has " << this_boundary_node_ids.size()
2247 14 : << " nodes on boundary `"
2248 497 : << this->get_boundary_info().get_sideset_name(this_mesh_boundary_id)
2249 14 : << "' (" << this_mesh_boundary_id << ").\n"
2250 14 : << "Other mesh has " << other_boundary_node_ids.size()
2251 14 : << " nodes on boundary `"
2252 497 : << other_mesh->get_boundary_info().get_sideset_name(other_mesh_boundary_id)
2253 14 : << "' (" << other_mesh_boundary_id << ").\n";
2254 :
2255 497 : if (h_min_updated)
2256 : {
2257 28 : libMesh::out << "Minimum edge length on both surfaces is " << h_min << ".\n";
2258 : }
2259 : else
2260 : {
2261 0 : libMesh::out << "No minimum edge length determined on specified surfaces." << std::endl;
2262 : }
2263 : }
2264 :
2265 : // At this point, if h_min==0 it means that there were at least two coincident
2266 : // nodes on the surfaces being stitched, and we don't currently support that case.
2267 : // (It might be possible to support, but getting it exactly right would be tricky
2268 : // and probably not worth the extra complications to the "normal" case.)
2269 1775 : libmesh_error_msg_if(h_min < std::numeric_limits<Real>::epsilon(),
2270 : "Coincident nodes detected on source and/or target "
2271 : "surface, stitching meshes is not possible.");
2272 :
2273 : // We require nanoflann for the "binary search" (really kd-tree)
2274 : // option to work. If it's not available, turn that option off,
2275 : // warn the user, and fall back on the N^2 search algorithm.
2276 : if (use_binary_search)
2277 : {
2278 : #ifndef LIBMESH_HAVE_NANOFLANN
2279 : use_binary_search = false;
2280 : libmesh_warning("The use_binary_search option in the "
2281 : "UnstructuredMesh stitching algorithms requires nanoflann "
2282 : "support. Falling back on N^2 search algorithm.");
2283 : #endif
2284 : }
2285 :
2286 1775 : if (!this_boundary_node_ids.empty())
2287 : {
2288 1775 : if (use_binary_search)
2289 : {
2290 : #ifdef LIBMESH_HAVE_NANOFLANN
2291 : typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, VectorOfNodesAdaptor>,
2292 : VectorOfNodesAdaptor, 3, std::size_t> kd_tree_t;
2293 :
2294 : // Create the dataset needed to build the kd tree with nanoflann
2295 0 : std::vector<std::pair<Point, dof_id_type>> this_mesh_nodes(this_boundary_node_ids.size());
2296 :
2297 0 : for (auto [it, ctr] = std::make_tuple(this_boundary_node_ids.begin(), 0u);
2298 0 : it != this_boundary_node_ids.end(); ++it, ++ctr)
2299 : {
2300 0 : this_mesh_nodes[ctr].first = this->point(*it);
2301 0 : this_mesh_nodes[ctr].second = *it;
2302 : }
2303 :
2304 0 : VectorOfNodesAdaptor vec_nodes_adaptor(this_mesh_nodes);
2305 :
2306 0 : kd_tree_t this_kd_tree(3, vec_nodes_adaptor, 10);
2307 0 : this_kd_tree.buildIndex();
2308 :
2309 : // Storage for nearest neighbor in the loop below
2310 : std::size_t ret_index;
2311 : Real ret_dist_sqr;
2312 :
2313 : // Loop over other mesh. For each node, find its nearest neighbor in this mesh, and fill in the maps.
2314 0 : for (const auto & node_id : other_boundary_node_ids)
2315 : {
2316 0 : const auto & p = other_mesh->point(node_id);
2317 0 : const Real query_pt[] = {p(0), p(1), p(2)};
2318 0 : this_kd_tree.knnSearch(&query_pt[0], 1, &ret_index, &ret_dist_sqr);
2319 :
2320 : // TODO: here we should use the user's specified tolerance
2321 : // and the previously determined value of h_min in the
2322 : // distance comparison, not just TOLERANCE^2.
2323 0 : if (ret_dist_sqr < TOLERANCE*TOLERANCE)
2324 : {
2325 0 : node_to_node_map[this_mesh_nodes[ret_index].second] = node_id;
2326 0 : other_to_this_node_map[node_id] = this_mesh_nodes[ret_index].second;
2327 : }
2328 : }
2329 :
2330 : // If the two maps don't have the same size, it means one
2331 : // node in this mesh is the nearest neighbor of several
2332 : // nodes in other mesh. Since the stitching is ambiguous in
2333 : // this case, we throw an error.
2334 0 : libmesh_error_msg_if(node_to_node_map.size() != other_to_this_node_map.size(),
2335 : "Error: Found multiple matching nodes in stitch_meshes");
2336 : #endif
2337 : }
2338 : else // !use_binary_search
2339 : {
2340 : // In the unlikely event that two meshes composed entirely of
2341 : // NodeElems are being stitched together, we will not have
2342 : // selected a valid h_min value yet, and the distance
2343 : // comparison below will be true for essentially any two
2344 : // nodes. In this case we simply fall back on an absolute
2345 : // distance check.
2346 1775 : if (!h_min_updated)
2347 : {
2348 : libmesh_warning("No valid h_min value was found, falling back on "
2349 : "absolute distance check in the N^2 search algorithm.");
2350 0 : h_min = 1.;
2351 : }
2352 :
2353 : // Otherwise, use a simple N^2 search to find the closest matching points. This can be helpful
2354 : // in the case that we have tolerance issues which cause mismatch between the two surfaces
2355 : // that are being stitched.
2356 32944 : for (const auto & this_node_id : this_boundary_node_ids)
2357 : {
2358 31169 : Node & this_node = this->node_ref(this_node_id);
2359 :
2360 878 : bool found_matching_nodes = false;
2361 :
2362 928538 : for (const auto & other_node_id : other_boundary_node_ids)
2363 : {
2364 897369 : const Node & other_node = other_mesh->node_ref(other_node_id);
2365 :
2366 872091 : Real node_distance = (this_node - other_node).norm();
2367 :
2368 897369 : if (node_distance < tol*h_min)
2369 : {
2370 : // Make sure we didn't already find a matching node!
2371 31169 : libmesh_error_msg_if(found_matching_nodes,
2372 : "Error: Found multiple matching nodes in stitch_meshes");
2373 :
2374 31169 : node_to_node_map[this_node_id] = other_node_id;
2375 31169 : other_to_this_node_map[other_node_id] = this_node_id;
2376 :
2377 878 : found_matching_nodes = true;
2378 : }
2379 : }
2380 : }
2381 : }
2382 : }
2383 :
2384 : // Build up the node_to_elems_map, using only one loop over other_mesh
2385 42830 : for (auto & el : other_mesh->element_ptr_range())
2386 : {
2387 : // For each node on the element, find the corresponding node
2388 : // on "this" Mesh, 'this_node_id', if it exists, and push
2389 : // the current element ID back onto node_to_elems_map[this_node_id].
2390 : // For that we will use the reverse mapping we created at
2391 : // the same time as the forward mapping.
2392 294865 : for (auto & n : el->node_ref_range())
2393 281780 : if (const auto it = other_to_this_node_map.find(/*other_node_id=*/n.id());
2394 7720 : it != other_to_this_node_map.end())
2395 50481 : node_to_elems_map[/*this_node_id=*/it->second].push_back( el->id() );
2396 1675 : }
2397 :
2398 1775 : if (verbose)
2399 : {
2400 14 : libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
2401 28 : << "Found " << node_to_node_map.size()
2402 14 : << " matching nodes.\n"
2403 14 : << std::endl;
2404 : }
2405 :
2406 1775 : if (enforce_all_nodes_match_on_boundaries)
2407 : {
2408 2 : std::size_t n_matching_nodes = node_to_node_map.size();
2409 2 : std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2410 2 : std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2411 71 : libmesh_error_msg_if((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes),
2412 : "Error: We expected the number of nodes to match.");
2413 : }
2414 :
2415 1775 : if (merge_boundary_nodes_all_or_nothing)
2416 : {
2417 2 : std::size_t n_matching_nodes = node_to_node_map.size();
2418 2 : std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2419 2 : std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2420 71 : if ((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes))
2421 : {
2422 0 : if (verbose)
2423 : {
2424 : libMesh::out << "Skipping node merging in "
2425 : "UnstructuredMesh::stitch_meshes because not "
2426 0 : "all boundary nodes were matched."
2427 0 : << std::endl;
2428 : }
2429 0 : node_to_node_map.clear();
2430 0 : other_to_this_node_map.clear();
2431 0 : node_to_elems_map.clear();
2432 : }
2433 100 : }
2434 3350 : }
2435 : else
2436 : {
2437 0 : if (verbose)
2438 : {
2439 0 : libMesh::out << "Skip node merging in UnstructuredMesh::stitch_meshes:" << std::endl;
2440 : }
2441 : }
2442 :
2443 1775 : dof_id_type node_delta = this->max_node_id();
2444 1775 : dof_id_type elem_delta = this->max_elem_id();
2445 :
2446 : unique_id_type unique_delta =
2447 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2448 1775 : this->parallel_max_unique_id();
2449 : #else
2450 : 0;
2451 : #endif
2452 :
2453 : // If other_mesh != nullptr, then we have to do a bunch of work
2454 : // in order to copy it to this mesh
2455 1775 : if (this!=other_mesh)
2456 : {
2457 88 : LOG_SCOPE("stitch_meshes copying", "UnstructuredMesh");
2458 :
2459 : #ifdef LIBMESH_ENABLE_PERIODIC
2460 : // Copy disjoint neighbor boundary pairs (PeriodicBoundary objects)
2461 : // from `other_mesh` to `this` mesh
2462 1562 : if (other_db && !other_db->empty())
2463 : {
2464 213 : for (const auto & [bdy_id, pb_ptr] : *other_db)
2465 : {
2466 4 : const auto & pb = *pb_ptr;
2467 142 : const boundary_id_type a = pb.myboundary;
2468 142 : const boundary_id_type b = pb.pairedboundary;
2469 :
2470 142 : if (this_db)
2471 : {
2472 : // Skip if identical pair already exists
2473 142 : if (const auto * existing_pb = this_db->boundary(a))
2474 71 : if ((existing_pb->myboundary == a && existing_pb->pairedboundary == b) ||
2475 0 : (existing_pb->myboundary == b && existing_pb->pairedboundary == a))
2476 69 : continue;
2477 :
2478 : // If both boundary ids exist on this mesh but aren't paired here, refuse to create a new pair
2479 2 : const auto & bdy_ids = this->get_boundary_info().get_boundary_ids();
2480 4 : const bool a_exists = bdy_ids.count(a);
2481 4 : const bool b_exists = bdy_ids.count(b);
2482 : // If a and b already exist on `this`, we should be screaming and dying
2483 : // unless they already have PeriodicBoundary objects connecting them too
2484 71 : if (a_exists && b_exists && !this_db->boundary(a))
2485 0 : libmesh_error_msg("Conflict: boundaries " << a << " and " << b
2486 : << " already exist on this mesh but are not paired.");
2487 : }
2488 :
2489 140 : this->add_disjoint_neighbor_boundary_pairs(a, b, pb.get_corresponding_pos(Point(0.0,0.0,0.0)));
2490 : }
2491 : }
2492 : #endif // LIBMESH_ENABLE_PERIODIC
2493 :
2494 :
2495 : // Increment the node_to_node_map and node_to_elems_map
2496 : // to account for id offsets
2497 32305 : for (auto & pr : node_to_node_map)
2498 30743 : pr.second += node_delta;
2499 :
2500 32305 : for (auto & pr : node_to_elems_map)
2501 80798 : for (auto & entry : pr.second)
2502 50055 : entry += elem_delta;
2503 :
2504 : // We run into problems when the libMesh subdomain standard (the
2505 : // id defines the subdomain; the name was an afterthought) and
2506 : // the MOOSE standard (the name defines the subdomain; the id
2507 : // might be autogenerated) clash.
2508 : //
2509 : // Subdomain ids with the same name in both meshes are surely
2510 : // meant to represent the same subdomain. We can just merge
2511 : // them.
2512 : //
2513 : // Subdomain ids which don't have a name in either mesh are
2514 : // almost surely meant to represent the same subdomain. We'll
2515 : // just merge them.
2516 : //
2517 : // Subdomain ids with different names in different meshes, or
2518 : // names with different ids in different meshes, are trickier.
2519 : // For backwards compatibility we default to the old "just copy
2520 : // all the subdomain ids over" behavior, but if requested we'll
2521 : // remap any ids that appear to be clear conflicts, and we'll
2522 : // scream and die if we see any ids that are ambiguous due to
2523 : // being named in one mesh but not the other.
2524 88 : std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
2525 1562 : if (remap_subdomain_ids)
2526 : {
2527 4 : const auto & this_map = this->get_subdomain_name_map();
2528 4 : const auto & other_map = other_mesh->get_subdomain_name_map();
2529 8 : std::unordered_map<std::string, subdomain_id_type> other_map_reversed;
2530 284 : for (auto & [sid, sname] : other_map)
2531 4 : other_map_reversed.emplace(sname, sid);
2532 :
2533 8 : std::unordered_map<std::string, subdomain_id_type> this_map_reversed;
2534 213 : for (auto & [sid, sname] : this_map)
2535 2 : this_map_reversed.emplace(sname, sid);
2536 :
2537 : // We don't require either mesh to be prepared, but that
2538 : // means we need to check for subdomains manually.
2539 284 : auto get_subdomains = [](const MeshBase & mesh) {
2540 8 : std::set<subdomain_id_type> all_subdomains;
2541 4976 : for (auto & el : mesh.element_ptr_range())
2542 2540 : all_subdomains.insert(el->subdomain_id());
2543 284 : return all_subdomains;
2544 : };
2545 :
2546 146 : const auto this_subdomains = get_subdomains(*this);
2547 146 : const auto other_subdomains = get_subdomains(*other_mesh);
2548 :
2549 213 : for (auto & [sid, sname] : this_map)
2550 : {
2551 : // The same name with the same id means we're fine. The
2552 : // same name with another id means we remap their id to
2553 : // ours
2554 2 : if (const auto other_reverse_it = other_map_reversed.find(sname);
2555 71 : other_reverse_it != other_map_reversed.end() && other_reverse_it->second != sid)
2556 71 : id_remapping[other_reverse_it->second] = sid;
2557 :
2558 : // The same id with a different name, we'll get to
2559 : // later. The same id without any name means we don't
2560 : // know what the user wants.
2561 2 : if (other_subdomains.count(sid) && !other_map.count(sid))
2562 0 : libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2563 : << sid << " but not subdomain name " << sname);
2564 : }
2565 :
2566 142 : subdomain_id_type next_free_id = 0;
2567 : // We might try to stitch empty meshes ...
2568 142 : if (!this_subdomains.empty())
2569 142 : next_free_id = *this_subdomains.rbegin() + 1;
2570 142 : if (!other_subdomains.empty())
2571 142 : next_free_id =
2572 142 : std::max(next_free_id,
2573 : cast_int<subdomain_id_type>
2574 215 : (*other_subdomains.rbegin() + 1));
2575 :
2576 213 : for (auto & [sid, sname] : other_map)
2577 : {
2578 : // At this point we've figured out any remapping
2579 : // necessary for an sname that we share. And we don't
2580 : // need to remap any sid we don't share.
2581 8 : if (!this_map_reversed.count(sname))
2582 : {
2583 : // But if we don't have this sname and we do have this
2584 : // sid then we can't just merge into that.
2585 2 : if (this_subdomains.count(sid))
2586 : {
2587 : // If we have this sid with no name, we don't
2588 : // know what the user wants.
2589 2 : if (!this_map.count(sid))
2590 211 : libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2591 : << sid << " but under subdomain name " << sname);
2592 :
2593 : // We have this sid under a different name, so
2594 : // we just need to give the other elements a new
2595 : // id.
2596 :
2597 : // Users might have done crazy things with id
2598 : // choice so let's make sure they didn't get too
2599 : // crazy.
2600 0 : libmesh_error_msg_if ((!this_subdomains.empty() &&
2601 : next_free_id < *this_subdomains.rbegin()) ||
2602 : (!other_subdomains.empty() &&
2603 : next_free_id < *other_subdomains.rbegin()),
2604 : "Subdomain id overflow");
2605 :
2606 0 : id_remapping[sid] = next_free_id++;
2607 0 : this->set_subdomain_name(next_free_id, sname);
2608 : }
2609 : // If we don't have this subdomain id, well, we're
2610 : // about to, so we should have its name too.
2611 : else
2612 0 : this->set_subdomain_name(sid, sname);
2613 : }
2614 : }
2615 : }
2616 :
2617 : // Copy mesh data. If we skip the call to find_neighbors(), the lists
2618 : // of neighbors will be copied verbatim from the other mesh
2619 1491 : this->copy_nodes_and_elements(*other_mesh, skip_find_neighbors,
2620 : elem_delta, node_delta,
2621 84 : unique_delta, &id_remapping);
2622 :
2623 : // Copy BoundaryInfo from other_mesh too. We do this via the
2624 : // list APIs rather than element-by-element for speed.
2625 42 : BoundaryInfo & boundary = this->get_boundary_info();
2626 42 : const BoundaryInfo & other_boundary = other_mesh->get_boundary_info();
2627 :
2628 158869 : for (const auto & t : other_boundary.build_node_list())
2629 157336 : boundary.add_node(std::get<0>(t) + node_delta,
2630 157336 : std::get<1>(t));
2631 :
2632 42855 : for (const auto & t : other_boundary.build_side_list())
2633 42486 : boundary.add_side(std::get<0>(t) + elem_delta,
2634 41322 : std::get<1>(t),
2635 41322 : std::get<2>(t));
2636 :
2637 1533 : for (const auto & t : other_boundary.build_edge_list())
2638 0 : boundary.add_edge(std::get<0>(t) + elem_delta,
2639 0 : std::get<1>(t),
2640 0 : std::get<2>(t));
2641 :
2642 1533 : for (const auto & t : other_boundary.build_shellface_list())
2643 0 : boundary.add_shellface(std::get<0>(t) + elem_delta,
2644 0 : std::get<1>(t),
2645 0 : std::get<2>(t));
2646 :
2647 42 : const auto & other_ns_id_to_name = other_boundary.get_nodeset_name_map();
2648 42 : auto & ns_id_to_name = boundary.set_nodeset_name_map();
2649 1491 : ns_id_to_name.insert(other_ns_id_to_name.begin(), other_ns_id_to_name.end());
2650 :
2651 42 : const auto & other_ss_id_to_name = other_boundary.get_sideset_name_map();
2652 42 : auto & ss_id_to_name = boundary.set_sideset_name_map();
2653 1491 : ss_id_to_name.insert(other_ss_id_to_name.begin(), other_ss_id_to_name.end());
2654 :
2655 42 : const auto & other_es_id_to_name = other_boundary.get_edgeset_name_map();
2656 42 : auto & es_id_to_name = boundary.set_edgeset_name_map();
2657 1491 : es_id_to_name.insert(other_es_id_to_name.begin(), other_es_id_to_name.end());
2658 :
2659 : // Merge other_mesh's elemset information with ours. Throw an
2660 : // error if this and other_mesh have overlapping elemset codes
2661 : // that refer to different elemset ids.
2662 1533 : std::vector<dof_id_type> this_elemset_codes = this->get_elemset_codes();
2663 84 : MeshBase::elemset_type this_id_set_to_fill, other_id_set_to_fill;
2664 1817 : for (const auto & elemset_code : other_mesh->get_elemset_codes())
2665 : {
2666 : // Get the elemset ids for this elemset_code on other_mesh
2667 284 : other_mesh->get_elemsets(elemset_code, other_id_set_to_fill);
2668 :
2669 : // Check that this elemset code does not already exist
2670 : // in this mesh, or if it does, that it has the same elemset
2671 : // ids associated with it.
2672 : //
2673 : // Note: get_elemset_codes() is guaranteed to return a
2674 : // sorted vector, so we can binary search in it.
2675 268 : auto it = Utility::binary_find(this_elemset_codes.begin(),
2676 : this_elemset_codes.end(),
2677 16 : elemset_code);
2678 :
2679 284 : if (it != this_elemset_codes.end())
2680 : {
2681 : // This mesh has the same elemset code. Does it refer to
2682 : // the same elemset ids?
2683 0 : this->get_elemsets(elemset_code, this_id_set_to_fill);
2684 :
2685 : // Throw an error if they don't match, otherwise we
2686 : // don't need to do anything
2687 0 : libmesh_error_msg_if(other_id_set_to_fill != this_id_set_to_fill,
2688 : "Attempted to stitch together meshes with conflicting elemset codes.");
2689 : }
2690 : else
2691 : {
2692 : // Add other_mesh's elemset code to this mesh
2693 560 : this->add_elemset_code(elemset_code, other_id_set_to_fill);
2694 : }
2695 : }
2696 :
2697 : } // end if (other_mesh)
2698 :
2699 : // Finally, we need to "merge" the overlapping nodes
2700 : // We do this by iterating over node_to_elems_map and updating
2701 : // the elements so that they "point" to the nodes that came
2702 : // from this mesh, rather than from other_mesh.
2703 : // Then we iterate over node_to_node_map and delete the
2704 : // duplicate nodes that came from other_mesh.
2705 :
2706 : {
2707 96 : LOG_SCOPE("stitch_meshes node updates", "UnstructuredMesh");
2708 :
2709 : // Container to catch boundary IDs passed back from BoundaryInfo.
2710 96 : std::vector<boundary_id_type> bc_ids;
2711 :
2712 32234 : for (const auto & [target_node_id, elem_vec] : node_to_elems_map)
2713 : {
2714 30530 : dof_id_type other_node_id = node_to_node_map[target_node_id];
2715 30530 : Node & target_node = this->node_ref(target_node_id);
2716 :
2717 1720 : std::size_t n_elems = elem_vec.size();
2718 79875 : for (std::size_t i=0; i<n_elems; i++)
2719 : {
2720 49345 : dof_id_type elem_id = elem_vec[i];
2721 49345 : Elem * el = this->elem_ptr(elem_id);
2722 :
2723 : // find the local node index that we want to update
2724 47955 : unsigned int local_node_index = el->local_node(other_node_id);
2725 1390 : libmesh_assert_not_equal_to(local_node_index, libMesh::invalid_uint);
2726 :
2727 : // We also need to copy over the nodeset info here,
2728 : // because the node will get deleted below
2729 50735 : this->get_boundary_info().boundary_ids(el->node_ptr(local_node_index), bc_ids);
2730 49345 : el->set_node(local_node_index, &target_node);
2731 49345 : this->get_boundary_info().add_node(&target_node, bc_ids);
2732 : }
2733 : }
2734 : }
2735 :
2736 : {
2737 96 : LOG_SCOPE("stitch_meshes node deletion", "UnstructuredMesh");
2738 32234 : for (const auto & [other_node_id, this_node_id] : node_to_node_map)
2739 : {
2740 : // In the case that this==other_mesh, the two nodes might be the same (e.g. if
2741 : // we're stitching a "sliver"), hence we need to skip node deletion in that case.
2742 30530 : if ((this == other_mesh) && (this_node_id == other_node_id))
2743 0 : continue;
2744 :
2745 30530 : this->delete_node( this->node_ptr(this_node_id) );
2746 : }
2747 : }
2748 :
2749 : // If find_neighbors() wasn't called in prepare_for_use(), we need to
2750 : // manually loop once more over all elements adjacent to the stitched boundary
2751 : // and fix their lists of neighbors.
2752 : // This is done according to the following steps:
2753 : // 1. Loop over all copied elements adjacent to the boundary using node_to_elems_map (trying to avoid duplicates)
2754 : // 2. Look at all their sides with a nullptr neighbor and update them using side_to_elem_map if necessary
2755 : // 3. Update the corresponding side in side_to_elem_map as well
2756 1704 : if (skip_find_neighbors)
2757 : {
2758 96 : LOG_SCOPE("stitch_meshes neighbor fixes", "UnstructuredMesh");
2759 :
2760 : // Pull objects out of the loop to reduce heap operations
2761 1704 : std::unique_ptr<const Elem> my_side, their_side;
2762 :
2763 96 : std::set<dof_id_type> fixed_elems;
2764 32234 : for (const auto & pr : node_to_elems_map)
2765 : {
2766 1720 : std::size_t n_elems = pr.second.size();
2767 79875 : for (std::size_t i=0; i<n_elems; i++)
2768 : {
2769 50735 : dof_id_type elem_id = pr.second[i];
2770 1390 : if (!fixed_elems.count(elem_id))
2771 : {
2772 10508 : Elem * el = this->elem_ptr(elem_id);
2773 10212 : fixed_elems.insert(elem_id);
2774 57226 : for (auto s : el->side_index_range())
2775 : {
2776 46718 : bool has_real_neighbor = (el->neighbor_ptr(s) != nullptr);
2777 47570 : bool has_disdjoint_neighbor = is_valid_disjoint_pair_to_stitch &&
2778 852 : (this->get_boundary_info().has_boundary_id(el, s, this_mesh_boundary_id)
2779 852 : || this->get_boundary_info().has_boundary_id(el, s, other_mesh_boundary_id));
2780 :
2781 46718 : if (!has_real_neighbor || has_disdjoint_neighbor)
2782 : {
2783 21016 : key_type key = el->low_order_key(s);
2784 592 : auto bounds = side_to_elem_map.equal_range(key);
2785 :
2786 21016 : if (bounds.first != bounds.second)
2787 : {
2788 : // Get the side for this element
2789 7952 : el->side_ptr(my_side, s);
2790 :
2791 : // Look at all the entries with an equivalent key
2792 7952 : while (bounds.first != bounds.second)
2793 : {
2794 : // Get the potential element
2795 7952 : Elem * neighbor = const_cast<Elem *>(bounds.first->second.first);
2796 :
2797 : // Get the side for the neighboring element
2798 7952 : const unsigned int ns = bounds.first->second.second;
2799 7952 : neighbor->side_ptr(their_side, ns);
2800 : //libmesh_assert(my_side.get());
2801 : //libmesh_assert(their_side.get());
2802 :
2803 : // If found a match with my side
2804 : //
2805 : // We need special tests here for 1D:
2806 : // since parents and children have an equal
2807 : // side (i.e. a node), we need to check
2808 : // ns != ms, and we also check level() to
2809 : // avoid setting our neighbor pointer to
2810 : // any of our neighbor's descendants
2811 15680 : if ((*my_side == *their_side) &&
2812 15904 : (el->level() == neighbor->level()) &&
2813 7952 : ((el->dim() != 1) || (ns != s)))
2814 : {
2815 : // So share a side. Is this a mixed pair
2816 : // of subactive and active/ancestor
2817 : // elements?
2818 : // If not, then we're neighbors.
2819 : // If so, then the subactive's neighbor is
2820 :
2821 8176 : if (el->subactive() ==
2822 7952 : neighbor->subactive())
2823 : {
2824 : // an element is only subactive if it has
2825 : // been coarsened but not deleted
2826 448 : el->set_neighbor (s,neighbor);
2827 448 : neighbor->set_neighbor(ns,el);
2828 : }
2829 0 : else if (el->subactive())
2830 : {
2831 0 : el->set_neighbor(s,neighbor);
2832 : }
2833 0 : else if (neighbor->subactive())
2834 : {
2835 0 : neighbor->set_neighbor(ns,el);
2836 : }
2837 : // It's OK to invalidate the
2838 : // bounds.first iterator here,
2839 : // as we are immediately going
2840 : // to break out of this while
2841 : // loop. bounds.first will
2842 : // therefore not be used for
2843 : // anything else.
2844 224 : side_to_elem_map.erase (bounds.first);
2845 224 : break;
2846 : }
2847 :
2848 0 : ++bounds.first;
2849 : }
2850 : }
2851 : }
2852 : }
2853 : }
2854 : }
2855 : }
2856 1608 : }
2857 :
2858 : #ifdef LIBMESH_ENABLE_PERIODIC
2859 : // Remove only the disjoint pair that was actually stitched.
2860 : // Safe because `is_valid_disjoint_pair_to_stitch` is true
2861 : // only if this exact (a,b) pair exists in the registry.
2862 : // Other disjoint pairs remain untouched.
2863 1704 : if (is_valid_disjoint_pair_to_stitch)
2864 213 : this->remove_disjoint_boundary_pair(this_mesh_boundary_id, other_mesh_boundary_id);
2865 : #endif
2866 :
2867 1704 : if (prepare_after_stitching)
2868 : {
2869 : // We set our new neighbor pointers already
2870 92 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
2871 46 : this->allow_find_neighbors(!skip_find_neighbors);
2872 :
2873 : // We haven't newly remoted any elements
2874 92 : const bool old_allow_remote_element_removal = this->allow_remote_element_removal();
2875 46 : this->allow_remote_element_removal(false);
2876 :
2877 1633 : this->prepare_for_use();
2878 :
2879 46 : this->allow_find_neighbors(old_allow_find_neighbors);
2880 46 : this->allow_remote_element_removal(old_allow_remote_element_removal);
2881 : }
2882 :
2883 : // After the stitching, we may want to clear boundary IDs from element
2884 : // faces that are now internal to the mesh
2885 1704 : if (clear_stitched_boundary_ids)
2886 : {
2887 92 : LOG_SCOPE("stitch_meshes clear bcids", "UnstructuredMesh");
2888 :
2889 1633 : this->get_boundary_info().clear_stitched_boundary_side_ids(
2890 : this_mesh_boundary_id, other_mesh_boundary_id, /*clear_nodeset_data=*/true);
2891 : }
2892 :
2893 : // Return the number of nodes which were merged.
2894 1752 : return node_to_node_map.size();
2895 1742 : }
2896 :
2897 :
2898 : } // namespace libMesh
|
