Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/boundary_info.h"
22 : #include "libmesh/ghosting_functor.h"
23 : #include "libmesh/ghost_point_neighbors.h"
24 : #include "libmesh/unstructured_mesh.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/elem_range.h"
28 : #include "libmesh/mesh_tools.h" // For n_levels
29 : #include "libmesh/parallel.h"
30 : #include "libmesh/remote_elem.h"
31 : #include "libmesh/namebased_io.h"
32 : #include "libmesh/partitioner.h"
33 : #include "libmesh/enum_order.h"
34 : #include "libmesh/mesh_communication.h"
35 : #include "libmesh/enum_to_string.h"
36 : #include "libmesh/mesh_serializer.h"
37 : #include "libmesh/utility.h"
38 :
39 : #ifdef LIBMESH_HAVE_NANOFLANN
40 : #include "libmesh/nanoflann.hpp"
41 : #endif
42 :
43 : // C++ includes
44 : #include <algorithm> // std::all_of
45 : #include <atomic>
46 : #include <fstream>
47 : #include <iomanip>
48 : #include <map>
49 : #include <sstream>
50 : #include <unordered_map>
51 :
52 : // for disjoint neighbors
53 : #include "libmesh/periodic_boundaries.h"
54 : #include "libmesh/periodic_boundary.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 34384812 : 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 34384812 : hi_elem.n_second_order_adjacent_vertices(hon);
74 :
75 34384812 : std::vector<dof_id_type> adjacent_vertices_ids(n_adjacent_vertices);
76 :
77 115627833 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
78 83401987 : adjacent_vertices_ids[v] =
79 81243021 : 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 34384812 : 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 36194916 : return adj_vertices_to_ho_nodes.try_emplace(adjacent_vertices_ids, nullptr).first;
93 : }
94 :
95 3656873 : 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 94994 : libmesh_assert_equal_to (lo_elem.n_vertices(), hi_elem->n_vertices());
106 :
107 189988 : const processor_id_type my_pid = mesh.processor_id();
108 3656873 : 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 3656873 : const unsigned int hon_begin = lo_elem.n_nodes();
119 3656873 : const unsigned int hon_end = hi_elem->n_nodes();
120 :
121 37756122 : for (unsigned int hon=hon_begin; hon<hon_end; hon++)
122 : {
123 34099249 : auto pos = map_hi_order_node(hon, *hi_elem, adj_vertices_to_ho_nodes);
124 :
125 : // no, not added yet
126 34099249 : if (!pos->second)
127 : {
128 321164 : 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 321164 : Point new_location = 0;
138 43546056 : for (dof_id_type vertex_id : adjacent_vertices_ids)
139 31347154 : new_location += mesh.point(vertex_id);
140 :
141 12520066 : 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 12198902 : (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 12198902 : unique_id_type new_unique_id = max_unique_id +
176 12520066 : max_new_nodes_per_elem * lo_elem.id() +
177 12198902 : hon - hon_begin;
178 :
179 321164 : 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 12198902 : pos->second = hi_node;
189 :
190 12198902 : hi_elem->set_node(hon, hi_node);
191 : }
192 : // yes, already added.
193 : else
194 : {
195 570126 : Node * hi_node = pos->second;
196 570126 : libmesh_assert(hi_node);
197 570126 : libmesh_assert_equal_to(mesh.node_ptr(hi_node->id()), hi_node);
198 :
199 21900347 : 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 21900347 : 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 21900347 : if (!mesh.is_replicated() &&
213 21900347 : hi_node->processor_id() != my_pid &&
214 : chosen_pid == my_pid)
215 4130 : mesh.own_node(*hi_node);
216 :
217 21900347 : 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 18720234 : for (auto s : lo_elem.side_index_range())
228 : {
229 15063361 : Elem * neigh = lo_elem.neighbor_ptr(s);
230 15063361 : if (!neigh)
231 14174856 : continue;
232 :
233 522311 : if (neigh != remote_elem)
234 : {
235 : // We don't support AMR even outside our own range yet.
236 30186 : libmesh_assert_equal_to (neigh->level(), 0);
237 :
238 488414 : const unsigned int ns = neigh->which_neighbor_am_i(&lo_elem);
239 30186 : libmesh_assert_not_equal_to(ns, libMesh::invalid_uint);
240 :
241 60372 : neigh->set_neighbor(ns, hi_elem.get());
242 : }
243 :
244 60692 : 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 3656873 : Elem * interior_p = lo_elem.interior_parent();
254 3656873 : if (interior_p)
255 0 : hi_elem->set_interior_parent(interior_p);
256 :
257 3656873 : if (auto parent_exterior_it = exterior_children_of.find(interior_p);
258 94994 : 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 3751867 : if (auto exterior_it = exterior_children_of.find(&lo_elem);
274 94994 : exterior_it != exterior_children_of.end())
275 : {
276 3656873 : 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 94994 : mesh.get_boundary_info().copy_boundary_ids
295 3656873 : (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 189988 : hi_elem->set_id(lo_elem.id());
303 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
304 189988 : hi_elem->set_unique_id(lo_elem.unique_id());
305 : #endif
306 :
307 3656873 : const unsigned int nei = lo_elem.n_extra_integers();
308 3656873 : hi_elem->add_extra_integers(nei);
309 3657192 : for (unsigned int i=0; i != nei; ++i)
310 319 : hi_elem->set_extra_integer(i, lo_elem.get_extra_integer(i));
311 :
312 3561879 : hi_elem->inherit_data_from(lo_elem);
313 :
314 3751867 : mesh.insert_elem(std::move(hi_elem));
315 3656873 : }
316 :
317 :
318 : template <typename ElemTypeConverter>
319 : void
320 41998 : 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 1198 : 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 1198 : 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 1198 : libmesh_assert(mesh.comm().verify(mesh.n_elem()));
341 1198 : 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 41998 : 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 103107 : auto is_higher_order = [&elem_type_converter](const Elem * elem) {
360 40197 : ElemType old_type = elem->type();
361 2078 : ElemType new_type = elem_type_converter(old_type);
362 40197 : return old_type == new_type;
363 : };
364 :
365 41998 : bool already_higher_order =
366 82798 : std::all_of(range.begin(), range.end(), is_higher_order);
367 :
368 : // Check with other processors and possibly return early
369 41998 : mesh.comm().min(already_higher_order);
370 41998 : 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 2164 : 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 2164 : 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 38014 : 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 38014 : dof_id_type n_unpartitioned_elem = 0,
415 38014 : 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 7310888 : auto track_if_necessary = [&adj_vertices_to_ho_nodes,
432 : &exterior_children_of,
433 193852 : &elem_type_converter](Elem * elem) {
434 3694226 : if (elem && elem != remote_elem)
435 : {
436 3694226 : 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 3694226 : const ElemType old_type = elem->type();
444 131474 : const ElemType new_type = elem_type_converter(old_type);
445 3694226 : if (old_type != new_type)
446 3674602 : 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 114900 : if (range.begin() == mesh.elements_begin() &&
453 139346 : range.end() == mesh.elements_end())
454 : {
455 7148224 : for (auto & elem : range)
456 3651881 : 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 7401906 : for (auto & elem : mesh.element_ptr_range())
479 3843643 : if (auto exterior_map_it = exterior_children_of.find(elem->interior_parent());
480 99044 : 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 7162710 : for (auto & lo_elem : range)
490 : {
491 : // Now we can skip the elements in the range that are already
492 : // higher-order.
493 3657364 : const ElemType old_type = lo_elem->type();
494 129102 : const ElemType new_type = elem_type_converter(old_type);
495 :
496 3657364 : if (old_type == new_type)
497 491 : continue;
498 :
499 : // this does _not_ work for refined elements
500 94994 : libmesh_assert_equal_to (lo_elem->level(), 0);
501 :
502 3656873 : if (lo_elem->processor_id() == DofObject::invalid_processor_id)
503 3555003 : ++n_unpartitioned_elem;
504 : else
505 101870 : ++n_partitioned_elem;
506 :
507 : /*
508 : * Build the higher-order equivalent; add to
509 : * the new_elements list.
510 : */
511 3656873 : auto ho_elem = Elem::build (new_type);
512 :
513 94994 : 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 18879982 : for (unsigned int v=0, lnn=lo_elem->n_nodes(); v < lnn; v++)
520 15628789 : ho_elem->set_node(v, lo_elem->node_ptr(v));
521 :
522 3846861 : 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 1082 : adj_vertices_to_ho_nodes.clear();
532 :
533 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
534 38014 : const unique_id_type new_max_unique_id = max_unique_id +
535 38014 : max_new_nodes_per_elem * mesh.n_elem();
536 38014 : 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 38014 : if (!mesh.is_replicated())
546 : {
547 32866 : dof_id_type max_unpartitioned_elem = n_unpartitioned_elem;
548 32866 : mesh.comm().max(max_unpartitioned_elem);
549 32866 : 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 51204 : if (!mesh.comm().verify(n_unpartitioned_elem) ||
555 51226 : !mesh.comm().verify(n_partitioned_elem) ||
556 25613 : !mesh.is_serial())
557 0 : libmesh_not_implemented();
558 : }
559 : else
560 : {
561 7253 : 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 2164 : const bool old_find_neighbors = mesh.allow_find_neighbors();
570 38014 : if (mesh.is_prepared())
571 260 : mesh.allow_find_neighbors(false);
572 38014 : mesh.prepare_for_use();
573 1082 : 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 333137 : UnstructuredMesh::UnstructuredMesh (const Parallel::Communicator & comm_in,
629 333137 : unsigned char d) :
630 333137 : MeshBase (comm_in,d)
631 : {
632 9834 : libmesh_assert (libMesh::initialized());
633 333137 : }
634 :
635 :
636 :
637 36742 : UnstructuredMesh::UnstructuredMesh (const MeshBase & other_mesh) :
638 36742 : MeshBase (other_mesh)
639 : {
640 13154 : libmesh_assert (libMesh::initialized());
641 36742 : }
642 :
643 :
644 :
645 38517 : 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 : const bool skip_preparation)
657 : {
658 26408 : LOG_SCOPE("copy_nodes_and_elements()", "UnstructuredMesh");
659 :
660 : // If we're asked to skip all preparation, we should be skipping
661 : // find_neighbors specifically.
662 13204 : libmesh_assert(!skip_preparation || skip_find_neighbors);
663 :
664 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
665 52559 : extra_int_maps = this->merge_extra_integer_names(other_mesh);
666 :
667 38517 : const unsigned int n_old_node_ints = extra_int_maps.second.size(),
668 38517 : n_new_node_ints = _node_integer_names.size(),
669 38517 : n_old_elem_ints = extra_int_maps.first.size(),
670 38517 : n_new_elem_ints = _elem_integer_names.size();
671 :
672 : // If we are partitioned into fewer parts than the incoming mesh has
673 : // processors to handle, then we need to "wrap" the other Mesh's
674 : // processor ids to fit within our range. This can happen, for
675 : // example, while stitching meshes with small numbers of elements in
676 : // parallel...
677 14042 : bool wrap_proc_ids = (this->n_processors() <
678 14042 : other_mesh.n_partitions());
679 :
680 : // We're assuming the other mesh has proper element number ordering,
681 : // so that we add parents before their children, and that the other
682 : // mesh is consistently partitioned. We're not assuming that node
683 : // proc ids are topologically consistent, so we don't just
684 : // libmesh_assert_valid_procids.
685 : #ifdef DEBUG
686 13204 : MeshTools::libmesh_assert_valid_amr_elem_ids(other_mesh);
687 13204 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(other_mesh);
688 : #endif
689 :
690 : //Copy in Nodes
691 : {
692 : //Preallocate Memory if necessary
693 38517 : this->reserve_nodes(other_mesh.n_nodes());
694 :
695 10173396 : for (const auto & oldn : other_mesh.node_ptr_range())
696 : {
697 : processor_id_type added_pid = cast_int<processor_id_type>
698 6483342 : (wrap_proc_ids ? oldn->processor_id() % this->n_processors() : oldn->processor_id());
699 :
700 : // Add new nodes in old node Point locations
701 : Node * newn =
702 14504412 : this->add_point(*oldn,
703 6483342 : oldn->id() + node_id_offset,
704 2857118 : added_pid);
705 :
706 6483342 : newn->add_extra_integers(n_new_node_ints);
707 6751350 : for (unsigned int i = 0; i != n_old_node_ints; ++i)
708 279956 : newn->set_extra_integer(extra_int_maps.second[i],
709 268008 : oldn->get_extra_integer(i));
710 :
711 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
712 6483342 : newn->set_unique_id(oldn->unique_id() + unique_id_offset);
713 : #endif
714 24475 : }
715 : }
716 :
717 : //Copy in Elements
718 : {
719 : //Preallocate Memory if necessary
720 38517 : this->reserve_elem(other_mesh.n_elem());
721 :
722 : // Declare a map linking old and new elements, needed to copy the neighbor lists
723 : typedef std::unordered_map<const Elem *, Elem *> map_type;
724 26408 : map_type old_elems_to_new_elems, ip_map;
725 :
726 : // Loop over the elements
727 20923338 : for (const auto & old : other_mesh.element_ptr_range())
728 : {
729 : // Build a new element
730 11729941 : Elem * newparent = old->parent() ?
731 1054665 : this->elem_ptr(old->parent()->id() + element_id_offset) :
732 1663204 : nullptr;
733 13080755 : auto el = old->disconnected_clone();
734 1975594 : el->set_parent(newparent);
735 :
736 11417551 : subdomain_id_type sbd_id = old->subdomain_id();
737 11417551 : if (id_remapping)
738 : {
739 538 : auto remapping_it = id_remapping->find(sbd_id);
740 19099 : if (remapping_it != id_remapping->end())
741 568 : sbd_id = remapping_it->second;
742 : }
743 11417551 : el->subdomain_id() = sbd_id;
744 :
745 : // Hold off on trying to set the interior parent because we may actually
746 : // add lower dimensional elements before their interior parents
747 11417551 : if (old->interior_parent())
748 3008 : ip_map[old] = el.get();
749 :
750 : #ifdef LIBMESH_ENABLE_AMR
751 11417551 : if (old->has_children())
752 1388409 : for (unsigned int c = 0, nc = old->n_children(); c != nc; ++c)
753 1147414 : if (old->child_ptr(c) == remote_elem)
754 68325 : el->add_child(const_cast<RemoteElem *>(remote_elem), c);
755 :
756 : //Create the parent's child pointers if necessary
757 11417551 : if (newparent)
758 : {
759 1079069 : unsigned int oldc = old->parent()->which_child_am_i(old);
760 1054665 : newparent->add_child(el.get(), oldc);
761 : }
762 :
763 : // Copy the refinement flags
764 11417551 : el->set_refinement_flag(old->refinement_flag());
765 :
766 : // Use hack_p_level since we may not have sibling elements
767 : // added yet
768 1975594 : el->hack_p_level(old->p_level());
769 :
770 1975594 : el->set_p_refinement_flag(old->p_refinement_flag());
771 : #endif // #ifdef LIBMESH_ENABLE_AMR
772 :
773 : //Assign all the nodes
774 64040941 : for (auto i : el->node_index_range())
775 63760308 : el->set_node(i,
776 52623390 : this->node_ptr(old->node_id(i) + node_id_offset));
777 :
778 : // And start it off with the same processor id (mod _n_parts).
779 11417551 : el->processor_id() = cast_int<processor_id_type>
780 11417551 : (wrap_proc_ids ? old->processor_id() % this->n_processors() : old->processor_id());
781 :
782 : // Give it the same element and unique ids
783 11417551 : el->set_id(old->id() + element_id_offset);
784 :
785 11417551 : el->add_extra_integers(n_new_elem_ints);
786 11440052 : for (unsigned int i = 0; i != n_old_elem_ints; ++i)
787 28721 : el->set_extra_integer(extra_int_maps.first[i],
788 22501 : old->get_extra_integer(i));
789 :
790 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
791 11417551 : el->set_unique_id(old->unique_id() + unique_id_offset);
792 : #endif
793 :
794 : //Hold onto it
795 11417551 : if (!skip_find_neighbors)
796 : {
797 220630 : for (auto s : old->side_index_range())
798 182136 : if (old->neighbor_ptr(s) == remote_elem)
799 256 : el->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
800 46942 : this->add_elem(std::move(el));
801 : }
802 : else
803 : {
804 11684407 : Elem * new_el = this->add_elem(std::move(el));
805 11373425 : old_elems_to_new_elems[old] = new_el;
806 : }
807 9466432 : }
808 :
809 : // If the other_mesh had some interior parents, we may need to
810 : // copy those pointers (if they're to elements in a third mesh),
811 : // or create new equivalent pointers (if they're to elements we
812 : // just copied), or scream and die (if the other mesh had interior
813 : // parents from a third mesh but we already have interior parents
814 : // that aren't to that same third mesh.
815 38517 : if (!ip_map.empty())
816 : {
817 298 : std::atomic<bool> existing_interior_parents{false};
818 :
819 : Threads::parallel_for
820 298 : (this->element_stored_range(),
821 610 : [&existing_interior_parents](const ElemRange & range)
822 : {
823 10582 : for (Elem * elem : range)
824 10284 : if (elem->interior_parent())
825 : {
826 0 : existing_interior_parents = true;
827 0 : break;
828 : }
829 298 : });
830 :
831 : MeshBase * other_interior_mesh =
832 142 : const_cast<MeshBase *>(&other_mesh.interior_mesh());
833 :
834 : // If we don't already have interior parents, then we can just
835 : // use whatever interior_mesh we need for the incoming
836 : // elements.
837 298 : if (!existing_interior_parents)
838 : {
839 298 : if (other_interior_mesh == &other_mesh)
840 62 : this->set_interior_mesh(*this);
841 : else
842 80 : this->set_interior_mesh(*other_interior_mesh);
843 : }
844 :
845 298 : if (other_interior_mesh == &other_mesh &&
846 218 : _interior_mesh == this)
847 1538 : for (auto & elem_pair : ip_map)
848 1320 : elem_pair.second->set_interior_parent(
849 1382 : this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset));
850 80 : else if (other_interior_mesh == _interior_mesh)
851 1768 : for (auto & elem_pair : ip_map)
852 : {
853 1688 : Elem * ip = const_cast<Elem *>(elem_pair.first->interior_parent());
854 1688 : libmesh_assert(ip == remote_elem ||
855 : ip == other_interior_mesh->elem_ptr(ip->id()));
856 1688 : elem_pair.second->set_interior_parent(ip);
857 : }
858 : else
859 0 : libmesh_error_msg("Cannot copy boundary elements between meshes with different interior meshes");
860 : }
861 :
862 : // Loop (again) over the elements to fill in the neighbors
863 38517 : if (skip_find_neighbors)
864 : {
865 38233 : old_elems_to_new_elems[remote_elem] = const_cast<RemoteElem*>(remote_elem);
866 :
867 20837342 : for (const auto & old_elem : other_mesh.element_ptr_range())
868 : {
869 11373425 : Elem * new_elem = old_elems_to_new_elems[old_elem];
870 57096052 : for (auto s : old_elem->side_index_range())
871 : {
872 46646885 : const Elem * old_neighbor = old_elem->neighbor_ptr(s);
873 45411645 : Elem * new_neighbor = old_elems_to_new_elems[old_neighbor];
874 7919292 : new_elem->set_neighbor(s, new_neighbor);
875 : }
876 24207 : }
877 : }
878 : }
879 :
880 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
881 : // We set the unique ids of nodes after adding them to the mesh such that our value of
882 : // _next_unique_id may be wrong. So we amend that here
883 38517 : this->set_next_unique_id(other_mesh.parallel_max_unique_id() + unique_id_offset + 1);
884 : #endif
885 :
886 : // Finally, partially prepare the new Mesh for use, if that isn't
887 : // being skipped.
888 : // Even the default behavior here is for backwards compatibility,
889 : // and we don't want to prepare everything.
890 :
891 38517 : if (!skip_preparation)
892 : {
893 : // Keep the same numbering and partitioning and distribution
894 : // status for now, but save our original policies to restore
895 : // later.
896 100 : const bool allowed_renumbering = this->allow_renumbering();
897 100 : const bool allowed_find_neighbors = this->allow_find_neighbors();
898 100 : const bool allowed_elem_removal = this->allow_remote_element_removal();
899 100 : const bool allowed_detect_detect_interior_parents = this->allow_detect_interior_parents();
900 50 : this->allow_renumbering(false);
901 50 : this->allow_remote_element_removal(false);
902 50 : this->allow_find_neighbors(!skip_find_neighbors);
903 100 : this->allow_detect_interior_parents(other_mesh.allow_detect_interior_parents());
904 :
905 : // We should generally be able to skip *all* partitioning here
906 : // because we're only adding one already-consistent mesh to
907 : // another.
908 100 : const bool skipped_partitioning = this->skip_partitioning();
909 50 : this->skip_partitioning(true);
910 :
911 50 : const Preparation old_preparation = this->preparation();
912 1775 : this->prepare_for_use();
913 :
914 : //But in the long term, don't change our policies.
915 50 : this->allow_find_neighbors(allowed_find_neighbors);
916 50 : this->allow_renumbering(allowed_renumbering);
917 50 : this->allow_remote_element_removal(allowed_elem_removal);
918 50 : this->skip_partitioning(skipped_partitioning);
919 50 : this->allow_detect_interior_parents(allowed_detect_detect_interior_parents);
920 :
921 : // That prepare_for_use() call marked us as prepared, but we
922 : // specifically avoided some important preparation, so we might not
923 : // actually be prepared now.
924 1775 : if (skip_find_neighbors)
925 42 : this->unset_has_neighbor_ptrs();
926 :
927 50 : const Preparation other_preparation = other_mesh.preparation();
928 1775 : if (!old_preparation.is_partitioned ||
929 1775 : !other_preparation.is_partitioned)
930 0 : this->unset_is_partitioned();
931 1775 : if (!old_preparation.has_removed_orphaned_nodes ||
932 1775 : !other_preparation.has_removed_orphaned_nodes)
933 0 : this->unset_has_removed_orphaned_nodes();
934 1775 : if (!old_preparation.has_removed_remote_elements ||
935 1775 : !other_preparation.has_removed_remote_elements)
936 0 : this->unset_has_removed_remote_elements();
937 : }
938 :
939 : // In general we've just invalidated just about everything, and we'd
940 : // like to unset_is_prepared(), but specific use cases might know a
941 : // priori that they're still partitioned well, or that they've
942 : // copied in a disjoint mesh component and don't need new neighbor
943 : // pointers, or that they're not adding anything that would change
944 : // cached subdomain/element/boundary sets, etc., so we'll rely on
945 : // users of the "advanced" skip_preparation option to also set what
946 : // preparation they still need.
947 :
948 : // else
949 : // this->unset_is_prepared();
950 38517 : }
951 :
952 :
953 :
954 369879 : UnstructuredMesh::~UnstructuredMesh ()
955 : {
956 : // this->clear (); // Nothing to clear at this level
957 :
958 22988 : libmesh_exceptionless_assert (!libMesh::closed());
959 369879 : }
960 :
961 :
962 :
963 :
964 :
965 618160 : void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
966 : const bool reset_current_list,
967 : const bool assert_valid,
968 : const bool check_non_remote)
969 : {
970 : // We might actually want to run this on an empty mesh
971 : // (e.g. the boundary mesh for a nonexistent bcid!)
972 : // libmesh_assert_not_equal_to (this->n_nodes(), 0);
973 : // libmesh_assert_not_equal_to (this->n_elem(), 0);
974 :
975 : // This function must be run on all processors at once
976 25932 : parallel_object_only();
977 :
978 25932 : LOG_SCOPE("find_neighbors()", "Mesh");
979 :
980 : //TODO:[BSK] This should be removed later?!
981 618160 : if (reset_current_list)
982 : Threads::parallel_for
983 549613 : (this->element_stored_range(),
984 1578385 : [reset_remote_elements](const ElemRange & range)
985 : {
986 71394813 : for (Elem * e : range)
987 354428884 : for (auto s : e->side_index_range())
988 289262890 : if (e->neighbor_ptr(s) != remote_elem || reset_remote_elements)
989 11126548 : e->set_neighbor(s, nullptr);
990 550119 : });
991 :
992 : // Find neighboring elements by first finding elements
993 : // with identical side keys and then check to see if they
994 : // are neighbors
995 : {
996 : // data structures -- Use the hash_multimap if available
997 : typedef dof_id_type key_type;
998 : typedef std::pair<Elem *, unsigned char> val_type;
999 : typedef std::unordered_multimap<key_type, val_type> map_type;
1000 :
1001 : // A map from side keys to corresponding elements & side numbers
1002 51864 : map_type side_to_elem_map;
1003 :
1004 : // Pull objects out of the loop to reduce heap operations
1005 618160 : std::unique_ptr<Elem> my_side, their_side;
1006 :
1007 150557850 : for (const auto & element : this->element_ptr_range())
1008 : {
1009 387233293 : for (auto ms : element->side_index_range())
1010 : {
1011 440093220 : next_side:
1012 : // If we haven't yet found a neighbor on this side, try.
1013 : // Even if we think our neighbor is remote, that
1014 : // information may be out of date.
1015 : //
1016 : // If we're only checking remote neighbors, after a
1017 : // redistribution, then we'll skip the non-remote ones
1018 453564381 : if ((element->neighbor_ptr(ms) == nullptr && check_non_remote) ||
1019 158604005 : element->neighbor_ptr(ms) == remote_elem)
1020 : {
1021 : // Get the key for the side of this element. Use the
1022 : // low_order_key so we can find neighbors in
1023 : // mixed-order meshes if necessary.
1024 285344085 : const dof_id_type key = element->low_order_key(ms);
1025 :
1026 : // Look for elements that have an identical side key
1027 11146034 : auto bounds = side_to_elem_map.equal_range(key);
1028 :
1029 : // May be multiple keys, check all the possible
1030 : // elements which _might_ be neighbors.
1031 285344085 : if (bounds.first != bounds.second)
1032 : {
1033 : // Get the side for this element
1034 131152723 : element->side_ptr(my_side, ms);
1035 :
1036 : // Look at all the entries with an equivalent key
1037 131360750 : while (bounds.first != bounds.second)
1038 : {
1039 : // Get the potential element
1040 131271206 : Elem * neighbor = bounds.first->second.first;
1041 :
1042 : // Get the side for the neighboring element
1043 131271206 : const unsigned int ns = bounds.first->second.second;
1044 131271206 : neighbor->side_ptr(their_side, ns);
1045 : //libmesh_assert(my_side.get());
1046 : //libmesh_assert(their_side.get());
1047 :
1048 : // If found a match with my side
1049 : //
1050 : // In 1D, since parents and children have an
1051 : // equal side (i.e. a node) we need to check
1052 : // for matching level() to avoid setting our
1053 : // neighbor pointer to any of our neighbor's
1054 : // descendants.
1055 262542412 : if ((*my_side == *their_side) &&
1056 131271206 : (element->level() == neighbor->level()))
1057 : {
1058 : // So share a side. Is this a mixed pair
1059 : // of subactive and active/ancestor
1060 : // elements?
1061 : // If not, then we're neighbors.
1062 : // If so, then the subactive's neighbor is
1063 :
1064 133652677 : if (element->subactive() ==
1065 131063179 : neighbor->subactive())
1066 : {
1067 : // an element is only subactive if it has
1068 : // been coarsened but not deleted
1069 7645841 : element->set_neighbor (ms,neighbor);
1070 130882242 : neighbor->set_neighbor(ns,element);
1071 : }
1072 180937 : else if (element->subactive())
1073 : {
1074 7644 : element->set_neighbor(ms,neighbor);
1075 : }
1076 99221 : else if (neighbor->subactive())
1077 : {
1078 14052 : neighbor->set_neighbor(ns,element);
1079 : }
1080 5078039 : side_to_elem_map.erase (bounds.first);
1081 :
1082 : // get out of this nested crap
1083 131063179 : goto next_side;
1084 : }
1085 :
1086 5860 : ++bounds.first;
1087 : }
1088 : }
1089 :
1090 : // didn't find a match...
1091 : // Build the map entry for this element
1092 : side_to_elem_map.emplace
1093 154280906 : (key, std::make_pair(element, cast_int<unsigned char>(ms)));
1094 : }
1095 : }
1096 578830 : }
1097 578830 : }
1098 :
1099 : #ifdef LIBMESH_ENABLE_PERIODIC
1100 : // Get the disjoint neighbor boundary pairs object (from periodic BCs)
1101 618160 : auto * db = this->get_disjoint_neighbor_boundary_pairs();
1102 :
1103 618160 : if (db)
1104 : {
1105 : // Obtain a point locator
1106 1085 : std::unique_ptr<PointLocatorBase> point_locator = this->sub_point_locator();
1107 :
1108 7202 : for (const auto & element : this->element_ptr_range())
1109 : {
1110 13322 : for (auto ms : element->side_index_range())
1111 : {
1112 : // Skip if this side already has a valid neighbor (including remote neighbors)
1113 10912 : if (element->neighbor_ptr(ms) != nullptr &&
1114 2047 : element->neighbor_ptr(ms) != remote_elem)
1115 2003 : continue;
1116 :
1117 32779 : for (const auto & [id, boundary_ptr] : *db)
1118 : {
1119 24210 : if (!this->get_boundary_info().has_boundary_id(element, ms, id))
1120 22747 : continue;
1121 :
1122 : unsigned int neigh_side;
1123 : const Elem * neigh =
1124 1491 : db->neighbor(id, *point_locator, element, ms, &neigh_side);
1125 :
1126 1463 : if (neigh && neigh != remote_elem && neigh != element)
1127 : {
1128 1463 : auto neigh_changeable = this->elem_ptr(neigh->id());
1129 1463 : element->set_neighbor(ms, neigh_changeable);
1130 1463 : neigh_changeable->set_neighbor(neigh_side, element);
1131 : }
1132 : }
1133 : }
1134 985 : }
1135 985 : }
1136 : #endif // LIBMESH_ENABLE_PERIODIC
1137 :
1138 : #ifdef LIBMESH_ENABLE_AMR
1139 :
1140 : /**
1141 : * Here we look at all of the child elements which
1142 : * don't already have valid neighbors.
1143 : *
1144 : * If a child element has a nullptr neighbor it is
1145 : * either because it is on the boundary or because
1146 : * its neighbor is at a different level. In the
1147 : * latter case we must get the neighbor from the
1148 : * parent.
1149 : *
1150 : * If a child element has a remote_elem neighbor
1151 : * on a boundary it shares with its parent, that
1152 : * info may have become out-dated through coarsening
1153 : * of the neighbor's parent. In this case, if the
1154 : * parent's neighbor is active then the child should
1155 : * share it.
1156 : *
1157 : * Furthermore, that neighbor better be active,
1158 : * otherwise we missed a child somewhere.
1159 : *
1160 : *
1161 : * We also need to look through children ordered by increasing
1162 : * refinement level in order to add new interior_parent() links in
1163 : * boundary elements which have just been generated by refinement,
1164 : * and fix links in boundary elements whose previous
1165 : * interior_parent() has just been coarsened away.
1166 : */
1167 618160 : const unsigned int n_levels = MeshTools::n_levels(*this);
1168 897176 : for (unsigned int level = 1; level < n_levels; ++level)
1169 : {
1170 2194256 : for (auto & current_elem : as_range(level_elements_begin(level),
1171 101913478 : level_elements_end(level)))
1172 : {
1173 1645192 : libmesh_assert(current_elem);
1174 50823296 : Elem * parent = current_elem->parent();
1175 1645192 : libmesh_assert(parent);
1176 50823296 : const unsigned int my_child_num = parent->which_child_am_i(current_elem);
1177 :
1178 252916484 : for (auto s : current_elem->side_index_range())
1179 : {
1180 210513968 : if (current_elem->neighbor_ptr(s) == nullptr ||
1181 193620418 : (current_elem->neighbor_ptr(s) == remote_elem &&
1182 1858887 : parent->is_child_on_side(my_child_num, s)))
1183 : {
1184 698312 : Elem * neigh = parent->neighbor_ptr(s);
1185 :
1186 : // If neigh was refined and had non-subactive children
1187 : // made remote earlier, then our current elem should
1188 : // actually have one of those remote children as a
1189 : // neighbor
1190 15074141 : if (neigh &&
1191 4476743 : (neigh->ancestor() ||
1192 : // If neigh has subactive children which should have
1193 : // matched as neighbors of the current element but
1194 : // did not, then those likewise must be remote
1195 : // children.
1196 4289740 : (current_elem->subactive() && neigh->has_children() &&
1197 186 : (neigh->level()+1) == current_elem->level())))
1198 : {
1199 : #ifdef DEBUG
1200 : // Let's make sure that "had children made remote"
1201 : // situation is actually the case
1202 0 : libmesh_assert(neigh->has_children());
1203 0 : bool neigh_has_remote_children = false;
1204 0 : for (auto & child : neigh->child_ref_range())
1205 0 : if (&child == remote_elem)
1206 0 : neigh_has_remote_children = true;
1207 0 : libmesh_assert(neigh_has_remote_children);
1208 :
1209 : // And let's double-check that we don't have
1210 : // a remote_elem neighboring an active local element
1211 0 : if (current_elem->active())
1212 0 : libmesh_assert_not_equal_to (current_elem->processor_id(),
1213 : this->processor_id());
1214 : #endif // DEBUG
1215 78681 : neigh = const_cast<RemoteElem *>(remote_elem);
1216 : }
1217 : // If neigh and current_elem are more than one level
1218 : // apart, figuring out whether we have a remote
1219 : // neighbor here becomes much harder.
1220 10643211 : else if (neigh && (current_elem->subactive() &&
1221 15800 : neigh->has_children()))
1222 : {
1223 : // Find the deepest descendant of neigh which
1224 : // we could consider for a neighbor. If we run
1225 : // out of neigh children, then that's our
1226 : // neighbor. If we find a potential neighbor
1227 : // with remote_children and we don't find any
1228 : // potential neighbors among its non-remote
1229 : // children, then our neighbor must be remote.
1230 0 : while (neigh != remote_elem &&
1231 0 : neigh->has_children())
1232 : {
1233 0 : bool found_neigh = false;
1234 0 : for (unsigned int c = 0, nc = neigh->n_children();
1235 0 : !found_neigh && c != nc; ++c)
1236 : {
1237 0 : Elem * child = neigh->child_ptr(c);
1238 0 : if (child == remote_elem)
1239 0 : continue;
1240 0 : for (auto ncn : child->neighbor_ptr_range())
1241 : {
1242 0 : if (ncn != remote_elem &&
1243 0 : ncn->is_ancestor_of(current_elem))
1244 : {
1245 0 : neigh = ncn;
1246 0 : found_neigh = true;
1247 0 : break;
1248 : }
1249 : }
1250 : }
1251 0 : if (!found_neigh)
1252 0 : neigh = const_cast<RemoteElem *>(remote_elem);
1253 : }
1254 : }
1255 10706092 : current_elem->set_neighbor(s, neigh);
1256 : #ifdef DEBUG
1257 463912 : if (neigh != nullptr && neigh != remote_elem)
1258 : // We ignore subactive elements here because
1259 : // we don't care about neighbors of subactive element.
1260 214304 : if ((!neigh->active()) && (!current_elem->subactive()))
1261 : {
1262 0 : libMesh::err << "On processor " << this->processor_id()
1263 0 : << std::endl;
1264 0 : libMesh::err << "Bad element ID = " << current_elem->id()
1265 0 : << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
1266 0 : libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
1267 0 : << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
1268 0 : libMesh::err << "Bad element size = " << current_elem->hmin()
1269 0 : << ", Bad neighbor size = " << neigh->hmin() << std::endl;
1270 0 : libMesh::err << "Bad element center = " << current_elem->vertex_average()
1271 0 : << ", Bad neighbor center = " << neigh->vertex_average() << std::endl;
1272 0 : libMesh::err << "ERROR: "
1273 0 : << (current_elem->active()?"Active":"Ancestor")
1274 0 : << " Element at level "
1275 0 : << current_elem->level() << std::endl;
1276 0 : libMesh::err << "with "
1277 0 : << (parent->active()?"active":
1278 0 : (parent->subactive()?"subactive":"ancestor"))
1279 0 : << " parent share "
1280 0 : << (neigh->subactive()?"subactive":"ancestor")
1281 0 : << " neighbor at level " << neigh->level()
1282 0 : << std::endl;
1283 0 : NameBasedIO(*this).write ("bad_mesh.gmv");
1284 0 : libmesh_error_msg("Problematic mesh written to bad_mesh.gmv.");
1285 : }
1286 : #endif // DEBUG
1287 : }
1288 : }
1289 :
1290 : // We can skip to the next element if we're full-dimension
1291 : // and therefore don't have any interior parents
1292 50823296 : if (current_elem->dim() >= LIBMESH_DIM)
1293 8453681 : continue;
1294 :
1295 : // We have no interior parents unless we can find one later
1296 43403479 : current_elem->set_interior_parent(nullptr);
1297 :
1298 43403479 : Elem * pip = parent->interior_parent();
1299 :
1300 43403479 : if (!pip)
1301 42752618 : continue;
1302 :
1303 : // If there's no interior_parent children, whether due to a
1304 : // remote element or a non-conformity, then there's no
1305 : // children to search.
1306 25419 : if (pip == remote_elem || pip->active())
1307 : {
1308 2052 : current_elem->set_interior_parent(pip);
1309 2052 : continue;
1310 : }
1311 :
1312 : // For node comparisons we'll need a sensible tolerance
1313 23367 : Real node_tolerance = current_elem->hmin() * TOLERANCE;
1314 :
1315 : // Otherwise our interior_parent should be a child of our
1316 : // parent's interior_parent.
1317 77739 : for (auto & child : pip->child_ref_range())
1318 : {
1319 : // If we have a remote_elem, that might be our
1320 : // interior_parent. We'll set it provisionally now and
1321 : // keep trying to find something better.
1322 76601 : if (&child == remote_elem)
1323 : {
1324 : current_elem->set_interior_parent
1325 4788 : (const_cast<RemoteElem *>(remote_elem));
1326 4788 : continue;
1327 : }
1328 :
1329 5860 : bool child_contains_our_nodes = true;
1330 145312 : for (auto & n : current_elem->node_ref_range())
1331 : {
1332 10176 : bool child_contains_this_node = false;
1333 795528 : for (auto & cn : child.node_ref_range())
1334 745944 : if (cn.absolute_fuzzy_equals
1335 714744 : (n, node_tolerance))
1336 : {
1337 6268 : child_contains_this_node = true;
1338 6268 : break;
1339 : }
1340 120075 : if (!child_contains_this_node)
1341 : {
1342 3908 : child_contains_our_nodes = false;
1343 3908 : break;
1344 : }
1345 : }
1346 71813 : if (child_contains_our_nodes)
1347 : {
1348 22229 : current_elem->set_interior_parent(&child);
1349 1952 : break;
1350 : }
1351 : }
1352 :
1353 : // We should have found *some* interior_parent at this
1354 : // point, whether semilocal or remote.
1355 1952 : libmesh_assert(current_elem->interior_parent());
1356 265498 : }
1357 : }
1358 : #endif // AMR
1359 :
1360 : #ifdef DEBUG
1361 25932 : if (assert_valid)
1362 : {
1363 25796 : MeshTools::libmesh_assert_valid_neighbors(*this,
1364 25796 : !reset_remote_elements);
1365 25796 : MeshTools::libmesh_assert_valid_amr_interior_parents(*this);
1366 : }
1367 : #else
1368 : libmesh_ignore(assert_valid);
1369 : #endif
1370 :
1371 618160 : this->_preparation.has_neighbor_ptrs = true;
1372 618160 : }
1373 :
1374 :
1375 :
1376 5634 : void UnstructuredMesh::read (const std::string & name,
1377 : void *,
1378 : bool skip_renumber_nodes_and_elements,
1379 : bool skip_find_neighbors,
1380 : bool skip_detect_interior_parents)
1381 : {
1382 : // Set the skip_renumber_nodes_and_elements flag on all processors
1383 : // if necessary.
1384 : // This ensures that renumber_nodes_and_elements is *not* called
1385 : // during prepare_for_use() for certain types of mesh files.
1386 : // This is required in cases where there is an associated solution
1387 : // file which expects a certain ordering of the nodes.
1388 5634 : if (Utility::ends_with(name, ".gmv"))
1389 0 : this->allow_renumbering(false);
1390 :
1391 5634 : NameBasedIO(*this).read(name);
1392 :
1393 5634 : if (skip_renumber_nodes_and_elements)
1394 : {
1395 : // Use MeshBase::allow_renumbering() yourself instead.
1396 : libmesh_deprecated();
1397 0 : this->allow_renumbering(false);
1398 : }
1399 :
1400 : // Done reading the mesh. Now prepare it for use.
1401 348 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
1402 348 : const bool old_allow_detect_interior_parents = this->allow_detect_interior_parents();
1403 :
1404 174 : this->allow_find_neighbors(!skip_find_neighbors);
1405 174 : this->allow_detect_interior_parents(!skip_detect_interior_parents);
1406 :
1407 5634 : this->prepare_for_use();
1408 :
1409 174 : this->allow_find_neighbors(old_allow_find_neighbors);
1410 174 : this->allow_detect_interior_parents(old_allow_detect_interior_parents);
1411 5634 : }
1412 :
1413 :
1414 :
1415 2938 : void UnstructuredMesh::write (const std::string & name) const
1416 : {
1417 598 : LOG_SCOPE("write()", "Mesh");
1418 :
1419 3536 : NameBasedIO(*this).write(name);
1420 2938 : }
1421 :
1422 :
1423 :
1424 0 : void UnstructuredMesh::write (const std::string & name,
1425 : const std::vector<Number> & v,
1426 : const std::vector<std::string> & vn) const
1427 : {
1428 0 : LOG_SCOPE("write()", "Mesh");
1429 :
1430 0 : NameBasedIO(*this).write_nodal_data(name, v, vn);
1431 0 : }
1432 :
1433 :
1434 :
1435 :
1436 :
1437 0 : void UnstructuredMesh::create_pid_mesh(UnstructuredMesh & pid_mesh,
1438 : const processor_id_type pid) const
1439 : {
1440 :
1441 : // Issue a warning if the number the number of processors
1442 : // currently available is less that that requested for
1443 : // partitioning. This is not necessarily an error since
1444 : // you may run on one processor and still partition the
1445 : // mesh into several partitions.
1446 : #ifdef DEBUG
1447 0 : if (this->n_processors() < pid)
1448 : {
1449 0 : libMesh::out << "WARNING: You are creating a "
1450 0 : << "mesh for a processor id (="
1451 0 : << pid
1452 0 : << ") greater than "
1453 0 : << "the number of processors available for "
1454 0 : << "the calculation. (="
1455 0 : << this->n_processors()
1456 0 : << ")."
1457 0 : << std::endl;
1458 : }
1459 : #endif
1460 :
1461 0 : this->create_submesh (pid_mesh,
1462 0 : this->active_pid_elements_begin(pid),
1463 0 : this->active_pid_elements_end(pid));
1464 0 : }
1465 :
1466 :
1467 :
1468 :
1469 :
1470 :
1471 :
1472 0 : void UnstructuredMesh::create_submesh (UnstructuredMesh & new_mesh,
1473 : const const_element_iterator & it,
1474 : const const_element_iterator & it_end) const
1475 : {
1476 : // Just in case the subdomain_mesh already has some information
1477 : // in it, get rid of it.
1478 0 : new_mesh.clear();
1479 :
1480 : // If we're not serial, our submesh isn't either.
1481 : // There are no remote elements to delete on an empty mesh, but
1482 : // calling the method to do so marks the mesh as parallel.
1483 0 : if (!this->is_serial())
1484 0 : new_mesh.delete_remote_elements();
1485 :
1486 : // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
1487 : // This may happen if the user accidentally passes the original mesh into
1488 : // this function! We will check this by making sure we did not just
1489 : // clear ourself.
1490 0 : libmesh_assert_not_equal_to (this->n_nodes(), 0);
1491 0 : libmesh_assert_not_equal_to (this->n_elem(), 0);
1492 :
1493 : // Container to catch boundary IDs handed back by BoundaryInfo
1494 0 : std::vector<boundary_id_type> bc_ids;
1495 :
1496 : // Put any extra integers on the new mesh too
1497 0 : new_mesh.merge_extra_integer_names(*this);
1498 0 : const unsigned int n_node_ints = _node_integer_names.size();
1499 :
1500 0 : for (const auto & old_elem : as_range(it, it_end))
1501 : {
1502 : // Add an equivalent element type to the new_mesh.
1503 : // disconnected_clone() copies ids, extra element integers, etc.
1504 0 : auto uelem = old_elem->disconnected_clone();
1505 0 : Elem * new_elem = new_mesh.add_elem(std::move(uelem));
1506 0 : libmesh_assert(new_elem);
1507 :
1508 : // Loop over the nodes on this element.
1509 0 : for (auto n : old_elem->node_index_range())
1510 : {
1511 0 : const dof_id_type this_node_id = old_elem->node_id(n);
1512 :
1513 : // Add this node to the new mesh if it's not there already
1514 0 : if (!new_mesh.query_node_ptr(this_node_id))
1515 : {
1516 : Node * newn =
1517 0 : new_mesh.add_point (old_elem->point(n),
1518 : this_node_id,
1519 0 : old_elem->node_ptr(n)->processor_id());
1520 :
1521 0 : newn->add_extra_integers(n_node_ints);
1522 0 : for (unsigned int i = 0; i != n_node_ints; ++i)
1523 0 : newn->set_extra_integer(i, old_elem->node_ptr(n)->get_extra_integer(i));
1524 :
1525 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1526 0 : newn->set_unique_id(old_elem->node_ptr(n)->unique_id());
1527 : #endif
1528 : }
1529 :
1530 : // Define this element's connectivity on the new mesh
1531 0 : new_elem->set_node(n, new_mesh.node_ptr(this_node_id));
1532 : }
1533 :
1534 : // Maybe add boundary conditions for this element
1535 0 : for (auto s : old_elem->side_index_range())
1536 : {
1537 0 : this->get_boundary_info().boundary_ids(old_elem, s, bc_ids);
1538 0 : new_mesh.get_boundary_info().add_side (new_elem, s, bc_ids);
1539 : }
1540 0 : } // end loop over elements
1541 :
1542 : // Prepare the new_mesh for use
1543 0 : new_mesh.prepare_for_use();
1544 0 : }
1545 :
1546 :
1547 :
1548 : #ifdef LIBMESH_ENABLE_AMR
1549 25826 : bool UnstructuredMesh::contract ()
1550 : {
1551 858 : LOG_SCOPE ("contract()", "Mesh");
1552 :
1553 : // Flag indicating if this call actually changes the mesh
1554 858 : bool mesh_changed = false;
1555 :
1556 : #ifdef DEBUG
1557 545312 : for (const auto & elem : this->element_ptr_range())
1558 544454 : libmesh_assert(elem->active() || elem->subactive() || elem->ancestor());
1559 : #endif
1560 :
1561 : // Loop over the elements.
1562 20049984 : for (auto & elem : this->element_ptr_range())
1563 : {
1564 : // Delete all the subactive ones
1565 10544049 : if (elem->subactive())
1566 : {
1567 : // No level-0 element should be subactive.
1568 : // Note that we CAN'T test elem->level(), as that
1569 : // touches elem->parent()->dim(), and elem->parent()
1570 : // might have already been deleted!
1571 71776 : libmesh_assert(elem->parent());
1572 :
1573 : // Delete the element
1574 : // This just sets a pointer to nullptr, and doesn't
1575 : // invalidate any iterators
1576 1498518 : this->delete_elem(elem);
1577 :
1578 : // the mesh has certainly changed
1579 71776 : mesh_changed = true;
1580 : }
1581 : else
1582 : {
1583 : // Compress all the active ones
1584 472678 : if (elem->active())
1585 6776769 : elem->contract();
1586 : else
1587 113362 : libmesh_assert (elem->ancestor());
1588 : }
1589 24110 : }
1590 :
1591 : // Strip any newly-created nullptr voids out of the element array
1592 25826 : this->renumber_nodes_and_elements();
1593 :
1594 : // FIXME: Need to understand why deleting subactive children
1595 : // invalidates the point locator. For now we will clear it explicitly
1596 25826 : this->clear_point_locator();
1597 :
1598 : // Allow our GhostingFunctor objects to reinit if necessary.
1599 27788 : for (auto & gf : as_range(this->ghosting_functors_begin(),
1600 116257 : this->ghosting_functors_end()))
1601 : {
1602 2820 : libmesh_assert(gf);
1603 86753 : gf->mesh_reinit();
1604 : }
1605 :
1606 26684 : return mesh_changed;
1607 : }
1608 : #endif // #ifdef LIBMESH_ENABLE_AMR
1609 :
1610 :
1611 :
1612 11259 : void UnstructuredMesh::all_first_order ()
1613 : {
1614 784 : LOG_SCOPE("all_first_order()", "Mesh");
1615 :
1616 : /**
1617 : * Prepare to identify (and then delete) a bunch of no-longer-used nodes.
1618 : */
1619 11651 : std::vector<bool> node_touched_by_me(this->max_node_id(), false);
1620 :
1621 : // Loop over the high-ordered elements.
1622 : // First make sure they _are_ indeed high-order, and then replace
1623 : // them with an equivalent first-order element.
1624 481878 : for (auto & so_elem : element_ptr_range())
1625 : {
1626 25848 : libmesh_assert(so_elem);
1627 :
1628 : /*
1629 : * build the first-order equivalent, add to
1630 : * the new_elements list.
1631 : */
1632 : auto lo_elem = Elem::build
1633 : (Elem::first_order_equivalent_type
1634 281572 : (so_elem->type()), so_elem->parent());
1635 :
1636 255724 : const unsigned short n_sides = so_elem->n_sides();
1637 :
1638 1220978 : for (unsigned short s=0; s != n_sides; ++s)
1639 1063058 : if (so_elem->neighbor_ptr(s) == remote_elem)
1640 0 : lo_elem->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1641 :
1642 : #ifdef LIBMESH_ENABLE_AMR
1643 : /*
1644 : * Reset the parent links of any child elements
1645 : */
1646 255724 : if (so_elem->has_children())
1647 376697 : for (unsigned int c = 0, nc = so_elem->n_children(); c != nc; ++c)
1648 : {
1649 295596 : Elem * child = so_elem->child_ptr(c);
1650 295596 : if (child != remote_elem)
1651 48208 : child->set_parent(lo_elem.get());
1652 295596 : lo_elem->add_child(child, c);
1653 : }
1654 :
1655 : /*
1656 : * Reset the child link of any parent element
1657 : */
1658 281572 : if (so_elem->parent())
1659 : {
1660 : unsigned int c =
1661 231063 : so_elem->parent()->which_child_am_i(so_elem);
1662 255167 : lo_elem->parent()->replace_child(lo_elem.get(), c);
1663 : }
1664 :
1665 : /*
1666 : * Copy as much data to the new element as makes sense
1667 : */
1668 281572 : lo_elem->set_p_level(so_elem->p_level());
1669 255724 : lo_elem->set_refinement_flag(so_elem->refinement_flag());
1670 51696 : lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
1671 : #endif
1672 :
1673 25848 : libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
1674 :
1675 : /*
1676 : * By definition the vertices of the linear and
1677 : * second order element are identically numbered.
1678 : * transfer these.
1679 : */
1680 1229918 : for (unsigned int v=0, snv=so_elem->n_vertices(); v < snv; v++)
1681 : {
1682 1072478 : lo_elem->set_node(v, so_elem->node_ptr(v));
1683 294852 : node_touched_by_me[lo_elem->node_id(v)] = true;
1684 : }
1685 :
1686 : /*
1687 : * find_neighbors relies on remote_elem neighbor links being
1688 : * properly maintained.
1689 : */
1690 1220978 : for (unsigned short s=0; s != n_sides; s++)
1691 : {
1692 1063058 : if (so_elem->neighbor_ptr(s) == remote_elem)
1693 0 : lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1694 : }
1695 :
1696 : /**
1697 : * If the second order element had any boundary conditions they
1698 : * should be transferred to the first-order element. The old
1699 : * boundary conditions will be removed from the BoundaryInfo
1700 : * data structure by insert_elem.
1701 : */
1702 25848 : this->get_boundary_info().copy_boundary_ids
1703 255724 : (this->get_boundary_info(), so_elem, lo_elem.get());
1704 :
1705 : /*
1706 : * The new first-order element is ready.
1707 : * Inserting it into the mesh will replace and delete
1708 : * the second-order element.
1709 : */
1710 255724 : lo_elem->set_id(so_elem->id());
1711 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1712 51696 : lo_elem->set_unique_id(so_elem->unique_id());
1713 : #endif
1714 :
1715 255724 : const unsigned int nei = so_elem->n_extra_integers();
1716 255724 : lo_elem->add_extra_integers(nei);
1717 258551 : for (unsigned int i=0; i != nei; ++i)
1718 2827 : lo_elem->set_extra_integer(i, so_elem->get_extra_integer(i));
1719 :
1720 255724 : lo_elem->inherit_data_from(*so_elem);
1721 :
1722 307420 : this->insert_elem(std::move(lo_elem));
1723 214503 : }
1724 :
1725 : // Deleting nodes does not invalidate iterators, so this is safe.
1726 1694776 : for (const auto & node : this->node_ptr_range())
1727 1011185 : if (!node_touched_by_me[node->id()])
1728 628888 : this->delete_node(node);
1729 :
1730 : // If crazy people applied boundary info to non-vertices and then
1731 : // deleted those non-vertices, we should make sure their boundary id
1732 : // caches are correct.
1733 11259 : this->get_boundary_info().regenerate_id_sets();
1734 :
1735 : // On hanging nodes that used to also be second order nodes, we
1736 : // might now have an invalid nodal processor_id()
1737 11259 : Partitioner::set_node_processor_ids(*this);
1738 :
1739 : // delete or renumber nodes if desired
1740 11259 : this->prepare_for_use();
1741 11259 : }
1742 :
1743 :
1744 :
1745 : void
1746 21631 : UnstructuredMesh::all_second_order_range (const SimpleRange<element_iterator> & range,
1747 : const bool full_ordered)
1748 : {
1749 628 : LOG_SCOPE("all_second_order_range()", "Mesh");
1750 :
1751 : /*
1752 : * The maximum number of new second order nodes we might be adding,
1753 : * for use when picking unique unique_id values later. This variable
1754 : * is not used unless unique ids are enabled.
1755 : */
1756 : unsigned int max_new_nodes_per_elem;
1757 :
1758 : /*
1759 : * For speed-up of the \p add_point() method, we
1760 : * can reserve memory. Guess the number of additional
1761 : * nodes based on the element spatial dimensions and the
1762 : * total number of nodes in the mesh as an upper bound.
1763 : */
1764 21631 : switch (this->mesh_dimension())
1765 : {
1766 923 : case 1:
1767 : /*
1768 : * in 1D, there can only be order-increase from Edge2
1769 : * to Edge3. Something like 1/2 of n_nodes() have
1770 : * to be added
1771 : */
1772 26 : max_new_nodes_per_elem = 3 - 2;
1773 1846 : this->reserve_nodes(static_cast<unsigned int>
1774 923 : (1.5*static_cast<double>(this->n_nodes())));
1775 897 : break;
1776 :
1777 3657 : case 2:
1778 : /*
1779 : * in 2D, either refine from Tri3 to Tri6 (double the nodes)
1780 : * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
1781 : */
1782 118 : max_new_nodes_per_elem = 9 - 4;
1783 7314 : this->reserve_nodes(static_cast<unsigned int>
1784 3657 : (2*static_cast<double>(this->n_nodes())));
1785 3539 : break;
1786 :
1787 :
1788 17051 : case 3:
1789 : /*
1790 : * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
1791 : * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
1792 : * quite some nodes, and since we do not want to overburden the memory by
1793 : * a too conservative guess, use the lower bound
1794 : */
1795 484 : max_new_nodes_per_elem = 27 - 8;
1796 34102 : this->reserve_nodes(static_cast<unsigned int>
1797 17051 : (2.5*static_cast<double>(this->n_nodes())));
1798 16567 : break;
1799 :
1800 0 : default:
1801 : // Hm?
1802 0 : libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1803 : }
1804 :
1805 : // All the real work is done in the helper function
1806 21631 : all_increased_order_range(*this, range, max_new_nodes_per_elem,
1807 89761 : [full_ordered](ElemType t) {
1808 2137930 : return Elem::second_order_equivalent_type(t, full_ordered);
1809 : });
1810 21631 : }
1811 :
1812 :
1813 :
1814 20367 : void UnstructuredMesh::all_complete_order_range(const SimpleRange<element_iterator> & range)
1815 : {
1816 570 : LOG_SCOPE("all_complete_order()", "Mesh");
1817 :
1818 : /*
1819 : * The maximum number of new higher-order nodes we might be adding,
1820 : * for use when picking unique unique_id values later. This variable
1821 : * is not used unless unique ids are enabled.
1822 : */
1823 : unsigned int max_new_nodes_per_elem;
1824 :
1825 : /*
1826 : * for speed-up of the \p add_point() method, we
1827 : * can reserve memory. Guess the number of additional
1828 : * nodes based on the element spatial dimensions and the
1829 : * total number of nodes in the mesh as an upper bound.
1830 : */
1831 20367 : switch (this->mesh_dimension())
1832 : {
1833 0 : case 1:
1834 : /*
1835 : * in 1D, there can only be order-increase from Edge2
1836 : * to Edge3. Something like 1/2 of n_nodes() have
1837 : * to be added
1838 : */
1839 0 : max_new_nodes_per_elem = 3 - 2;
1840 0 : this->reserve_nodes(static_cast<unsigned int>
1841 0 : (1.5*static_cast<double>(this->n_nodes())));
1842 0 : break;
1843 :
1844 1704 : case 2:
1845 : /*
1846 : * in 2D, we typically refine from Tri6 to Tri7 (1.1667 times
1847 : * the nodes) but might refine from Quad4 to Quad9
1848 : * (2.25 times the nodes)
1849 : */
1850 48 : max_new_nodes_per_elem = 9 - 4;
1851 3408 : this->reserve_nodes(static_cast<unsigned int>
1852 1704 : (2*static_cast<double>(this->n_nodes())));
1853 1656 : break;
1854 :
1855 :
1856 18663 : case 3:
1857 : /*
1858 : * in 3D, we typically refine from Tet10 to Tet14 (factor = 1.4)
1859 : * but may go Hex8 to Hex27 (something > 3). Since in 3D there
1860 : * _are_ already quite some nodes, and since we do not want to
1861 : * overburden the memory by a too conservative guess, use a
1862 : * moderate bound
1863 : */
1864 522 : max_new_nodes_per_elem = 27 - 8;
1865 37326 : this->reserve_nodes(static_cast<unsigned int>
1866 18663 : (2.5*static_cast<double>(this->n_nodes())));
1867 18141 : break;
1868 :
1869 0 : default:
1870 : // Hm?
1871 0 : libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1872 : }
1873 :
1874 : // All the real work is done in the helper function
1875 20367 : all_increased_order_range(*this, range, max_new_nodes_per_elem,
1876 143109 : [](ElemType t) {
1877 5253857 : return Elem::complete_order_equivalent_type(t);
1878 : });
1879 20367 : }
1880 :
1881 :
1882 : std::size_t
1883 1633 : UnstructuredMesh::stitch_meshes (const MeshBase & other_mesh,
1884 : boundary_id_type this_mesh_boundary_id,
1885 : boundary_id_type other_mesh_boundary_id,
1886 : Real tol,
1887 : bool clear_stitched_boundary_ids,
1888 : bool verbose,
1889 : bool use_binary_search,
1890 : bool enforce_all_nodes_match_on_boundaries,
1891 : bool merge_boundary_nodes_all_or_nothing,
1892 : bool remap_subdomain_ids,
1893 : bool prepare_after_stitching)
1894 : {
1895 92 : LOG_SCOPE("stitch_meshes()", "UnstructuredMesh");
1896 1633 : return stitching_helper(&other_mesh,
1897 : this_mesh_boundary_id,
1898 : other_mesh_boundary_id,
1899 : tol,
1900 : clear_stitched_boundary_ids,
1901 : verbose,
1902 : use_binary_search,
1903 : enforce_all_nodes_match_on_boundaries,
1904 : true,
1905 : merge_boundary_nodes_all_or_nothing,
1906 : remap_subdomain_ids,
1907 1533 : prepare_after_stitching);
1908 : }
1909 :
1910 :
1911 : std::size_t
1912 213 : UnstructuredMesh::stitch_surfaces (boundary_id_type boundary_id_1,
1913 : boundary_id_type boundary_id_2,
1914 : Real tol,
1915 : bool clear_stitched_boundary_ids,
1916 : bool verbose,
1917 : bool use_binary_search,
1918 : bool enforce_all_nodes_match_on_boundaries,
1919 : bool merge_boundary_nodes_all_or_nothing,
1920 : bool prepare_after_stitching)
1921 :
1922 : {
1923 213 : return stitching_helper(nullptr,
1924 : boundary_id_1,
1925 : boundary_id_2,
1926 : tol,
1927 : clear_stitched_boundary_ids,
1928 : verbose,
1929 : use_binary_search,
1930 : enforce_all_nodes_match_on_boundaries,
1931 : /* skip_find_neighbors = */ true,
1932 : merge_boundary_nodes_all_or_nothing,
1933 : /* remap_subdomain_ids = */ false,
1934 213 : prepare_after_stitching);
1935 : }
1936 :
1937 :
1938 : std::size_t
1939 1846 : UnstructuredMesh::stitching_helper (const MeshBase * other_mesh,
1940 : boundary_id_type this_mesh_boundary_id,
1941 : boundary_id_type other_mesh_boundary_id,
1942 : Real tol,
1943 : bool clear_stitched_boundary_ids,
1944 : bool verbose,
1945 : bool use_binary_search,
1946 : bool enforce_all_nodes_match_on_boundaries,
1947 : bool skip_find_neighbors,
1948 : bool merge_boundary_nodes_all_or_nothing,
1949 : bool remap_subdomain_ids,
1950 : bool prepare_after_stitching)
1951 : {
1952 : #ifdef DEBUG
1953 : // We rely on neighbor links here
1954 52 : MeshTools::libmesh_assert_valid_neighbors(*this);
1955 : #endif
1956 :
1957 52 : bool is_valid_disjoint_pair_to_stitch = false;
1958 :
1959 : #ifdef LIBMESH_ENABLE_PERIODIC
1960 1846 : auto * this_db = this->get_disjoint_neighbor_boundary_pairs();
1961 1846 : auto * other_db = (other_mesh ? other_mesh->get_disjoint_neighbor_boundary_pairs() : nullptr);
1962 : const bool have_disc_bdys =
1963 1846 : (this_db && !this_db->empty()) || (other_db && !other_db->empty());
1964 :
1965 52 : if (have_disc_bdys)
1966 : {
1967 10 : const boundary_id_type a = this_mesh_boundary_id;
1968 10 : const boundary_id_type b = other_mesh_boundary_id;
1969 :
1970 40 : auto get_pb = [](const PeriodicBoundaries * db, boundary_id_type id)
1971 : {
1972 539 : return db ? db->boundary(id) : nullptr;
1973 : };
1974 :
1975 : // this mesh
1976 355 : const auto * pb_this_a = get_pb(this_db, a);
1977 355 : const auto * pb_this_b = get_pb(this_db, b);
1978 10 : const bool in_this =
1979 355 : (pb_this_a && pb_this_a->pairedboundary == b) ||
1980 0 : (pb_this_b && pb_this_b->pairedboundary == a);
1981 :
1982 : // other mesh
1983 345 : const auto * pb_other_b = get_pb(other_db, b);
1984 10 : const auto * pb_other_a = get_pb(other_db, a);
1985 10 : const bool in_other =
1986 355 : (pb_other_b && pb_other_b->pairedboundary == a) ||
1987 0 : (pb_other_a && pb_other_a->pairedboundary == b);
1988 :
1989 : // Conflict conditions:
1990 : // Case 1: On "this" mesh, a or b exist but are not paired,
1991 : // while the other mesh pairs them.
1992 355 : if (!in_this && (pb_this_a || pb_this_b) && in_other)
1993 146 : libmesh_error_msg("Disjoint neighbor boundary pairing mismatch: on 'this' mesh, "
1994 : "boundary (" << a << " or " << b
1995 : << ") exists but is not paired; on 'other' mesh the pair is present.");
1996 :
1997 : // Case 2: On "other" mesh, a or b exist but are not paired,
1998 : // while this mesh pairs them.
1999 284 : if (!in_other && (pb_other_a || pb_other_b) && in_this)
2000 0 : libmesh_error_msg("Disjoint neighbor boundary pairing mismatch: on 'other' mesh, "
2001 : "boundary (" << a << " or " << b
2002 : << ") exists but is not paired; on 'this' mesh the pair is present.");
2003 :
2004 : // Legal conditions: either side has a correct pairing
2005 284 : if (in_this || in_other)
2006 6 : is_valid_disjoint_pair_to_stitch = true;
2007 : }
2008 : #endif // LIBMESH_ENABLE_PERIODIC
2009 :
2010 : // We can't even afford any unset neighbor links here.
2011 1775 : if (!this->is_prepared())
2012 71 : this->find_neighbors();
2013 :
2014 : // FIXME: make distributed mesh support efficient.
2015 : // Yes, we currently suck.
2016 1875 : MeshSerializer serialize(*this);
2017 :
2018 : // *Badly*.
2019 1725 : std::unique_ptr<MeshSerializer> serialize_other;
2020 1775 : if (other_mesh)
2021 : serialize_other = std::make_unique<MeshSerializer>
2022 3036 : (*const_cast<MeshBase *>(other_mesh));
2023 :
2024 102 : std::map<dof_id_type, dof_id_type> node_to_node_map, other_to_this_node_map; // The second is the inverse map of the first
2025 100 : std::map<dof_id_type, std::vector<dof_id_type>> node_to_elems_map;
2026 :
2027 : typedef dof_id_type key_type;
2028 : typedef std::pair<const Elem *, unsigned char> val_type;
2029 : typedef std::pair<key_type, val_type> key_val_pair;
2030 : typedef std::unordered_multimap<key_type, val_type> map_type;
2031 : // Mapping between all side keys in this mesh and elements+side numbers relevant to the boundary in this mesh as well.
2032 100 : map_type side_to_elem_map;
2033 :
2034 : // If there is only one mesh (i.e. other_mesh == nullptr), then loop over this mesh twice
2035 1775 : if (!other_mesh)
2036 : {
2037 6 : other_mesh = this;
2038 : }
2039 :
2040 1775 : if ((this_mesh_boundary_id != BoundaryInfo::invalid_id) &&
2041 50 : (other_mesh_boundary_id != BoundaryInfo::invalid_id))
2042 : {
2043 100 : LOG_SCOPE("stitch_meshes node merging", "UnstructuredMesh");
2044 :
2045 : // While finding nodes on the boundary, also find the minimum edge length
2046 : // of all faces on both boundaries. This will later be used in relative
2047 : // distance checks when stitching nodes.
2048 1775 : Real h_min = std::numeric_limits<Real>::max();
2049 50 : bool h_min_updated = false;
2050 :
2051 : // Loop below fills in these sets for the two meshes.
2052 100 : std::set<dof_id_type> this_boundary_node_ids, other_boundary_node_ids;
2053 :
2054 : // Pull objects out of the loop to reduce heap operations
2055 1775 : std::unique_ptr<const Elem> side;
2056 :
2057 : {
2058 : // Make temporary fixed-size arrays for loop
2059 1775 : boundary_id_type id_array[2] = {this_mesh_boundary_id, other_mesh_boundary_id};
2060 1775 : std::set<dof_id_type> * set_array[2] = {&this_boundary_node_ids, &other_boundary_node_ids};
2061 1775 : const MeshBase * mesh_array[2] = {this, other_mesh};
2062 :
2063 5325 : for (unsigned i=0; i<2; ++i)
2064 : {
2065 : // First we deal with node boundary IDs. We only enter
2066 : // this loop if we have at least one nodeset. Note that we
2067 : // do not attempt to make an h_min determination here.
2068 : // The h_min determination is done while looping over the
2069 : // Elems and checking their sides and edges for boundary
2070 : // information, below.
2071 3650 : if (mesh_array[i]->get_boundary_info().n_nodeset_conds() > 0)
2072 : {
2073 : // build_node_list() returns a vector of (node-id, bc-id) tuples
2074 324270 : for (const auto & t : mesh_array[i]->get_boundary_info().build_node_list())
2075 : {
2076 321204 : boundary_id_type node_bc_id = std::get<1>(t);
2077 321204 : if (node_bc_id == id_array[i])
2078 : {
2079 61415 : dof_id_type this_node_id = std::get<0>(t);
2080 61415 : set_array[i]->insert( this_node_id );
2081 : }
2082 : }
2083 : }
2084 :
2085 : // Container to catch boundary IDs passed back from BoundaryInfo.
2086 200 : std::vector<boundary_id_type> bc_ids;
2087 :
2088 : // Pointers to boundary NodeElems encountered while looping over the entire Mesh
2089 : // and checking side and edge boundary ids. The Nodes associated with NodeElems
2090 : // may be in a boundary nodeset, but not connected to any other Elems. In this
2091 : // case, we also consider the "minimum node separation distance" amongst all
2092 : // NodeElems when determining the relevant h_min value for this mesh.
2093 200 : std::vector<const Elem *> boundary_node_elems;
2094 :
2095 86902 : for (auto & el : mesh_array[i]->element_ptr_range())
2096 : {
2097 : // Now check whether elem has a face on the specified boundary
2098 267337 : for (auto side_id : el->side_index_range())
2099 : {
2100 : bool should_stitch_this_side =
2101 231450 : (el->neighbor_ptr(side_id) == nullptr) ||
2102 1420 : (is_valid_disjoint_pair_to_stitch &&
2103 1460 : mesh_array[i]->get_boundary_info().has_boundary_id(el, side_id, id_array[i]));
2104 :
2105 6340 : if (should_stitch_this_side)
2106 : {
2107 : // Get *all* boundary IDs on this side, not just the first one!
2108 99718 : mesh_array[i]->get_boundary_info().boundary_ids (el, side_id, bc_ids);
2109 :
2110 99718 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2111 : {
2112 16472 : el->build_side_ptr(side, side_id);
2113 112535 : for (auto & n : side->node_ref_range())
2114 96063 : set_array[i]->insert(n.id());
2115 :
2116 16472 : h_min = std::min(h_min, side->hmin());
2117 464 : h_min_updated = true;
2118 :
2119 : // This side is on the boundary, add its information to side_to_elem
2120 16472 : if (skip_find_neighbors && (i==0))
2121 : {
2122 8236 : key_type key = el->low_order_key(side_id);
2123 232 : val_type val;
2124 8236 : val.first = el;
2125 232 : val.second = cast_int<unsigned char>(side_id);
2126 :
2127 8236 : key_val_pair kvp;
2128 8236 : kvp.first = key;
2129 232 : kvp.second = val;
2130 232 : side_to_elem_map.insert (kvp);
2131 : }
2132 : }
2133 :
2134 : // Also, check the edges on this side. We don't have to worry about
2135 : // updating neighbor info in this case since elements don't store
2136 : // neighbor info on edges.
2137 1211294 : for (auto edge_id : el->edge_index_range())
2138 : {
2139 1111576 : if (el->is_edge_on_side(edge_id, side_id))
2140 : {
2141 : // Get *all* boundary IDs on this edge, not just the first one!
2142 368490 : mesh_array[i]->get_boundary_info().edge_boundary_ids (el, edge_id, bc_ids);
2143 :
2144 368490 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2145 : {
2146 0 : std::unique_ptr<const Elem> edge (el->build_edge_ptr(edge_id));
2147 0 : for (auto & n : edge->node_ref_range())
2148 0 : set_array[i]->insert( n.id() );
2149 :
2150 0 : h_min = std::min(h_min, edge->hmin());
2151 0 : h_min_updated = true;
2152 0 : }
2153 : }
2154 : } // end for (edge_id)
2155 : } // end if (should_stitch_this_side)
2156 : } // end for (side_id)
2157 :
2158 : // Alternatively, is this a boundary NodeElem? If so,
2159 : // add it to a list of NodeElems that will later be
2160 : // used to set h_min based on the minimum node
2161 : // separation distance between all pairs of boundary
2162 : // NodeElems.
2163 41109 : if (el->type() == NODEELEM)
2164 : {
2165 2700 : mesh_array[i]->get_boundary_info().boundary_ids(el->node_ptr(0), bc_ids);
2166 2628 : if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
2167 : {
2168 2556 : boundary_node_elems.push_back(el);
2169 :
2170 : // Debugging:
2171 : // libMesh::out << "Elem " << el->id() << " is a NodeElem on boundary " << id_array[i] << std::endl;
2172 : }
2173 : } // end if (el->type() == NODEELEM)
2174 3350 : } // end for (el)
2175 :
2176 : // Compute the minimum node separation distance amongst
2177 : // all boundary NodeElem pairs.
2178 : {
2179 200 : const auto N = boundary_node_elems.size();
2180 6106 : for (auto node_elem_i : make_range(N))
2181 47286 : for (auto node_elem_j : make_range(node_elem_i+1, N))
2182 : {
2183 : Real node_sep =
2184 47250 : (boundary_node_elems[node_elem_i]->point(0) - boundary_node_elems[node_elem_j]->point(0)).norm();
2185 :
2186 : // We only want to consider non-coincident
2187 : // boundary NodeElem pairs when determining the
2188 : // minimum node separation distance.
2189 44730 : if (node_sep > 0.)
2190 : {
2191 44730 : h_min = std::min(h_min, node_sep);
2192 1260 : h_min_updated = true;
2193 : }
2194 : } // end for (node_elem_j)
2195 : } // end minimum NodeElem separation scope
2196 : } // end for (i)
2197 : } // end scope
2198 :
2199 1775 : if (verbose)
2200 : {
2201 14 : libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
2202 14 : << "This mesh has " << this_boundary_node_ids.size()
2203 14 : << " nodes on boundary `"
2204 497 : << this->get_boundary_info().get_sideset_name(this_mesh_boundary_id)
2205 14 : << "' (" << this_mesh_boundary_id << ").\n"
2206 14 : << "Other mesh has " << other_boundary_node_ids.size()
2207 14 : << " nodes on boundary `"
2208 497 : << other_mesh->get_boundary_info().get_sideset_name(other_mesh_boundary_id)
2209 14 : << "' (" << other_mesh_boundary_id << ").\n";
2210 :
2211 497 : if (h_min_updated)
2212 : {
2213 28 : libMesh::out << "Minimum edge length on both surfaces is " << h_min << ".\n";
2214 : }
2215 : else
2216 : {
2217 0 : libMesh::out << "No minimum edge length determined on specified surfaces." << std::endl;
2218 : }
2219 : }
2220 :
2221 : // At this point, if h_min==0 it means that there were at least two coincident
2222 : // nodes on the surfaces being stitched, and we don't currently support that case.
2223 : // (It might be possible to support, but getting it exactly right would be tricky
2224 : // and probably not worth the extra complications to the "normal" case.)
2225 1775 : libmesh_error_msg_if(h_min < std::numeric_limits<Real>::epsilon(),
2226 : "Coincident nodes detected on source and/or target "
2227 : "surface, stitching meshes is not possible.");
2228 :
2229 : // We require nanoflann for the "binary search" (really kd-tree)
2230 : // option to work. If it's not available, turn that option off,
2231 : // warn the user, and fall back on the N^2 search algorithm.
2232 : if (use_binary_search)
2233 : {
2234 : #ifndef LIBMESH_HAVE_NANOFLANN
2235 : use_binary_search = false;
2236 : libmesh_warning("The use_binary_search option in the "
2237 : "UnstructuredMesh stitching algorithms requires nanoflann "
2238 : "support. Falling back on N^2 search algorithm.");
2239 : #endif
2240 : }
2241 :
2242 1775 : if (!this_boundary_node_ids.empty())
2243 : {
2244 1775 : if (use_binary_search)
2245 : {
2246 : #ifdef LIBMESH_HAVE_NANOFLANN
2247 : typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, VectorOfNodesAdaptor>,
2248 : VectorOfNodesAdaptor, 3, std::size_t> kd_tree_t;
2249 :
2250 : // Create the dataset needed to build the kd tree with nanoflann
2251 0 : std::vector<std::pair<Point, dof_id_type>> this_mesh_nodes(this_boundary_node_ids.size());
2252 :
2253 0 : for (auto [it, ctr] = std::make_tuple(this_boundary_node_ids.begin(), 0u);
2254 0 : it != this_boundary_node_ids.end(); ++it, ++ctr)
2255 : {
2256 0 : this_mesh_nodes[ctr].first = this->point(*it);
2257 0 : this_mesh_nodes[ctr].second = *it;
2258 : }
2259 :
2260 0 : VectorOfNodesAdaptor vec_nodes_adaptor(this_mesh_nodes);
2261 :
2262 0 : kd_tree_t this_kd_tree(3, vec_nodes_adaptor, 10);
2263 0 : this_kd_tree.buildIndex();
2264 :
2265 : // Storage for nearest neighbor in the loop below
2266 : std::size_t ret_index;
2267 : Real ret_dist_sqr;
2268 :
2269 : // Loop over other mesh. For each node, find its nearest neighbor in this mesh, and fill in the maps.
2270 0 : for (const auto & node_id : other_boundary_node_ids)
2271 : {
2272 0 : const auto & p = other_mesh->point(node_id);
2273 0 : const Real query_pt[] = {p(0), p(1), p(2)};
2274 0 : this_kd_tree.knnSearch(&query_pt[0], 1, &ret_index, &ret_dist_sqr);
2275 :
2276 : // TODO: here we should use the user's specified tolerance
2277 : // and the previously determined value of h_min in the
2278 : // distance comparison, not just TOLERANCE^2.
2279 0 : if (ret_dist_sqr < TOLERANCE*TOLERANCE)
2280 : {
2281 0 : node_to_node_map[this_mesh_nodes[ret_index].second] = node_id;
2282 0 : other_to_this_node_map[node_id] = this_mesh_nodes[ret_index].second;
2283 : }
2284 : }
2285 :
2286 : // If the two maps don't have the same size, it means one
2287 : // node in this mesh is the nearest neighbor of several
2288 : // nodes in other mesh. Since the stitching is ambiguous in
2289 : // this case, we throw an error.
2290 0 : libmesh_error_msg_if(node_to_node_map.size() != other_to_this_node_map.size(),
2291 : "Error: Found multiple matching nodes in stitch_meshes");
2292 : #endif
2293 : }
2294 : else // !use_binary_search
2295 : {
2296 : // In the unlikely event that two meshes composed entirely of
2297 : // NodeElems are being stitched together, we will not have
2298 : // selected a valid h_min value yet, and the distance
2299 : // comparison below will be true for essentially any two
2300 : // nodes. In this case we simply fall back on an absolute
2301 : // distance check.
2302 1775 : if (!h_min_updated)
2303 : {
2304 : libmesh_warning("No valid h_min value was found, falling back on "
2305 : "absolute distance check in the N^2 search algorithm.");
2306 0 : h_min = 1.;
2307 : }
2308 :
2309 : // Otherwise, use a simple N^2 search to find the closest matching points. This can be helpful
2310 : // in the case that we have tolerance issues which cause mismatch between the two surfaces
2311 : // that are being stitched.
2312 32944 : for (const auto & this_node_id : this_boundary_node_ids)
2313 : {
2314 31169 : Node & this_node = this->node_ref(this_node_id);
2315 :
2316 878 : bool found_matching_nodes = false;
2317 :
2318 928538 : for (const auto & other_node_id : other_boundary_node_ids)
2319 : {
2320 897369 : const Node & other_node = other_mesh->node_ref(other_node_id);
2321 :
2322 872091 : Real node_distance = (this_node - other_node).norm();
2323 :
2324 897369 : if (node_distance < tol*h_min)
2325 : {
2326 : // Make sure we didn't already find a matching node!
2327 31169 : libmesh_error_msg_if(found_matching_nodes,
2328 : "Error: Found multiple matching nodes in stitch_meshes");
2329 :
2330 31169 : node_to_node_map[this_node_id] = other_node_id;
2331 31169 : other_to_this_node_map[other_node_id] = this_node_id;
2332 :
2333 878 : found_matching_nodes = true;
2334 : }
2335 : }
2336 : }
2337 : }
2338 : }
2339 :
2340 : // Build up the node_to_elems_map, using only one loop over other_mesh
2341 42830 : for (auto & el : other_mesh->element_ptr_range())
2342 : {
2343 : // For each node on the element, find the corresponding node
2344 : // on "this" Mesh, 'this_node_id', if it exists, and push
2345 : // the current element ID back onto node_to_elems_map[this_node_id].
2346 : // For that we will use the reverse mapping we created at
2347 : // the same time as the forward mapping.
2348 294865 : for (auto & n : el->node_ref_range())
2349 281780 : if (const auto it = other_to_this_node_map.find(/*other_node_id=*/n.id());
2350 7720 : it != other_to_this_node_map.end())
2351 50481 : node_to_elems_map[/*this_node_id=*/it->second].push_back( el->id() );
2352 1675 : }
2353 :
2354 1775 : if (verbose)
2355 : {
2356 14 : libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
2357 28 : << "Found " << node_to_node_map.size()
2358 14 : << " matching nodes.\n"
2359 14 : << std::endl;
2360 : }
2361 :
2362 1775 : if (enforce_all_nodes_match_on_boundaries)
2363 : {
2364 2 : std::size_t n_matching_nodes = node_to_node_map.size();
2365 2 : std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2366 2 : std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2367 71 : libmesh_error_msg_if((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes),
2368 : "Error: We expected the number of nodes to match.");
2369 : }
2370 :
2371 1775 : if (merge_boundary_nodes_all_or_nothing)
2372 : {
2373 2 : std::size_t n_matching_nodes = node_to_node_map.size();
2374 2 : std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2375 2 : std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2376 71 : if ((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes))
2377 : {
2378 0 : if (verbose)
2379 : {
2380 : libMesh::out << "Skipping node merging in "
2381 : "UnstructuredMesh::stitch_meshes because not "
2382 0 : "all boundary nodes were matched."
2383 0 : << std::endl;
2384 : }
2385 0 : node_to_node_map.clear();
2386 0 : other_to_this_node_map.clear();
2387 0 : node_to_elems_map.clear();
2388 : }
2389 100 : }
2390 3350 : }
2391 : else
2392 : {
2393 0 : if (verbose)
2394 : {
2395 0 : libMesh::out << "Skip node merging in UnstructuredMesh::stitch_meshes:" << std::endl;
2396 : }
2397 : }
2398 :
2399 1775 : dof_id_type node_delta = this->max_node_id();
2400 1775 : dof_id_type elem_delta = this->max_elem_id();
2401 :
2402 : unique_id_type unique_delta =
2403 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2404 1775 : this->parallel_max_unique_id();
2405 : #else
2406 : 0;
2407 : #endif
2408 :
2409 : // If other_mesh != nullptr, then we have to do a bunch of work
2410 : // in order to copy it to this mesh
2411 1775 : if (this!=other_mesh)
2412 : {
2413 88 : LOG_SCOPE("stitch_meshes copying", "UnstructuredMesh");
2414 :
2415 : #ifdef LIBMESH_ENABLE_PERIODIC
2416 : // Copy disjoint neighbor boundary pairs (PeriodicBoundary objects)
2417 : // from `other_mesh` to `this` mesh
2418 1562 : if (other_db && !other_db->empty())
2419 : {
2420 213 : for (const auto & [bdy_id, pb_ptr] : *other_db)
2421 : {
2422 4 : const auto & pb = *pb_ptr;
2423 142 : const boundary_id_type a = pb.myboundary;
2424 142 : const boundary_id_type b = pb.pairedboundary;
2425 :
2426 142 : if (this_db)
2427 : {
2428 : // Skip if identical pair already exists
2429 142 : if (const auto * existing_pb = this_db->boundary(a))
2430 71 : if ((existing_pb->myboundary == a && existing_pb->pairedboundary == b) ||
2431 0 : (existing_pb->myboundary == b && existing_pb->pairedboundary == a))
2432 69 : continue;
2433 :
2434 : // If both boundary ids exist on this mesh but aren't paired here, refuse to create a new pair
2435 2 : const auto & bdy_ids = this->get_boundary_info().get_boundary_ids();
2436 4 : const bool a_exists = bdy_ids.count(a);
2437 4 : const bool b_exists = bdy_ids.count(b);
2438 : // If a and b already exist on `this`, we should be screaming and dying
2439 : // unless they already have PeriodicBoundary objects connecting them too
2440 71 : if (a_exists && b_exists && !this_db->boundary(a))
2441 0 : libmesh_error_msg("Conflict: boundaries " << a << " and " << b
2442 : << " already exist on this mesh but are not paired.");
2443 : }
2444 :
2445 140 : this->add_disjoint_neighbor_boundary_pairs(a, b, pb.get_corresponding_pos(Point(0.0,0.0,0.0)));
2446 : }
2447 : }
2448 : #endif // LIBMESH_ENABLE_PERIODIC
2449 :
2450 :
2451 : // Increment the node_to_node_map and node_to_elems_map
2452 : // to account for id offsets
2453 32305 : for (auto & pr : node_to_node_map)
2454 30743 : pr.second += node_delta;
2455 :
2456 32305 : for (auto & pr : node_to_elems_map)
2457 80798 : for (auto & entry : pr.second)
2458 50055 : entry += elem_delta;
2459 :
2460 : // We run into problems when the libMesh subdomain standard (the
2461 : // id defines the subdomain; the name was an afterthought) and
2462 : // the MOOSE standard (the name defines the subdomain; the id
2463 : // might be autogenerated) clash.
2464 : //
2465 : // Subdomain ids with the same name in both meshes are surely
2466 : // meant to represent the same subdomain. We can just merge
2467 : // them.
2468 : //
2469 : // Subdomain ids which don't have a name in either mesh are
2470 : // almost surely meant to represent the same subdomain. We'll
2471 : // just merge them.
2472 : //
2473 : // Subdomain ids with different names in different meshes, or
2474 : // names with different ids in different meshes, are trickier.
2475 : // For backwards compatibility we default to the old "just copy
2476 : // all the subdomain ids over" behavior, but if requested we'll
2477 : // remap any ids that appear to be clear conflicts, and we'll
2478 : // scream and die if we see any ids that are ambiguous due to
2479 : // being named in one mesh but not the other.
2480 88 : std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
2481 1562 : if (remap_subdomain_ids)
2482 : {
2483 4 : const auto & this_map = this->get_subdomain_name_map();
2484 4 : const auto & other_map = other_mesh->get_subdomain_name_map();
2485 8 : std::unordered_map<std::string, subdomain_id_type> other_map_reversed;
2486 284 : for (auto & [sid, sname] : other_map)
2487 4 : other_map_reversed.emplace(sname, sid);
2488 :
2489 8 : std::unordered_map<std::string, subdomain_id_type> this_map_reversed;
2490 213 : for (auto & [sid, sname] : this_map)
2491 2 : this_map_reversed.emplace(sname, sid);
2492 :
2493 : // We don't require either mesh to be prepared, but that
2494 : // means we need to check for subdomains manually.
2495 284 : auto get_subdomains = [](const MeshBase & mesh) {
2496 8 : std::set<subdomain_id_type> all_subdomains;
2497 4976 : for (auto & el : mesh.element_ptr_range())
2498 2540 : all_subdomains.insert(el->subdomain_id());
2499 284 : return all_subdomains;
2500 : };
2501 :
2502 146 : const auto this_subdomains = get_subdomains(*this);
2503 146 : const auto other_subdomains = get_subdomains(*other_mesh);
2504 :
2505 213 : for (auto & [sid, sname] : this_map)
2506 : {
2507 : // The same name with the same id means we're fine. The
2508 : // same name with another id means we remap their id to
2509 : // ours
2510 2 : if (const auto other_reverse_it = other_map_reversed.find(sname);
2511 71 : other_reverse_it != other_map_reversed.end() && other_reverse_it->second != sid)
2512 71 : id_remapping[other_reverse_it->second] = sid;
2513 :
2514 : // The same id with a different name, we'll get to
2515 : // later. The same id without any name means we don't
2516 : // know what the user wants.
2517 2 : if (other_subdomains.count(sid) && !other_map.count(sid))
2518 0 : libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2519 : << sid << " but not subdomain name " << sname);
2520 : }
2521 :
2522 142 : subdomain_id_type next_free_id = 0;
2523 : // We might try to stitch empty meshes ...
2524 142 : if (!this_subdomains.empty())
2525 142 : next_free_id = *this_subdomains.rbegin() + 1;
2526 142 : if (!other_subdomains.empty())
2527 142 : next_free_id =
2528 142 : std::max(next_free_id,
2529 : cast_int<subdomain_id_type>
2530 215 : (*other_subdomains.rbegin() + 1));
2531 :
2532 213 : for (auto & [sid, sname] : other_map)
2533 : {
2534 : // At this point we've figured out any remapping
2535 : // necessary for an sname that we share. And we don't
2536 : // need to remap any sid we don't share.
2537 8 : if (!this_map_reversed.count(sname))
2538 : {
2539 : // But if we don't have this sname and we do have this
2540 : // sid then we can't just merge into that.
2541 2 : if (this_subdomains.count(sid))
2542 : {
2543 : // If we have this sid with no name, we don't
2544 : // know what the user wants.
2545 2 : if (!this_map.count(sid))
2546 211 : libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2547 : << sid << " but under subdomain name " << sname);
2548 :
2549 : // We have this sid under a different name, so
2550 : // we just need to give the other elements a new
2551 : // id.
2552 :
2553 : // Users might have done crazy things with id
2554 : // choice so let's make sure they didn't get too
2555 : // crazy.
2556 0 : libmesh_error_msg_if ((!this_subdomains.empty() &&
2557 : next_free_id < *this_subdomains.rbegin()) ||
2558 : (!other_subdomains.empty() &&
2559 : next_free_id < *other_subdomains.rbegin()),
2560 : "Subdomain id overflow");
2561 :
2562 0 : id_remapping[sid] = next_free_id++;
2563 0 : this->subdomain_name(next_free_id) = sname;
2564 : }
2565 : // If we don't have this subdomain id, well, we're
2566 : // about to, so we should have its name too.
2567 : else
2568 0 : this->subdomain_name(sid) = sname;
2569 : }
2570 : }
2571 : }
2572 :
2573 : // Copy mesh data. If we skip the call to find_neighbors(), the lists
2574 : // of neighbors will be copied verbatim from the other mesh
2575 1491 : this->copy_nodes_and_elements(*other_mesh, skip_find_neighbors,
2576 : elem_delta, node_delta,
2577 84 : unique_delta, &id_remapping);
2578 :
2579 : // Copy BoundaryInfo from other_mesh too. We do this via the
2580 : // list APIs rather than element-by-element for speed.
2581 42 : BoundaryInfo & boundary = this->get_boundary_info();
2582 42 : const BoundaryInfo & other_boundary = other_mesh->get_boundary_info();
2583 :
2584 158869 : for (const auto & t : other_boundary.build_node_list())
2585 157336 : boundary.add_node(std::get<0>(t) + node_delta,
2586 157336 : std::get<1>(t));
2587 :
2588 42855 : for (const auto & t : other_boundary.build_side_list())
2589 42486 : boundary.add_side(std::get<0>(t) + elem_delta,
2590 41322 : std::get<1>(t),
2591 41322 : std::get<2>(t));
2592 :
2593 1533 : for (const auto & t : other_boundary.build_edge_list())
2594 0 : boundary.add_edge(std::get<0>(t) + elem_delta,
2595 0 : std::get<1>(t),
2596 0 : std::get<2>(t));
2597 :
2598 1533 : for (const auto & t : other_boundary.build_shellface_list())
2599 0 : boundary.add_shellface(std::get<0>(t) + elem_delta,
2600 0 : std::get<1>(t),
2601 0 : std::get<2>(t));
2602 :
2603 42 : const auto & other_ns_id_to_name = other_boundary.get_nodeset_name_map();
2604 42 : auto & ns_id_to_name = boundary.set_nodeset_name_map();
2605 1491 : ns_id_to_name.insert(other_ns_id_to_name.begin(), other_ns_id_to_name.end());
2606 :
2607 42 : const auto & other_ss_id_to_name = other_boundary.get_sideset_name_map();
2608 42 : auto & ss_id_to_name = boundary.set_sideset_name_map();
2609 1491 : ss_id_to_name.insert(other_ss_id_to_name.begin(), other_ss_id_to_name.end());
2610 :
2611 42 : const auto & other_es_id_to_name = other_boundary.get_edgeset_name_map();
2612 42 : auto & es_id_to_name = boundary.set_edgeset_name_map();
2613 1491 : es_id_to_name.insert(other_es_id_to_name.begin(), other_es_id_to_name.end());
2614 :
2615 : // Merge other_mesh's elemset information with ours. Throw an
2616 : // error if this and other_mesh have overlapping elemset codes
2617 : // that refer to different elemset ids.
2618 1533 : std::vector<dof_id_type> this_elemset_codes = this->get_elemset_codes();
2619 84 : MeshBase::elemset_type this_id_set_to_fill, other_id_set_to_fill;
2620 1817 : for (const auto & elemset_code : other_mesh->get_elemset_codes())
2621 : {
2622 : // Get the elemset ids for this elemset_code on other_mesh
2623 284 : other_mesh->get_elemsets(elemset_code, other_id_set_to_fill);
2624 :
2625 : // Check that this elemset code does not already exist
2626 : // in this mesh, or if it does, that it has the same elemset
2627 : // ids associated with it.
2628 : //
2629 : // Note: get_elemset_codes() is guaranteed to return a
2630 : // sorted vector, so we can binary search in it.
2631 268 : auto it = Utility::binary_find(this_elemset_codes.begin(),
2632 : this_elemset_codes.end(),
2633 16 : elemset_code);
2634 :
2635 284 : if (it != this_elemset_codes.end())
2636 : {
2637 : // This mesh has the same elemset code. Does it refer to
2638 : // the same elemset ids?
2639 0 : this->get_elemsets(elemset_code, this_id_set_to_fill);
2640 :
2641 : // Throw an error if they don't match, otherwise we
2642 : // don't need to do anything
2643 0 : libmesh_error_msg_if(other_id_set_to_fill != this_id_set_to_fill,
2644 : "Attempted to stitch together meshes with conflicting elemset codes.");
2645 : }
2646 : else
2647 : {
2648 : // Add other_mesh's elemset code to this mesh
2649 560 : this->add_elemset_code(elemset_code, other_id_set_to_fill);
2650 : }
2651 : }
2652 :
2653 : } // end if (other_mesh)
2654 :
2655 : // Finally, we need to "merge" the overlapping nodes
2656 : // We do this by iterating over node_to_elems_map and updating
2657 : // the elements so that they "point" to the nodes that came
2658 : // from this mesh, rather than from other_mesh.
2659 : // Then we iterate over node_to_node_map and delete the
2660 : // duplicate nodes that came from other_mesh.
2661 :
2662 : {
2663 96 : LOG_SCOPE("stitch_meshes node updates", "UnstructuredMesh");
2664 :
2665 : // Container to catch boundary IDs passed back from BoundaryInfo.
2666 96 : std::vector<boundary_id_type> bc_ids;
2667 :
2668 32234 : for (const auto & [target_node_id, elem_vec] : node_to_elems_map)
2669 : {
2670 30530 : dof_id_type other_node_id = node_to_node_map[target_node_id];
2671 30530 : Node & target_node = this->node_ref(target_node_id);
2672 :
2673 1720 : std::size_t n_elems = elem_vec.size();
2674 79875 : for (std::size_t i=0; i<n_elems; i++)
2675 : {
2676 49345 : dof_id_type elem_id = elem_vec[i];
2677 49345 : Elem * el = this->elem_ptr(elem_id);
2678 :
2679 : // find the local node index that we want to update
2680 47955 : unsigned int local_node_index = el->local_node(other_node_id);
2681 1390 : libmesh_assert_not_equal_to(local_node_index, libMesh::invalid_uint);
2682 :
2683 : // We also need to copy over the nodeset info here,
2684 : // because the node will get deleted below
2685 50735 : this->get_boundary_info().boundary_ids(el->node_ptr(local_node_index), bc_ids);
2686 49345 : el->set_node(local_node_index, &target_node);
2687 49345 : this->get_boundary_info().add_node(&target_node, bc_ids);
2688 : }
2689 : }
2690 : }
2691 :
2692 : {
2693 96 : LOG_SCOPE("stitch_meshes node deletion", "UnstructuredMesh");
2694 32234 : for (const auto & [other_node_id, this_node_id] : node_to_node_map)
2695 : {
2696 : // In the case that this==other_mesh, the two nodes might be the same (e.g. if
2697 : // we're stitching a "sliver"), hence we need to skip node deletion in that case.
2698 30530 : if ((this == other_mesh) && (this_node_id == other_node_id))
2699 0 : continue;
2700 :
2701 30530 : this->delete_node( this->node_ptr(this_node_id) );
2702 : }
2703 : }
2704 :
2705 : // If find_neighbors() wasn't called in prepare_for_use(), we need to
2706 : // manually loop once more over all elements adjacent to the stitched boundary
2707 : // and fix their lists of neighbors.
2708 : // This is done according to the following steps:
2709 : // 1. Loop over all copied elements adjacent to the boundary using node_to_elems_map (trying to avoid duplicates)
2710 : // 2. Look at all their sides with a nullptr neighbor and update them using side_to_elem_map if necessary
2711 : // 3. Update the corresponding side in side_to_elem_map as well
2712 1704 : if (skip_find_neighbors)
2713 : {
2714 96 : LOG_SCOPE("stitch_meshes neighbor fixes", "UnstructuredMesh");
2715 :
2716 : // Pull objects out of the loop to reduce heap operations
2717 1704 : std::unique_ptr<const Elem> my_side, their_side;
2718 :
2719 96 : std::set<dof_id_type> fixed_elems;
2720 32234 : for (const auto & pr : node_to_elems_map)
2721 : {
2722 1720 : std::size_t n_elems = pr.second.size();
2723 79875 : for (std::size_t i=0; i<n_elems; i++)
2724 : {
2725 50735 : dof_id_type elem_id = pr.second[i];
2726 1390 : if (!fixed_elems.count(elem_id))
2727 : {
2728 10508 : Elem * el = this->elem_ptr(elem_id);
2729 10212 : fixed_elems.insert(elem_id);
2730 57226 : for (auto s : el->side_index_range())
2731 : {
2732 46718 : bool has_real_neighbor = (el->neighbor_ptr(s) != nullptr);
2733 47570 : bool has_disdjoint_neighbor = is_valid_disjoint_pair_to_stitch &&
2734 852 : (this->get_boundary_info().has_boundary_id(el, s, this_mesh_boundary_id)
2735 852 : || this->get_boundary_info().has_boundary_id(el, s, other_mesh_boundary_id));
2736 :
2737 46718 : if (!has_real_neighbor || has_disdjoint_neighbor)
2738 : {
2739 21016 : key_type key = el->low_order_key(s);
2740 592 : auto bounds = side_to_elem_map.equal_range(key);
2741 :
2742 21016 : if (bounds.first != bounds.second)
2743 : {
2744 : // Get the side for this element
2745 7952 : el->side_ptr(my_side, s);
2746 :
2747 : // Look at all the entries with an equivalent key
2748 7952 : while (bounds.first != bounds.second)
2749 : {
2750 : // Get the potential element
2751 7952 : Elem * neighbor = const_cast<Elem *>(bounds.first->second.first);
2752 :
2753 : // Get the side for the neighboring element
2754 7952 : const unsigned int ns = bounds.first->second.second;
2755 7952 : neighbor->side_ptr(their_side, ns);
2756 : //libmesh_assert(my_side.get());
2757 : //libmesh_assert(their_side.get());
2758 :
2759 : // If found a match with my side
2760 : //
2761 : // We need special tests here for 1D:
2762 : // since parents and children have an equal
2763 : // side (i.e. a node), we need to check
2764 : // ns != ms, and we also check level() to
2765 : // avoid setting our neighbor pointer to
2766 : // any of our neighbor's descendants
2767 15680 : if ((*my_side == *their_side) &&
2768 15904 : (el->level() == neighbor->level()) &&
2769 7952 : ((el->dim() != 1) || (ns != s)))
2770 : {
2771 : // So share a side. Is this a mixed pair
2772 : // of subactive and active/ancestor
2773 : // elements?
2774 : // If not, then we're neighbors.
2775 : // If so, then the subactive's neighbor is
2776 :
2777 8176 : if (el->subactive() ==
2778 7952 : neighbor->subactive())
2779 : {
2780 : // an element is only subactive if it has
2781 : // been coarsened but not deleted
2782 448 : el->set_neighbor (s,neighbor);
2783 448 : neighbor->set_neighbor(ns,el);
2784 : }
2785 0 : else if (el->subactive())
2786 : {
2787 0 : el->set_neighbor(s,neighbor);
2788 : }
2789 0 : else if (neighbor->subactive())
2790 : {
2791 0 : neighbor->set_neighbor(ns,el);
2792 : }
2793 : // It's OK to invalidate the
2794 : // bounds.first iterator here,
2795 : // as we are immediately going
2796 : // to break out of this while
2797 : // loop. bounds.first will
2798 : // therefore not be used for
2799 : // anything else.
2800 224 : side_to_elem_map.erase (bounds.first);
2801 224 : break;
2802 : }
2803 :
2804 0 : ++bounds.first;
2805 : }
2806 : }
2807 : }
2808 : }
2809 : }
2810 : }
2811 : }
2812 1608 : }
2813 :
2814 : #ifdef LIBMESH_ENABLE_PERIODIC
2815 : // Remove only the disjoint pair that was actually stitched.
2816 : // Safe because `is_valid_disjoint_pair_to_stitch` is true
2817 : // only if this exact (a,b) pair exists in the registry.
2818 : // Other disjoint pairs remain untouched.
2819 1704 : if (is_valid_disjoint_pair_to_stitch)
2820 213 : this->remove_disjoint_boundary_pair(this_mesh_boundary_id, other_mesh_boundary_id);
2821 : #endif
2822 :
2823 1704 : if (prepare_after_stitching)
2824 : {
2825 : // We set our new neighbor pointers already
2826 92 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
2827 46 : this->allow_find_neighbors(!skip_find_neighbors);
2828 :
2829 : // We haven't newly remoted any elements
2830 92 : const bool old_allow_remote_element_removal = this->allow_remote_element_removal();
2831 46 : this->allow_remote_element_removal(false);
2832 :
2833 1633 : this->prepare_for_use();
2834 :
2835 46 : this->allow_find_neighbors(old_allow_find_neighbors);
2836 46 : this->allow_remote_element_removal(old_allow_remote_element_removal);
2837 : }
2838 :
2839 : // After the stitching, we may want to clear boundary IDs from element
2840 : // faces that are now internal to the mesh
2841 1704 : if (clear_stitched_boundary_ids)
2842 : {
2843 92 : LOG_SCOPE("stitch_meshes clear bcids", "UnstructuredMesh");
2844 :
2845 1633 : this->get_boundary_info().clear_stitched_boundary_side_ids(
2846 : this_mesh_boundary_id, other_mesh_boundary_id, /*clear_nodeset_data=*/true);
2847 : }
2848 :
2849 : // Return the number of nodes which were merged.
2850 1752 : return node_to_node_map.size();
2851 1742 : }
2852 :
2853 :
2854 : } // namespace libMesh
|
