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 : #ifndef LIBMESH_MESH_TOOLS_H
21 : #define LIBMESH_MESH_TOOLS_H
22 :
23 : // Local Includes
24 : #include "libmesh/libmesh.h"
25 : #include "libmesh/bounding_box.h"
26 : #include "libmesh/id_types.h"
27 : #include "libmesh/mesh_base.h"
28 :
29 : // C++ Includes
30 : #include <limits>
31 : #include <set>
32 : #include <unordered_set>
33 : #include <unordered_map>
34 : #include <vector>
35 :
36 : namespace libMesh
37 : {
38 :
39 : // forward declarations
40 : class Sphere;
41 : class Elem;
42 : enum ElemType : int;
43 :
44 : /**
45 : * Utility functions for operations on a \p Mesh object. Here is where
46 : * useful functions for interfacing with a \p Mesh should be defined.
47 : * In general this namespace should be used to prevent the \p Mesh class
48 : * from becoming too cluttered.
49 : *
50 : * \author Benjamin S. Kirk
51 : * \date 2004
52 : */
53 : namespace MeshTools
54 : {
55 :
56 : /**
57 : * \returns The sum over all the elements of the number
58 : * of nodes per element.
59 : *
60 : * This can be useful for partitioning hybrid meshes. A feasible load
61 : * balancing scheme is to keep the weight per processor as uniform as
62 : * possible.
63 : */
64 : dof_id_type total_weight (const MeshBase & mesh);
65 :
66 : /**
67 : * \returns The sum over all the elements on processor \p pid
68 : * of nodes per element.
69 : *
70 : * This can be useful for partitioning hybrid meshes. A feasible load
71 : * balancing scheme is to keep the weight per processor as uniform as
72 : * possible.
73 : */
74 : dof_id_type weight (const MeshBase & mesh,
75 : const processor_id_type pid);
76 :
77 : inline
78 : dof_id_type weight (const MeshBase & mesh)
79 : { return MeshTools::weight (mesh, mesh.processor_id()); }
80 :
81 : /**
82 : * After calling this function the input vector \p nodes_to_elem_map
83 : * will contain the node to element connectivity. That is to say
84 : * \p nodes_to_elem_map[i][j] is the global number of \f$ j^{th} \f$
85 : * element connected to node \p i.
86 : */
87 : void build_nodes_to_elem_map (const MeshBase & mesh,
88 : std::vector<std::vector<dof_id_type>> & nodes_to_elem_map);
89 :
90 : /**
91 : * The same, except element pointers are returned instead of indices.
92 : */
93 : void build_nodes_to_elem_map (const MeshBase & mesh,
94 : std::vector<std::vector<const Elem *>> & nodes_to_elem_map);
95 :
96 : /**
97 : * After calling this function the input map \p nodes_to_elem_map
98 : * will contain the node to element connectivity. That is to say
99 : * \p nodes_to_elem_map[i][j] is the global number of \f$ j^{th} \f$
100 : * element connected to node \p i.
101 : */
102 : void build_nodes_to_elem_map (const MeshBase & mesh,
103 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & nodes_to_elem_map);
104 :
105 : /**
106 : * The same, except element pointers are returned instead of indices.
107 : */
108 : void build_nodes_to_elem_map (const MeshBase & mesh,
109 : std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map);
110 :
111 :
112 : // /**
113 : // * Calling this function on a 2D mesh will convert all the elements
114 : // * to triangles. \p QUAD4s will be converted to \p TRI3s, \p QUAD8s
115 : // * and \p QUAD9s will be converted to \p TRI6s.
116 : // */
117 : // void all_tri (MeshBase & mesh);
118 :
119 : /**
120 : * Returns a std::set containing Node IDs for all of the boundary nodes
121 : */
122 : std::unordered_set<dof_id_type> find_boundary_nodes(const MeshBase & mesh);
123 :
124 : /**
125 : * Returns a std::set containing Node IDs for all of the block boundary nodes
126 : *
127 : * A "block boundary node" is a node that is connected to elements from 2
128 : * or more blocks
129 : */
130 : std::unordered_set<dof_id_type> find_block_boundary_nodes(const MeshBase & mesh);
131 :
132 : /**
133 : * \returns A BoundingBox that bounds the mesh.
134 : */
135 : libMesh::BoundingBox
136 : create_bounding_box (const MeshBase & mesh);
137 :
138 : /**
139 : * \returns A bounding sphere for \p mesh instead of a bounding box.
140 : */
141 : Sphere
142 : bounding_sphere (const MeshBase & mesh);
143 :
144 : /**
145 : * \returns Two points defining a cartesian box that bounds the
146 : * nodes of the mesh.
147 : *
148 : * In the case of curved elements, this box might *not* bound the
149 : * elements of the mesh.
150 : */
151 : libMesh::BoundingBox
152 : create_nodal_bounding_box (const MeshBase & mesh);
153 :
154 : /**
155 : * \returns Two points defining a cartesian box that bounds the
156 : * elements belonging to the local processor.
157 : *
158 : * Unlike the other bounding box creation functions, this does *not*
159 : * need to be run in parallel, because this is the only function we
160 : * can guarantee can be resolved with only local information.
161 : */
162 : libMesh::BoundingBox
163 : create_local_bounding_box (const MeshBase & mesh);
164 :
165 : /**
166 : * \returns A BoundingBox that bounds the elements belonging to
167 : * processor pid.
168 : */
169 : libMesh::BoundingBox
170 : create_processor_bounding_box (const MeshBase & mesh,
171 : const processor_id_type pid);
172 :
173 : /**
174 : * \returns A processor bounding sphere instead of a processor bounding box.
175 : */
176 : Sphere
177 : processor_bounding_sphere (const MeshBase & mesh,
178 : const processor_id_type pid);
179 :
180 : /**
181 : * \returns A BoundingBox that bounds the elements belonging to
182 : * subdomain sid.
183 : */
184 : libMesh::BoundingBox
185 : create_subdomain_bounding_box (const MeshBase & mesh,
186 : const subdomain_id_type sid);
187 :
188 : /**
189 : * \returns A subdomain bounding sphere instead of a subdomain bounding box.
190 : */
191 : Sphere
192 : subdomain_bounding_sphere (const MeshBase & mesh,
193 : const subdomain_id_type sid);
194 :
195 :
196 : /**
197 : * Fills in a vector of all element types in the mesh. Implemented
198 : * in terms of element_iterators.
199 : */
200 : void elem_types (const MeshBase & mesh,
201 : std::vector<ElemType> & et);
202 :
203 : /**
204 : * \returns The number of elements of type \p type.
205 : *
206 : * Implemented in terms of type_element_iterators.
207 : */
208 : dof_id_type n_elem_of_type (const MeshBase & mesh,
209 : const ElemType type);
210 :
211 : /**
212 : * \returns The number of active elements of type \p type.
213 : *
214 : * Implemented in terms of active_type_element_iterators.
215 : */
216 : dof_id_type n_active_elem_of_type (const MeshBase & mesh,
217 : const ElemType type);
218 :
219 : /**
220 : * \returns The number of elements of type \p type at the specified
221 : * refinement level.
222 : *
223 : * \todo Replace all of the n_xxx_elem() functions like this with
224 : * a single function which takes a range of iterators and computes the
225 : * std::distance between them.
226 : */
227 : dof_id_type n_non_subactive_elem_of_type_at_level(const MeshBase & mesh,
228 : const ElemType type,
229 : const unsigned int level);
230 :
231 : /**
232 : * \returns The number of levels of refinement in the mesh.
233 : *
234 : * Implemented by looping over all the local elements and
235 : * unpartitioned elements and finding the maximum level, then summing
236 : * in parallel.
237 : */
238 : unsigned int n_levels(const MeshBase & mesh);
239 :
240 : /**
241 : * \returns The number of levels of refinement in the local mesh.
242 : *
243 : * Implemented by looping over all the local elements and finding the
244 : * maximum level.
245 : */
246 : unsigned int n_local_levels(const MeshBase & mesh);
247 :
248 : /**
249 : * \returns The number of levels of refinement in the active mesh.
250 : *
251 : * Implemented by looping over all the active local elements and finding
252 : * the maximum level, then taking the max in parallel.
253 : */
254 : unsigned int n_active_levels(const MeshBase & mesh);
255 :
256 : /**
257 : * \returns The number of levels of refinement in the active local mesh.
258 : *
259 : * Implemented by looping over all the active local elements and finding
260 : * the maximum level.
261 : */
262 : unsigned int n_active_local_levels(const MeshBase & mesh);
263 :
264 : /**
265 : * \returns The number of p-levels of refinement in the mesh.
266 : *
267 : * Implemented by looping over all the local elements and finding the
268 : * maximum p-level, then summing in parallel.
269 : */
270 : unsigned int n_p_levels (const MeshBase & mesh);
271 :
272 : /**
273 : * \returns The number of levels of refinement in the mesh, even if that
274 : * mesh is not currently properly distributed or properly serialized.
275 : *
276 : * Implemented by looping over all elements and finding the maximum
277 : * level, then summing in parallel. This is much slower than
278 : * n_levels() but will return correct values even when the mesh is in
279 : * an inconsistent parallel state.
280 : */
281 : unsigned int paranoid_n_levels(const MeshBase & mesh);
282 :
283 : /**
284 : * \returns The number of topologically connected components in the
285 : * mesh, where a local connection is defined as two nodes connected to
286 : * the same element, two elements connected to the same node, and/or
287 : * two nodes connected by a constraint row; thus a connection can be
288 : * made at a single vertex, not just at an element side.
289 : *
290 : * To count sufficiently weak constraint row coefficients as a
291 : * lack-of-connection, set a constraint_tol greater than 0.
292 : */
293 : dof_id_type n_connected_components(const MeshBase & mesh,
294 : Real constraint_tol = 0);
295 :
296 : /**
297 : * Builds a set of node IDs for nodes which belong to non-subactive
298 : * elements. Non-subactive elements are those which are either active
299 : * or inactive. This is useful for determining which nodes should be
300 : * written to a data file, and is used by the XDA mesh writing methods.
301 : */
302 : void get_not_subactive_node_ids(const MeshBase & mesh,
303 : std::set<dof_id_type> & not_subactive_node_ids);
304 :
305 : /**
306 : * Count up the number of elements of a specific type
307 : * (as defined by an iterator range).
308 : */
309 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
310 : const MeshBase::const_element_iterator & end);
311 :
312 :
313 : /**
314 : * Count up the number of nodes of a specific type
315 : * (as defined by an iterator range).
316 : */
317 : dof_id_type n_nodes (const MeshBase::const_node_iterator & begin,
318 : const MeshBase::const_node_iterator & end);
319 :
320 :
321 : /**
322 : * Find the total volume of a mesh (interpreting that as area for \p
323 : * dim = 2, or total arc length for \p dim = 1, or number of NodeElem
324 : * in the mesh for \p dim = 0). By default the dimension chosen will
325 : * be the mesh_dimension(), which if not manually overridden will be
326 : * the largest dimension of elements in the mesh.
327 : */
328 : Real volume (const MeshBase & mesh,
329 : unsigned int dim=libMesh::invalid_uint);
330 :
331 :
332 : /**
333 : * Given a mesh and a node in the mesh, the vector will be filled with
334 : * every node directly attached to the given one.
335 : */
336 : void find_nodal_neighbors(const MeshBase & mesh,
337 : const Node & n,
338 : const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
339 : std::vector<const Node *> & neighbors);
340 :
341 : /**
342 : * Given a mesh and a node in the mesh, the vector will be filled with
343 : * every node directly attached to the given one.
344 : */
345 : void find_nodal_neighbors(const MeshBase & mesh,
346 : const Node & n,
347 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
348 : std::vector<const Node *> & neighbors);
349 :
350 : /**
351 : * Given a mesh hanging_nodes will be filled with an associative array keyed off the
352 : * global id of all the hanging nodes in the mesh. It will hold an array of the
353 : * parents of the node (meaning the two nodes to either side of it that make up
354 : * the side the hanging node is on.
355 : */
356 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
357 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes);
358 :
359 : /**
360 : * Changes the processor ids on each node so be the same as the id of the
361 : * lowest element touching that node.
362 : *
363 : * This corrects "orphaned" processor ids that may occur from element
364 : * coarsening.
365 : *
366 : * On a distributed mesh, this function must be called in parallel
367 : * to sync everyone's corrected processor ids on ghost nodes.
368 : */
369 : void correct_node_proc_ids(MeshBase &);
370 :
371 : /**
372 : * Remove spline node (for IsoGeometric Analysis meshes) elements
373 : * and nodes and constraints from the mesh. This should be done
374 : * after the mesh is read but before the constraints are used by
375 : * System initialization. The result is a mesh and solution space
376 : * with lower required continuity and more unconstrained degrees of
377 : * freedom, but with fewer total degrees of freedom and far fewer
378 : * constraint equations.
379 : */
380 : void clear_spline_nodes(MeshBase &);
381 :
382 : /**
383 : * A function for testing whether a mesh's cached is_prepared() setting
384 : * is not a false positive. If the mesh is marked as not prepared, or
385 : * if preparing the already-partitioned mesh (without any
386 : * repartitioning or renumbering) does not change it, then we return
387 : * true. If the mesh believes it is prepared but prepare_for_use()
388 : * would change it, we return false.
389 : */
390 : bool valid_is_prepared (const MeshBase & mesh);
391 :
392 : ///@{
393 :
394 : /**
395 : * The following functions, only available in builds with NDEBUG
396 : * undefined, are for asserting internal consistency that we hope
397 : * should never be broken in opt
398 : */
399 :
400 : #ifndef NDEBUG
401 :
402 : /**
403 : * A function for testing that all DofObjects within a mesh
404 : * have the same n_systems count
405 : */
406 : void libmesh_assert_equal_n_systems (const MeshBase & mesh);
407 :
408 : /**
409 : * A function for testing that all non-recently-created DofObjects
410 : * within a mesh have old_dof_object data. This is not expected to
411 : * be true at all points within a simulation code.
412 : */
413 : void libmesh_assert_old_dof_objects (const MeshBase & mesh);
414 :
415 : /**
416 : * A function for walking across the mesh to try and ferret out
417 : * invalidated or misassigned pointers
418 : */
419 : void libmesh_assert_valid_node_pointers (const MeshBase & mesh);
420 :
421 : /**
422 : * A function for verifying that active local elements' neighbors
423 : * are never remote elements
424 : */
425 : void libmesh_assert_valid_remote_elems (const MeshBase & mesh);
426 :
427 : /**
428 : * A function for verifying that ids and processor assignment of elements
429 : * are correctly sorted (monotone increasing)
430 : */
431 : void libmesh_assert_valid_elem_ids (const MeshBase & mesh);
432 :
433 : /**
434 : * A function for verifying that ids of elements are correctly
435 : * sorted for AMR (parents have lower ids than children)
436 : */
437 : void libmesh_assert_valid_amr_elem_ids (const MeshBase & mesh);
438 :
439 : /**
440 : * A function for verifying that any interior_parent pointers on
441 : * elements are consistent with AMR (parents' interior_parents are
442 : * interior_parents' parents)
443 : */
444 : void libmesh_assert_valid_amr_interior_parents (const MeshBase & mesh);
445 :
446 : /**
447 : * A function for verifying that degree of freedom indexes are
448 : * contiguous on each processor, as is required by libMesh numeric
449 : * classes.
450 : *
451 : * Verify a particular system by specifying that system's number.
452 : */
453 : void libmesh_assert_contiguous_dof_ids (const MeshBase & mesh,
454 : unsigned int sysnum);
455 :
456 : /**
457 : * A function for verifying that processor assignment is
458 : * topologically consistent on nodes (each node part of an active
459 : * element on its processor) or elements (each parent has the
460 : * processor id of one of its children).
461 : */
462 : template <typename DofObjectSubclass>
463 : void libmesh_assert_topology_consistent_procids (const MeshBase & mesh);
464 :
465 : /**
466 : * A function for verifying that processor assignment of nodes matches
467 : * the heuristic specified in Node::choose_processor_id()
468 : */
469 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh);
470 :
471 : /**
472 : * A function for verifying that elements on this processor have
473 : * valid descendants and consistent active flags.
474 : */
475 : void libmesh_assert_valid_refinement_tree (const MeshBase & mesh);
476 :
477 : #endif // !NDEBUG
478 :
479 : ///@}
480 :
481 : ///@{
482 :
483 : /**
484 : * The following functions, only available in builds with DEBUG
485 : * defined, typically have surprisingly slow asymptotic behavior, and
486 : * so are unsuitable for use even with METHOD=devel
487 : */
488 :
489 : #ifdef DEBUG
490 :
491 : /**
492 : * A function for verifying that an element has been cut off
493 : * from the rest of the mesh
494 : */
495 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
496 : const Elem * bad_elem);
497 :
498 : /**
499 : * A function for testing that node locations match across processors.
500 : */
501 : void libmesh_assert_equal_points (const MeshBase & mesh);
502 :
503 : /**
504 : * A function for testing that element nodal connectivities match
505 : * across processors.
506 : */
507 : void libmesh_assert_equal_connectivity (const MeshBase & mesh);
508 :
509 : /**
510 : * A function for verifying that all nodes are connected to at least
511 : * one element.
512 : *
513 : * This will fail in the most general case. When DistributedMesh and
514 : * NodeConstraints are enabled, we expect the possibility that a
515 : * processor will be given remote nodes to satisfy node constraints
516 : * without also being given the remote elements connected to those
517 : * nodes.
518 : */
519 : void libmesh_assert_connected_nodes (const MeshBase & mesh);
520 :
521 : /**
522 : * A function for verifying that all mesh constraint rows express
523 : * relations between nodes and elements that are semilocal (local or
524 : * ghosted) to the current processor's portion of the mesh.
525 : */
526 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh);
527 :
528 : /**
529 : * A function for verifying that boundary condition ids match
530 : * across processors.
531 : */
532 : void libmesh_assert_valid_boundary_ids (const MeshBase & mesh);
533 :
534 : /**
535 : * A function for verifying that degree of freedom indexing matches
536 : * across processors.
537 : *
538 : * Verify a particular system by specifying that system's number, or
539 : * verify all systems at once by leaving \p sysnum unspecified.
540 : */
541 : void libmesh_assert_valid_dof_ids (const MeshBase & mesh,
542 : unsigned int sysnum = libMesh::invalid_uint);
543 :
544 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
545 : /**
546 : * A function for verifying that unique ids match across processors.
547 : *
548 : * FIXME: we ought to check for uniqueness too.
549 : */
550 : void libmesh_assert_valid_unique_ids (const MeshBase & mesh);
551 : #endif
552 :
553 : /**
554 : * A function for verifying that distribution of dof objects is
555 : * parallel consistent (every processor can see every node or element
556 : * it owns)
557 : */
558 : void libmesh_assert_consistent_distributed(const MeshBase & mesh);
559 :
560 : /**
561 : * A function for verifying that distribution of nodes is parallel
562 : * consistent (every processor can see every node it owns) even before
563 : * node ids have been made consistent
564 : */
565 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh);
566 :
567 : /**
568 : * A function for verifying that processor assignment is parallel
569 : * consistent (every processor agrees on the processor id of each node
570 : * it can see) even on nodes which have not yet received consistent
571 : * DofObject::id(), using element topology to identify matching nodes.
572 : */
573 : void libmesh_assert_parallel_consistent_new_node_procids (const MeshBase & mesh);
574 :
575 : /**
576 : * A function for verifying that processor assignment is parallel
577 : * consistent (every processor agrees on the processor id of each dof
578 : * object it can see)
579 : */
580 : template <typename DofObjectSubclass>
581 : void libmesh_assert_parallel_consistent_procids (const MeshBase & mesh);
582 :
583 : /**
584 : * A function for verifying that processor assignment is
585 : * both parallel and topologically consistent.
586 : */
587 : template <typename DofObjectSubclass>
588 44158 : void libmesh_assert_valid_procids (const MeshBase & mesh)
589 : {
590 44158 : libmesh_assert_parallel_consistent_procids<DofObjectSubclass>(mesh);
591 44158 : libmesh_assert_topology_consistent_procids<DofObjectSubclass>(mesh);
592 44158 : }
593 :
594 : /**
595 : * A function for verifying that refinement flags on elements
596 : * are consistent between processors
597 : */
598 : void libmesh_assert_valid_refinement_flags (const MeshBase & mesh);
599 :
600 : /**
601 : * A function for verifying that neighbor connectivity is correct (each
602 : * element is a neighbor of or descendant of a neighbor of its neighbors)
603 : * and consistent (each neighbor link goes to either the same neighbor
604 : * or to a RemoteElem on each processor)
605 : *
606 : * If assert_valid_remote_elems is set to false, then no error will be
607 : * thrown for neighbor links where a remote_elem should exist but a nullptr
608 : * exists instead.
609 : */
610 : void libmesh_assert_valid_neighbors (const MeshBase & mesh,
611 : bool assert_valid_remote_elems=true);
612 :
613 : #endif // DEBUG
614 :
615 : ///@}
616 :
617 : // There is no reason for users to call functions in the MeshTools::Private namespace.
618 : namespace Private {
619 : /**
620 : * There is no reason for a user to ever call this function.
621 : *
622 : * This function determines partition-agnostic global indices for all
623 : * nodes and elements in the mesh.
624 : *
625 : * \note After this function is called, the mesh will likely be in an
626 : * inconsistent state, i.e. \p mesh.nodes(i)->id() != i in the nodes
627 : * container. Direct node/element access via the \p mesh.node(n) or
628 : * \p mesh.elem(e) functions will likely fail. The original numbering
629 : * can (and should) be restored with a subsequent call to \p
630 : * fix_node_and_element_numbering().
631 : *
632 : */
633 : void globally_renumber_nodes_and_elements (MeshBase &);
634 : } // end namespace Private
635 :
636 : } // end namespace MeshTools
637 :
638 : } // namespace libMesh
639 :
640 :
641 : #endif // LIBMESH_MESH_TOOLS_H
|
