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