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