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