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/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 48994 : ReplicatedMesh::ReplicatedMesh (const Parallel::Communicator & comm_in,
43 48994 : unsigned char d) :
44 : UnstructuredMesh (comm_in,d),
45 48994 : _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 48994 : _next_unique_id = 0;
51 : #endif
52 :
53 57588 : const std::string default_partitioner = "metis";
54 : const std::string my_partitioner =
55 : libMesh::command_line_value("--default-partitioner",
56 97988 : default_partitioner);
57 : _partitioner = Partitioner::build
58 89394 : (Utility::string_to_enum<PartitionerType>(my_partitioner));
59 48994 : }
60 :
61 :
62 1414 : bool ReplicatedMesh::subclass_locally_equals(const MeshBase & other_mesh_base) const
63 : {
64 770 : const ReplicatedMesh * rep_mesh_ptr =
65 1414 : dynamic_cast<const ReplicatedMesh *>(&other_mesh_base);
66 1414 : if (!rep_mesh_ptr)
67 0 : return false;
68 770 : const ReplicatedMesh & other_mesh = *rep_mesh_ptr;
69 :
70 3598 : if (_n_nodes != other_mesh._n_nodes ||
71 1414 : _n_elem != other_mesh._n_elem ||
72 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
73 4242 : _next_unique_id != other_mesh._next_unique_id ||
74 : #endif
75 1414 : !this->nodes_and_elements_equal(other_mesh))
76 78 : return false;
77 :
78 766 : return true;
79 : }
80 :
81 :
82 62334 : ReplicatedMesh::~ReplicatedMesh ()
83 : {
84 52777 : this->ReplicatedMesh::clear(); // Free nodes and elements
85 62334 : }
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 3712 : ReplicatedMesh::ReplicatedMesh (const ReplicatedMesh & other_mesh) :
92 3712 : ReplicatedMesh(static_cast<const MeshBase&>(other_mesh))
93 : {
94 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
95 3712 : this->_next_unique_id = other_mesh._next_unique_id;
96 : #endif
97 3712 : }
98 :
99 :
100 3783 : ReplicatedMesh::ReplicatedMesh (const MeshBase & other_mesh) :
101 : UnstructuredMesh (other_mesh),
102 3783 : _n_nodes(0), _n_elem(0) // copy_* will increment this
103 : {
104 3783 : this->copy_nodes_and_elements(other_mesh, true);
105 :
106 1372 : this->allow_find_neighbors(other_mesh.allow_find_neighbors());
107 1372 : this->allow_renumbering(other_mesh.allow_renumbering());
108 1372 : this->allow_remote_element_removal(other_mesh.allow_remote_element_removal());
109 1372 : this->skip_partitioning(other_mesh.skip_partitioning());
110 :
111 : // The prepare_for_use() in copy_nodes_and_elements() is going to be
112 : // tricky to remove without breaking backwards compatibility, but it
113 : // updates some things we want to just copy.
114 3783 : this->copy_cached_data(other_mesh);
115 :
116 3783 : this->copy_constraint_rows(other_mesh);
117 :
118 3783 : this->_is_prepared = other_mesh.is_prepared();
119 :
120 694 : auto & this_boundary_info = this->get_boundary_info();
121 694 : const auto & other_boundary_info = other_mesh.get_boundary_info();
122 :
123 3783 : this_boundary_info = other_boundary_info;
124 :
125 694 : this->set_subdomain_name_map() = other_mesh.get_subdomain_name_map();
126 :
127 : // If other_mesh is distributed, then we've got parts of it on each
128 : // processor but we're not replicated yet; fix that.
129 3783 : if (!other_mesh.is_serial())
130 71 : MeshCommunication().allgather(*this);
131 3783 : }
132 :
133 156 : ReplicatedMesh & ReplicatedMesh::operator= (ReplicatedMesh && other_mesh)
134 : {
135 8 : LOG_SCOPE("operator=(&&)", "ReplicatedMesh");
136 :
137 : // Move assign as an UnstructuredMesh
138 8 : this->UnstructuredMesh::operator=(std::move(other_mesh));
139 :
140 : // Nodes and elements belong to ReplicatedMesh and have to be
141 : // moved before we can move arbitrary GhostingFunctor, Partitioner,
142 : // etc. subclasses.
143 156 : this->move_nodes_and_elements(std::move(other_mesh));
144 :
145 : // Handle those remaining moves.
146 156 : this->post_dofobject_moves(std::move(other_mesh));
147 :
148 164 : return *this;
149 : }
150 :
151 156 : MeshBase & ReplicatedMesh::assign(MeshBase && other_mesh)
152 : {
153 156 : *this = std::move(cast_ref<ReplicatedMesh&>(other_mesh));
154 :
155 156 : return *this;
156 : }
157 :
158 156 : void ReplicatedMesh::move_nodes_and_elements(MeshBase && other_meshbase)
159 : {
160 8 : ReplicatedMesh & other_mesh = cast_ref<ReplicatedMesh&>(other_meshbase);
161 :
162 156 : this->_nodes = std::move(other_mesh._nodes);
163 156 : this->_n_nodes = other_mesh.n_nodes();
164 :
165 156 : this->_elements = std::move(other_mesh._elements);
166 156 : this->_n_elem = other_mesh.n_elem();
167 156 : }
168 :
169 :
170 15323614 : const Point & ReplicatedMesh::point (const dof_id_type i) const
171 : {
172 15323614 : return this->node_ref(i);
173 : }
174 :
175 :
176 :
177 :
178 100853687 : const Node * ReplicatedMesh::node_ptr (const dof_id_type i) const
179 : {
180 5035399 : libmesh_assert_less (i, this->max_node_id());
181 5035399 : libmesh_assert(_nodes[i]);
182 5035399 : libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
183 :
184 105851186 : return _nodes[i];
185 : }
186 :
187 :
188 :
189 :
190 102814018 : Node * ReplicatedMesh::node_ptr (const dof_id_type i)
191 : {
192 12422848 : libmesh_assert_less (i, this->max_node_id());
193 12422848 : libmesh_assert(_nodes[i]);
194 12422848 : libmesh_assert_equal_to (_nodes[i]->id(), i); // This will change soon
195 :
196 114640131 : return _nodes[i];
197 : }
198 :
199 :
200 :
201 :
202 15892260 : const Node * ReplicatedMesh::query_node_ptr (const dof_id_type i) const
203 : {
204 15892260 : if (i >= this->max_node_id())
205 2 : return nullptr;
206 15765053 : libmesh_assert (_nodes[i] == nullptr ||
207 : _nodes[i]->id() == i); // This will change soon
208 :
209 15930368 : return _nodes[i];
210 : }
211 :
212 :
213 :
214 :
215 1565962 : Node * ReplicatedMesh::query_node_ptr (const dof_id_type i)
216 : {
217 1565962 : if (i >= this->max_node_id())
218 59038 : return nullptr;
219 100713 : libmesh_assert (_nodes[i] == nullptr ||
220 : _nodes[i]->id() == i); // This will change soon
221 :
222 981341 : return _nodes[i];
223 : }
224 :
225 :
226 :
227 :
228 2305519 : const Elem * ReplicatedMesh::elem_ptr (const dof_id_type i) const
229 : {
230 497460 : libmesh_assert_less (i, this->max_elem_id());
231 497460 : libmesh_assert(_elements[i]);
232 497460 : libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
233 :
234 2703868 : return _elements[i];
235 : }
236 :
237 :
238 :
239 :
240 77624096 : Elem * ReplicatedMesh::elem_ptr (const dof_id_type i)
241 : {
242 4853499 : libmesh_assert_less (i, this->max_elem_id());
243 4853499 : libmesh_assert(_elements[i]);
244 4853499 : libmesh_assert_equal_to (_elements[i]->id(), i); // This will change soon
245 :
246 81845915 : return _elements[i];
247 : }
248 :
249 :
250 :
251 :
252 15515357 : const Elem * ReplicatedMesh::query_elem_ptr (const dof_id_type i) const
253 : {
254 15515357 : if (i >= this->max_elem_id())
255 4 : return nullptr;
256 15470536 : libmesh_assert (_elements[i] == nullptr ||
257 : _elements[i]->id() == i); // This will change soon
258 :
259 15529816 : return _elements[i];
260 : }
261 :
262 :
263 :
264 :
265 1414850 : Elem * ReplicatedMesh::query_elem_ptr (const dof_id_type i)
266 : {
267 1414850 : if (i >= this->max_elem_id())
268 42891 : return nullptr;
269 188466 : libmesh_assert (_elements[i] == nullptr ||
270 : _elements[i]->id() == i); // This will change soon
271 :
272 765367 : return _elements[i];
273 : }
274 :
275 :
276 :
277 :
278 15273193 : Elem * ReplicatedMesh::add_elem (Elem * e)
279 : {
280 1114768 : libmesh_assert(e);
281 :
282 : // We no longer merely append elements with ReplicatedMesh
283 :
284 : // If the user requests a valid id that doesn't correspond to an
285 : // existing element, let's give them that id, resizing the elements
286 : // container if necessary.
287 15273193 : if (!e->valid_id())
288 793044 : e->set_id (cast_int<dof_id_type>(_elements.size()));
289 :
290 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
291 15273193 : if (!e->valid_unique_id())
292 4581042 : e->set_unique_id(_next_unique_id++);
293 : else
294 11525956 : _next_unique_id = std::max(_next_unique_id, e->unique_id()+1);
295 : #endif
296 :
297 2218716 : const dof_id_type id = e->id();
298 :
299 16377141 : if (id < _elements.size())
300 : {
301 : // This should *almost* never happen, but we rely on it when
302 : // using allgather to replicate a not-yet-actually-replicated
303 : // ReplicatedMesh under construction in parallel.
304 144798 : if (e == _elements[id])
305 0 : return e;
306 :
307 : // Overwriting existing elements is still probably a mistake.
308 12996 : libmesh_assert(!_elements[id]);
309 : }
310 : else
311 : {
312 15128395 : _elements.resize(id+1, nullptr);
313 : }
314 :
315 15273193 : ++_n_elem;
316 15273193 : _elements[id] = e;
317 :
318 : // Make sure any new element is given space for any extra integers
319 : // we've requested
320 15273193 : e->add_extra_integers(_elem_integer_names.size(),
321 15273193 : _elem_integer_default_values);
322 :
323 : // And set mapping type and data on any new element
324 2218716 : e->set_mapping_type(this->default_mapping_type());
325 2218716 : e->set_mapping_data(this->default_mapping_data());
326 :
327 15273193 : return e;
328 : }
329 :
330 14465708 : Elem * ReplicatedMesh::add_elem (std::unique_ptr<Elem> e)
331 : {
332 : // The mesh now takes ownership of the Elem. Eventually the guts of
333 : // add_elem() will get moved to a private helper function, and
334 : // calling add_elem() directly will be deprecated.
335 14465708 : return add_elem(e.release());
336 : }
337 :
338 :
339 :
340 554231 : Elem * ReplicatedMesh::insert_elem (Elem * e)
341 : {
342 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
343 554231 : if (!e->valid_unique_id())
344 18176 : e->set_unique_id(_next_unique_id++);
345 : else
346 536055 : _next_unique_id = std::max(_next_unique_id, e->unique_id()+1);
347 : #endif
348 :
349 239108 : dof_id_type eid = e->id();
350 119554 : libmesh_assert_less (eid, _elements.size());
351 554231 : Elem * oldelem = _elements[eid];
352 :
353 554231 : if (oldelem)
354 : {
355 119554 : libmesh_assert_equal_to (oldelem->id(), eid);
356 554231 : this->delete_elem(oldelem);
357 : }
358 :
359 554231 : ++_n_elem;
360 554231 : _elements[eid] = e;
361 :
362 : // Make sure any new element is given space for any extra integers
363 : // we've requested
364 554231 : e->add_extra_integers(_elem_integer_names.size(),
365 554231 : _elem_integer_default_values);
366 :
367 : // And set mapping type and data on any new element
368 239108 : e->set_mapping_type(this->default_mapping_type());
369 239108 : e->set_mapping_data(this->default_mapping_data());
370 :
371 554231 : return e;
372 : }
373 :
374 554231 : Elem * ReplicatedMesh::insert_elem (std::unique_ptr<Elem> e)
375 : {
376 : // The mesh now takes ownership of the Elem. Eventually the guts of
377 : // insert_elem(Elem*) will get moved to a private helper function, and
378 : // calling insert_elem(Elem*) directly will be deprecated.
379 554231 : return insert_elem(e.release());
380 : }
381 :
382 :
383 :
384 2007902 : void ReplicatedMesh::delete_elem(Elem * e)
385 : {
386 217195 : libmesh_assert(e);
387 :
388 : // Initialize an iterator to eventually point to the element we want to delete
389 217195 : std::vector<Elem *>::iterator pos = _elements.end();
390 :
391 : // In many cases, e->id() gives us a clue as to where e
392 : // is located in the _elements vector. Try that first
393 : // before trying the O(n_elem) search.
394 217195 : libmesh_assert_less (e->id(), _elements.size());
395 :
396 2225097 : if (_elements[e->id()] == e)
397 : {
398 : // We found it!
399 1790707 : pos = _elements.begin();
400 217195 : std::advance(pos, e->id());
401 : }
402 :
403 : else
404 : {
405 : // This search is O(n_elem)
406 0 : pos = std::find (_elements.begin(),
407 : _elements.end(),
408 0 : e);
409 : }
410 :
411 : // Huh? Element not in the vector?
412 217195 : libmesh_assert (pos != _elements.end());
413 :
414 : // Remove the element from the BoundaryInfo object
415 2007902 : this->get_boundary_info().remove(e);
416 :
417 : // delete the element
418 2007902 : --_n_elem;
419 2007902 : delete e;
420 :
421 : // explicitly zero the pointer
422 2007902 : *pos = nullptr;
423 2007902 : }
424 :
425 :
426 :
427 0 : void ReplicatedMesh::renumber_elem(const dof_id_type old_id,
428 : const dof_id_type new_id)
429 : {
430 : // This could be a no-op
431 0 : if (old_id == new_id)
432 0 : return;
433 :
434 : // This doesn't get used in serial yet
435 0 : Elem * el = _elements[old_id];
436 0 : libmesh_assert (el);
437 :
438 0 : if (new_id >= _elements.size())
439 0 : _elements.resize(new_id+1, nullptr);
440 :
441 0 : el->set_id(new_id);
442 0 : libmesh_assert (!_elements[new_id]);
443 0 : _elements[new_id] = el;
444 0 : _elements[old_id] = nullptr;
445 : }
446 :
447 :
448 :
449 10843485 : Node * ReplicatedMesh::add_point (const Point & p,
450 : const dof_id_type id,
451 : const processor_id_type proc_id)
452 : {
453 1751405 : Node * n = nullptr;
454 :
455 : // If the user requests a valid id, either
456 : // provide the existing node or resize the container
457 : // to fit the new node.
458 10843485 : if (id != DofObject::invalid_id)
459 7674858 : if (id < _nodes.size())
460 10309 : n = _nodes[id];
461 : else
462 6747872 : _nodes.resize(id+1);
463 : else
464 4085304 : _nodes.push_back (static_cast<Node *>(nullptr));
465 :
466 : // if the node already exists, then assign new (x,y,z) values
467 1760864 : if (n)
468 0 : *n = p;
469 : // otherwise build a new node, put it in the right spot, and return
470 : // a valid pointer.
471 : else
472 : {
473 20770289 : n = Node::build(p, (id == DofObject::invalid_id) ?
474 6669297 : cast_int<dof_id_type>(_nodes.size()-1) : id).release();
475 10843485 : n->processor_id() = proc_id;
476 :
477 10843485 : n->add_extra_integers(_node_integer_names.size(),
478 10843485 : _node_integer_default_values);
479 :
480 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
481 10843485 : if (!n->valid_unique_id())
482 10843485 : n->set_unique_id(_next_unique_id++);
483 : else
484 0 : _next_unique_id = std::max(_next_unique_id, n->unique_id()+1);
485 : #endif
486 :
487 10843485 : ++_n_nodes;
488 10843485 : if (id == DofObject::invalid_id)
489 4085304 : _nodes.back() = n;
490 : else
491 7674858 : _nodes[id] = n;
492 : }
493 :
494 : // better not pass back a nullptr.
495 1751405 : libmesh_assert (n);
496 :
497 10843485 : return n;
498 : }
499 :
500 :
501 :
502 697761 : Node * ReplicatedMesh::add_node (Node * n)
503 : {
504 62046 : libmesh_assert(n);
505 :
506 : // If the user requests a valid id, either set the existing
507 : // container entry or resize the container to fit the new node.
508 697761 : if (n->valid_id())
509 : {
510 62046 : const dof_id_type id = n->id();
511 759807 : if (id < _nodes.size())
512 6710 : libmesh_assert(!_nodes[id]);
513 : else
514 625722 : _nodes.resize(id+1); // default nullptr
515 :
516 759807 : _nodes[id] = n;
517 : }
518 : else
519 : {
520 0 : n->set_id (cast_int<dof_id_type>(_nodes.size()));
521 0 : _nodes.push_back(n);
522 : }
523 :
524 697761 : ++_n_nodes;
525 :
526 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
527 697761 : if (!n->valid_unique_id())
528 0 : n->set_unique_id(_next_unique_id++);
529 : else
530 1158512 : _next_unique_id = std::max(_next_unique_id, n->unique_id()+1);
531 : #endif
532 :
533 697761 : n->add_extra_integers(_node_integer_names.size(),
534 697761 : _node_integer_default_values);
535 :
536 697761 : return n;
537 : }
538 :
539 697761 : Node * ReplicatedMesh::add_node (std::unique_ptr<Node> n)
540 : {
541 : // The mesh now takes ownership of the Node. Eventually the guts of
542 : // add_node() will get moved to a private helper function, and
543 : // calling add_node() directly will be deprecated.
544 697761 : return add_node(n.release());
545 : }
546 :
547 : #ifdef LIBMESH_ENABLE_DEPRECATED
548 :
549 0 : Node * ReplicatedMesh::insert_node(Node * n)
550 : {
551 : libmesh_deprecated();
552 0 : libmesh_error_msg_if(!n, "Error, attempting to insert nullptr node.");
553 0 : libmesh_error_msg_if(n->id() == DofObject::invalid_id, "Error, cannot insert node with invalid id.");
554 :
555 0 : if (n->id() < _nodes.size())
556 : {
557 : // Trying to add an existing node is a no-op. This lets us use
558 : // a straightforward MeshCommunication::allgather() to fix a
559 : // not-actually-replicated-yet ReplicatedMesh, as occurs when
560 : // reading Nemesis files.
561 0 : if (n->valid_id() && _nodes[n->id()] == n)
562 0 : return n;
563 :
564 : // Don't allow anything else to replace an existing Node.
565 0 : libmesh_error_msg_if(_nodes[ n->id() ] != nullptr,
566 : "Error, cannot insert new node on top of existing node.");
567 : }
568 : else
569 : {
570 : // Allocate just enough space to store the new node. This will
571 : // cause highly non-ideal memory allocation behavior if called
572 : // repeatedly...
573 0 : _nodes.resize(n->id() + 1);
574 : }
575 :
576 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
577 0 : if (!n->valid_unique_id())
578 0 : n->set_unique_id(_next_unique_id++);
579 : else
580 0 : _next_unique_id = std::max(_next_unique_id, n->unique_id()+1);
581 : #endif
582 :
583 0 : n->add_extra_integers(_node_integer_names.size(),
584 0 : _node_integer_default_values);
585 :
586 : // We have enough space and this spot isn't already occupied by
587 : // another node, so go ahead and add it.
588 0 : ++_n_nodes;
589 0 : _nodes[ n->id() ] = n;
590 :
591 : // If we made it this far, we just inserted the node the user handed
592 : // us, so we can give it right back.
593 0 : return n;
594 : }
595 :
596 0 : Node * ReplicatedMesh::insert_node(std::unique_ptr<Node> n)
597 : {
598 : libmesh_deprecated();
599 : // The mesh now takes ownership of the Node. Eventually the guts of
600 : // insert_node(Node*) will get moved to a private helper function, and
601 : // calling insert_node(Node*) directly will be deprecated.
602 0 : return insert_node(n.release());
603 : }
604 :
605 : #endif
606 :
607 250701 : void ReplicatedMesh::delete_node(Node * n)
608 : {
609 63997 : libmesh_assert(n);
610 63997 : libmesh_assert_less (n->id(), _nodes.size());
611 :
612 : // Initialize an iterator to eventually point to the element we want
613 : // to delete
614 63997 : std::vector<Node *>::iterator pos;
615 :
616 : // In many cases, e->id() gives us a clue as to where e
617 : // is located in the _elements vector. Try that first
618 : // before trying the O(n_elem) search.
619 314698 : if (_nodes[n->id()] == n)
620 : {
621 186704 : pos = _nodes.begin();
622 63997 : std::advance(pos, n->id());
623 : }
624 : else
625 : {
626 0 : pos = std::find (_nodes.begin(),
627 : _nodes.end(),
628 0 : n);
629 : }
630 :
631 : // Huh? Node not in the vector?
632 63997 : libmesh_assert (pos != _nodes.end());
633 :
634 : // Delete the node from the BoundaryInfo object
635 250701 : this->get_boundary_info().remove(n);
636 127994 : _constraint_rows.erase(n);
637 :
638 : // delete the node
639 250701 : --_n_nodes;
640 437405 : delete n;
641 :
642 : // explicitly zero the pointer
643 250701 : *pos = nullptr;
644 250701 : }
645 :
646 :
647 :
648 0 : void ReplicatedMesh::renumber_node(const dof_id_type old_id,
649 : const dof_id_type new_id)
650 : {
651 : // This could be a no-op
652 0 : if (old_id == new_id)
653 0 : return;
654 :
655 : // This doesn't get used in serial yet
656 0 : Node * nd = _nodes[old_id];
657 0 : libmesh_assert (nd);
658 :
659 0 : if (new_id >= _nodes.size())
660 0 : _nodes.resize(new_id+1, nullptr);
661 :
662 0 : nd->set_id(new_id);
663 0 : libmesh_assert (!_nodes[new_id]);
664 0 : _nodes[new_id] = nd;
665 0 : _nodes[old_id] = nullptr;
666 : }
667 :
668 :
669 :
670 94746 : void ReplicatedMesh::clear ()
671 : {
672 : // Call parent clear function
673 94746 : MeshBase::clear();
674 :
675 : // Clear our elements and nodes
676 : // There is no need to remove them from
677 : // the BoundaryInfo data structure since we
678 : // already cleared it.
679 94746 : this->ReplicatedMesh::clear_elems();
680 :
681 10330200 : for (auto & node : _nodes)
682 18835884 : delete node;
683 :
684 94746 : _n_nodes = 0;
685 18024 : _nodes.clear();
686 94746 : }
687 :
688 :
689 :
690 95398 : void ReplicatedMesh::clear_elems ()
691 : {
692 13994839 : for (auto & elem : _elements)
693 13899441 : delete elem;
694 :
695 95398 : _n_elem = 0;
696 18212 : _elements.clear();
697 95398 : }
698 :
699 :
700 :
701 153816 : void ReplicatedMesh::update_parallel_id_counts()
702 : {
703 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
704 153816 : _next_unique_id = this->parallel_max_unique_id();
705 : #endif
706 153816 : }
707 :
708 :
709 :
710 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
711 168048 : unique_id_type ReplicatedMesh::parallel_max_unique_id() const
712 : {
713 : // This function must be run on all processors at once
714 27622 : parallel_object_only();
715 :
716 168048 : unique_id_type max_local = _next_unique_id;
717 168048 : this->comm().max(max_local);
718 168048 : return max_local;
719 : }
720 :
721 :
722 :
723 10376 : void ReplicatedMesh::set_next_unique_id(unique_id_type id)
724 : {
725 10376 : _next_unique_id = id;
726 10376 : }
727 : #endif
728 :
729 :
730 :
731 124263 : void ReplicatedMesh::renumber_nodes_and_elements ()
732 : {
733 44448 : LOG_SCOPE("renumber_nodes_and_elem()", "Mesh");
734 :
735 : // node and element id counters
736 22224 : dof_id_type next_free_elem = 0;
737 22224 : dof_id_type next_free_node = 0;
738 :
739 : // Will hold the set of nodes that are currently connected to elements
740 44448 : std::unordered_set<Node *> connected_nodes;
741 :
742 : // Loop over the elements. Note that there may
743 : // be nullptrs in the _elements vector from the coarsening
744 : // process. Pack the elements in to a contiguous array
745 : // and then trim any excess.
746 : {
747 22224 : std::vector<Elem *>::iterator in = _elements.begin();
748 22224 : std::vector<Elem *>::iterator out_iter = _elements.begin();
749 22224 : const std::vector<Elem *>::iterator end = _elements.end();
750 :
751 34595300 : for (; in != end; ++in)
752 34471037 : if (*in != nullptr)
753 : {
754 2875984 : Elem * el = *in;
755 :
756 33188219 : *out_iter = *in;
757 2875984 : ++out_iter;
758 :
759 : // Increment the element counter
760 33188219 : el->set_id (next_free_elem++);
761 :
762 33188219 : if (_skip_renumber_nodes_and_elements)
763 : {
764 : // Add this elements nodes to the connected list
765 2051315 : for (auto & n : el->node_ref_range())
766 1846096 : connected_nodes.insert(&n);
767 : }
768 : else // We DO want node renumbering
769 : {
770 : // Loop over this element's nodes. Number them,
771 : // if they have not been numbered already. Also,
772 : // position them in the _nodes vector so that they
773 : // are packed contiguously from the beginning.
774 191792174 : for (auto & n : el->node_ref_range())
775 158809174 : if (n.id() == next_free_node) // don't need to process
776 25867116 : next_free_node++; // [(src == dst) below]
777 :
778 132942058 : else if (n.id() > next_free_node) // need to process
779 : {
780 : // The source and destination indices
781 : // for this node
782 1061994 : const dof_id_type src_idx = n.id();
783 7304460 : const dof_id_type dst_idx = next_free_node++;
784 :
785 : // ensure we want to swap a valid nodes
786 1061994 : libmesh_assert(_nodes[src_idx]);
787 :
788 : // Swap the source and destination nodes
789 2644164 : std::swap(_nodes[src_idx],
790 2644164 : _nodes[dst_idx] );
791 :
792 : // Set proper indices where that makes sense
793 7304460 : if (_nodes[src_idx] != nullptr)
794 1048668 : _nodes[src_idx]->set_id (src_idx);
795 1061994 : _nodes[dst_idx]->set_id (dst_idx);
796 : }
797 : }
798 : }
799 :
800 : // Erase any additional storage. These elements have been
801 : // copied into nullptr voids by the procedure above, and are
802 : // thus repeated and unnecessary.
803 124263 : _elements.erase (out_iter, end);
804 : }
805 :
806 :
807 124263 : if (_skip_renumber_nodes_and_elements)
808 : {
809 : // Loop over the nodes. Note that there may
810 : // be nullptrs in the _nodes vector from the coarsening
811 : // process. Pack the nodes in to a contiguous array
812 : // and then trim any excess.
813 :
814 56 : std::vector<Node *>::iterator in = _nodes.begin();
815 56 : std::vector<Node *>::iterator out_iter = _nodes.begin();
816 56 : const std::vector<Node *>::iterator end = _nodes.end();
817 :
818 901831 : for (; in != end; ++in)
819 901446 : if (*in != nullptr)
820 : {
821 : // This is a reference so that if we change the pointer it will change in the vector
822 200046 : Node * & nd = *in;
823 :
824 : // If this node is still connected to an elem, put it in the list
825 400092 : if (connected_nodes.count(nd))
826 : {
827 649782 : *out_iter = nd;
828 128142 : ++out_iter;
829 :
830 : // Increment the node counter
831 649782 : nd->set_id (next_free_node++);
832 : }
833 : else // This node is orphaned, delete it!
834 : {
835 251664 : this->get_boundary_info().remove (nd);
836 143808 : _constraint_rows.erase(nd);
837 :
838 : // delete the node
839 251664 : --_n_nodes;
840 431424 : delete nd;
841 251664 : nd = nullptr;
842 : }
843 : }
844 :
845 : // Erase any additional storage. Whatever was
846 385 : _nodes.erase (out_iter, end);
847 : }
848 : else // We really DO want node renumbering
849 : {
850 : // Any nodes in the vector >= _nodes[next_free_node]
851 : // are not connected to any elements and may be deleted
852 : // if desired.
853 :
854 : // Now, delete the unused nodes
855 : {
856 101710 : std::vector<Node *>::iterator nd = _nodes.begin();
857 22168 : const std::vector<Node *>::iterator end = _nodes.end();
858 :
859 22168 : std::advance (nd, next_free_node);
860 :
861 1193307 : for (auto & node : as_range(nd, end))
862 : {
863 : // Mesh modification code might have already deleted some
864 : // nodes
865 1069429 : if (node == nullptr)
866 164880 : continue;
867 :
868 : // remove any boundary information associated with
869 : // this node
870 840651 : this->get_boundary_info().remove (node);
871 159500 : _constraint_rows.erase(node);
872 :
873 : // delete the node
874 840651 : --_n_nodes;
875 1601552 : delete node;
876 840651 : node = nullptr;
877 : }
878 :
879 123878 : _nodes.erase (nd, end);
880 : }
881 : }
882 :
883 22224 : libmesh_assert_equal_to (next_free_elem, _elements.size());
884 22224 : libmesh_assert_equal_to (next_free_node, _nodes.size());
885 :
886 124263 : this->update_parallel_id_counts();
887 124263 : }
888 :
889 :
890 :
891 2080 : void ReplicatedMesh::fix_broken_node_and_element_numbering ()
892 : {
893 : // Nodes first
894 2948210 : for (auto n : index_range(_nodes))
895 2946130 : if (this->_nodes[n] != nullptr)
896 2946130 : this->_nodes[n]->set_id() = cast_int<dof_id_type>(n);
897 :
898 : // Elements next
899 3585538 : for (auto e : index_range(_elements))
900 3583458 : if (this->_elements[e] != nullptr)
901 3583458 : this->_elements[e]->set_id() = cast_int<dof_id_type>(e);
902 2080 : }
903 :
904 :
905 86474 : dof_id_type ReplicatedMesh::n_active_elem () const
906 : {
907 157732 : return static_cast<dof_id_type>(std::distance (this->active_elements_begin(),
908 244206 : this->active_elements_end()));
909 : }
910 :
911 : std::vector<dof_id_type>
912 142 : ReplicatedMesh::get_disconnected_subdomains(std::vector<subdomain_id_type> * subdomain_ids) const
913 : {
914 : // find number of disconnected subdomains
915 4 : std::vector<dof_id_type> representative_elem_ids;
916 :
917 : // use subdomain_ids as markers for all elements to indicate if the elements
918 : // have been visited. Note: here subdomain ID is unrelated with element
919 : // subdomain_id().
920 8 : std::vector<subdomain_id_type> subdomains;
921 142 : if (!subdomain_ids)
922 0 : subdomain_ids = &subdomains;
923 4 : subdomain_ids->clear();
924 142 : subdomain_ids->resize(max_elem_id() + 1, Elem::invalid_subdomain_id);
925 :
926 : // counter of disconnected subdomains
927 4 : subdomain_id_type subdomain_counter = 0;
928 :
929 : // a stack for visiting elements, make its capacity sufficiently large to avoid
930 : // memory allocation and deallocation when the vector size changes
931 8 : std::vector<const Elem *> list;
932 142 : list.reserve(n_elem());
933 :
934 : // counter of visited elements
935 4 : dof_id_type visited = 0;
936 142 : dof_id_type n_active = n_active_elem();
937 4 : do
938 : {
939 1672 : for (const auto & elem : active_element_ptr_range())
940 876 : if ((*subdomain_ids)[elem->id()] == Elem::invalid_subdomain_id)
941 : {
942 284 : list.push_back(elem);
943 284 : (*subdomain_ids)[elem->id()] = subdomain_counter;
944 284 : break;
945 268 : }
946 : // we should be able to find a seed here
947 8 : libmesh_assert(list.size() > 0);
948 :
949 284 : dof_id_type min_id = std::numeric_limits<dof_id_type>::max();
950 1988 : while (list.size() > 0)
951 : {
952 : // pop up an element
953 1704 : const Elem * elem = list.back(); list.pop_back(); ++visited;
954 :
955 3084 : min_id = std::min(elem->id(), min_id);
956 :
957 8520 : for (auto s : elem->side_index_range())
958 : {
959 6816 : const Elem * neighbor = elem->neighbor_ptr(s);
960 6816 : if (neighbor != nullptr && (*subdomain_ids)[neighbor->id()] == Elem::invalid_subdomain_id)
961 : {
962 : // neighbor must be active
963 40 : libmesh_assert(neighbor->active());
964 1420 : list.push_back(neighbor);
965 1460 : (*subdomain_ids)[neighbor->id()] = subdomain_counter;
966 : }
967 : }
968 : }
969 :
970 284 : representative_elem_ids.push_back(min_id);
971 284 : subdomain_counter++;
972 : }
973 284 : while (visited != n_active);
974 :
975 146 : return representative_elem_ids;
976 : }
977 :
978 : std::unordered_map<dof_id_type, std::vector<std::vector<Point>>>
979 142 : ReplicatedMesh::get_boundary_points() const
980 : {
981 142 : libmesh_error_msg_if(mesh_dimension() != 2,
982 : "Error: get_boundary_points only works for 2D now");
983 :
984 : // find number of disconnected subdomains
985 : // subdomains will hold the IDs of disconnected subdomains for all elements.
986 8 : std::vector<subdomain_id_type> subdomains;
987 146 : std::vector<dof_id_type> elem_ids = get_disconnected_subdomains(&subdomains);
988 :
989 4 : std::unordered_map<dof_id_type, std::vector<std::vector<Point>>> boundary_points;
990 :
991 : // get all boundary sides that are to be erased later during visiting
992 : // use a comparison functor to avoid run-time randomness due to pointers
993 : struct boundary_side_compare
994 : {
995 1584 : bool operator()(const std::pair<const Elem *, unsigned int> & lhs,
996 : const std::pair<const Elem *, unsigned int> & rhs) const
997 : {
998 42340 : if (lhs.first->id() < rhs.first->id())
999 328 : return true;
1000 29298 : else if (lhs.first->id() == rhs.first->id())
1001 : {
1002 14668 : if (lhs.second < rhs.second)
1003 112 : return true;
1004 : }
1005 1144 : return false;
1006 : }
1007 : };
1008 8 : std::set<std::pair<const Elem *, unsigned int>, boundary_side_compare> boundary_elements;
1009 3592 : for (const auto & elem : active_element_ptr_range())
1010 8568 : for (auto s : elem->side_index_range())
1011 7008 : if (elem->neighbor_ptr(s) == nullptr)
1012 3542 : boundary_elements.insert(std::pair<const Elem *, unsigned int>(elem, s));
1013 :
1014 568 : while (!boundary_elements.empty())
1015 : {
1016 : // get the first entry as the seed
1017 426 : const Elem * eseed = boundary_elements.begin()->first;
1018 426 : unsigned int sseed = boundary_elements.begin()->second;
1019 :
1020 : // get the subdomain ID that these boundary sides attached to
1021 438 : subdomain_id_type subdomain_id = subdomains[eseed->id()];
1022 :
1023 : // start visiting the mesh to find all boundary nodes with the seed
1024 24 : std::vector<Point> bpoints;
1025 426 : const Elem * elem = eseed;
1026 12 : unsigned int s = sseed;
1027 438 : std::vector<unsigned int> local_side_nodes = elem->nodes_on_side(s);
1028 : while (true)
1029 : {
1030 96 : std::pair<const Elem *, unsigned int> side(elem, s);
1031 96 : libmesh_assert(boundary_elements.count(side));
1032 96 : boundary_elements.erase(side);
1033 :
1034 : // push all nodes on the side except the node on the other end of the side (index 1)
1035 11928 : for (auto i : index_range(local_side_nodes))
1036 8520 : if (i != 1)
1037 5400 : bpoints.push_back(*static_cast<const Point *>(elem->node_ptr(local_side_nodes[i])));
1038 :
1039 : // use the last node to find next element and side
1040 3408 : const Node * node = elem->node_ptr(local_side_nodes[1]);
1041 96 : std::set<const Elem *> neighbors;
1042 3408 : elem->find_point_neighbors(*node, neighbors);
1043 :
1044 : // if only one neighbor is found (itself), this node is a cornor node on boundary
1045 3408 : if (neighbors.size() != 1)
1046 64 : neighbors.erase(elem);
1047 :
1048 : // find the connecting side
1049 96 : bool found = false;
1050 3687 : for (const auto & neighbor : neighbors)
1051 : {
1052 10019 : for (auto ss : neighbor->side_index_range())
1053 9908 : if (neighbor->neighbor_ptr(ss) == nullptr && !(elem == neighbor && s == ss))
1054 : {
1055 4960 : local_side_nodes = neighbor->nodes_on_side(ss);
1056 : // we expect the starting point of the side to be the same as the end of the previous side
1057 5100 : if (neighbor->node_ptr(local_side_nodes[0]) == node)
1058 : {
1059 3408 : elem = neighbor;
1060 96 : s = ss;
1061 96 : found = true;
1062 96 : break;
1063 : }
1064 1596 : else if (neighbor->node_ptr(local_side_nodes[1]) == node)
1065 : {
1066 0 : elem = neighbor;
1067 0 : s = ss;
1068 0 : found = true;
1069 : // flip nodes in local_side_nodes because the side is in an opposite direction
1070 0 : auto temp(local_side_nodes);
1071 0 : local_side_nodes[0] = temp[1];
1072 0 : local_side_nodes[1] = temp[0];
1073 0 : for (unsigned int i = 2; i < temp.size(); ++i)
1074 0 : local_side_nodes[temp.size() + 1 - i] = temp[i];
1075 0 : break;
1076 : }
1077 : }
1078 105 : if (found)
1079 96 : break;
1080 : }
1081 :
1082 3408 : libmesh_error_msg_if(!found, "ERROR: mesh topology error on visiting boundary sides");
1083 :
1084 : // exit if we reach the starting point
1085 3408 : if (elem == eseed && s == sseed)
1086 12 : break;
1087 84 : }
1088 426 : boundary_points[elem_ids[subdomain_id]].push_back(bpoints);
1089 : }
1090 :
1091 146 : return boundary_points;
1092 : }
1093 :
1094 : } // namespace libMesh
|
