Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/replicated_mesh.h"
22 :
23 : #include "libmesh/boundary_info.h"
24 : #include "libmesh/elem.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/mesh_communication.h"
27 : #include "libmesh/parallel_implementation.h"
28 : #include "libmesh/partitioner.h"
29 : #include "libmesh/point.h"
30 : #include "libmesh/string_to_enum.h"
31 : #include "libmesh/utility.h"
32 :
33 : // C++ includes
34 : #include <unordered_map>
35 : #include <unordered_set>
36 :
37 : namespace libMesh
38 : {
39 :
40 : // ------------------------------------------------------------
41 : // ReplicatedMesh class member functions
42 76957 : ReplicatedMesh::ReplicatedMesh (const Parallel::Communicator & comm_in,
43 76957 : unsigned char d) :
44 : UnstructuredMesh (comm_in,d),
45 76957 : _n_nodes(0), _n_elem(0)
46 : {
47 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
48 : // In serial we just need to reset the next unique id to zero
49 : // here in the constructor.
50 76957 : _next_unique_id = 0;
51 : #endif
52 :
53 90530 : const std::string default_partitioner = "metis";
54 : const std::string my_partitioner =
55 : libMesh::command_line_value("--default-partitioner",
56 153914 : default_partitioner);
57 : _partitioner = Partitioner::build
58 140341 : (Utility::string_to_enum<PartitionerType>(my_partitioner));
59 76957 : }
60 :
61 :
62 2393 : bool ReplicatedMesh::subclass_locally_equals(const MeshBase & other_mesh_base) const
63 : {
64 1036 : const ReplicatedMesh * rep_mesh_ptr =
65 2393 : dynamic_cast<const ReplicatedMesh *>(&other_mesh_base);
66 2393 : if (!rep_mesh_ptr)
67 0 : return false;
68 1036 : const ReplicatedMesh & other_mesh = *rep_mesh_ptr;
69 :
70 5822 : if (_n_nodes != other_mesh._n_nodes ||
71 2393 : _n_elem != other_mesh._n_elem ||
72 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
73 7179 : _next_unique_id != other_mesh._next_unique_id ||
74 : #endif
75 2393 : !this->nodes_and_elements_equal(other_mesh))
76 78 : return false;
77 :
78 1032 : return true;
79 : }
80 :
81 :
82 91450 : ReplicatedMesh::~ReplicatedMesh ()
83 : {
84 81139 : this->ReplicatedMesh::clear(); // Free nodes and elements
85 91450 : }
86 :
87 :
88 : // This might be specialized later, but right now it's just here to
89 : // make sure the compiler doesn't give us a default (non-deep) copy
90 : // constructor instead.
91 4111 : ReplicatedMesh::ReplicatedMesh (const ReplicatedMesh & other_mesh) :
92 4111 : ReplicatedMesh(static_cast<const MeshBase&>(other_mesh))
93 : {
94 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
95 4111 : this->_next_unique_id = other_mesh._next_unique_id;
96 : #endif
97 4111 : }
98 :
99 :
100 4182 : ReplicatedMesh::ReplicatedMesh (const MeshBase & other_mesh) :
101 : UnstructuredMesh (other_mesh),
102 4182 : _n_nodes(0), _n_elem(0) // copy_* will increment this
103 : {
104 : // Just copy, skipping preparation
105 4182 : this->copy_nodes_and_elements(other_mesh, true, 0, 0, 0, nullptr, true);
106 :
107 1420 : this->allow_find_neighbors(other_mesh.allow_find_neighbors());
108 1420 : this->allow_detect_interior_parents(other_mesh.allow_detect_interior_parents());
109 1420 : this->allow_renumbering(other_mesh.allow_renumbering());
110 1420 : this->allow_remote_element_removal(other_mesh.allow_remote_element_removal());
111 1420 : this->skip_partitioning(other_mesh.skip_partitioning());
112 :
113 4182 : this->copy_constraint_rows(other_mesh);
114 :
115 722 : auto & this_boundary_info = this->get_boundary_info();
116 722 : const auto & other_boundary_info = other_mesh.get_boundary_info();
117 :
118 4182 : this_boundary_info = other_boundary_info;
119 :
120 722 : this->set_subdomain_name_map() = other_mesh.get_subdomain_name_map();
121 :
122 4182 : this->_preparation = other_mesh.preparation();
123 :
124 : // If other_mesh is distributed, then we've got parts of it on each
125 : // processor but we're not replicated yet; fix that.
126 4182 : if (!other_mesh.is_serial())
127 71 : MeshCommunication().allgather(*this);
128 4182 : }
129 :
130 156 : ReplicatedMesh & ReplicatedMesh::operator= (ReplicatedMesh && other_mesh)
131 : {
132 16 : LOG_SCOPE("operator=(&&)", "ReplicatedMesh");
133 :
134 : // Move assign as an UnstructuredMesh
135 8 : this->UnstructuredMesh::operator=(std::move(other_mesh));
136 :
137 : // Nodes and elements belong to ReplicatedMesh and have to be
138 : // moved before we can move arbitrary GhostingFunctor, Partitioner,
139 : // etc. subclasses.
140 156 : this->move_nodes_and_elements(std::move(other_mesh));
141 :
142 : // Handle those remaining moves.
143 156 : this->post_dofobject_moves(std::move(other_mesh));
144 :
145 164 : return *this;
146 : }
147 :
148 156 : MeshBase & ReplicatedMesh::assign(MeshBase && other_mesh)
149 : {
150 156 : *this = std::move(cast_ref<ReplicatedMesh&>(other_mesh));
151 :
152 156 : return *this;
153 : }
154 :
155 156 : void ReplicatedMesh::move_nodes_and_elements(MeshBase && other_meshbase)
156 : {
157 8 : ReplicatedMesh & other_mesh = cast_ref<ReplicatedMesh&>(other_meshbase);
158 :
159 156 : this->_nodes = std::move(other_mesh._nodes);
160 156 : this->_n_nodes = other_mesh.n_nodes();
161 :
162 156 : this->_elements = std::move(other_mesh._elements);
163 156 : this->_n_elem = other_mesh.n_elem();
164 156 : }
165 :
166 :
167 16049743 : const Point & ReplicatedMesh::point (const dof_id_type i) const
168 : {
169 16049743 : return this->node_ref(i);
170 : }
171 :
172 :
173 :
174 :
175 101625346 : const Node * ReplicatedMesh::node_ptr (const dof_id_type i) const
176 : {
177 5274200 : libmesh_assert_less (i, this->max_node_id());
178 5274200 : libmesh_assert(_nodes[i]);
179 5274200 : libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
180 :
181 106861652 : return _nodes[i];
182 : }
183 :
184 :
185 :
186 :
187 104223830 : Node * ReplicatedMesh::node_ptr (const dof_id_type i)
188 : {
189 10154960 : libmesh_assert_less (i, this->max_node_id());
190 10154960 : libmesh_assert(_nodes[i]);
191 10154960 : libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
192 :
193 113730430 : return _nodes[i];
194 : }
195 :
196 :
197 :
198 :
199 17557946 : const Node * ReplicatedMesh::query_node_ptr (const dof_id_type i) const
200 : {
201 17557946 : if (i >= this->max_node_id())
202 0 : return nullptr;
203 17307941 : libmesh_assert (_nodes[i] == nullptr ||
204 : _nodes[i]->id() == i); // This will change soon
205 :
206 17643128 : return _nodes[i];
207 : }
208 :
209 :
210 :
211 :
212 2818155 : Node * ReplicatedMesh::query_node_ptr (const dof_id_type i)
213 : {
214 2818155 : if (i >= this->max_node_id())
215 87338 : return nullptr;
216 224566 : libmesh_assert (_nodes[i] == nullptr ||
217 : _nodes[i]->id() == i); // This will change soon
218 :
219 2162231 : return _nodes[i];
220 : }
221 :
222 :
223 :
224 :
225 3203770 : const Elem * ReplicatedMesh::elem_ptr (const dof_id_type i) const
226 : {
227 517714 : libmesh_assert_less (i, this->max_elem_id());
228 517714 : libmesh_assert(_elements[i]);
229 517714 : libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
230 :
231 3625055 : return _elements[i];
232 : }
233 :
234 :
235 :
236 :
237 92881946 : Elem * ReplicatedMesh::elem_ptr (const dof_id_type i)
238 : {
239 5266263 : libmesh_assert_less (i, this->max_elem_id());
240 5266263 : libmesh_assert(_elements[i]);
241 5266263 : libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
242 :
243 97389201 : return _elements[i];
244 : }
245 :
246 :
247 :
248 :
249 15411598 : const Elem * ReplicatedMesh::query_elem_ptr (const dof_id_type i) const
250 : {
251 15411598 : if (i >= this->max_elem_id())
252 0 : return nullptr;
253 15327258 : libmesh_assert (_elements[i] == nullptr ||
254 : _elements[i]->id() == i); // This will change soon
255 :
256 15441182 : return _elements[i];
257 : }
258 :
259 :
260 :
261 :
262 1476394 : Elem * ReplicatedMesh::query_elem_ptr (const dof_id_type i)
263 : {
264 1476394 : if (i >= this->max_elem_id())
265 44895 : return nullptr;
266 190638 : libmesh_assert (_elements[i] == nullptr ||
267 : _elements[i]->id() == i); // This will change soon
268 :
269 823071 : return _elements[i];
270 : }
271 :
272 :
273 :
274 :
275 15479538 : Elem * ReplicatedMesh::add_elem (Elem * e)
276 : {
277 1150946 : libmesh_assert(e);
278 :
279 : // We no longer merely append elements with ReplicatedMesh
280 :
281 : // If the user requests a valid id that doesn't correspond to an
282 : // existing element, let's give them that id, resizing the elements
283 : // container if necessary.
284 15479538 : if (!e->valid_id())
285 824884 : e->set_id (cast_int<dof_id_type>(_elements.size()));
286 :
287 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
288 15479538 : if (!e->valid_unique_id())
289 5239062 : e->set_unique_id(_next_unique_id++);
290 : else
291 11099693 : _next_unique_id = std::max(_next_unique_id, e->unique_id()+1);
292 : #endif
293 :
294 2291004 : const dof_id_type id = e->id();
295 :
296 16619596 : if (id < _elements.size())
297 : {
298 : // This should *almost* never happen, but we rely on it when
299 : // using allgather to replicate a not-yet-actually-replicated
300 : // ReplicatedMesh under construction in parallel.
301 169654 : if (e == _elements[id])
302 0 : return e;
303 :
304 : // Overwriting existing elements is still probably a mistake.
305 17233 : libmesh_assert(!_elements[id]);
306 : }
307 : else
308 : {
309 15309884 : _elements.resize(id+1, nullptr);
310 : }
311 :
312 15479538 : ++_n_elem;
313 15479538 : _elements[id] = e;
314 :
315 : // We actually added a new element. Some of our caches might still
316 : // be valid, but we should clear the ones which definitely are not.
317 15479538 : this->clear_point_locator();
318 15479538 : this->clear_stored_ranges();
319 :
320 : // Make sure any new element is given space for any extra integers
321 : // we've requested
322 15479538 : e->add_extra_integers(_elem_integer_names.size(),
323 15479538 : _elem_integer_default_values);
324 :
325 : // And set mapping type and data on any new element
326 2291004 : e->set_mapping_type(this->default_mapping_type());
327 2291004 : e->set_mapping_data(this->default_mapping_data());
328 :
329 15479538 : return e;
330 : }
331 :
332 14665225 : Elem * ReplicatedMesh::add_elem (std::unique_ptr<Elem> e)
333 : {
334 : // The mesh now takes ownership of the Elem. Eventually the guts of
335 : // add_elem() will get moved to a private helper function, and
336 : // calling add_elem() directly will be deprecated.
337 14665225 : return add_elem(e.release());
338 : }
339 :
340 :
341 :
342 630413 : Elem * ReplicatedMesh::insert_elem (Elem * e)
343 : {
344 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
345 630413 : if (!e->valid_unique_id())
346 18176 : e->set_unique_id(_next_unique_id++);
347 : else
348 612237 : _next_unique_id = std::max(_next_unique_id, e->unique_id()+1);
349 : #endif
350 :
351 262692 : dof_id_type eid = e->id();
352 131346 : libmesh_assert_less (eid, _elements.size());
353 630413 : Elem * oldelem = _elements[eid];
354 :
355 630413 : if (oldelem)
356 : {
357 131346 : libmesh_assert_equal_to (oldelem->id(), eid);
358 630413 : this->delete_elem(oldelem);
359 : }
360 :
361 630413 : ++_n_elem;
362 630413 : _elements[eid] = e;
363 :
364 : // We actually added a new element. Some of our caches might still
365 : // be valid, but we should clear the ones which definitely are not.
366 630413 : this->clear_point_locator();
367 630413 : this->clear_stored_ranges();
368 :
369 : // Make sure any new element is given space for any extra integers
370 : // we've requested
371 630413 : e->add_extra_integers(_elem_integer_names.size(),
372 630413 : _elem_integer_default_values);
373 :
374 : // And set mapping type and data on any new element
375 262692 : e->set_mapping_type(this->default_mapping_type());
376 262692 : e->set_mapping_data(this->default_mapping_data());
377 :
378 630413 : return e;
379 : }
380 :
381 630413 : Elem * ReplicatedMesh::insert_elem (std::unique_ptr<Elem> e)
382 : {
383 : // The mesh now takes ownership of the Elem. Eventually the guts of
384 : // insert_elem(Elem*) will get moved to a private helper function, and
385 : // calling insert_elem(Elem*) directly will be deprecated.
386 630413 : return insert_elem(e.release());
387 : }
388 :
389 :
390 :
391 2450698 : void ReplicatedMesh::delete_elem(Elem * e)
392 : {
393 244175 : libmesh_assert(e);
394 :
395 : // Initialize an iterator to eventually point to the element we want to delete
396 244175 : std::vector<Elem *>::iterator pos = _elements.end();
397 :
398 : // In many cases, e->id() gives us a clue as to where e
399 : // is located in the _elements vector. Try that first
400 : // before trying the O(n_elem) search.
401 244175 : libmesh_assert_less (e->id(), _elements.size());
402 :
403 2694873 : if (_elements[e->id()] == e)
404 : {
405 : // We found it!
406 2206523 : pos = _elements.begin();
407 244175 : std::advance(pos, e->id());
408 : }
409 :
410 : else
411 : {
412 : // This search is O(n_elem)
413 0 : pos = std::find (_elements.begin(),
414 : _elements.end(),
415 0 : e);
416 : }
417 :
418 : // Huh? Element not in the vector?
419 244175 : libmesh_assert (pos != _elements.end());
420 :
421 : // Remove the element from the BoundaryInfo object
422 2450698 : this->get_boundary_info().remove(e);
423 :
424 : // delete the element
425 2450698 : --_n_elem;
426 2450698 : delete e;
427 :
428 : // explicitly zero the pointer
429 2450698 : *pos = nullptr;
430 :
431 : // Some of our caches might still be valid, but we should clear the
432 : // ones which definitely are not.
433 2450698 : this->clear_point_locator();
434 2450698 : this->clear_stored_ranges();
435 2450698 : }
436 :
437 :
438 :
439 21420 : void ReplicatedMesh::renumber_elem(const dof_id_type old_id,
440 : const dof_id_type new_id)
441 : {
442 : // This could be a no-op
443 21420 : if (old_id == new_id)
444 0 : return;
445 :
446 : // This doesn't get used in serial yet
447 21420 : Elem * el = _elements[old_id];
448 7140 : libmesh_assert (el);
449 :
450 28560 : if (new_id >= _elements.size())
451 288 : _elements.resize(new_id+1, nullptr);
452 :
453 7140 : el->set_id(new_id);
454 7140 : libmesh_assert (!_elements[new_id]);
455 28560 : _elements[new_id] = el;
456 21420 : _elements[old_id] = nullptr;
457 :
458 : // Should we delete any caches here? Our point locator indexes by
459 : // element pointer and should be fine with an id change. Our stored
460 : // ranges are no longer sorted, which is *probably* fine, but let's
461 : // just be safe.
462 21420 : this->clear_stored_ranges();
463 : }
464 :
465 :
466 :
467 11673471 : Node * ReplicatedMesh::add_point (const Point & p,
468 : const dof_id_type id,
469 : const processor_id_type proc_id)
470 : {
471 1830270 : Node * n = nullptr;
472 :
473 : // If the user requests a valid id, either
474 : // provide the existing node or resize the container
475 : // to fit the new node.
476 11673471 : if (id != DofObject::invalid_id)
477 7727330 : if (id < _nodes.size())
478 14932 : n = _nodes[id];
479 : else
480 6779612 : _nodes.resize(id+1);
481 : else
482 4878927 : _nodes.push_back (static_cast<Node *>(nullptr));
483 :
484 : // if the node already exists, then assign new (x,y,z) values
485 1843039 : if (n)
486 2 : *n = p;
487 : // otherwise build a new node, put it in the right spot, and return
488 : // a valid pointer.
489 : else
490 : {
491 22414016 : n = Node::build(p, (id == DofObject::invalid_id) ?
492 7604423 : cast_int<dof_id_type>(_nodes.size()-1) : id).release();
493 11673400 : n->processor_id() = proc_id;
494 :
495 : // A new node invalidates any cached node ranges
496 11673400 : this->clear_stored_ranges();
497 :
498 11673400 : n->add_extra_integers(_node_integer_names.size(),
499 11673400 : _node_integer_default_values);
500 :
501 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
502 11673400 : if (!n->valid_unique_id())
503 11673400 : n->set_unique_id(_next_unique_id++);
504 : else
505 0 : _next_unique_id = std::max(_next_unique_id, n->unique_id()+1);
506 : #endif
507 :
508 11673400 : ++_n_nodes;
509 11673400 : if (id == DofObject::invalid_id)
510 4878927 : _nodes.back() = n;
511 : else
512 7727257 : _nodes[id] = n;
513 : }
514 :
515 : // better not pass back a nullptr.
516 1830270 : libmesh_assert (n);
517 :
518 11673471 : return n;
519 : }
520 :
521 :
522 :
523 948677 : Node * ReplicatedMesh::add_node (Node * n)
524 : {
525 96824 : libmesh_assert(n);
526 :
527 : // If the user requests a valid id, either set the existing
528 : // container entry or resize the container to fit the new node.
529 948677 : if (n->valid_id())
530 : {
531 96734 : const dof_id_type id = n->id();
532 1044895 : if (id < _nodes.size())
533 15080 : libmesh_assert(!_nodes[id]);
534 : else
535 751637 : _nodes.resize(id+1); // default nullptr
536 :
537 1044895 : _nodes[id] = n;
538 : }
539 : else
540 : {
541 180 : n->set_id (cast_int<dof_id_type>(_nodes.size()));
542 516 : _nodes.push_back(n);
543 : }
544 :
545 948677 : ++_n_nodes;
546 :
547 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
548 948677 : if (!n->valid_unique_id())
549 524 : n->set_unique_id(_next_unique_id++);
550 : else
551 1449250 : _next_unique_id = std::max(_next_unique_id, n->unique_id()+1);
552 : #endif
553 :
554 : // A new node invalidates any cached node ranges
555 948677 : this->clear_stored_ranges();
556 :
557 948677 : n->add_extra_integers(_node_integer_names.size(),
558 948677 : _node_integer_default_values);
559 :
560 948677 : return n;
561 : }
562 :
563 948677 : Node * ReplicatedMesh::add_node (std::unique_ptr<Node> n)
564 : {
565 : // The mesh now takes ownership of the Node. Eventually the guts of
566 : // add_node() will get moved to a private helper function, and
567 : // calling add_node() directly will be deprecated.
568 948677 : return add_node(n.release());
569 : }
570 :
571 250912 : void ReplicatedMesh::delete_node(Node * n)
572 : {
573 63787 : libmesh_assert(n);
574 63787 : libmesh_assert_less (n->id(), _nodes.size());
575 :
576 : // Initialize an iterator to eventually point to the element we want
577 : // to delete
578 63787 : std::vector<Node *>::iterator pos;
579 :
580 : // In many cases, e->id() gives us a clue as to where e
581 : // is located in the _elements vector. Try that first
582 : // before trying the O(n_elem) search.
583 314699 : if (_nodes[n->id()] == n)
584 : {
585 187125 : pos = _nodes.begin();
586 63787 : std::advance(pos, n->id());
587 : }
588 : else
589 : {
590 0 : pos = std::find (_nodes.begin(),
591 : _nodes.end(),
592 0 : n);
593 : }
594 :
595 : // Huh? Node not in the vector?
596 63787 : libmesh_assert (pos != _nodes.end());
597 :
598 : // Delete the node from the BoundaryInfo object
599 250912 : this->get_boundary_info().remove(n);
600 127574 : _constraint_rows.erase(n);
601 :
602 : // delete the node
603 250912 : --_n_nodes;
604 438037 : delete n;
605 :
606 : // explicitly zero the pointer
607 250912 : *pos = nullptr;
608 :
609 : // We've invalidated any ranges that contained this node
610 250912 : this->clear_stored_ranges();
611 250912 : }
612 :
613 :
614 :
615 70092 : void ReplicatedMesh::renumber_node(const dof_id_type old_id,
616 : const dof_id_type new_id)
617 : {
618 : // This could be a no-op
619 70092 : if (old_id == new_id)
620 0 : return;
621 :
622 : // This doesn't get used in serial yet
623 70092 : Node * nd = _nodes[old_id];
624 23364 : libmesh_assert (nd);
625 :
626 93456 : if (new_id >= _nodes.size())
627 288 : _nodes.resize(new_id+1, nullptr);
628 :
629 23364 : nd->set_id(new_id);
630 23364 : libmesh_assert (!_nodes[new_id]);
631 93456 : _nodes[new_id] = nd;
632 70092 : _nodes[old_id] = nullptr;
633 : }
634 :
635 :
636 :
637 140082 : void ReplicatedMesh::clear ()
638 : {
639 : // Call parent clear function
640 140082 : MeshBase::clear();
641 :
642 : // Clear our elements and nodes
643 : // There is no need to remove them from
644 : // the BoundaryInfo data structure since we
645 : // already cleared it.
646 140082 : this->ReplicatedMesh::clear_elems();
647 :
648 10908522 : for (auto & node : _nodes)
649 19820013 : delete node;
650 :
651 140082 : _n_nodes = 0;
652 27243 : _nodes.clear();
653 140082 : }
654 :
655 :
656 :
657 140794 : void ReplicatedMesh::clear_elems ()
658 : {
659 14018437 : for (auto & elem : _elements)
660 13877643 : delete elem;
661 :
662 140794 : _n_elem = 0;
663 27449 : _elements.clear();
664 140794 : }
665 :
666 :
667 :
668 255372 : void ReplicatedMesh::update_parallel_id_counts()
669 : {
670 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
671 255372 : _next_unique_id = this->parallel_max_unique_id();
672 : #endif
673 :
674 : // Implicitly get max_elem_id/max_node_id by trimming the vectors
675 510744 : auto trim_vec = [](auto & dof_vec) {
676 632990 : auto last_non_null = std::find_if(dof_vec.rbegin(), dof_vec.rend(), [](DofObject * d) { return (d != nullptr); });
677 510744 : dof_vec.resize(last_non_null.base()-dof_vec.begin());
678 510744 : };
679 :
680 255372 : trim_vec(this->_nodes);
681 255372 : trim_vec(this->_elements);
682 :
683 255372 : _n_constraint_rows = this->count_constraint_rows();
684 :
685 255372 : this->_preparation.has_synched_id_counts = true;
686 255372 : }
687 :
688 :
689 :
690 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
691 273037 : unique_id_type ReplicatedMesh::parallel_max_unique_id() const
692 : {
693 : // This function must be run on all processors at once
694 49503 : parallel_object_only();
695 :
696 273037 : unique_id_type max_local = _next_unique_id;
697 273037 : this->comm().max(max_local);
698 273037 : return max_local;
699 : }
700 :
701 :
702 :
703 21362 : void ReplicatedMesh::set_next_unique_id(unique_id_type id)
704 : {
705 21362 : _next_unique_id = id;
706 21362 : }
707 : #endif
708 :
709 :
710 :
711 206938 : void ReplicatedMesh::renumber_nodes_and_elements ()
712 : {
713 124848 : LOG_SCOPE("renumber_nodes_and_elem()", "Mesh");
714 :
715 : // node and element id counters
716 41616 : dof_id_type next_free_elem = 0;
717 41616 : dof_id_type next_free_node = 0;
718 :
719 : // Will hold the set of nodes that are currently connected to elements
720 83232 : std::unordered_set<Node *> connected_nodes;
721 :
722 : // Loop over the elements. Note that there may
723 : // be nullptrs in the _elements vector from the coarsening
724 : // process. Pack the elements in to a contiguous array
725 : // and then trim any excess.
726 : {
727 41616 : std::vector<Elem *>::iterator in = _elements.begin();
728 41616 : std::vector<Elem *>::iterator out_iter = _elements.begin();
729 41616 : const std::vector<Elem *>::iterator end = _elements.end();
730 :
731 40063655 : for (; in != end; ++in)
732 39856717 : if (*in != nullptr)
733 : {
734 3136572 : Elem * el = *in;
735 :
736 38399549 : *out_iter = *in;
737 3136572 : ++out_iter;
738 :
739 : // Increment the element counter
740 38399549 : el->set_id (next_free_elem++);
741 :
742 38399549 : if (_skip_renumber_nodes_and_elements)
743 : {
744 : // Add this elements nodes to the connected list
745 2145255 : for (auto & n : el->node_ref_range())
746 1934776 : connected_nodes.insert(&n);
747 : }
748 : else // We DO want node renumbering
749 : {
750 : // Loop over this element's nodes. Number them,
751 : // if they have not been numbered already. Also,
752 : // position them in the _nodes vector so that they
753 : // are packed contiguously from the beginning.
754 236845979 : for (auto & n : el->node_ref_range())
755 198656909 : if (n.id() == next_free_node) // don't need to process
756 37726765 : next_free_node++; // [(src == dst) below]
757 :
758 160930144 : else if (n.id() > next_free_node) // need to process
759 : {
760 : // The source and destination indices
761 : // for this node
762 1209428 : const dof_id_type src_idx = n.id();
763 9057502 : const dof_id_type dst_idx = next_free_node++;
764 :
765 : // ensure we want to swap a valid nodes
766 1209428 : libmesh_assert(_nodes[src_idx]);
767 :
768 : // Swap the source and destination nodes
769 2991213 : std::swap(_nodes[src_idx],
770 2991213 : _nodes[dst_idx] );
771 :
772 : // Set proper indices where that makes sense
773 9057502 : if (_nodes[src_idx] != nullptr)
774 1194670 : _nodes[src_idx]->set_id (src_idx);
775 1209428 : _nodes[dst_idx]->set_id (dst_idx);
776 : }
777 : }
778 : }
779 :
780 : // Erase any additional storage. These elements have been
781 : // copied into nullptr voids by the procedure above, and are
782 : // thus repeated and unnecessary.
783 206938 : _elements.erase (out_iter, end);
784 : }
785 :
786 :
787 206938 : if (_skip_renumber_nodes_and_elements)
788 : {
789 : // Loop over the nodes. Note that there may
790 : // be nullptrs in the _nodes vector from the coarsening
791 : // process. Pack the nodes in to a contiguous array
792 : // and then trim any excess.
793 :
794 60 : std::vector<Node *>::iterator in = _nodes.begin();
795 60 : std::vector<Node *>::iterator out_iter = _nodes.begin();
796 60 : const std::vector<Node *>::iterator end = _nodes.end();
797 :
798 925936 : for (; in != end; ++in)
799 925541 : if (*in != nullptr)
800 : {
801 : // This is a reference so that if we change the pointer it will change in the vector
802 209684 : Node * & nd = *in;
803 :
804 : // If this node is still connected to an elem, put it in the list
805 419368 : if (connected_nodes.count(nd))
806 : {
807 673877 : *out_iter = nd;
808 137780 : ++out_iter;
809 :
810 : // Increment the node counter
811 673877 : nd->set_id (next_free_node++);
812 : }
813 : else // This node is orphaned, delete it!
814 : {
815 251664 : this->get_boundary_info().remove (nd);
816 143808 : _constraint_rows.erase(nd);
817 :
818 : // delete the node
819 251664 : --_n_nodes;
820 431424 : delete nd;
821 251664 : nd = nullptr;
822 : }
823 : }
824 :
825 : // Erase any additional storage. Whatever was
826 395 : _nodes.erase (out_iter, end);
827 : }
828 : else // We really DO want node renumbering
829 : {
830 : // Any nodes in the vector >= _nodes[next_free_node]
831 : // are not connected to any elements and may be deleted
832 : // if desired.
833 :
834 : // Now, delete the unused nodes
835 : {
836 164987 : std::vector<Node *>::iterator nd = _nodes.begin();
837 41556 : const std::vector<Node *>::iterator end = _nodes.end();
838 :
839 41556 : std::advance (nd, next_free_node);
840 :
841 1735589 : for (auto & node : as_range(nd, end))
842 : {
843 : // Mesh modification code might have already deleted some
844 : // nodes
845 1529046 : if (node == nullptr)
846 104431 : continue;
847 :
848 : // remove any boundary information associated with
849 : // this node
850 1389113 : this->get_boundary_info().remove (node);
851 225172 : _constraint_rows.erase(node);
852 :
853 : // delete the node
854 1389113 : --_n_nodes;
855 2665640 : delete node;
856 1389113 : node = nullptr;
857 : }
858 :
859 206543 : _nodes.erase (nd, end);
860 : }
861 : }
862 :
863 : // We've probably done *some* node removal and/or renumbering; let's
864 : // clear caches to be safe.
865 206938 : this->clear_stored_ranges();
866 :
867 206938 : this->_preparation.has_removed_orphaned_nodes = true;
868 :
869 41616 : libmesh_assert_equal_to (next_free_elem, _elements.size());
870 41616 : libmesh_assert_equal_to (next_free_node, _nodes.size());
871 :
872 206938 : this->update_parallel_id_counts();
873 206938 : }
874 :
875 :
876 :
877 2220 : void ReplicatedMesh::fix_broken_node_and_element_numbering ()
878 : {
879 : // Nodes first
880 3180890 : for (auto n : index_range(_nodes))
881 3178670 : if (this->_nodes[n] != nullptr)
882 3178670 : this->_nodes[n]->set_id() = cast_int<dof_id_type>(n);
883 :
884 : // Elements next
885 3664918 : for (auto e : index_range(_elements))
886 3662698 : if (this->_elements[e] != nullptr)
887 3662698 : this->_elements[e]->set_id() = cast_int<dof_id_type>(e);
888 2220 : }
889 :
890 :
891 121297 : dof_id_type ReplicatedMesh::n_active_elem () const
892 : {
893 221707 : return static_cast<dof_id_type>(std::distance (this->active_elements_begin(),
894 343004 : this->active_elements_end()));
895 : }
896 :
897 : std::vector<dof_id_type>
898 142 : ReplicatedMesh::get_disconnected_subdomains(std::vector<subdomain_id_type> * subdomain_ids) const
899 : {
900 : // find number of disconnected subdomains
901 4 : std::vector<dof_id_type> representative_elem_ids;
902 :
903 : // use subdomain_ids as markers for all elements to indicate if the elements
904 : // have been visited. Note: here subdomain ID is unrelated with element
905 : // subdomain_id().
906 8 : std::vector<subdomain_id_type> subdomains;
907 142 : if (!subdomain_ids)
908 0 : subdomain_ids = &subdomains;
909 4 : subdomain_ids->clear();
910 142 : subdomain_ids->resize(max_elem_id() + 1, Elem::invalid_subdomain_id);
911 :
912 : // counter of disconnected subdomains
913 4 : subdomain_id_type subdomain_counter = 0;
914 :
915 : // a stack for visiting elements, make its capacity sufficiently large to avoid
916 : // memory allocation and deallocation when the vector size changes
917 8 : std::vector<const Elem *> list;
918 142 : list.reserve(n_elem());
919 :
920 : // counter of visited elements
921 4 : dof_id_type visited = 0;
922 142 : dof_id_type n_active = n_active_elem();
923 4 : do
924 : {
925 1672 : for (const auto & elem : active_element_ptr_range())
926 876 : if ((*subdomain_ids)[elem->id()] == Elem::invalid_subdomain_id)
927 : {
928 284 : list.push_back(elem);
929 284 : (*subdomain_ids)[elem->id()] = subdomain_counter;
930 284 : break;
931 268 : }
932 : // we should be able to find a seed here
933 8 : libmesh_assert(list.size() > 0);
934 :
935 284 : dof_id_type min_id = std::numeric_limits<dof_id_type>::max();
936 1988 : while (list.size() > 0)
937 : {
938 : // pop up an element
939 1704 : const Elem * elem = list.back(); list.pop_back(); ++visited;
940 :
941 3084 : min_id = std::min(elem->id(), min_id);
942 :
943 8520 : for (auto s : elem->side_index_range())
944 : {
945 6816 : const Elem * neighbor = elem->neighbor_ptr(s);
946 6816 : if (neighbor != nullptr && (*subdomain_ids)[neighbor->id()] == Elem::invalid_subdomain_id)
947 : {
948 : // neighbor must be active
949 40 : libmesh_assert(neighbor->active());
950 1420 : list.push_back(neighbor);
951 1460 : (*subdomain_ids)[neighbor->id()] = subdomain_counter;
952 : }
953 : }
954 : }
955 :
956 284 : representative_elem_ids.push_back(min_id);
957 284 : subdomain_counter++;
958 : }
959 284 : while (visited != n_active);
960 :
961 146 : return representative_elem_ids;
962 : }
963 :
964 : std::unordered_map<dof_id_type, std::vector<std::vector<Point>>>
965 142 : ReplicatedMesh::get_boundary_points() const
966 : {
967 142 : libmesh_error_msg_if(mesh_dimension() != 2,
968 : "Error: get_boundary_points only works for 2D now");
969 :
970 : // find number of disconnected subdomains
971 : // subdomains will hold the IDs of disconnected subdomains for all elements.
972 8 : std::vector<subdomain_id_type> subdomains;
973 146 : std::vector<dof_id_type> elem_ids = get_disconnected_subdomains(&subdomains);
974 :
975 4 : std::unordered_map<dof_id_type, std::vector<std::vector<Point>>> boundary_points;
976 :
977 : // get all boundary sides that are to be erased later during visiting
978 : // use a comparison functor to avoid run-time randomness due to pointers
979 : struct boundary_side_compare
980 : {
981 1584 : bool operator()(const std::pair<const Elem *, unsigned int> & lhs,
982 : const std::pair<const Elem *, unsigned int> & rhs) const
983 : {
984 42340 : if (lhs.first->id() < rhs.first->id())
985 328 : return true;
986 29298 : else if (lhs.first->id() == rhs.first->id())
987 : {
988 14668 : if (lhs.second < rhs.second)
989 112 : return true;
990 : }
991 1144 : return false;
992 : }
993 : };
994 8 : std::set<std::pair<const Elem *, unsigned int>, boundary_side_compare> boundary_elements;
995 3592 : for (const auto & elem : active_element_ptr_range())
996 8568 : for (auto s : elem->side_index_range())
997 7008 : if (elem->neighbor_ptr(s) == nullptr)
998 3542 : boundary_elements.insert(std::pair<const Elem *, unsigned int>(elem, s));
999 :
1000 568 : while (!boundary_elements.empty())
1001 : {
1002 : // get the first entry as the seed
1003 426 : const Elem * eseed = boundary_elements.begin()->first;
1004 426 : unsigned int sseed = boundary_elements.begin()->second;
1005 :
1006 : // get the subdomain ID that these boundary sides attached to
1007 438 : subdomain_id_type subdomain_id = subdomains[eseed->id()];
1008 :
1009 : // start visiting the mesh to find all boundary nodes with the seed
1010 24 : std::vector<Point> bpoints;
1011 426 : const Elem * elem = eseed;
1012 12 : unsigned int s = sseed;
1013 438 : std::vector<unsigned int> local_side_nodes = elem->nodes_on_side(s);
1014 : while (true)
1015 : {
1016 96 : std::pair<const Elem *, unsigned int> side(elem, s);
1017 96 : libmesh_assert(boundary_elements.count(side));
1018 96 : boundary_elements.erase(side);
1019 :
1020 : // push all nodes on the side except the node on the other end of the side (index 1)
1021 11928 : for (auto i : index_range(local_side_nodes))
1022 8520 : if (i != 1)
1023 5400 : bpoints.push_back(*static_cast<const Point *>(elem->node_ptr(local_side_nodes[i])));
1024 :
1025 : // use the last node to find next element and side
1026 3408 : const Node * node = elem->node_ptr(local_side_nodes[1]);
1027 96 : std::set<const Elem *> neighbors;
1028 3408 : elem->find_point_neighbors(*node, neighbors);
1029 :
1030 : // if only one neighbor is found (itself), this node is a cornor node on boundary
1031 3408 : if (neighbors.size() != 1)
1032 64 : neighbors.erase(elem);
1033 :
1034 : // find the connecting side
1035 96 : bool found = false;
1036 3681 : for (const auto & neighbor : neighbors)
1037 : {
1038 9990 : for (auto ss : neighbor->side_index_range())
1039 9888 : if (neighbor->neighbor_ptr(ss) == nullptr && !(elem == neighbor && s == ss))
1040 : {
1041 4948 : local_side_nodes = neighbor->nodes_on_side(ss);
1042 : // we expect the starting point of the side to be the same as the end of the previous side
1043 5090 : if (neighbor->node_ptr(local_side_nodes[0]) == node)
1044 : {
1045 3408 : elem = neighbor;
1046 96 : s = ss;
1047 96 : found = true;
1048 96 : break;
1049 : }
1050 1586 : else if (neighbor->node_ptr(local_side_nodes[1]) == node)
1051 : {
1052 0 : elem = neighbor;
1053 0 : s = ss;
1054 0 : found = true;
1055 : // flip nodes in local_side_nodes because the side is in an opposite direction
1056 0 : auto temp(local_side_nodes);
1057 0 : local_side_nodes[0] = temp[1];
1058 0 : local_side_nodes[1] = temp[0];
1059 0 : for (unsigned int i = 2; i < temp.size(); ++i)
1060 0 : local_side_nodes[temp.size() + 1 - i] = temp[i];
1061 0 : break;
1062 : }
1063 : }
1064 104 : if (found)
1065 96 : break;
1066 : }
1067 :
1068 3408 : libmesh_error_msg_if(!found, "ERROR: mesh topology error on visiting boundary sides");
1069 :
1070 : // exit if we reach the starting point
1071 3408 : if (elem == eseed && s == sseed)
1072 12 : break;
1073 84 : }
1074 426 : boundary_points[elem_ids[subdomain_id]].push_back(bpoints);
1075 : }
1076 :
1077 146 : return boundary_points;
1078 : }
1079 :
1080 : } // namespace libMesh
|
