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