1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute 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 // library configuration
21 #include "libmesh/libmesh_config.h"
22
23 // C++ includes
24 #include <algorithm> // for std::min
25 #include <map> // for std::multimap
26 #include <sstream> // for std::ostringstream
27 #include <unordered_map>
28
29 // Local includes
30 #include "libmesh/boundary_info.h"
31 #include "libmesh/libmesh_logging.h"
32 #include "libmesh/elem.h"
33 #include "libmesh/ghost_point_neighbors.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/mesh_communication.h"
36 #include "libmesh/mesh_tools.h"
37 #include "libmesh/parallel.h"
38 #include "libmesh/partitioner.h"
39 #include "libmesh/point_locator_base.h"
40 #include "libmesh/threads.h"
41 #include "libmesh/enum_elem_type.h"
42 #include "libmesh/enum_point_locator_type.h"
43 #include "libmesh/auto_ptr.h" // libmesh_make_unique
44
45 namespace libMesh
46 {
47
48
49
50 // ------------------------------------------------------------
51 // MeshBase class member functions
MeshBase(const Parallel::Communicator & comm_in,unsigned char d)52 MeshBase::MeshBase (const Parallel::Communicator & comm_in,
53 unsigned char d) :
54 ParallelObject (comm_in),
55 boundary_info (new BoundaryInfo(*this)),
56 _n_parts (1),
57 _default_mapping_type(LAGRANGE_MAP),
58 _default_mapping_data(0),
59 _is_prepared (false),
60 _point_locator (),
61 _count_lower_dim_elems_in_point_locator(true),
62 _partitioner (),
63 #ifdef LIBMESH_ENABLE_UNIQUE_ID
64 _next_unique_id(DofObject::invalid_unique_id),
65 #endif
66 _skip_noncritical_partitioning(false),
67 _skip_all_partitioning(libMesh::on_command_line("--skip-partitioning")),
68 _skip_renumber_nodes_and_elements(false),
69 _skip_find_neighbors(false),
70 _allow_remote_element_removal(true),
71 _spatial_dimension(d),
72 _default_ghosting(libmesh_make_unique<GhostPointNeighbors>(*this)),
73 _point_locator_close_to_point_tol(0.)
74 {
75 _elem_dims.insert(d);
76 _ghosting_functors.insert(_default_ghosting.get());
77 libmesh_assert_less_equal (LIBMESH_DIM, 3);
78 libmesh_assert_greater_equal (LIBMESH_DIM, d);
79 libmesh_assert (libMesh::initialized());
80 }
81
82
83
MeshBase(const MeshBase & other_mesh)84 MeshBase::MeshBase (const MeshBase & other_mesh) :
85 ParallelObject (other_mesh),
86 boundary_info (new BoundaryInfo(*this)),
87 _n_parts (other_mesh._n_parts),
88 _default_mapping_type(other_mesh._default_mapping_type),
89 _default_mapping_data(other_mesh._default_mapping_data),
90 _is_prepared (other_mesh._is_prepared),
91 _point_locator (),
92 _count_lower_dim_elems_in_point_locator(other_mesh._count_lower_dim_elems_in_point_locator),
93 _partitioner (),
94 #ifdef LIBMESH_ENABLE_UNIQUE_ID
95 _next_unique_id(other_mesh._next_unique_id),
96 #endif
97 _skip_noncritical_partitioning(false),
98 _skip_all_partitioning(libMesh::on_command_line("--skip-partitioning")),
99 _skip_renumber_nodes_and_elements(other_mesh._skip_renumber_nodes_and_elements),
100 _skip_find_neighbors(other_mesh._skip_find_neighbors),
101 _allow_remote_element_removal(true),
102 _elem_dims(other_mesh._elem_dims),
103 _spatial_dimension(other_mesh._spatial_dimension),
104 _default_ghosting(libmesh_make_unique<GhostPointNeighbors>(*this)),
105 _point_locator_close_to_point_tol(other_mesh._point_locator_close_to_point_tol)
106 {
107 for (const auto & gf : other_mesh._ghosting_functors )
108 {
109 std::shared_ptr<GhostingFunctor> clone_gf = gf->clone();
110 // Some subclasses of GhostingFunctor might not override the
111 // clone function yet. If this is the case, GhostingFunctor will
112 // return nullptr by default. The clone function should be overridden
113 // in all derived classes. This following code ("else") is written
114 // for API upgrade. That will allow users gradually to update their code.
115 // Once the API upgrade is done, we will come back and delete "else."
116 if (clone_gf)
117 {
118 clone_gf->set_mesh(this);
119 add_ghosting_functor(clone_gf);
120 }
121 else
122 {
123 libmesh_deprecated();
124 add_ghosting_functor(*gf);
125 }
126 }
127
128 // Make sure we don't accidentally delete the other mesh's default
129 // ghosting functor; we'll use our own if that's needed.
130 if (other_mesh._ghosting_functors.count(other_mesh._default_ghosting.get()))
131 {
132 _ghosting_functors.erase(other_mesh._default_ghosting.get());
133 _ghosting_functors.insert(_default_ghosting.get());
134 }
135
136 if (other_mesh._partitioner.get())
137 {
138 _partitioner = other_mesh._partitioner->clone();
139 }
140 }
141
~MeshBase()142 MeshBase::~MeshBase()
143 {
144 this->clear();
145
146 libmesh_exceptionless_assert (!libMesh::closed());
147 }
148
149
150
mesh_dimension()151 unsigned int MeshBase::mesh_dimension() const
152 {
153 if (!_elem_dims.empty())
154 return cast_int<unsigned int>(*_elem_dims.rbegin());
155 return 0;
156 }
157
158
159
set_elem_dimensions(const std::set<unsigned char> & elem_dims)160 void MeshBase::set_elem_dimensions(const std::set<unsigned char> & elem_dims)
161 {
162 #ifdef DEBUG
163 // In debug mode, we call cache_elem_dims() and then make sure
164 // the result actually agrees with what the user specified.
165 parallel_object_only();
166
167 this->cache_elem_dims();
168 libmesh_assert_msg(_elem_dims == elem_dims, \
169 "Specified element dimensions does not match true element dimensions!");
170 #endif
171
172 _elem_dims = elem_dims;
173 }
174
spatial_dimension()175 unsigned int MeshBase::spatial_dimension () const
176 {
177 return cast_int<unsigned int>(_spatial_dimension);
178 }
179
180
181
set_spatial_dimension(unsigned char d)182 void MeshBase::set_spatial_dimension(unsigned char d)
183 {
184 // The user can set the _spatial_dimension however they wish,
185 // libMesh will only *increase* the spatial dimension, however,
186 // never decrease it.
187 _spatial_dimension = d;
188 }
189
190
191
add_elem_integer(const std::string & name,bool allocate_data,dof_id_type default_value)192 unsigned int MeshBase::add_elem_integer(const std::string & name,
193 bool allocate_data,
194 dof_id_type default_value)
195 {
196 for (auto i : index_range(_elem_integer_names))
197 if (_elem_integer_names[i] == name)
198 {
199 libmesh_assert_less(i, _elem_integer_default_values.size());
200 _elem_integer_default_values[i] = default_value;
201 return i;
202 }
203
204 libmesh_assert_equal_to(_elem_integer_names.size(),
205 _elem_integer_default_values.size());
206 _elem_integer_names.push_back(name);
207 _elem_integer_default_values.push_back(default_value);
208 if (allocate_data)
209 this->size_elem_extra_integers();
210 return _elem_integer_names.size()-1;
211 }
212
213
214
add_elem_integers(const std::vector<std::string> & names,bool allocate_data,const std::vector<dof_id_type> * default_values)215 std::vector<unsigned int> MeshBase::add_elem_integers(const std::vector<std::string> & names,
216 bool allocate_data,
217 const std::vector<dof_id_type> * default_values)
218 {
219 libmesh_assert(!default_values || default_values->size() == names.size());
220 libmesh_assert_equal_to(_elem_integer_names.size(), _elem_integer_default_values.size());
221
222 std::unordered_map<std::string, std::size_t> name_indices;
223 for (auto i : index_range(_elem_integer_names))
224 name_indices[_elem_integer_names[i]] = i;
225
226 std::vector<unsigned int> returnval(names.size());
227
228 bool added_an_integer = false;
229 for (auto i : index_range(names))
230 {
231 const std::string & name = names[i];
232 auto it = name_indices.find(name);
233 if (it != name_indices.end())
234 {
235 returnval[i] = it->second;
236 _elem_integer_default_values[it->second] =
237 default_values ? (*default_values)[i] : DofObject::invalid_id;
238 }
239 else
240 {
241 returnval[i] = _elem_integer_names.size();
242 name_indices[name] = returnval[i];
243 _elem_integer_names.push_back(name);
244 _elem_integer_default_values.push_back
245 (default_values ? (*default_values)[i] : DofObject::invalid_id);
246 added_an_integer = true;
247 }
248 }
249
250 if (allocate_data && added_an_integer)
251 this->size_elem_extra_integers();
252
253 return returnval;
254 }
255
256
257
get_elem_integer_index(const std::string & name)258 unsigned int MeshBase::get_elem_integer_index(const std::string & name) const
259 {
260 for (auto i : index_range(_elem_integer_names))
261 if (_elem_integer_names[i] == name)
262 return i;
263
264 libmesh_error_msg("Unknown elem integer " << name);
265 return libMesh::invalid_uint;
266 }
267
268
269
has_elem_integer(const std::string & name)270 bool MeshBase::has_elem_integer(const std::string & name) const
271 {
272 for (auto & entry : _elem_integer_names)
273 if (entry == name)
274 return true;
275
276 return false;
277 }
278
279
280
add_node_integer(const std::string & name,bool allocate_data,dof_id_type default_value)281 unsigned int MeshBase::add_node_integer(const std::string & name,
282 bool allocate_data,
283 dof_id_type default_value)
284 {
285 for (auto i : index_range(_node_integer_names))
286 if (_node_integer_names[i] == name)
287 {
288 libmesh_assert_less(i, _node_integer_default_values.size());
289 _node_integer_default_values[i] = default_value;
290 return i;
291 }
292
293 libmesh_assert_equal_to(_node_integer_names.size(),
294 _node_integer_default_values.size());
295 _node_integer_names.push_back(name);
296 _node_integer_default_values.push_back(default_value);
297 if (allocate_data)
298 this->size_node_extra_integers();
299 return _node_integer_names.size()-1;
300 }
301
302
303
add_node_integers(const std::vector<std::string> & names,bool allocate_data,const std::vector<dof_id_type> * default_values)304 std::vector<unsigned int> MeshBase::add_node_integers(const std::vector<std::string> & names,
305 bool allocate_data,
306 const std::vector<dof_id_type> * default_values)
307 {
308 libmesh_assert(!default_values || default_values->size() == names.size());
309 libmesh_assert_equal_to(_node_integer_names.size(), _node_integer_default_values.size());
310
311 std::unordered_map<std::string, std::size_t> name_indices;
312 for (auto i : index_range(_node_integer_names))
313 name_indices[_node_integer_names[i]] = i;
314
315 std::vector<unsigned int> returnval(names.size());
316
317 bool added_an_integer = false;
318 for (auto i : index_range(names))
319 {
320 const std::string & name = names[i];
321 auto it = name_indices.find(name);
322 if (it != name_indices.end())
323 {
324 returnval[i] = it->second;
325 _node_integer_default_values[it->second] =
326 default_values ? (*default_values)[i] : DofObject::invalid_id;
327 }
328 else
329 {
330 returnval[i] = _node_integer_names.size();
331 name_indices[name] = returnval[i];
332 _node_integer_names.push_back(name);
333 _node_integer_default_values.push_back
334 (default_values ? (*default_values)[i] : DofObject::invalid_id);
335 added_an_integer = true;
336 }
337 }
338
339 if (allocate_data && added_an_integer)
340 this->size_node_extra_integers();
341
342 return returnval;
343 }
344
345
346
get_node_integer_index(const std::string & name)347 unsigned int MeshBase::get_node_integer_index(const std::string & name) const
348 {
349 for (auto i : index_range(_node_integer_names))
350 if (_node_integer_names[i] == name)
351 return i;
352
353 libmesh_error_msg("Unknown node integer " << name);
354 return libMesh::invalid_uint;
355 }
356
357
358
has_node_integer(const std::string & name)359 bool MeshBase::has_node_integer(const std::string & name) const
360 {
361 for (auto & entry : _node_integer_names)
362 if (entry == name)
363 return true;
364
365 return false;
366 }
367
368
369
remove_orphaned_nodes()370 void MeshBase::remove_orphaned_nodes ()
371 {
372 LOG_SCOPE("remove_orphaned_nodes()", "MeshBase");
373
374 // Will hold the set of nodes that are currently connected to elements
375 std::unordered_set<Node *> connected_nodes;
376
377 // Loop over the elements. Find which nodes are connected to at
378 // least one of them.
379 for (const auto & element : this->element_ptr_range())
380 for (auto & n : element->node_ref_range())
381 connected_nodes.insert(&n);
382
383 for (const auto & node : this->node_ptr_range())
384 if (!connected_nodes.count(node))
385 this->delete_node(node);
386 }
387
388
389
prepare_for_use(const bool skip_renumber_nodes_and_elements,const bool skip_find_neighbors)390 void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
391 {
392 libmesh_deprecated();
393
394 // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
395 // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
396 // precedence
397 if (skip_renumber_nodes_and_elements)
398 this->allow_renumbering(false);
399
400 // We always accept the user's value for skip_find_neighbors, in contrast to skip_renumber
401 const bool old_allow_find_neighbors = this->allow_find_neighbors();
402 this->allow_find_neighbors(!skip_find_neighbors);
403
404 this->prepare_for_use();
405
406 this->allow_find_neighbors(old_allow_find_neighbors);
407 }
408
prepare_for_use(const bool skip_renumber_nodes_and_elements)409 void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements)
410 {
411 libmesh_deprecated();
412
413 // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
414 // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
415 // precedence
416 if (skip_renumber_nodes_and_elements)
417 this->allow_renumbering(false);
418
419 this->prepare_for_use();
420 }
421
prepare_for_use()422 void MeshBase::prepare_for_use ()
423 {
424 LOG_SCOPE("prepare_for_use()", "MeshBase");
425
426 parallel_object_only();
427
428 libmesh_assert(this->comm().verify(this->is_serial()));
429
430 // A distributed mesh may have processors with no elements (or
431 // processors with no elements of higher dimension, if we ever
432 // support mixed-dimension meshes), but we want consistent
433 // mesh_dimension anyways.
434 //
435 // cache_elem_dims() should get the elem_dimensions() and
436 // mesh_dimension() correct later, and we don't need it earlier.
437
438
439 // Renumber the nodes and elements so that they in contiguous
440 // blocks. By default, _skip_renumber_nodes_and_elements is false.
441 //
442 // Instances where you if prepare_for_use() should not renumber the nodes
443 // and elements include reading in e.g. an xda/r or gmv file. In
444 // this case, the ordering of the nodes may depend on an accompanying
445 // solution, and the node ordering cannot be changed.
446
447
448 // Mesh modification operations might not leave us with consistent
449 // id counts, or might leave us with orphaned nodes we're no longer
450 // using, but our partitioner might need that consistency and/or
451 // might be confused by orphaned nodes.
452 if (!_skip_renumber_nodes_and_elements)
453 this->renumber_nodes_and_elements();
454 else
455 {
456 this->remove_orphaned_nodes();
457 this->update_parallel_id_counts();
458 }
459
460 // Let all the elements find their neighbors
461 if (!_skip_find_neighbors)
462 this->find_neighbors();
463
464 // The user may have set boundary conditions. We require that the
465 // boundary conditions were set consistently. Because we examine
466 // neighbors when evaluating non-raw boundary condition IDs, this
467 // assert is only valid when our neighbor links are in place.
468 #ifdef DEBUG
469 MeshTools::libmesh_assert_valid_boundary_ids(*this);
470 #endif
471
472 // Search the mesh for all the dimensions of the elements
473 // and cache them.
474 this->cache_elem_dims();
475
476 // Search the mesh for elements that have a neighboring element
477 // of dim+1 and set that element as the interior parent
478 this->detect_interior_parents();
479
480 // Fix up node unique ids in case mesh generation code didn't take
481 // exceptional care to do so.
482 // MeshCommunication().make_node_unique_ids_parallel_consistent(*this);
483
484 // We're going to still require that mesh generation code gets
485 // element unique ids consistent.
486 #if defined(DEBUG) && defined(LIBMESH_ENABLE_UNIQUE_ID)
487 MeshTools::libmesh_assert_valid_unique_ids(*this);
488 #endif
489
490 // Reset our PointLocator. Any old locator is invalidated any time
491 // the elements in the underlying elements in the mesh have changed,
492 // so we clear it here.
493 this->clear_point_locator();
494
495 // Allow our GhostingFunctor objects to reinit if necessary.
496 // Do this before partitioning and redistributing, and before
497 // deleting remote elements.
498 for (auto & gf : _ghosting_functors)
499 {
500 libmesh_assert(gf);
501 gf->mesh_reinit();
502 }
503
504 // Partition the mesh unless *all* partitioning is to be skipped.
505 // If only noncritical partitioning is to be skipped, the
506 // partition() call will still check for orphaned nodes.
507 if (!skip_partitioning())
508 this->partition();
509
510 // If we're using DistributedMesh, we'll probably want it
511 // parallelized.
512 if (this->_allow_remote_element_removal)
513 this->delete_remote_elements();
514
515 if (!_skip_renumber_nodes_and_elements)
516 this->renumber_nodes_and_elements();
517
518 // The mesh is now prepared for use.
519 _is_prepared = true;
520
521 #if defined(DEBUG) && defined(LIBMESH_ENABLE_UNIQUE_ID)
522 MeshTools::libmesh_assert_valid_boundary_ids(*this);
523 MeshTools::libmesh_assert_valid_unique_ids(*this);
524 #endif
525 }
526
527
528
clear()529 void MeshBase::clear ()
530 {
531 // Reset the number of partitions
532 _n_parts = 1;
533
534 // Reset the _is_prepared flag
535 _is_prepared = false;
536
537 // Clear boundary information
538 if (boundary_info)
539 boundary_info->clear();
540
541 // Clear element dimensions
542 _elem_dims.clear();
543
544 // Clear our point locator.
545 this->clear_point_locator();
546 }
547
548
549
remove_ghosting_functor(GhostingFunctor & ghosting_functor)550 void MeshBase::remove_ghosting_functor(GhostingFunctor & ghosting_functor)
551 {
552 _ghosting_functors.erase(&ghosting_functor);
553
554 auto it = _shared_functors.find(&ghosting_functor);
555 if (it != _shared_functors.end())
556 _shared_functors.erase(it);
557 }
558
559
560
subdomain_ids(std::set<subdomain_id_type> & ids)561 void MeshBase::subdomain_ids (std::set<subdomain_id_type> & ids) const
562 {
563 // This requires an inspection on every processor
564 parallel_object_only();
565
566 ids.clear();
567
568 for (const auto & elem : this->active_local_element_ptr_range())
569 ids.insert(elem->subdomain_id());
570
571 // Some subdomains may only live on other processors
572 this->comm().set_union(ids);
573 }
574
575
576
n_subdomains()577 subdomain_id_type MeshBase::n_subdomains() const
578 {
579 // This requires an inspection on every processor
580 parallel_object_only();
581
582 std::set<subdomain_id_type> ids;
583
584 this->subdomain_ids (ids);
585
586 return cast_int<subdomain_id_type>(ids.size());
587 }
588
589
590
591
n_nodes_on_proc(const processor_id_type proc_id)592 dof_id_type MeshBase::n_nodes_on_proc (const processor_id_type proc_id) const
593 {
594 // We're either counting a processor's nodes or unpartitioned
595 // nodes
596 libmesh_assert (proc_id < this->n_processors() ||
597 proc_id == DofObject::invalid_processor_id);
598
599 return static_cast<dof_id_type>(std::distance (this->pid_nodes_begin(proc_id),
600 this->pid_nodes_end (proc_id)));
601 }
602
603
604
n_elem_on_proc(const processor_id_type proc_id)605 dof_id_type MeshBase::n_elem_on_proc (const processor_id_type proc_id) const
606 {
607 // We're either counting a processor's elements or unpartitioned
608 // elements
609 libmesh_assert (proc_id < this->n_processors() ||
610 proc_id == DofObject::invalid_processor_id);
611
612 return static_cast<dof_id_type>(std::distance (this->pid_elements_begin(proc_id),
613 this->pid_elements_end (proc_id)));
614 }
615
616
617
n_active_elem_on_proc(const processor_id_type proc_id)618 dof_id_type MeshBase::n_active_elem_on_proc (const processor_id_type proc_id) const
619 {
620 libmesh_assert_less (proc_id, this->n_processors());
621 return static_cast<dof_id_type>(std::distance (this->active_pid_elements_begin(proc_id),
622 this->active_pid_elements_end (proc_id)));
623 }
624
625
626
n_sub_elem()627 dof_id_type MeshBase::n_sub_elem () const
628 {
629 dof_id_type ne=0;
630
631 for (const auto & elem : this->element_ptr_range())
632 ne += elem->n_sub_elem();
633
634 return ne;
635 }
636
637
638
n_active_sub_elem()639 dof_id_type MeshBase::n_active_sub_elem () const
640 {
641 dof_id_type ne=0;
642
643 for (const auto & elem : this->active_element_ptr_range())
644 ne += elem->n_sub_elem();
645
646 return ne;
647 }
648
649
650
get_info()651 std::string MeshBase::get_info() const
652 {
653 std::ostringstream oss;
654
655 oss << " Mesh Information:" << '\n';
656
657 if (!_elem_dims.empty())
658 {
659 oss << " elem_dimensions()={";
660 std::copy(_elem_dims.begin(),
661 --_elem_dims.end(), // --end() is valid if the set is non-empty
662 std::ostream_iterator<unsigned int>(oss, ", "));
663 oss << cast_int<unsigned int>(*_elem_dims.rbegin());
664 oss << "}\n";
665 }
666
667 oss << " spatial_dimension()=" << this->spatial_dimension() << '\n'
668 << " n_nodes()=" << this->n_nodes() << '\n'
669 << " n_local_nodes()=" << this->n_local_nodes() << '\n'
670 << " n_elem()=" << this->n_elem() << '\n'
671 << " n_local_elem()=" << this->n_local_elem() << '\n'
672 #ifdef LIBMESH_ENABLE_AMR
673 << " n_active_elem()=" << this->n_active_elem() << '\n'
674 #endif
675 << " n_subdomains()=" << static_cast<std::size_t>(this->n_subdomains()) << '\n'
676 << " n_partitions()=" << static_cast<std::size_t>(this->n_partitions()) << '\n'
677 << " n_processors()=" << static_cast<std::size_t>(this->n_processors()) << '\n'
678 << " n_threads()=" << static_cast<std::size_t>(libMesh::n_threads()) << '\n'
679 << " processor_id()=" << static_cast<std::size_t>(this->processor_id()) << '\n';
680
681 return oss.str();
682 }
683
684
print_info(std::ostream & os)685 void MeshBase::print_info(std::ostream & os) const
686 {
687 os << this->get_info()
688 << std::endl;
689 }
690
691
692 std::ostream & operator << (std::ostream & os, const MeshBase & m)
693 {
694 m.print_info(os);
695 return os;
696 }
697
698
partition(const unsigned int n_parts)699 void MeshBase::partition (const unsigned int n_parts)
700 {
701 // If we get here and we have unpartitioned elements, we need that
702 // fixed.
703 if (this->n_unpartitioned_elem() > 0)
704 {
705 libmesh_assert (partitioner().get());
706 libmesh_assert (this->is_serial());
707 partitioner()->partition (*this, n_parts);
708 }
709 // A nullptr partitioner or a skip_partitioning(true) call or a
710 // skip_noncritical_partitioning(true) call means don't repartition;
711 // skip_noncritical_partitioning() checks all these.
712 else if (!skip_noncritical_partitioning())
713 {
714 partitioner()->partition (*this, n_parts);
715 }
716 else
717 {
718 // Adaptive coarsening may have "orphaned" nodes on processors
719 // whose elements no longer share them. We need to check for
720 // and possibly fix that.
721 MeshTools::correct_node_proc_ids(*this);
722
723 // Make sure locally cached partition count is correct
724 this->recalculate_n_partitions();
725
726 // Make sure any other locally cached data is correct
727 this->update_post_partitioning();
728 }
729 }
730
recalculate_n_partitions()731 unsigned int MeshBase::recalculate_n_partitions()
732 {
733 // This requires an inspection on every processor
734 parallel_object_only();
735
736 unsigned int max_proc_id=0;
737
738 for (const auto & elem : this->active_local_element_ptr_range())
739 max_proc_id = std::max(max_proc_id, static_cast<unsigned int>(elem->processor_id()));
740
741 // The number of partitions is one more than the max processor ID.
742 _n_parts = max_proc_id+1;
743
744 this->comm().max(_n_parts);
745
746 return _n_parts;
747 }
748
749
750
sub_point_locator()751 std::unique_ptr<PointLocatorBase> MeshBase::sub_point_locator () const
752 {
753 // If there's no master point locator, then we need one.
754 if (_point_locator.get() == nullptr)
755 {
756 // PointLocator construction may not be safe within threads
757 libmesh_assert(!Threads::in_threads);
758
759 // And it may require parallel communication
760 parallel_object_only();
761
762 _point_locator = PointLocatorBase::build(TREE_ELEMENTS, *this);
763
764 if (_point_locator_close_to_point_tol > 0.)
765 _point_locator->set_close_to_point_tol(_point_locator_close_to_point_tol);
766 }
767
768 // Otherwise there was a master point locator, and we can grab a
769 // sub-locator easily.
770 return PointLocatorBase::build(TREE_ELEMENTS, *this, _point_locator.get());
771 }
772
773
774
clear_point_locator()775 void MeshBase::clear_point_locator ()
776 {
777 _point_locator.reset(nullptr);
778 }
779
780
781
set_count_lower_dim_elems_in_point_locator(bool count_lower_dim_elems)782 void MeshBase::set_count_lower_dim_elems_in_point_locator(bool count_lower_dim_elems)
783 {
784 _count_lower_dim_elems_in_point_locator = count_lower_dim_elems;
785 }
786
787
788
get_count_lower_dim_elems_in_point_locator()789 bool MeshBase::get_count_lower_dim_elems_in_point_locator() const
790 {
791 return _count_lower_dim_elems_in_point_locator;
792 }
793
794
795
subdomain_name(subdomain_id_type id)796 std::string & MeshBase::subdomain_name(subdomain_id_type id)
797 {
798 return _block_id_to_name[id];
799 }
800
subdomain_name(subdomain_id_type id)801 const std::string & MeshBase::subdomain_name(subdomain_id_type id) const
802 {
803 // An empty string to return when no matching subdomain name is found
804 static const std::string empty;
805
806 std::map<subdomain_id_type, std::string>::const_iterator iter = _block_id_to_name.find(id);
807 if (iter == _block_id_to_name.end())
808 return empty;
809 else
810 return iter->second;
811 }
812
813
814
815
get_id_by_name(const std::string & name)816 subdomain_id_type MeshBase::get_id_by_name(const std::string & name) const
817 {
818 // Linear search over the map values.
819 std::map<subdomain_id_type, std::string>::const_iterator
820 iter = _block_id_to_name.begin(),
821 end_iter = _block_id_to_name.end();
822
823 for ( ; iter != end_iter; ++iter)
824 if (iter->second == name)
825 return iter->first;
826
827 // If we made it here without returning, we don't have a subdomain
828 // with the requested name, so return Elem::invalid_subdomain_id.
829 return Elem::invalid_subdomain_id;
830 }
831
cache_elem_dims()832 void MeshBase::cache_elem_dims()
833 {
834 // This requires an inspection on every processor
835 parallel_object_only();
836
837 // Need to clear _elem_dims first in case all elements of a
838 // particular dimension have been deleted.
839 _elem_dims.clear();
840
841 for (const auto & elem : this->active_element_ptr_range())
842 _elem_dims.insert(cast_int<unsigned char>(elem->dim()));
843
844 // Some different dimension elements may only live on other processors
845 this->comm().set_union(_elem_dims);
846
847 // If the largest element dimension found is larger than the current
848 // _spatial_dimension, increase _spatial_dimension.
849 unsigned int max_dim = this->mesh_dimension();
850 if (max_dim > _spatial_dimension)
851 _spatial_dimension = cast_int<unsigned char>(max_dim);
852
853 // _spatial_dimension may need to increase from 1->2 or 2->3 if the
854 // mesh is full of 1D elements but they are not x-aligned, or the
855 // mesh is full of 2D elements but they are not in the x-y plane.
856 // If the mesh is x-aligned or x-y planar, we will end up checking
857 // every node's coordinates and not breaking out of the loop
858 // early...
859 #if LIBMESH_DIM > 1
860 if (_spatial_dimension < 3)
861 {
862 for (const auto & node : this->node_ptr_range())
863 {
864 // Note: the exact floating point comparison is intentional,
865 // we don't want to get tripped up by tolerances.
866 if ((*node)(1) != 0.)
867 {
868 _spatial_dimension = 2;
869 #if LIBMESH_DIM == 2
870 // If libmesh is compiled in 2D mode, this is the
871 // largest spatial dimension possible so we can break
872 // out.
873 break;
874 #endif
875 }
876
877 #if LIBMESH_DIM > 2
878 if ((*node)(2) != 0.)
879 {
880 // Spatial dimension can't get any higher than this, so
881 // we can break out.
882 _spatial_dimension = 3;
883 break;
884 }
885 #endif
886 }
887 }
888 #endif // LIBMESH_DIM > 1
889 }
890
detect_interior_parents()891 void MeshBase::detect_interior_parents()
892 {
893 // This requires an inspection on every processor
894 parallel_object_only();
895
896 // Check if the mesh contains mixed dimensions. If so, then set interior parents, otherwise return.
897 if (this->elem_dimensions().size() == 1)
898 return;
899
900 //This map will be used to set interior parents
901 std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elem;
902
903 for (const auto & elem : this->active_element_ptr_range())
904 {
905 // Populating the node_to_elem map, same as MeshTools::build_nodes_to_elem_map
906 for (auto n : make_range(elem->n_vertices()))
907 {
908 libmesh_assert_less (elem->id(), this->max_elem_id());
909
910 node_to_elem[elem->node_id(n)].push_back(elem->id());
911 }
912 }
913
914 // Automatically set interior parents
915 for (const auto & element : this->element_ptr_range())
916 {
917 // Ignore an 3D element or an element that already has an interior parent
918 if (element->dim()>=LIBMESH_DIM || element->interior_parent())
919 continue;
920
921 // Start by generating a SET of elements that are dim+1 to the current
922 // element at each vertex of the current element, thus ignoring interior nodes.
923 // If one of the SET of elements is empty, then we will not have an interior parent
924 // since an interior parent must be connected to all vertices of the current element
925 std::vector<std::set<dof_id_type>> neighbors( element->n_vertices() );
926
927 bool found_interior_parents = false;
928
929 for (auto n : make_range(element->n_vertices()))
930 {
931 std::vector<dof_id_type> & element_ids = node_to_elem[element->node_id(n)];
932 for (const auto & eid : element_ids)
933 if (this->elem_ref(eid).dim() == element->dim()+1)
934 neighbors[n].insert(eid);
935
936 if (neighbors[n].size()>0)
937 {
938 found_interior_parents = true;
939 }
940 else
941 {
942 // We have found an empty set, no reason to continue
943 // Ensure we set this flag to false before the break since it could have
944 // been set to true for previous vertex
945 found_interior_parents = false;
946 break;
947 }
948 }
949
950 // If we have successfully generated a set of elements for each vertex, we will compare
951 // the set for vertex 0 will the sets for the vertices until we find a id that exists in
952 // all sets. If found, this is our an interior parent id. The interior parent id found
953 // will be the lowest element id if there is potential for multiple interior parents.
954 if (found_interior_parents)
955 {
956 std::set<dof_id_type> & neighbors_0 = neighbors[0];
957 for (const auto & interior_parent_id : neighbors_0)
958 {
959 found_interior_parents = false;
960 for (auto n : make_range(1u, element->n_vertices()))
961 {
962 if (neighbors[n].find(interior_parent_id)!=neighbors[n].end())
963 {
964 found_interior_parents=true;
965 }
966 else
967 {
968 found_interior_parents=false;
969 break;
970 }
971 }
972 if (found_interior_parents)
973 {
974 element->set_interior_parent(this->elem_ptr(interior_parent_id));
975 break;
976 }
977 }
978 }
979 }
980 }
981
982
983
set_point_locator_close_to_point_tol(Real val)984 void MeshBase::set_point_locator_close_to_point_tol(Real val)
985 {
986 _point_locator_close_to_point_tol = val;
987 if (_point_locator)
988 {
989 if (val > 0.)
990 _point_locator->set_close_to_point_tol(val);
991 else
992 _point_locator->unset_close_to_point_tol();
993 }
994 }
995
996
997
get_point_locator_close_to_point_tol()998 Real MeshBase::get_point_locator_close_to_point_tol() const
999 {
1000 return _point_locator_close_to_point_tol;
1001 }
1002
1003
1004
size_elem_extra_integers()1005 void MeshBase::size_elem_extra_integers()
1006 {
1007 const std::size_t new_size = _elem_integer_names.size();
1008 for (auto elem : this->element_ptr_range())
1009 elem->add_extra_integers(new_size, _elem_integer_default_values);
1010 }
1011
1012
1013
size_node_extra_integers()1014 void MeshBase::size_node_extra_integers()
1015 {
1016 const std::size_t new_size = _node_integer_names.size();
1017 for (auto node : this->node_ptr_range())
1018 node->add_extra_integers(new_size, _node_integer_default_values);
1019 }
1020
1021
1022 std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
merge_extra_integer_names(const MeshBase & other)1023 MeshBase::merge_extra_integer_names(const MeshBase & other)
1024 {
1025 std::pair<std::vector<unsigned int>, std::vector<unsigned int>> returnval;
1026 returnval.first = this->add_elem_integers(other._elem_integer_names, true, &other._elem_integer_default_values);
1027 returnval.second = this->add_node_integers(other._node_integer_names, true, &other._node_integer_default_values);
1028 return returnval;
1029 }
1030
1031
1032 } // namespace libMesh
1033