Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "MooseError.h"
11 : #include "MooseMesh.h"
12 : #include "Factory.h"
13 : #include "CacheChangedListsThread.h"
14 : #include "MooseUtils.h"
15 : #include "MooseApp.h"
16 : #include "RelationshipManager.h"
17 : #include "PointListAdaptor.h"
18 : #include "Executioner.h"
19 : #include "NonlinearSystemBase.h"
20 : #include "LinearSystem.h"
21 : #include "AuxiliarySystem.h"
22 : #include "Assembly.h"
23 : #include "SubProblem.h"
24 : #include "MooseVariableBase.h"
25 : #include "MooseMeshUtils.h"
26 : #include "MooseAppCoordTransform.h"
27 : #include "FEProblemBase.h"
28 :
29 : #include <utility>
30 :
31 : // libMesh
32 : #include "libmesh/bounding_box.h"
33 : #include "libmesh/boundary_info.h"
34 : #include "libmesh/mesh_tools.h"
35 : #include "libmesh/parallel.h"
36 : #include "libmesh/mesh_communication.h"
37 : #include "libmesh/periodic_boundary_base.h"
38 : #include "libmesh/fe_base.h"
39 : #include "libmesh/fe_interface.h"
40 : #include "libmesh/mesh_communication.h"
41 : #include "libmesh/mesh_tools.h"
42 : #include "libmesh/parallel.h"
43 : #include "libmesh/parallel_elem.h"
44 : #include "libmesh/parallel_node.h"
45 : #include "libmesh/parallel_ghost_sync.h"
46 : #include "libmesh/utility.h"
47 : #include "libmesh/remote_elem.h"
48 : #include "libmesh/linear_partitioner.h"
49 : #include "libmesh/centroid_partitioner.h"
50 : #include "libmesh/parmetis_partitioner.h"
51 : #include "libmesh/hilbert_sfc_partitioner.h"
52 : #include "libmesh/morton_sfc_partitioner.h"
53 : #include "libmesh/edge_edge2.h"
54 : #include "libmesh/checkpoint_io.h"
55 : #include "libmesh/mesh_refinement.h"
56 : #include "libmesh/quadrature.h"
57 : #include "libmesh/boundary_info.h"
58 : #include "libmesh/periodic_boundaries.h"
59 : #include "libmesh/quadrature_gauss.h"
60 : #include "libmesh/point_locator_base.h"
61 : #include "libmesh/default_coupling.h"
62 : #include "libmesh/ghost_point_neighbors.h"
63 : #include "libmesh/fe_type.h"
64 : #include "libmesh/enum_to_string.h"
65 : #include "libmesh/elem_side_builder.h"
66 :
67 : using namespace libMesh;
68 :
69 : // Make newer nanoflann API compatible with older nanoflann versions
70 : #if NANOFLANN_VERSION < 0x150
71 : namespace nanoflann
72 : {
73 : typedef SearchParams SearchParameters;
74 :
75 : template <typename T, typename U>
76 : using ResultItem = std::pair<T, U>;
77 : }
78 : #endif
79 :
80 : const std::array<bool, 3> MooseMesh::periodic_dim_default{false, false, false};
81 :
82 : InputParameters
83 204250 : MooseMesh::validParams()
84 : {
85 204250 : InputParameters params = MooseObject::validParams();
86 :
87 817000 : MooseEnum parallel_type("DEFAULT REPLICATED DISTRIBUTED", "DEFAULT");
88 1021250 : params.addParam<MooseEnum>("parallel_type",
89 : parallel_type,
90 : "DEFAULT: Use libMesh::ReplicatedMesh unless --distributed-mesh is "
91 : "specified on the command line "
92 : "REPLICATED: Always use libMesh::ReplicatedMesh "
93 : "DISTRIBUTED: Always use libMesh::DistributedMesh");
94 :
95 612750 : params.addParam<bool>(
96 : "allow_renumbering",
97 408500 : true,
98 : "If allow_renumbering=false, node and element numbers are kept fixed until deletion");
99 :
100 612750 : params.addParam<MooseEnum>(
101 : "partitioner",
102 408500 : partitioning(),
103 : "Specifies a mesh partitioner to use when splitting the mesh for a parallel computation.");
104 817000 : MooseEnum direction("x y z radial");
105 817000 : params.addParam<MooseEnum>("centroid_partitioner_direction",
106 : direction,
107 : "Specifies the sort direction if using the centroid partitioner. "
108 : "Available options: x, y, z, radial");
109 :
110 817000 : MooseEnum patch_update_strategy("never always auto iteration", "never");
111 817000 : params.addParam<MooseEnum>(
112 : "patch_update_strategy",
113 : patch_update_strategy,
114 : "How often to update the geometric search 'patch'. The default is to "
115 : "never update it (which is the most efficient but could be a problem "
116 : "with lots of relative motion). 'always' will update the patch for all "
117 : "secondary nodes at the beginning of every timestep which might be time "
118 : "consuming. 'auto' will attempt to determine at the start of which "
119 : "timesteps the patch for all secondary nodes needs to be updated automatically."
120 : "'iteration' updates the patch at every nonlinear iteration for a "
121 : "subset of secondary nodes for which penetration is not detected. If there "
122 : "can be substantial relative motion between the primary and secondary surfaces "
123 : "during the nonlinear iterations within a timestep, it is advisable to use "
124 : "'iteration' option to ensure accurate contact detection.");
125 :
126 : // Note: This parameter is named to match 'construct_side_list_from_node_list' in SetupMeshAction
127 612750 : params.addParam<bool>(
128 : "construct_node_list_from_side_list",
129 408500 : true,
130 : "Whether or not to generate nodesets from the sidesets (currently often required).");
131 612750 : params.addParam<bool>(
132 : "displace_node_list_by_side_list",
133 408500 : true,
134 : "Whether to renumber existing nodesets with ids matching sidesets that "
135 : "lack names matching sidesets, when constructing nodesets from sidesets via the default "
136 : "'construct_node_list_from_side_list' option, rather than to merge them with the sideset.");
137 612750 : params.addParam<unsigned int>(
138 408500 : "patch_size", 40, "The number of nodes to consider in the NearestNode neighborhood.");
139 817000 : params.addParam<unsigned int>("ghosting_patch_size",
140 : "The number of nearest neighbors considered "
141 : "for ghosting purposes when 'iteration' "
142 : "patch update strategy is used. Default is "
143 : "5 * patch_size.");
144 612750 : params.addParam<unsigned int>("max_leaf_size",
145 408500 : 10,
146 : "The maximum number of points in each leaf of the KDTree used in "
147 : "the nearest neighbor search. As the leaf size becomes larger,"
148 : "KDTree construction becomes faster but the nearest neighbor search"
149 : "becomes slower.");
150 :
151 612750 : params.addParam<bool>("build_all_side_lowerd_mesh",
152 408500 : false,
153 : "True to build the lower-dimensional mesh for all sides.");
154 :
155 612750 : params.addParam<bool>("skip_refine_when_use_split",
156 408500 : true,
157 : "True to skip uniform refinements when using a pre-split mesh.");
158 :
159 817000 : params.addParam<std::vector<SubdomainID>>(
160 : "add_subdomain_ids",
161 : "The listed subdomain ids will be assumed valid for the mesh. This permits setting up "
162 : "subdomain restrictions for subdomains initially containing no elements, which can occur, "
163 : "for example, in additive manufacturing simulations which dynamically add and remove "
164 : "elements. Names for this subdomains may be provided using add_subdomain_names. In this case "
165 : "this list and add_subdomain_names must contain the same number of items.");
166 817000 : params.addParam<std::vector<SubdomainName>>(
167 : "add_subdomain_names",
168 : "The listed subdomain names will be assumed valid for the mesh. This permits setting up "
169 : "subdomain restrictions for subdomains initially containing no elements, which can occur, "
170 : "for example, in additive manufacturing simulations which dynamically add and remove "
171 : "elements. IDs for this subdomains may be provided using add_subdomain_ids. Otherwise IDs "
172 : "are automatically assigned. In case add_subdomain_ids is set too, both lists must contain "
173 : "the same number of items.");
174 :
175 817000 : params.addParam<std::vector<BoundaryID>>(
176 : "add_sideset_ids",
177 : "The listed sideset ids will be assumed valid for the mesh. This permits setting up boundary "
178 : "restrictions for sidesets initially containing no sides. Names for this sidesets may be "
179 : "provided using add_sideset_names. In this case this list and add_sideset_names must contain "
180 : "the same number of items.");
181 817000 : params.addParam<std::vector<BoundaryName>>(
182 : "add_sideset_names",
183 : "The listed sideset names will be assumed valid for the mesh. This permits setting up "
184 : "boundary restrictions for sidesets initially containing no sides. Ids for this sidesets may "
185 : "be provided using add_sideset_ids. In this case this list and add_sideset_ids must contain "
186 : "the same number of items.");
187 :
188 817000 : params.addParam<std::vector<BoundaryID>>(
189 : "add_nodeset_ids",
190 : "The listed nodeset ids will be assumed valid for the mesh. This permits setting up boundary "
191 : "restrictions for node initially containing no sides. Names for this nodesets may be "
192 : "provided using add_nodeset_names. In this case this list and add_nodeset_names must contain "
193 : "the same number of items.");
194 612750 : params.addParam<std::vector<BoundaryName>>(
195 : "add_nodeset_names",
196 : "The listed nodeset names will be assumed valid for the mesh. This permits setting up "
197 : "boundary restrictions for nodesets initially containing no sides. Ids for this nodesets may "
198 : "be provided using add_nodesets_ids. In this case this list and add_nodesets_ids must "
199 : "contain the same number of items.");
200 :
201 204250 : params += MooseAppCoordTransform::validParams();
202 :
203 : // This indicates that the derived mesh type accepts a MeshGenerator, and should be set to true in
204 : // derived types that do so.
205 408500 : params.addPrivateParam<bool>("_mesh_generator_mesh", false);
206 :
207 : // Whether or not the mesh is pre split
208 612750 : params.addPrivateParam<bool>("_is_split", false);
209 :
210 408500 : params.registerBase("MooseMesh");
211 :
212 : // groups
213 817000 : params.addParamNamesToGroup("patch_update_strategy patch_size max_leaf_size", "Geometric search");
214 817000 : params.addParamNamesToGroup("add_subdomain_ids add_subdomain_names add_sideset_ids "
215 : "add_sideset_names add_nodeset_ids add_nodeset_names",
216 : "Pre-declaration of future mesh sub-entities");
217 817000 : params.addParamNamesToGroup("construct_node_list_from_side_list build_all_side_lowerd_mesh "
218 : "displace_node_list_by_side_list",
219 : "Automatic definition of mesh element sides entities");
220 612750 : params.addParamNamesToGroup("partitioner centroid_partitioner_direction", "Partitioning");
221 :
222 408500 : return params;
223 204250 : }
224 :
225 67216 : MooseMesh::MooseMesh(const InputParameters & parameters)
226 : : MooseObject(parameters),
227 : Restartable(this, "Mesh"),
228 : PerfGraphInterface(this),
229 67216 : _parallel_type(getParam<MooseEnum>("parallel_type").getEnum<MooseMesh::ParallelType>()),
230 67216 : _use_distributed_mesh(false),
231 67216 : _distribution_overridden(false),
232 67216 : _parallel_type_overridden(false),
233 67216 : _mesh(nullptr),
234 134432 : _partitioner_name(getParam<MooseEnum>("partitioner")),
235 67216 : _partitioner_overridden(false),
236 67216 : _custom_partitioner_requested(false),
237 67216 : _uniform_refine_level(0),
238 134432 : _skip_refine_when_use_split(getParam<bool>("skip_refine_when_use_split")),
239 67216 : _skip_deletion_repartition_after_refine(false),
240 67216 : _is_nemesis(false),
241 134432 : _patch_size(getParam<unsigned int>("patch_size")),
242 134432 : _ghosting_patch_size(isParamValid("ghosting_patch_size")
243 134432 : ? getParam<unsigned int>("ghosting_patch_size")
244 67216 : : 5 * _patch_size),
245 134432 : _max_leaf_size(getParam<unsigned int>("max_leaf_size")),
246 67216 : _patch_update_strategy(
247 134432 : getParam<MooseEnum>("patch_update_strategy").getEnum<Moose::PatchUpdateType>()),
248 67216 : _regular_orthogonal_mesh(false),
249 134432 : _is_split(getParam<bool>("_is_split")),
250 67216 : _allow_recovery(true),
251 134432 : _construct_node_list_from_side_list(getParam<bool>("construct_node_list_from_side_list")),
252 134432 : _displace_node_list_by_side_list(getParam<bool>("displace_node_list_by_side_list")),
253 67216 : _need_delete(false),
254 67216 : _allow_remote_element_removal(true),
255 67216 : _need_ghost_ghosted_boundaries(true),
256 67216 : _is_displaced(false),
257 67216 : _coord_sys(
258 134432 : declareRestartableData<std::map<SubdomainID, Moose::CoordinateSystemType>>("coord_sys")),
259 134432 : _rz_coord_axis(getParam<MooseEnum>("rz_coord_axis")),
260 67216 : _coord_system_set(false),
261 655497 : _doing_p_refinement(false)
262 : {
263 201648 : if (isParamValid("ghosting_patch_size") && (_patch_update_strategy != Moose::Iteration))
264 0 : mooseError("Ghosting patch size parameter has to be set in the mesh block "
265 : "only when 'iteration' patch update strategy is used.");
266 :
267 201648 : if (isParamValid("coord_block"))
268 : {
269 72 : if (isParamValid("block"))
270 0 : paramWarning("block",
271 : "You set both 'Mesh/block' and 'Mesh/coord_block'. The value of "
272 : "'Mesh/coord_block' will be used.");
273 :
274 72 : _provided_coord_blocks = getParam<std::vector<SubdomainName>>("coord_block");
275 : }
276 201576 : else if (isParamValid("block"))
277 765 : _provided_coord_blocks = getParam<std::vector<SubdomainName>>("block");
278 :
279 201648 : if (getParam<bool>("build_all_side_lowerd_mesh"))
280 : // Do not initially allow removal of remote elements
281 223 : allowRemoteElementRemoval(false);
282 :
283 67216 : determineUseDistributedMesh();
284 :
285 : #ifdef MOOSE_KOKKOS_ENABLED
286 50553 : if (_app.isKokkosAvailable())
287 50553 : _kokkos_mesh = std::make_unique<Moose::Kokkos::Mesh>(*this);
288 : #endif
289 67216 : }
290 :
291 2987 : MooseMesh::MooseMesh(const MooseMesh & other_mesh)
292 : : MooseObject(other_mesh._pars),
293 : Restartable(this, "Mesh"),
294 : PerfGraphInterface(this, "CopiedMesh"),
295 2987 : _built_from_other_mesh(true),
296 2987 : _parallel_type(other_mesh._parallel_type),
297 2987 : _use_distributed_mesh(other_mesh._use_distributed_mesh),
298 2987 : _distribution_overridden(other_mesh._distribution_overridden),
299 2987 : _parallel_type_overridden(other_mesh._parallel_type_overridden),
300 2987 : _mesh(other_mesh.getMesh().clone()),
301 2987 : _partitioner_name(other_mesh._partitioner_name),
302 2987 : _partitioner_overridden(other_mesh._partitioner_overridden),
303 2987 : _custom_partitioner_requested(other_mesh._custom_partitioner_requested),
304 2987 : _uniform_refine_level(other_mesh.uniformRefineLevel()),
305 2987 : _skip_refine_when_use_split(other_mesh._skip_refine_when_use_split),
306 2987 : _skip_deletion_repartition_after_refine(other_mesh._skip_deletion_repartition_after_refine),
307 2987 : _is_nemesis(other_mesh._is_nemesis),
308 2987 : _patch_size(other_mesh._patch_size),
309 2987 : _ghosting_patch_size(other_mesh._ghosting_patch_size),
310 2987 : _max_leaf_size(other_mesh._max_leaf_size),
311 2987 : _patch_update_strategy(other_mesh._patch_update_strategy),
312 2987 : _regular_orthogonal_mesh(false),
313 2987 : _is_split(other_mesh._is_split),
314 2987 : _lower_d_interior_blocks(other_mesh._lower_d_interior_blocks),
315 2987 : _lower_d_boundary_blocks(other_mesh._lower_d_boundary_blocks),
316 2987 : _allow_recovery(other_mesh._allow_recovery),
317 2987 : _construct_node_list_from_side_list(other_mesh._construct_node_list_from_side_list),
318 2987 : _displace_node_list_by_side_list(other_mesh._displace_node_list_by_side_list),
319 2987 : _need_delete(other_mesh._need_delete),
320 2987 : _allow_remote_element_removal(other_mesh._allow_remote_element_removal),
321 2987 : _need_ghost_ghosted_boundaries(other_mesh._need_ghost_ghosted_boundaries),
322 2987 : _coord_sys(other_mesh._coord_sys),
323 2987 : _rz_coord_axis(other_mesh._rz_coord_axis),
324 2987 : _subdomain_id_to_rz_coord_axis(other_mesh._subdomain_id_to_rz_coord_axis),
325 2987 : _coord_system_set(other_mesh._coord_system_set),
326 2987 : _provided_coord_blocks(other_mesh._provided_coord_blocks),
327 29110 : _doing_p_refinement(other_mesh._doing_p_refinement)
328 : {
329 2987 : _bounds.resize(other_mesh._bounds.size());
330 3296 : for (std::size_t i = 0; i < _bounds.size(); ++i)
331 : {
332 309 : _bounds[i].resize(other_mesh._bounds[i].size());
333 927 : for (std::size_t j = 0; j < _bounds[i].size(); ++j)
334 618 : _bounds[i][j] = other_mesh._bounds[i][j];
335 : }
336 :
337 2987 : updateCoordTransform();
338 :
339 : #ifdef MOOSE_KOKKOS_ENABLED
340 2227 : if (_app.isKokkosAvailable())
341 2227 : _kokkos_mesh = std::make_unique<Moose::Kokkos::Mesh>(*this);
342 : #endif
343 2987 : }
344 :
345 66236 : MooseMesh::~MooseMesh()
346 : {
347 66236 : freeBndNodes();
348 66236 : freeBndElems();
349 66236 : clearQuadratureNodes();
350 66236 : }
351 :
352 : void
353 221648 : MooseMesh::freeBndNodes()
354 : {
355 : // free memory
356 12616183 : for (auto & bnode : _bnd_nodes)
357 12394535 : delete bnode;
358 :
359 822622 : for (auto & it : _node_set_nodes)
360 600974 : it.second.clear();
361 :
362 221648 : _node_set_nodes.clear();
363 :
364 822773 : for (auto & it : _bnd_node_ids)
365 601125 : it.second.clear();
366 :
367 221648 : _bnd_node_ids.clear();
368 221648 : _bnd_node_range.reset();
369 221648 : }
370 :
371 : void
372 221648 : MooseMesh::freeBndElems()
373 : {
374 : // free memory
375 9638580 : for (auto & belem : _bnd_elems)
376 9416932 : delete belem;
377 :
378 801907 : for (auto & it : _bnd_elem_ids)
379 580259 : it.second.clear();
380 :
381 221648 : _bnd_elem_ids.clear();
382 221648 : _bnd_elem_range.reset();
383 221648 : }
384 :
385 : bool
386 137659 : MooseMesh::prepare(const MeshBase * const mesh_to_clone)
387 : {
388 688295 : TIME_SECTION("prepare", 2, "Preparing Mesh", true);
389 :
390 : parallel_object_only();
391 :
392 137659 : bool libmesh_mesh_prepared = false;
393 :
394 : mooseAssert(_mesh, "The MeshBase has not been constructed");
395 :
396 137659 : if (!dynamic_cast<DistributedMesh *>(&getMesh()) || _is_nemesis)
397 : // For whatever reason we do not want to allow renumbering here nor ever in the future?
398 114613 : getMesh().allow_renumbering(false);
399 :
400 137659 : if (mesh_to_clone)
401 : {
402 : mooseAssert(mesh_to_clone->is_prepared(),
403 : "The mesh we wish to clone from must already be prepared");
404 149 : _mesh = mesh_to_clone->clone();
405 149 : _moose_mesh_prepared = false;
406 : }
407 137510 : else if (!_mesh->is_prepared())
408 : {
409 18184 : _mesh->complete_preparation();
410 18184 : _moose_mesh_prepared = false;
411 18184 : libmesh_mesh_prepared = true;
412 : }
413 :
414 137659 : if (_moose_mesh_prepared)
415 68890 : return libmesh_mesh_prepared;
416 :
417 : // Collect (local) subdomain IDs
418 68769 : _mesh_subdomains.clear();
419 13349369 : for (const auto & elem : getMesh().element_ptr_range())
420 13349369 : _mesh_subdomains.insert(elem->subdomain_id());
421 :
422 68769 : bool need_subdomain_name_map_sync = false;
423 : // add explicitly requested subdomains
424 206575 : if (isParamValid("add_subdomain_ids") && !isParamValid("add_subdomain_names"))
425 : {
426 : // only subdomain ids are explicitly given
427 72 : const auto & add_subdomain_id = getParam<std::vector<SubdomainID>>("add_subdomain_ids");
428 36 : _mesh_subdomains.insert(add_subdomain_id.begin(), add_subdomain_id.end());
429 : }
430 206395 : else if (isParamValid("add_subdomain_ids") && isParamValid("add_subdomain_names"))
431 : {
432 : const auto add_subdomain =
433 392 : getParam<SubdomainID, SubdomainName>("add_subdomain_ids", "add_subdomain_names");
434 244 : for (const auto & [sub_id, sub_name] : add_subdomain)
435 : {
436 : // add subdomain id
437 146 : _mesh_subdomains.insert(sub_id);
438 : // set name of the subdomain just added
439 146 : setSubdomainName(sub_id, sub_name);
440 : }
441 98 : need_subdomain_name_map_sync = true;
442 98 : }
443 205905 : else if (isParamValid("add_subdomain_names"))
444 : {
445 : // the user has defined add_subdomain_names, but not add_subdomain_ids
446 24 : const auto & add_subdomain_names = getParam<std::vector<SubdomainName>>("add_subdomain_names");
447 :
448 : // to define subdomain ids, we need the largest subdomain id defined yet.
449 12 : subdomain_id_type offset = 0;
450 12 : if (!_mesh_subdomains.empty())
451 12 : offset = *_mesh_subdomains.rbegin();
452 :
453 : // add all subdomains (and auto-assign ids)
454 48 : for (const SubdomainName & sub_name : add_subdomain_names)
455 : {
456 : // to avoid two subdomains with the same ID (notably on recover)
457 36 : if (getSubdomainID(sub_name) != libMesh::Elem::invalid_subdomain_id)
458 3 : continue;
459 33 : const auto sub_id = ++offset;
460 : // add subdomain id
461 33 : _mesh_subdomains.insert(sub_id);
462 : // set name of the subdomain just added
463 33 : setSubdomainName(sub_id, sub_name);
464 : }
465 12 : need_subdomain_name_map_sync = true;
466 : }
467 68769 : if (need_subdomain_name_map_sync)
468 110 : _mesh->sync_subdomain_name_map();
469 :
470 : // Make sure nodesets have been generated
471 68769 : buildNodeListFromSideList();
472 :
473 : // Collect (local) boundary IDs
474 68769 : const std::set<BoundaryID> & local_bids = getMesh().get_boundary_info().get_boundary_ids();
475 68769 : _mesh_boundary_ids.insert(local_bids.begin(), local_bids.end());
476 :
477 : const std::set<BoundaryID> & local_node_bids =
478 68769 : getMesh().get_boundary_info().get_node_boundary_ids();
479 68769 : _mesh_nodeset_ids.insert(local_node_bids.begin(), local_node_bids.end());
480 :
481 : const std::set<BoundaryID> & local_side_bids =
482 68769 : getMesh().get_boundary_info().get_side_boundary_ids();
483 68769 : _mesh_sideset_ids.insert(local_side_bids.begin(), local_side_bids.end());
484 :
485 : // Add explicitly requested sidesets/nodesets
486 : // This is done *after* the side boundaries (e.g. "right", ...) have been generated.
487 137538 : auto add_sets = [this](const bool sidesets, auto & set_ids)
488 : {
489 137538 : const std::string type = sidesets ? "sideset" : "nodeset";
490 137538 : const std::string id_param = "add_" + type + "_ids";
491 137538 : const std::string name_param = "add_" + type + "_names";
492 :
493 137538 : if (isParamValid(id_param))
494 : {
495 54 : const auto & add_ids = getParam<std::vector<BoundaryID>>(id_param);
496 54 : _mesh_boundary_ids.insert(add_ids.begin(), add_ids.end());
497 54 : set_ids.insert(add_ids.begin(), add_ids.end());
498 54 : if (isParamValid(name_param))
499 : {
500 42 : const auto & add_names = getParam<std::vector<BoundaryName>>(name_param);
501 : mooseAssert(add_names.size() == add_ids.size(),
502 : "Id and name sets must be the same size when adding.");
503 114 : for (const auto i : index_range(add_ids))
504 72 : setBoundaryName(add_ids[i], add_names[i]);
505 : }
506 : }
507 137484 : else if (isParamValid(name_param))
508 : {
509 : // the user has defined names, but not ids
510 12 : const auto & add_names = getParam<std::vector<BoundaryName>>(name_param);
511 :
512 12 : auto & mesh_ids = sidesets ? _mesh_sideset_ids : _mesh_nodeset_ids;
513 :
514 : // to define ids, we need the largest id defined yet.
515 12 : boundary_id_type offset = 0;
516 12 : if (!mesh_ids.empty())
517 12 : offset = *mesh_ids.rbegin();
518 12 : if (!_mesh_boundary_ids.empty())
519 12 : offset = std::max(offset, *_mesh_boundary_ids.rbegin());
520 :
521 : // add all sidesets/nodesets (and auto-assign ids)
522 24 : for (const auto & name : add_names)
523 : {
524 : // to avoid two sets with the same ID (notably on recover)
525 12 : if (getBoundaryID(name) != Moose::INVALID_BOUNDARY_ID)
526 1 : continue;
527 11 : const auto id = ++offset;
528 : // add sideset id
529 11 : _mesh_boundary_ids.insert(id);
530 11 : set_ids.insert(id);
531 : // set name of the sideset just added
532 11 : setBoundaryName(id, name);
533 : }
534 : }
535 137538 : };
536 :
537 68769 : add_sets(true, _mesh_sideset_ids);
538 68769 : add_sets(false, _mesh_nodeset_ids);
539 :
540 : // Communicate subdomain and boundary IDs if this is a parallel mesh
541 68769 : if (!getMesh().is_serial())
542 : {
543 8937 : _communicator.set_union(_mesh_subdomains);
544 8937 : _communicator.set_union(_mesh_boundary_ids);
545 8937 : _communicator.set_union(_mesh_nodeset_ids);
546 8937 : _communicator.set_union(_mesh_sideset_ids);
547 : }
548 :
549 68769 : if (!_built_from_other_mesh)
550 : {
551 65956 : if (!_coord_system_set)
552 197640 : setCoordSystem(_provided_coord_blocks, getParam<MultiMooseEnum>("coord_type"));
553 228 : else if (_pars.isParamSetByUser("coord_type"))
554 0 : mooseError(
555 : "Trying to set coordinate system type information based on the user input file, but "
556 : "the coordinate system type information has already been set programmatically! "
557 : "Either remove your coordinate system type information from the input file, or contact "
558 : "your application developer");
559 : }
560 :
561 : // Set general axisymmetric axes if provided
562 275127 : if (isParamValid("rz_coord_blocks") && isParamValid("rz_coord_origins") &&
563 68820 : isParamValid("rz_coord_directions"))
564 : {
565 34 : const auto rz_coord_blocks = getParam<std::vector<SubdomainName>>("rz_coord_blocks");
566 34 : const auto rz_coord_origins = getParam<std::vector<Point>>("rz_coord_origins");
567 34 : const auto rz_coord_directions = getParam<std::vector<RealVectorValue>>("rz_coord_directions");
568 34 : if (rz_coord_origins.size() == rz_coord_blocks.size() &&
569 17 : rz_coord_directions.size() == rz_coord_blocks.size())
570 : {
571 17 : std::vector<std::pair<Point, RealVectorValue>> rz_coord_axes;
572 58 : for (unsigned int i = 0; i < rz_coord_origins.size(); ++i)
573 41 : rz_coord_axes.push_back(std::make_pair(rz_coord_origins[i], rz_coord_directions[i]));
574 :
575 17 : setGeneralAxisymmetricCoordAxes(rz_coord_blocks, rz_coord_axes);
576 :
577 51 : if (isParamSetByUser("rz_coord_axis"))
578 0 : mooseError("The parameter 'rz_coord_axis' may not be provided if 'rz_coord_blocks', "
579 : "'rz_coord_origins', and 'rz_coord_directions' are provided.");
580 17 : }
581 : else
582 0 : mooseError("The parameters 'rz_coord_blocks', 'rz_coord_origins', and "
583 : "'rz_coord_directions' must all have the same size.");
584 17 : }
585 481264 : else if (isParamValid("rz_coord_blocks") || isParamValid("rz_coord_origins") ||
586 275008 : isParamValid("rz_coord_directions"))
587 0 : mooseError("If any of the parameters 'rz_coord_blocks', 'rz_coord_origins', and "
588 : "'rz_coord_directions' are provided, then all must be provided.");
589 :
590 68769 : detectOrthogonalDimRanges();
591 :
592 68769 : update();
593 :
594 : // Check if there is subdomain name duplication for the same subdomain ID
595 68769 : checkDuplicateSubdomainNames();
596 :
597 68766 : _moose_mesh_prepared = true;
598 :
599 68766 : return libmesh_mesh_prepared;
600 137656 : }
601 :
602 : bool
603 155412 : MooseMesh::possiblyRebuildNodeToElemMap()
604 : {
605 : // *Rebuild* the node to element map. I emphasize rebuild because if it has not been built
606 : // previously we won't do anything
607 155412 : if (!_node_to_elem_map_built)
608 : {
609 : mooseAssert(_node_to_elem_map.empty(), "If it hasn't been built, it better well be empty");
610 147613 : return false;
611 : }
612 :
613 7799 : _node_to_elem_map.clear();
614 7799 : _node_to_elem_map_built = false;
615 7799 : internalNodeToElemMap();
616 7799 : return true;
617 : }
618 :
619 : void
620 155412 : MooseMesh::update()
621 : {
622 777060 : TIME_SECTION("update", 3, "Updating Mesh", true);
623 :
624 : // Rebuild the boundary conditions
625 155412 : buildNodeListFromSideList();
626 :
627 155412 : buildNodeList();
628 155412 : buildBndElemList();
629 155412 : cacheInfo();
630 155412 : buildElemIDInfo();
631 :
632 : // this will make moose mesh aware of p-refinement added by mesh generators including
633 : // a file mesh generator loading a restart checkpoint file
634 155412 : _max_p_level = 0;
635 155412 : _max_h_level = 0;
636 27780138 : for (const auto & elem : getMesh().active_local_element_ptr_range())
637 : {
638 27624726 : if (elem->p_level() > _max_p_level)
639 658 : _max_p_level = elem->p_level();
640 27624726 : if (elem->level() > _max_h_level)
641 25822 : _max_h_level = elem->level();
642 155412 : }
643 155412 : comm().max(_max_p_level);
644 155412 : comm().max(_max_h_level);
645 :
646 : // the flag might have been set by calling doingPRefinement(true)
647 155412 : _doing_p_refinement = _doing_p_refinement || (_max_p_level > 0);
648 :
649 155412 : computeMaxPerElemAndSide();
650 :
651 : #ifdef MOOSE_KOKKOS_ENABLED
652 129109 : if (_app.getExecutioner() && _app.feProblem().initialized() &&
653 12599 : _app.feProblem().hasKokkosObjects())
654 0 : _kokkos_mesh->update();
655 : #endif
656 :
657 155412 : _finite_volume_info_dirty = true;
658 :
659 155412 : possiblyRebuildNodeToElemMap();
660 155412 : }
661 :
662 : void
663 205 : MooseMesh::buildLowerDMesh()
664 : {
665 205 : auto & mesh = getMesh();
666 :
667 205 : if (!mesh.is_serial())
668 0 : mooseError(
669 : "Hybrid finite element method must use replicated mesh.\nCurrently lower-dimensional mesh "
670 : "does not support mesh re-partitioning and a debug assertion being hit related with "
671 : "neighbors of lower-dimensional element, with distributed mesh.");
672 :
673 : // Lower-D element build requires neighboring element information
674 205 : if (!mesh.is_prepared())
675 194 : mesh.find_neighbors();
676 :
677 : // maximum number of sides of all elements
678 205 : unsigned int max_n_sides = 0;
679 :
680 : // remove existing lower-d element first
681 205 : std::set<Elem *> deleteable_elems;
682 4347 : for (auto & elem : mesh.element_ptr_range())
683 4142 : if (_lower_d_interior_blocks.count(elem->subdomain_id()) ||
684 2071 : _lower_d_boundary_blocks.count(elem->subdomain_id()))
685 0 : deleteable_elems.insert(elem);
686 2071 : else if (elem->n_sides() > max_n_sides)
687 410 : max_n_sides = elem->n_sides();
688 :
689 205 : for (auto & elem : deleteable_elems)
690 0 : mesh.delete_elem(elem);
691 205 : for (const auto & id : _lower_d_interior_blocks)
692 0 : _mesh_subdomains.erase(id);
693 205 : for (const auto & id : _lower_d_boundary_blocks)
694 0 : _mesh_subdomains.erase(id);
695 205 : _lower_d_interior_blocks.clear();
696 205 : _lower_d_boundary_blocks.clear();
697 :
698 205 : mesh.comm().max(max_n_sides);
699 :
700 205 : deleteable_elems.clear();
701 :
702 : // get all side types
703 205 : std::set<int> interior_side_types;
704 205 : std::set<int> boundary_side_types;
705 4347 : for (const auto & elem : mesh.active_element_ptr_range())
706 11173 : for (const auto side : elem->side_index_range())
707 : {
708 9102 : Elem * neig = elem->neighbor_ptr(side);
709 9102 : std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
710 9102 : if (neig)
711 5956 : interior_side_types.insert(side_elem->type());
712 : else
713 3146 : boundary_side_types.insert(side_elem->type());
714 9307 : }
715 205 : mesh.comm().set_union(interior_side_types);
716 205 : mesh.comm().set_union(boundary_side_types);
717 :
718 : // assign block ids for different side types
719 205 : std::map<ElemType, SubdomainID> interior_block_ids;
720 205 : std::map<ElemType, SubdomainID> boundary_block_ids;
721 : // we assume this id is not used by the mesh
722 205 : auto id = libMesh::Elem::invalid_subdomain_id - 2;
723 424 : for (const auto & tpid : interior_side_types)
724 : {
725 219 : const auto type = ElemType(tpid);
726 219 : mesh.subdomain_name(id) = "INTERNAL_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type);
727 219 : interior_block_ids[type] = id;
728 219 : _lower_d_interior_blocks.insert(id);
729 219 : if (_mesh_subdomains.count(id) > 0)
730 0 : mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh");
731 219 : _mesh_subdomains.insert(id);
732 219 : --id;
733 : }
734 424 : for (const auto & tpid : boundary_side_types)
735 : {
736 219 : const auto type = ElemType(tpid);
737 219 : mesh.subdomain_name(id) = "BOUNDARY_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type);
738 219 : boundary_block_ids[type] = id;
739 219 : _lower_d_boundary_blocks.insert(id);
740 219 : if (_mesh_subdomains.count(id) > 0)
741 0 : mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh");
742 219 : _mesh_subdomains.insert(id);
743 219 : --id;
744 : }
745 :
746 205 : dof_id_type max_elem_id = mesh.max_elem_id();
747 205 : unique_id_type max_unique_id = mesh.parallel_max_unique_id();
748 :
749 205 : std::vector<Elem *> side_elems;
750 205 : _higher_d_elem_side_to_lower_d_elem.clear();
751 4347 : for (const auto & elem : mesh.active_element_ptr_range())
752 : {
753 : // skip existing lower-d elements
754 2071 : if (elem->interior_parent())
755 0 : continue;
756 :
757 11173 : for (const auto side : elem->side_index_range())
758 : {
759 9102 : Elem * neig = elem->neighbor_ptr(side);
760 :
761 9102 : bool build_side = false;
762 9102 : if (!neig)
763 3146 : build_side = true;
764 : else
765 : {
766 : mooseAssert(!neig->is_remote(), "We error if the mesh is not serial");
767 5956 : if (!neig->active())
768 0 : build_side = true;
769 5956 : else if (neig->level() == elem->level() && elem->id() < neig->id())
770 2978 : build_side = true;
771 : }
772 :
773 9102 : if (build_side)
774 : {
775 6124 : std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
776 :
777 : // The side will be added with the same processor id as the parent.
778 6124 : side_elem->processor_id() = elem->processor_id();
779 :
780 : // Add subdomain ID
781 6124 : if (neig)
782 2978 : side_elem->subdomain_id() = interior_block_ids.at(side_elem->type());
783 : else
784 3146 : side_elem->subdomain_id() = boundary_block_ids.at(side_elem->type());
785 :
786 : // set ids consistently across processors (these ids will be temporary)
787 6124 : side_elem->set_id(max_elem_id + elem->id() * max_n_sides + side);
788 6124 : side_elem->set_unique_id(max_unique_id + elem->id() * max_n_sides + side);
789 :
790 : // Also assign the side's interior parent, so it is always
791 : // easy to figure out the Elem we came from.
792 : // Note: the interior parent could be a ghost element.
793 6124 : side_elem->set_interior_parent(elem);
794 :
795 6124 : side_elems.push_back(side_elem.release());
796 :
797 : // add link between higher d element to lower d element
798 6124 : auto pair = std::make_pair(elem, side);
799 6124 : auto link = std::make_pair(pair, side_elems.back());
800 6124 : auto ilink = std::make_pair(side_elems.back(), side);
801 6124 : _lower_d_elem_to_higher_d_elem_side.insert(ilink);
802 6124 : _higher_d_elem_side_to_lower_d_elem.insert(link);
803 6124 : }
804 : }
805 205 : }
806 :
807 : // finally, add the lower-dimensional element to the mesh
808 : // Note: lower-d interior element will exist on a processor if its associated interior
809 : // parent exists on a processor whether or not being a ghost. Lower-d elements will
810 : // get its interior parent's processor id.
811 6329 : for (auto & elem : side_elems)
812 6124 : mesh.add_elem(elem);
813 :
814 : // we do all the stuff in prepare_for_use such as renumber_nodes_and_elements(),
815 : // update_parallel_id_counts(), cache_elem_dims(), etc. except partitioning here.
816 205 : const bool skip_partitioning_old = mesh.skip_partitioning();
817 205 : mesh.skip_partitioning(true);
818 : // Finding neighbors is ambiguous for lower-dimensional elements on interior faces
819 205 : mesh.allow_find_neighbors(false);
820 205 : mesh.prepare_for_use();
821 205 : mesh.skip_partitioning(skip_partitioning_old);
822 205 : }
823 :
824 : const Node &
825 0 : MooseMesh::node(const dof_id_type i) const
826 : {
827 0 : mooseDeprecated("MooseMesh::node() is deprecated, please use MooseMesh::nodeRef() instead");
828 0 : return nodeRef(i);
829 : }
830 :
831 : Node &
832 0 : MooseMesh::node(const dof_id_type i)
833 : {
834 0 : mooseDeprecated("MooseMesh::node() is deprecated, please use MooseMesh::nodeRef() instead");
835 0 : return nodeRef(i);
836 : }
837 :
838 : const Node &
839 42686230 : MooseMesh::nodeRef(const dof_id_type i) const
840 : {
841 42686230 : const auto node_ptr = queryNodePtr(i);
842 : mooseAssert(node_ptr, "Missing node");
843 42686230 : return *node_ptr;
844 : }
845 :
846 : Node &
847 24013498 : MooseMesh::nodeRef(const dof_id_type i)
848 : {
849 24013498 : return const_cast<Node &>(const_cast<const MooseMesh *>(this)->nodeRef(i));
850 : }
851 :
852 : const Node *
853 0 : MooseMesh::nodePtr(const dof_id_type i) const
854 : {
855 0 : return &nodeRef(i);
856 : }
857 :
858 : Node *
859 2096 : MooseMesh::nodePtr(const dof_id_type i)
860 : {
861 2096 : return &nodeRef(i);
862 : }
863 :
864 : const Node *
865 42688974 : MooseMesh::queryNodePtr(const dof_id_type i) const
866 : {
867 42688974 : if (i > getMesh().max_node_id())
868 : {
869 196773 : auto it = _quadrature_nodes.find(i);
870 196773 : if (it == _quadrature_nodes.end())
871 0 : return nullptr;
872 196773 : auto & node_ptr = it->second;
873 : mooseAssert(node_ptr, "Uninitialized quadrature node");
874 196773 : return node_ptr;
875 : }
876 :
877 42492201 : return getMesh().query_node_ptr(i);
878 : }
879 :
880 : Node *
881 2744 : MooseMesh::queryNodePtr(const dof_id_type i)
882 : {
883 2744 : return const_cast<Node *>(const_cast<const MooseMesh *>(this)->queryNodePtr(i));
884 : }
885 :
886 : void
887 84008 : MooseMesh::meshChanged()
888 : {
889 420040 : TIME_SECTION("meshChanged", 3, "Updating Because Mesh Changed");
890 :
891 84008 : update();
892 :
893 : // Delete all of the cached ranges
894 84008 : _active_node_range.reset();
895 84008 : _active_semilocal_node_range.reset();
896 84008 : _local_node_range.reset();
897 84008 : _bnd_node_range.reset();
898 84008 : _bnd_elem_range.reset();
899 :
900 : // Rebuild the ranges
901 84008 : getActiveLocalElementRange();
902 84008 : getActiveNodeRange();
903 84008 : getLocalNodeRange();
904 84008 : getBoundaryNodeRange();
905 84008 : getBoundaryElementRange();
906 :
907 : // Call the callback function onMeshChanged
908 84008 : onMeshChanged();
909 84008 : }
910 :
911 : void
912 84008 : MooseMesh::onMeshChanged()
913 : {
914 84008 : }
915 :
916 : void
917 208 : MooseMesh::cacheChangedLists()
918 : {
919 1040 : TIME_SECTION("cacheChangedLists", 5, "Caching Changed Lists");
920 :
921 208 : ConstElemRange elem_range(getMesh().local_elements_begin(), getMesh().local_elements_end(), 1);
922 208 : CacheChangedListsThread cclt(*this);
923 208 : Threads::parallel_reduce(elem_range, cclt);
924 :
925 208 : _coarsened_element_children.clear();
926 :
927 416 : _refined_elements = std::make_unique<ConstElemPointerRange>(cclt._refined_elements.begin(),
928 416 : cclt._refined_elements.end());
929 416 : _coarsened_elements = std::make_unique<ConstElemPointerRange>(cclt._coarsened_elements.begin(),
930 416 : cclt._coarsened_elements.end());
931 208 : _coarsened_element_children = cclt._coarsened_element_children;
932 208 : }
933 :
934 : ConstElemPointerRange *
935 208 : MooseMesh::refinedElementRange() const
936 : {
937 208 : return _refined_elements.get();
938 : }
939 :
940 : ConstElemPointerRange *
941 208 : MooseMesh::coarsenedElementRange() const
942 : {
943 208 : return _coarsened_elements.get();
944 : }
945 :
946 : const std::vector<const Elem *> &
947 2468 : MooseMesh::coarsenedElementChildren(const Elem * elem) const
948 : {
949 2468 : auto elem_to_child_pair = _coarsened_element_children.find(elem);
950 : mooseAssert(elem_to_child_pair != _coarsened_element_children.end(), "Missing element in map");
951 4936 : return elem_to_child_pair->second;
952 : }
953 :
954 : void
955 71870 : MooseMesh::updateActiveSemiLocalNodeRange(std::set<dof_id_type> & ghosted_elems)
956 : {
957 359350 : TIME_SECTION("updateActiveSemiLocalNodeRange", 5, "Updating ActiveSemiLocalNode Range");
958 :
959 71870 : _semilocal_node_list.clear();
960 :
961 : // First add the nodes connected to local elems
962 71870 : const ConstElemRange * active_local_elems = getActiveLocalElementRange();
963 13079648 : for (const auto & elem : *active_local_elems)
964 : {
965 82411462 : for (unsigned int n = 0; n < elem->n_nodes(); ++n)
966 : {
967 : // Since elem is const here but we require a non-const Node * to
968 : // store in the _semilocal_node_list (otherwise things like
969 : // UpdateDisplacedMeshThread don't work), we are using a
970 : // const_cast. A more long-term fix would be to have
971 : // getActiveLocalElementRange return a non-const ElemRange.
972 69403684 : Node * node = const_cast<Node *>(elem->node_ptr(n));
973 :
974 69403684 : _semilocal_node_list.insert(node);
975 : }
976 : }
977 :
978 : // Now add the nodes connected to ghosted_elems
979 117081 : for (const auto & ghost_elem_id : ghosted_elems)
980 : {
981 45211 : Elem * elem = getMesh().elem_ptr(ghost_elem_id);
982 250228 : for (unsigned int n = 0; n < elem->n_nodes(); n++)
983 : {
984 205017 : Node * node = elem->node_ptr(n);
985 :
986 205017 : _semilocal_node_list.insert(node);
987 : }
988 : }
989 :
990 : // Now create the actual range
991 143740 : _active_semilocal_node_range = std::make_unique<SemiLocalNodeRange>(_semilocal_node_list.begin(),
992 143740 : _semilocal_node_list.end());
993 71870 : }
994 :
995 : bool
996 26489 : MooseMesh::isSemiLocal(Node * const node) const
997 : {
998 26489 : return _semilocal_node_list.find(node) != _semilocal_node_list.end();
999 : }
1000 :
1001 : /**
1002 : * Helper class for sorting Boundary Nodes so that we always get the same
1003 : * order of application for boundary conditions.
1004 : */
1005 : class BndNodeCompare
1006 : {
1007 : public:
1008 155412 : BndNodeCompare() {}
1009 :
1010 121361976 : bool operator()(const BndNode * const & lhs, const BndNode * const & rhs)
1011 : {
1012 121361976 : if (lhs->_bnd_id < rhs->_bnd_id)
1013 22593622 : return true;
1014 :
1015 98768354 : if (lhs->_bnd_id > rhs->_bnd_id)
1016 10264838 : return false;
1017 :
1018 88503516 : if (lhs->_node->id() < rhs->_node->id())
1019 57194272 : return true;
1020 :
1021 31309244 : if (lhs->_node->id() > rhs->_node->id())
1022 31309244 : return false;
1023 :
1024 0 : return false;
1025 : }
1026 : };
1027 :
1028 : void
1029 155412 : MooseMesh::buildNodeList()
1030 : {
1031 777060 : TIME_SECTION("buildNodeList", 5, "Building Node List");
1032 :
1033 155412 : freeBndNodes();
1034 :
1035 155412 : auto bc_tuples = getMesh().get_boundary_info().build_node_list();
1036 :
1037 155412 : int n = bc_tuples.size();
1038 155412 : _bnd_nodes.clear();
1039 155412 : _bnd_nodes.reserve(n);
1040 12677648 : for (const auto & t : bc_tuples)
1041 : {
1042 12522236 : auto node_id = std::get<0>(t);
1043 12522236 : auto bc_id = std::get<1>(t);
1044 :
1045 12522236 : _bnd_nodes.push_back(new BndNode(getMesh().node_ptr(node_id), bc_id));
1046 12522236 : _node_set_nodes[bc_id].push_back(node_id);
1047 12522236 : _bnd_node_ids[bc_id].insert(node_id);
1048 : }
1049 :
1050 155412 : _bnd_nodes.reserve(_bnd_nodes.size() + _extra_bnd_nodes.size());
1051 155466 : for (unsigned int i = 0; i < _extra_bnd_nodes.size(); i++)
1052 : {
1053 54 : BndNode * bnode = new BndNode(_extra_bnd_nodes[i]._node, _extra_bnd_nodes[i]._bnd_id);
1054 54 : _bnd_nodes.push_back(bnode);
1055 54 : _bnd_node_ids[std::get<1>(bc_tuples[i])].insert(_extra_bnd_nodes[i]._node->id());
1056 : }
1057 :
1058 : // This sort is here so that boundary conditions are always applied in the same order
1059 155412 : std::sort(_bnd_nodes.begin(), _bnd_nodes.end(), BndNodeCompare());
1060 155412 : }
1061 :
1062 : void
1063 155412 : MooseMesh::computeMaxPerElemAndSide()
1064 : {
1065 155412 : auto & mesh = getMesh();
1066 :
1067 155412 : _max_sides_per_elem = 0;
1068 155412 : _max_nodes_per_elem = 0;
1069 155412 : _max_nodes_per_side = 0;
1070 :
1071 59598722 : for (auto & elem : as_range(mesh.local_elements_begin(), mesh.local_elements_end()))
1072 : {
1073 29721655 : _max_sides_per_elem = std::max(_max_sides_per_elem, elem->n_sides());
1074 29721655 : _max_nodes_per_elem = std::max(_max_nodes_per_elem, elem->n_nodes());
1075 :
1076 162573703 : for (unsigned int side = 0; side < elem->n_sides(); ++side)
1077 132852048 : _max_nodes_per_side = std::max(_max_nodes_per_side, elem->side_ptr(side)->n_nodes());
1078 155412 : }
1079 :
1080 155412 : mesh.comm().max(_max_sides_per_elem);
1081 155412 : mesh.comm().max(_max_nodes_per_elem);
1082 155412 : mesh.comm().max(_max_nodes_per_side);
1083 155412 : }
1084 :
1085 : void
1086 155412 : MooseMesh::buildElemIDInfo()
1087 : {
1088 155412 : unsigned int n = getMesh().n_elem_integers() + 1;
1089 :
1090 155412 : _block_id_mapping.clear();
1091 155412 : _max_ids.clear();
1092 155412 : _min_ids.clear();
1093 155412 : _id_identical_flag.clear();
1094 :
1095 155412 : _block_id_mapping.resize(n);
1096 155412 : _max_ids.resize(n, std::numeric_limits<dof_id_type>::min());
1097 155412 : _min_ids.resize(n, std::numeric_limits<dof_id_type>::max());
1098 310824 : _id_identical_flag.resize(n, std::vector<bool>(n, true));
1099 27780138 : for (const auto & elem : getMesh().active_local_element_ptr_range())
1100 59374603 : for (unsigned int i = 0; i < n; ++i)
1101 : {
1102 31749877 : auto id = (i == n - 1 ? elem->subdomain_id() : elem->get_extra_integer(i));
1103 31749877 : _block_id_mapping[i][elem->subdomain_id()].insert(id);
1104 31749877 : if (id > _max_ids[i])
1105 122193 : _max_ids[i] = id;
1106 31749877 : if (id < _min_ids[i])
1107 160213 : _min_ids[i] = id;
1108 76437984 : for (unsigned int j = 0; j < n; ++j)
1109 : {
1110 44688107 : auto idj = (j == n - 1 ? elem->subdomain_id() : elem->get_extra_integer(j));
1111 44688107 : if (i != j && _id_identical_flag[i][j] && id != idj)
1112 6790 : _id_identical_flag[i][j] = false;
1113 : }
1114 155412 : }
1115 :
1116 313285 : for (unsigned int i = 0; i < n; ++i)
1117 : {
1118 383989 : for (auto & blk : meshSubdomains())
1119 226116 : comm().set_union(_block_id_mapping[i][blk]);
1120 157873 : comm().min(_id_identical_flag[i]);
1121 : }
1122 155412 : comm().max(_max_ids);
1123 155412 : comm().min(_min_ids);
1124 155412 : }
1125 :
1126 : std::unordered_map<dof_id_type, std::set<dof_id_type>>
1127 11 : MooseMesh::getElemIDMapping(const std::string & from_id_name, const std::string & to_id_name) const
1128 : {
1129 11 : auto & mesh_base = getMesh();
1130 :
1131 11 : if (!mesh_base.has_elem_integer(from_id_name))
1132 0 : mooseError("Mesh does not have the element integer name '", from_id_name, "'");
1133 11 : if (!mesh_base.has_elem_integer(to_id_name))
1134 0 : mooseError("Mesh does not have the element integer name '", to_id_name, "'");
1135 :
1136 11 : const auto id1 = mesh_base.get_elem_integer_index(from_id_name);
1137 11 : const auto id2 = mesh_base.get_elem_integer_index(to_id_name);
1138 :
1139 11 : std::unordered_map<dof_id_type, std::set<dof_id_type>> id_map;
1140 33 : for (const auto id : getAllElemIDs(id1))
1141 33 : id_map[id] = std::set<dof_id_type>();
1142 :
1143 811 : for (const auto & elem : mesh_base.active_local_element_ptr_range())
1144 811 : id_map[elem->get_extra_integer(id1)].insert(elem->get_extra_integer(id2));
1145 :
1146 33 : for (auto & [id, ids] : id_map)
1147 : {
1148 22 : libmesh_ignore(id); // avoid overzealous gcc 9.4 unused var warning
1149 22 : comm().set_union(ids);
1150 : }
1151 :
1152 11 : return id_map;
1153 0 : }
1154 :
1155 : std::set<dof_id_type>
1156 50 : MooseMesh::getAllElemIDs(unsigned int elem_id_index) const
1157 : {
1158 50 : std::set<dof_id_type> unique_ids;
1159 139 : for (auto & pair : _block_id_mapping[elem_id_index])
1160 319 : for (auto & id : pair.second)
1161 230 : unique_ids.insert(id);
1162 50 : return unique_ids;
1163 0 : }
1164 :
1165 : std::set<dof_id_type>
1166 152 : MooseMesh::getElemIDsOnBlocks(unsigned int elem_id_index, const std::set<SubdomainID> & blks) const
1167 : {
1168 152 : std::set<dof_id_type> unique_ids;
1169 379 : for (auto & blk : blks)
1170 : {
1171 227 : auto it = _block_id_mapping[elem_id_index].find(blk);
1172 227 : if (it == _block_id_mapping[elem_id_index].end())
1173 0 : mooseError("Block ", blk, " is not available on the mesh");
1174 :
1175 532 : for (auto & mid : it->second)
1176 305 : unique_ids.insert(mid);
1177 : }
1178 152 : return unique_ids;
1179 0 : }
1180 :
1181 : void
1182 155412 : MooseMesh::buildBndElemList()
1183 : {
1184 777060 : TIME_SECTION("buildBndElemList", 5, "Building Boundary Elements List");
1185 :
1186 155412 : freeBndElems();
1187 :
1188 155412 : auto bc_tuples = getMesh().get_boundary_info().build_active_side_list();
1189 :
1190 155412 : int n = bc_tuples.size();
1191 155412 : _bnd_elems.clear();
1192 155412 : _bnd_elems.reserve(n);
1193 9686238 : for (const auto & t : bc_tuples)
1194 : {
1195 9530826 : auto elem_id = std::get<0>(t);
1196 9530826 : auto side_id = std::get<1>(t);
1197 9530826 : auto bc_id = std::get<2>(t);
1198 :
1199 9530826 : _bnd_elems.push_back(new BndElement(getMesh().elem_ptr(elem_id), side_id, bc_id));
1200 9530826 : _bnd_elem_ids[bc_id].insert(elem_id);
1201 : }
1202 155412 : }
1203 :
1204 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> &
1205 997114 : MooseMesh::internalNodeToElemMap()
1206 : {
1207 997114 : if (!_node_to_elem_map_built) // Guard the creation with a double checked lock
1208 : {
1209 11867 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1210 :
1211 11867 : if (!_node_to_elem_map_built)
1212 : {
1213 : // This is allowing the timing to be run even with threads
1214 : // This is safe because all threads will be waiting on this section when it runs
1215 : // NOTE: Do not copy this construction to other places without thinking REALLY hard about it
1216 : // The PerfGraph is NOT threadsafe and will cause all kinds of havok if care isn't taken
1217 11867 : auto in_threads = Threads::in_threads;
1218 11867 : Threads::in_threads = false;
1219 59335 : TIME_SECTION("nodeToElemMap", 5, "Building Node To Elem Map");
1220 11867 : Threads::in_threads = in_threads;
1221 :
1222 : mooseAssert(_node_to_elem_map.empty(), "Expected empty map before building");
1223 2984309 : for (const auto & elem : getMesh().active_element_ptr_range())
1224 18205478 : for (unsigned int n = 0; n < elem->n_nodes(); n++)
1225 15244903 : _node_to_elem_map[elem->node_id(n)].push_back(elem->id());
1226 :
1227 11867 : _node_to_elem_map_built = true; // MUST be set at the end for double-checked locking to work!
1228 11867 : }
1229 11867 : }
1230 997114 : return _node_to_elem_map;
1231 : }
1232 :
1233 : const std::unordered_map<dof_id_type, std::vector<dof_id_type>> &
1234 983958 : MooseMesh::nodeToElemMap()
1235 : {
1236 983958 : return internalNodeToElemMap();
1237 : }
1238 :
1239 : const ConstElemRange *
1240 10384864 : MooseMesh::getActiveLocalElementRange()
1241 : {
1242 10384864 : return &getMesh().active_local_element_stored_range();
1243 : }
1244 :
1245 : NodeRange *
1246 84072 : MooseMesh::getActiveNodeRange()
1247 : {
1248 84072 : if (!_active_node_range)
1249 : {
1250 252024 : TIME_SECTION("getActiveNodeRange", 5);
1251 :
1252 : _active_node_range =
1253 84008 : std::make_unique<NodeRange>(getMesh().active_nodes_begin(), getMesh().active_nodes_end());
1254 84008 : }
1255 :
1256 84072 : return _active_node_range.get();
1257 : }
1258 :
1259 : SemiLocalNodeRange *
1260 0 : MooseMesh::getActiveSemiLocalNodeRange() const
1261 : {
1262 : mooseAssert(_active_semilocal_node_range,
1263 : "_active_semilocal_node_range has not been created yet!");
1264 :
1265 0 : return _active_semilocal_node_range.get();
1266 : }
1267 :
1268 : ConstNodeRange *
1269 318671 : MooseMesh::getLocalNodeRange()
1270 : {
1271 318671 : if (!_local_node_range)
1272 : {
1273 252024 : TIME_SECTION("getLocalNodeRange", 5);
1274 :
1275 168016 : _local_node_range = std::make_unique<ConstNodeRange>(getMesh().local_nodes_begin(),
1276 252024 : getMesh().local_nodes_end());
1277 84008 : }
1278 :
1279 318671 : return _local_node_range.get();
1280 : }
1281 :
1282 : ConstBndNodeRange *
1283 3675138 : MooseMesh::getBoundaryNodeRange()
1284 : {
1285 3675138 : if (!_bnd_node_range)
1286 : {
1287 252540 : TIME_SECTION("getBoundaryNodeRange", 5);
1288 :
1289 84180 : _bnd_node_range = std::make_unique<ConstBndNodeRange>(bndNodesBegin(), bndNodesEnd());
1290 84180 : }
1291 :
1292 3675138 : return _bnd_node_range.get();
1293 : }
1294 :
1295 : ConstBndElemRange *
1296 189164 : MooseMesh::getBoundaryElementRange()
1297 : {
1298 189164 : if (!_bnd_elem_range)
1299 : {
1300 252024 : TIME_SECTION("getBoundaryElementRange", 5);
1301 :
1302 84008 : _bnd_elem_range = std::make_unique<ConstBndElemRange>(bndElemsBegin(), bndElemsEnd());
1303 84008 : }
1304 :
1305 189164 : return _bnd_elem_range.get();
1306 : }
1307 :
1308 : const std::unordered_map<boundary_id_type, std::unordered_set<dof_id_type>> &
1309 0 : MooseMesh::getBoundariesToElems() const
1310 : {
1311 0 : mooseDeprecated("MooseMesh::getBoundariesToElems is deprecated, "
1312 : "use MooseMesh::getBoundariesToActiveSemiLocalElemIds");
1313 0 : return getBoundariesToActiveSemiLocalElemIds();
1314 : }
1315 :
1316 : const std::unordered_map<boundary_id_type, std::unordered_set<dof_id_type>> &
1317 61 : MooseMesh::getBoundariesToActiveSemiLocalElemIds() const
1318 : {
1319 61 : return _bnd_elem_ids;
1320 : }
1321 :
1322 : std::unordered_set<dof_id_type>
1323 3671 : MooseMesh::getBoundaryActiveSemiLocalElemIds(BoundaryID bid) const
1324 : {
1325 : // The boundary to element map is computed on every mesh update
1326 3671 : const auto it = _bnd_elem_ids.find(bid);
1327 3671 : if (it == _bnd_elem_ids.end())
1328 : // Boundary is not local to this domain, return an empty set
1329 118 : return std::unordered_set<dof_id_type>{};
1330 3553 : return it->second;
1331 : }
1332 :
1333 : std::unordered_set<dof_id_type>
1334 0 : MooseMesh::getBoundaryActiveNeighborElemIds(BoundaryID bid) const
1335 : {
1336 : // Vector of boundary elems is updated every mesh update
1337 0 : std::unordered_set<dof_id_type> neighbor_elems;
1338 0 : for (const auto & bnd_elem : _bnd_elems)
1339 : {
1340 0 : const auto & [elem_ptr, elem_side, elem_bid] = *bnd_elem;
1341 0 : if (elem_bid == bid)
1342 : {
1343 0 : const auto * neighbor = elem_ptr->neighbor_ptr(elem_side);
1344 : // Dont add fully remote elements, ghosted is fine
1345 0 : if (neighbor && neighbor != libMesh::remote_elem)
1346 : {
1347 : // handle mesh refinement, only return active elements near the boundary
1348 0 : if (neighbor->active())
1349 0 : neighbor_elems.insert(neighbor->id());
1350 : else
1351 : {
1352 0 : std::vector<const Elem *> family;
1353 0 : neighbor->active_family_tree_by_neighbor(family, elem_ptr);
1354 0 : for (const auto & child_neighbor : family)
1355 0 : neighbor_elems.insert(child_neighbor->id());
1356 0 : }
1357 : }
1358 : }
1359 : }
1360 :
1361 0 : return neighbor_elems;
1362 0 : }
1363 :
1364 : bool
1365 0 : MooseMesh::isBoundaryFullyExternalToSubdomains(BoundaryID bid,
1366 : const std::set<SubdomainID> & blk_group) const
1367 : {
1368 : mooseAssert(_bnd_elem_range, "Boundary element range is not initialized");
1369 :
1370 : // Loop over all side elements of the mesh, select those on the boundary
1371 0 : for (const auto & bnd_elem : *_bnd_elem_range)
1372 : {
1373 0 : const auto & [elem_ptr, elem_side, elem_bid] = *bnd_elem;
1374 0 : if (elem_bid == bid)
1375 : {
1376 : // If an element is internal to the group of subdomain, check the neighbor
1377 0 : if (blk_group.find(elem_ptr->subdomain_id()) != blk_group.end())
1378 : {
1379 0 : const auto * const neighbor = elem_ptr->neighbor_ptr(elem_side);
1380 :
1381 : // If we did not ghost the neighbor, we cannot decide
1382 0 : if (neighbor == libMesh::remote_elem)
1383 0 : mooseError("Insufficient level of geometrical ghosting to determine "
1384 : "if a boundary is internal to the mesh");
1385 : // If the neighbor does not exist, then we are on the edge of the mesh
1386 0 : if (!neighbor)
1387 0 : continue;
1388 : // If the neighbor is also in the group of subdomain,
1389 : // then the boundary cuts the subdomains
1390 0 : if (blk_group.find(neighbor->subdomain_id()) != blk_group.end())
1391 0 : return false;
1392 : }
1393 : }
1394 : }
1395 0 : return true;
1396 : }
1397 :
1398 : void
1399 155412 : MooseMesh::cacheInfo()
1400 : {
1401 466236 : TIME_SECTION("cacheInfo", 3);
1402 :
1403 155412 : _sub_to_data.clear();
1404 155412 : _neighbor_subdomain_boundary_ids.clear();
1405 155412 : _block_node_list.clear();
1406 155412 : _higher_d_elem_side_to_lower_d_elem.clear();
1407 155412 : _lower_d_elem_to_higher_d_elem_side.clear();
1408 155412 : _lower_d_interior_blocks.clear();
1409 155412 : _lower_d_boundary_blocks.clear();
1410 :
1411 155412 : const auto & mesh = getMesh();
1412 :
1413 : // Cache higher and lowerD element information
1414 38342458 : for (const auto & elem : mesh.element_ptr_range())
1415 : {
1416 38187046 : const Elem * ip_elem = elem->interior_parent();
1417 :
1418 38187046 : if (ip_elem)
1419 : {
1420 102485 : unsigned int ip_side = ip_elem->which_side_am_i(elem);
1421 :
1422 : // For some grid sequencing tests: ip_side == libMesh::invalid_uint
1423 102485 : if (ip_side != libMesh::invalid_uint)
1424 : {
1425 102325 : auto pair = std::make_pair(ip_elem, ip_side);
1426 102325 : _higher_d_elem_side_to_lower_d_elem.insert(
1427 102325 : std::pair<std::pair<const Elem *, unsigned short int>, const Elem *>(pair, elem));
1428 102325 : _lower_d_elem_to_higher_d_elem_side.insert(
1429 102325 : std::pair<const Elem *, unsigned short int>(elem, ip_side));
1430 :
1431 102325 : auto id = elem->subdomain_id();
1432 102325 : if (ip_elem->neighbor_ptr(ip_side))
1433 : {
1434 6680 : if (mesh.subdomain_name(id).find("INTERNAL_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos)
1435 6580 : _lower_d_interior_blocks.insert(id);
1436 : }
1437 : else
1438 : {
1439 95645 : if (mesh.subdomain_name(id).find("BOUNDARY_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos)
1440 6890 : _lower_d_boundary_blocks.insert(id);
1441 : }
1442 : }
1443 : }
1444 :
1445 244170802 : for (unsigned int nd = 0; nd < elem->n_nodes(); ++nd)
1446 : {
1447 205983756 : const Node & node = *elem->node_ptr(nd);
1448 205983756 : _block_node_list[node.id()].insert(elem->subdomain_id());
1449 : }
1450 155412 : }
1451 155412 : _communicator.set_union(_lower_d_interior_blocks);
1452 155412 : _communicator.set_union(_lower_d_boundary_blocks);
1453 :
1454 : // Cache the boundaries next to each subdomain
1455 27780138 : for (const auto & elem : mesh.active_local_element_ptr_range())
1456 : {
1457 27624726 : SubdomainID subdomain_id = elem->subdomain_id();
1458 27624726 : auto & sub_data = _sub_to_data[subdomain_id];
1459 27624726 : const auto elem_boundary_ids = getBoundaryIDs(elem);
1460 152028370 : for (unsigned int side = 0; side < elem->n_sides(); side++)
1461 : {
1462 124403644 : const auto & boundary_ids = elem_boundary_ids[side];
1463 124403644 : sub_data.boundary_ids.insert(boundary_ids.begin(), boundary_ids.end());
1464 :
1465 124403644 : const Elem * neig = elem->neighbor_ptr(side);
1466 124403644 : if (neig)
1467 : {
1468 117183227 : _neighbor_subdomain_boundary_ids[neig->subdomain_id()].insert(boundary_ids.begin(),
1469 : boundary_ids.end());
1470 117183227 : SubdomainID neighbor_subdomain_id = neig->subdomain_id();
1471 117183227 : if (neighbor_subdomain_id != subdomain_id)
1472 1820280 : sub_data.neighbor_subs.insert(neighbor_subdomain_id);
1473 : }
1474 : }
1475 27780138 : }
1476 :
1477 375280 : for (const auto blk_id : _mesh_subdomains)
1478 : {
1479 219868 : auto & sub_data = _sub_to_data[blk_id];
1480 219868 : _communicator.set_union(sub_data.neighbor_subs);
1481 219868 : _communicator.set_union(sub_data.boundary_ids);
1482 219868 : _communicator.set_union(_neighbor_subdomain_boundary_ids[blk_id]);
1483 : }
1484 155412 : }
1485 :
1486 : const std::set<SubdomainID> &
1487 96073733 : MooseMesh::getNodeBlockIds(const Node & node) const
1488 : {
1489 96073733 : auto it = _block_node_list.find(node.id());
1490 :
1491 96073733 : if (it == _block_node_list.end())
1492 0 : mooseError("Unable to find node: ", node.id(), " in any block list.");
1493 :
1494 192147466 : return it->second;
1495 : }
1496 :
1497 : MooseMesh::face_info_iterator
1498 164703 : MooseMesh::ownedFaceInfoBegin()
1499 : {
1500 : return face_info_iterator(
1501 164703 : _face_info.begin(),
1502 164703 : _face_info.end(),
1503 329406 : libMesh::Predicates::pid<std::vector<const FaceInfo *>::iterator>(this->processor_id()));
1504 : }
1505 :
1506 : MooseMesh::face_info_iterator
1507 164703 : MooseMesh::ownedFaceInfoEnd()
1508 : {
1509 : return face_info_iterator(
1510 164703 : _face_info.end(),
1511 164703 : _face_info.end(),
1512 329406 : libMesh::Predicates::pid<std::vector<const FaceInfo *>::iterator>(this->processor_id()));
1513 : }
1514 :
1515 : MooseMesh::elem_info_iterator
1516 84250 : MooseMesh::ownedElemInfoBegin()
1517 : {
1518 84250 : return elem_info_iterator(_elem_info.begin(),
1519 84250 : _elem_info.end(),
1520 168500 : Predicates::NotNull<std::vector<const ElemInfo *>::iterator>());
1521 : }
1522 :
1523 : MooseMesh::elem_info_iterator
1524 84250 : MooseMesh::ownedElemInfoEnd()
1525 : {
1526 84250 : return elem_info_iterator(_elem_info.end(),
1527 84250 : _elem_info.end(),
1528 168500 : Predicates::NotNull<std::vector<const ElemInfo *>::iterator>());
1529 : }
1530 :
1531 : // default begin() accessor
1532 : MooseMesh::bnd_node_iterator
1533 86802 : MooseMesh::bndNodesBegin()
1534 : {
1535 86802 : Predicates::NotNull<bnd_node_iterator_imp> p;
1536 173604 : return bnd_node_iterator(_bnd_nodes.begin(), _bnd_nodes.end(), p);
1537 86802 : }
1538 :
1539 : // default end() accessor
1540 : MooseMesh::bnd_node_iterator
1541 86802 : MooseMesh::bndNodesEnd()
1542 : {
1543 86802 : Predicates::NotNull<bnd_node_iterator_imp> p;
1544 173604 : return bnd_node_iterator(_bnd_nodes.end(), _bnd_nodes.end(), p);
1545 86802 : }
1546 :
1547 : // default begin() accessor
1548 : MooseMesh::bnd_elem_iterator
1549 84162 : MooseMesh::bndElemsBegin()
1550 : {
1551 84162 : Predicates::NotNull<bnd_elem_iterator_imp> p;
1552 168324 : return bnd_elem_iterator(_bnd_elems.begin(), _bnd_elems.end(), p);
1553 84162 : }
1554 :
1555 : // default end() accessor
1556 : MooseMesh::bnd_elem_iterator
1557 84162 : MooseMesh::bndElemsEnd()
1558 : {
1559 84162 : Predicates::NotNull<bnd_elem_iterator_imp> p;
1560 168324 : return bnd_elem_iterator(_bnd_elems.end(), _bnd_elems.end(), p);
1561 84162 : }
1562 :
1563 : const Node *
1564 0 : MooseMesh::addUniqueNode(const Point & p, Real tol)
1565 : {
1566 : /**
1567 : * Looping through the mesh nodes each time we add a point is very slow. To speed things
1568 : * up we keep a local data structure
1569 : */
1570 0 : if (getMesh().n_nodes() != _node_map.size())
1571 : {
1572 0 : _node_map.clear();
1573 0 : _node_map.reserve(getMesh().n_nodes());
1574 0 : for (const auto & node : getMesh().node_ptr_range())
1575 0 : _node_map.push_back(node);
1576 : }
1577 :
1578 0 : Node * node = nullptr;
1579 0 : for (unsigned int i = 0; i < _node_map.size(); ++i)
1580 : {
1581 0 : if (p.relative_fuzzy_equals(*_node_map[i], tol))
1582 : {
1583 0 : node = _node_map[i];
1584 0 : break;
1585 : }
1586 : }
1587 0 : if (node == nullptr)
1588 : {
1589 0 : node = getMesh().add_node(new Node(p));
1590 0 : _node_map.push_back(node);
1591 : }
1592 :
1593 : mooseAssert(node != nullptr, "Node is NULL");
1594 0 : return node;
1595 : }
1596 :
1597 : Node *
1598 5357 : MooseMesh::addQuadratureNode(const Elem * elem,
1599 : const unsigned short int side,
1600 : const unsigned int qp,
1601 : BoundaryID bid,
1602 : const Point & point)
1603 : {
1604 : Node * qnode;
1605 :
1606 5357 : if (_elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side].find(qp) ==
1607 10714 : _elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side].end())
1608 : {
1609 : // Create a new node id starting from the max node id and counting down. This will be the least
1610 : // likely to collide with an existing node id.
1611 : // Note that we are using numeric_limits<unsigned>::max even
1612 : // though max_id is stored as a dof_id_type. I tried this with
1613 : // numeric_limits<dof_id_type>::max and it broke several tests in
1614 : // MOOSE. So, this is some kind of a magic number that we will
1615 : // just continue to use...
1616 5357 : dof_id_type max_id = std::numeric_limits<unsigned int>::max() - 100;
1617 5357 : dof_id_type new_id = max_id - _quadrature_nodes.size();
1618 :
1619 5357 : if (new_id <= getMesh().max_node_id())
1620 0 : mooseError("Quadrature node id collides with existing node id!");
1621 :
1622 5357 : qnode = new Node(point, new_id);
1623 :
1624 : // Keep track of this new node in two different ways for easy lookup
1625 5357 : _quadrature_nodes[new_id] = qnode;
1626 5357 : _elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side][qp] = qnode;
1627 :
1628 5357 : if (elem->active())
1629 5357 : internalNodeToElemMap()[new_id].push_back(elem->id());
1630 : }
1631 : else
1632 0 : qnode = _elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side][qp];
1633 :
1634 5357 : BndNode * bnode = new BndNode(qnode, bid);
1635 5357 : _bnd_nodes.push_back(bnode);
1636 5357 : _bnd_node_ids[bid].insert(qnode->id());
1637 :
1638 5357 : _extra_bnd_nodes.push_back(*bnode);
1639 :
1640 : // Do this so the range will be regenerated next time it is accessed
1641 5357 : _bnd_node_range.reset();
1642 :
1643 5357 : return qnode;
1644 : }
1645 :
1646 : Node *
1647 137880 : MooseMesh::getQuadratureNode(const Elem * elem,
1648 : const unsigned short int side,
1649 : const unsigned int qp)
1650 : {
1651 : mooseAssert(_elem_to_side_to_qp_to_quadrature_nodes.find(elem->id()) !=
1652 : _elem_to_side_to_qp_to_quadrature_nodes.end(),
1653 : "Elem has no quadrature nodes!");
1654 : mooseAssert(_elem_to_side_to_qp_to_quadrature_nodes[elem->id()].find(side) !=
1655 : _elem_to_side_to_qp_to_quadrature_nodes[elem->id()].end(),
1656 : "Side has no quadrature nodes!");
1657 : mooseAssert(_elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side].find(qp) !=
1658 : _elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side].end(),
1659 : "qp not found on side!");
1660 :
1661 137880 : return _elem_to_side_to_qp_to_quadrature_nodes[elem->id()][side][qp];
1662 : }
1663 :
1664 : void
1665 74400 : MooseMesh::clearQuadratureNodes()
1666 : {
1667 : // Delete all the quadrature nodes
1668 79745 : for (auto & it : _quadrature_nodes)
1669 5345 : delete it.second;
1670 :
1671 74400 : _quadrature_nodes.clear();
1672 74400 : _elem_to_side_to_qp_to_quadrature_nodes.clear();
1673 74400 : _extra_bnd_nodes.clear();
1674 :
1675 : // NOTE: this does not clear them from the nodeToElem map
1676 74400 : }
1677 :
1678 : BoundaryID
1679 543878 : MooseMesh::getBoundaryID(const BoundaryName & boundary_name) const
1680 : {
1681 543878 : if (boundary_name == "ANY_BOUNDARY_ID")
1682 0 : mooseError("Please use getBoundaryIDs() when passing \"ANY_BOUNDARY_ID\"");
1683 :
1684 543878 : return MooseMeshUtils::getBoundaryID(boundary_name, getMesh());
1685 : }
1686 :
1687 : const Elem *
1688 1577895527 : MooseMesh::getLowerDElem(const Elem * elem, unsigned short int side) const
1689 : {
1690 1577895527 : auto it = _higher_d_elem_side_to_lower_d_elem.find(std::make_pair(elem, side));
1691 :
1692 1577895527 : if (it != _higher_d_elem_side_to_lower_d_elem.end())
1693 280705 : return it->second;
1694 : else
1695 1577614822 : return nullptr;
1696 : }
1697 :
1698 : unsigned int
1699 260 : MooseMesh::getHigherDSide(const Elem * elem) const
1700 : {
1701 260 : auto it = _lower_d_elem_to_higher_d_elem_side.find(elem);
1702 :
1703 260 : if (it != _lower_d_elem_to_higher_d_elem_side.end())
1704 260 : return it->second;
1705 : else
1706 0 : return libMesh::invalid_uint;
1707 : }
1708 :
1709 : std::vector<BoundaryID>
1710 117217 : MooseMesh::getBoundaryIDs(const std::vector<BoundaryName> & boundary_name,
1711 : bool generate_unknown) const
1712 : {
1713 : return MooseMeshUtils::getBoundaryIDs(
1714 117217 : getMesh(), boundary_name, generate_unknown, _mesh_boundary_ids);
1715 : }
1716 :
1717 : SubdomainID
1718 529217 : MooseMesh::getSubdomainID(const SubdomainName & subdomain_name) const
1719 : {
1720 529217 : return MooseMeshUtils::getSubdomainID(subdomain_name, getMesh());
1721 : }
1722 :
1723 : std::vector<SubdomainID>
1724 240364 : MooseMesh::getSubdomainIDs(const std::vector<SubdomainName> & subdomain_name) const
1725 : {
1726 240364 : return MooseMeshUtils::getSubdomainIDs(getMesh(), subdomain_name);
1727 : }
1728 :
1729 : std::set<SubdomainID>
1730 0 : MooseMesh::getSubdomainIDs(const std::set<SubdomainName> & subdomain_name) const
1731 : {
1732 0 : return MooseMeshUtils::getSubdomainIDs(getMesh(), subdomain_name);
1733 : }
1734 :
1735 : void
1736 253 : MooseMesh::setSubdomainName(SubdomainID subdomain_id, const SubdomainName & name)
1737 : {
1738 : mooseAssert(name != "ANY_BLOCK_ID", "Cannot set subdomain name to 'ANY_BLOCK_ID'");
1739 253 : getMesh().subdomain_name(subdomain_id) = name;
1740 253 : }
1741 :
1742 : void
1743 0 : MooseMesh::setSubdomainName(MeshBase & mesh, SubdomainID subdomain_id, const SubdomainName & name)
1744 : {
1745 : mooseAssert(name != "ANY_BLOCK_ID", "Cannot set subdomain name to 'ANY_BLOCK_ID'");
1746 0 : mesh.subdomain_name(subdomain_id) = name;
1747 0 : }
1748 :
1749 : const std::string &
1750 4313676 : MooseMesh::getSubdomainName(SubdomainID subdomain_id) const
1751 : {
1752 4313676 : return getMesh().subdomain_name(subdomain_id);
1753 : }
1754 :
1755 : std::vector<SubdomainName>
1756 71 : MooseMesh::getSubdomainNames(const std::vector<SubdomainID> & subdomain_ids) const
1757 : {
1758 71 : std::vector<SubdomainName> names(subdomain_ids.size());
1759 :
1760 142 : for (unsigned int i = 0; i < subdomain_ids.size(); i++)
1761 71 : names[i] = getSubdomainName(subdomain_ids[i]);
1762 :
1763 71 : return names;
1764 0 : }
1765 :
1766 : void
1767 110 : MooseMesh::setBoundaryName(BoundaryID boundary_id, BoundaryName name)
1768 : {
1769 110 : BoundaryInfo & boundary_info = getMesh().get_boundary_info();
1770 :
1771 : // We need to figure out if this boundary is a sideset or nodeset
1772 110 : if (boundary_info.get_side_boundary_ids().count(boundary_id))
1773 30 : boundary_info.sideset_name(boundary_id) = name;
1774 : else
1775 80 : boundary_info.nodeset_name(boundary_id) = name;
1776 110 : }
1777 :
1778 : const std::string &
1779 7607654 : MooseMesh::getBoundaryName(const BoundaryID boundary_id) const
1780 : {
1781 7607654 : const BoundaryInfo & boundary_info = getMesh().get_boundary_info();
1782 :
1783 : // We need to figure out if this boundary is a sideset or nodeset
1784 7607654 : if (boundary_info.get_side_boundary_ids().count(boundary_id))
1785 7487284 : return boundary_info.get_sideset_name(boundary_id);
1786 : else
1787 120370 : return boundary_info.get_nodeset_name(boundary_id);
1788 : }
1789 :
1790 : std::string
1791 27 : MooseMesh::getBoundaryString(BoundaryID boundary_id) const
1792 : {
1793 27 : const auto name = getBoundaryName(boundary_id);
1794 54 : return name.size() ? name : std::to_string(boundary_id);
1795 27 : }
1796 :
1797 : // specialization for PointListAdaptor<MooseMesh::PeriodicNodeInfo>
1798 : template <>
1799 : inline const Point &
1800 173430 : PointListAdaptor<MooseMesh::PeriodicNodeInfo>::getPoint(
1801 : const MooseMesh::PeriodicNodeInfo & item) const
1802 : {
1803 173430 : return *(item.first);
1804 : }
1805 :
1806 : void
1807 27 : MooseMesh::buildPeriodicNodeMap(std::multimap<dof_id_type, dof_id_type> & periodic_node_map,
1808 : unsigned int var_number,
1809 : libMesh::PeriodicBoundaries * pbs) const
1810 : {
1811 81 : TIME_SECTION("buildPeriodicNodeMap", 5);
1812 :
1813 : // clear existing map
1814 27 : periodic_node_map.clear();
1815 :
1816 : // get periodic nodes
1817 27 : std::vector<PeriodicNodeInfo> periodic_nodes;
1818 1575 : for (const auto & t : getMesh().get_boundary_info().build_node_list())
1819 : {
1820 : // unfortunately libMesh does not give us a pointer, so we have to look it up ourselves
1821 1548 : auto node = _mesh->node_ptr(std::get<0>(t));
1822 : mooseAssert(node != nullptr,
1823 : "libMesh::BoundaryInfo::build_node_list() returned an ID for a non-existing node");
1824 1548 : auto bc_id = std::get<1>(t);
1825 1548 : periodic_nodes.emplace_back(node, bc_id);
1826 27 : }
1827 :
1828 : // sort by boundary id
1829 27 : std::sort(periodic_nodes.begin(),
1830 : periodic_nodes.end(),
1831 8658 : [](const PeriodicNodeInfo & a, const PeriodicNodeInfo & b) -> bool
1832 8658 : { return a.second > b.second; });
1833 :
1834 : // build kd-tree
1835 : using KDTreeType = nanoflann::KDTreeSingleIndexAdaptor<
1836 : nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<PeriodicNodeInfo>, Real, std::size_t>,
1837 : PointListAdaptor<PeriodicNodeInfo>,
1838 : LIBMESH_DIM,
1839 : std::size_t>;
1840 27 : const unsigned int max_leaf_size = 20; // slightly affects runtime
1841 : auto point_list =
1842 27 : PointListAdaptor<PeriodicNodeInfo>(periodic_nodes.begin(), periodic_nodes.end());
1843 : auto kd_tree = std::make_unique<KDTreeType>(
1844 27 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
1845 : mooseAssert(kd_tree != nullptr, "KDTree was not properly initialized.");
1846 27 : kd_tree->buildIndex();
1847 :
1848 : // data structures for kd-tree search
1849 27 : nanoflann::SearchParameters search_params;
1850 27 : std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
1851 :
1852 : // iterate over periodic nodes (boundary ids are in contiguous blocks)
1853 27 : libMesh::PeriodicBoundaryBase * periodic = nullptr;
1854 27 : BoundaryID current_bc_id = BoundaryInfo::invalid_id;
1855 1575 : for (auto & pair : periodic_nodes)
1856 : {
1857 : // entering a new block of boundary IDs
1858 1548 : if (pair.second != current_bc_id)
1859 : {
1860 108 : current_bc_id = pair.second;
1861 108 : periodic = pbs->boundary(current_bc_id);
1862 108 : if (periodic && !periodic->is_my_variable(var_number))
1863 0 : periodic = nullptr;
1864 : }
1865 :
1866 : // variable is not periodic at this node, skip
1867 1548 : if (!periodic)
1868 0 : continue;
1869 :
1870 : // clear result buffer
1871 1548 : ret_matches.clear();
1872 :
1873 : // id of the current node
1874 1548 : const auto id = pair.first->id();
1875 :
1876 : // position where we expect a periodic partner for the current node and boundary
1877 1548 : Point search_point = periodic->get_corresponding_pos(*pair.first);
1878 :
1879 : // search at the expected point
1880 1548 : kd_tree->radiusSearch(&(search_point)(0), libMesh::TOLERANCE, ret_matches, search_params);
1881 4248 : for (auto & match_pair : ret_matches)
1882 : {
1883 2700 : const auto & match = periodic_nodes[match_pair.first];
1884 : // add matched node if the boundary id is the corresponding id in the periodic pair
1885 2700 : if (match.second == periodic->pairedboundary)
1886 1548 : periodic_node_map.emplace(id, match.first->id());
1887 : }
1888 : }
1889 27 : }
1890 :
1891 : void
1892 0 : MooseMesh::buildPeriodicNodeSets(std::map<BoundaryID, std::set<dof_id_type>> & periodic_node_sets,
1893 : unsigned int var_number,
1894 : libMesh::PeriodicBoundaries * pbs) const
1895 : {
1896 0 : TIME_SECTION("buildPeriodicNodeSets", 5);
1897 :
1898 0 : periodic_node_sets.clear();
1899 :
1900 : // Loop over all the boundary nodes adding the periodic nodes to the appropriate set
1901 0 : for (const auto & t : getMesh().get_boundary_info().build_node_list())
1902 : {
1903 0 : auto node_id = std::get<0>(t);
1904 0 : auto bc_id = std::get<1>(t);
1905 :
1906 : // Is this current node on a known periodic boundary?
1907 0 : if (periodic_node_sets.find(bc_id) != periodic_node_sets.end())
1908 0 : periodic_node_sets[bc_id].insert(node_id);
1909 : else // This still might be a periodic node but we just haven't seen this boundary_id yet
1910 : {
1911 0 : const libMesh::PeriodicBoundaryBase * periodic = pbs->boundary(bc_id);
1912 0 : if (periodic && periodic->is_my_variable(var_number))
1913 0 : periodic_node_sets[bc_id].insert(node_id);
1914 : }
1915 0 : }
1916 0 : }
1917 :
1918 : bool
1919 68871 : MooseMesh::detectOrthogonalDimRanges(Real tol)
1920 : {
1921 206613 : TIME_SECTION("detectOrthogonalDimRanges", 5);
1922 :
1923 68871 : if (_regular_orthogonal_mesh)
1924 34541 : return true;
1925 :
1926 34330 : std::vector<Real> min(3, std::numeric_limits<Real>::max());
1927 34330 : std::vector<Real> max(3, std::numeric_limits<Real>::min());
1928 34330 : unsigned int dim = getMesh().mesh_dimension();
1929 :
1930 : // Find the bounding box of our mesh
1931 10260478 : for (const auto & node : getMesh().node_ptr_range())
1932 : // Check all coordinates, we don't know if this mesh might be lying in a higher dim even if the
1933 : // mesh dimension is lower.
1934 40904592 : for (const auto i : make_range(Moose::dim))
1935 : {
1936 30678444 : if ((*node)(i) < min[i])
1937 210285 : min[i] = (*node)(i);
1938 30678444 : if ((*node)(i) > max[i])
1939 469737 : max[i] = (*node)(i);
1940 34330 : }
1941 :
1942 34330 : this->comm().max(max);
1943 34330 : this->comm().min(min);
1944 :
1945 34330 : _extreme_nodes.resize(8); // 2^LIBMESH_DIM
1946 : // Now make sure that there are actual nodes at all of the extremes
1947 34330 : std::vector<bool> extreme_matches(8, false);
1948 34330 : std::vector<unsigned int> comp_map(3);
1949 10260478 : for (const auto & node : getMesh().node_ptr_range())
1950 : {
1951 : // See if the current node is located at one of the extremes
1952 10226148 : unsigned int coord_match = 0;
1953 :
1954 40904592 : for (const auto i : make_range(Moose::dim))
1955 : {
1956 30678444 : if (std::abs((*node)(i)-min[i]) < tol)
1957 : {
1958 6353225 : comp_map[i] = MIN;
1959 6353225 : ++coord_match;
1960 : }
1961 24325219 : else if (std::abs((*node)(i)-max[i]) < tol)
1962 : {
1963 1394759 : comp_map[i] = MAX;
1964 1394759 : ++coord_match;
1965 : }
1966 : }
1967 :
1968 10226148 : if (coord_match == LIBMESH_DIM) // Found a coordinate at one of the extremes
1969 : {
1970 124889 : _extreme_nodes[comp_map[X] * 4 + comp_map[Y] * 2 + comp_map[Z]] = node;
1971 124889 : extreme_matches[comp_map[X] * 4 + comp_map[Y] * 2 + comp_map[Z]] = true;
1972 : }
1973 34330 : }
1974 :
1975 : // See if we matched all of the extremes for the mesh dimension
1976 34330 : this->comm().max(extreme_matches);
1977 34330 : if (std::count(extreme_matches.begin(), extreme_matches.end(), true) == (1 << dim))
1978 29939 : _regular_orthogonal_mesh = true;
1979 :
1980 : // Set the bounds
1981 34330 : _bounds.resize(LIBMESH_DIM);
1982 137320 : for (const auto i : make_range(Moose::dim))
1983 : {
1984 102990 : _bounds[i].resize(2);
1985 102990 : _bounds[i][MIN] = min[i];
1986 102990 : _bounds[i][MAX] = max[i];
1987 : }
1988 :
1989 34330 : return _regular_orthogonal_mesh;
1990 68871 : }
1991 :
1992 : void
1993 538 : MooseMesh::detectPairedSidesets()
1994 : {
1995 1614 : TIME_SECTION("detectPairedSidesets", 5);
1996 :
1997 538 : _paired_boundary = std::vector<std::pair<BoundaryID, BoundaryID>>();
1998 :
1999 : // Loop over level-0 elements (since boundary condition information
2000 : // is only directly stored for them) and find sidesets with normals
2001 : // that point in the -x, +x, -y, +y, and -z, +z direction. If there
2002 : // is a unique sideset id for each direction, then the paired
2003 : // sidesets consist of (-x,+x), (-y,+y), (-z,+z). If there are
2004 : // multiple sideset ids for a given direction, then we can't pick a
2005 : // single pair for that direction. In that case, we'll just return
2006 : // as was done in the original algorithm.
2007 :
2008 : // we need to test all element dimensions from dim down to 1
2009 538 : const unsigned int mesh_dim = getMesh().mesh_dimension();
2010 :
2011 : // Helper for iterating through unit dimensions (0=x, 1=y, 2=z)
2012 : static constexpr std::array<std::size_t, 3> unit_dims{0, 1, 2};
2013 : // Helper for mapping from unit dim -> name
2014 1614 : static const std::array<std::string, 3> unit_dim_names{"x", "y", "z"};
2015 :
2016 : // Boundary id sets for elements of different dimensions
2017 : // First index: side dimension; 0=1D, 1=2D, 2=3D
2018 : // Second index: unit dimension; 0=x, 1=y, 2=z
2019 : // Third index: false for minus, true for plus
2020 16678 : std::array<std::array<std::array<std::set<BoundaryID>, 2>, 3>, 3> ids{};
2021 :
2022 : // Build quadrature needed to evaluate side normals
2023 538 : std::array<std::unique_ptr<FEBase>, 3> fe_faces{};
2024 538 : std::array<std::unique_ptr<libMesh::QGauss>, 3> qfaces{};
2025 1596 : for (const auto side_dim : make_range(mesh_dim))
2026 : {
2027 : // Face is assumed to be flat, therefore normal is assumed to be
2028 : // constant over the face, therefore only compute it at 1 qp.
2029 1058 : qfaces[side_dim] = std::unique_ptr<libMesh::QGauss>(new libMesh::QGauss(side_dim, CONSTANT));
2030 :
2031 : // A first-order Lagrange FE for the face.
2032 1058 : fe_faces[side_dim] = FEBase::build(side_dim + 1, FEType(FIRST, libMesh::LAGRANGE));
2033 1058 : fe_faces[side_dim]->attach_quadrature_rule(qfaces[side_dim].get());
2034 1058 : fe_faces[side_dim]->get_normals();
2035 : }
2036 :
2037 : // Get boundary IDs for each dimension that are in the unit normal
2038 538 : const auto & boundary_info = getMesh().get_boundary_info();
2039 : // Temporary for evaluating boundary_ids
2040 538 : std::vector<boundary_id_type> face_ids;
2041 : // The side dimensions we've come across, so that we only report
2042 : // warnings for side dimensions that we have
2043 538 : std::set<unsigned int> side_dims;
2044 : // Normal dimensions that we found that were nonzero; lets us
2045 : // skip warnings for dimensions that we don't have
2046 538 : std::array<bool, 3> nonzero_dims = periodic_dim_default;
2047 307171 : for (auto & elem : as_range(getMesh().level_elements_begin(0), getMesh().level_elements_end(0)))
2048 : {
2049 : // If not on the boundary, nothing to do
2050 306633 : if (!elem->on_boundary())
2051 259116 : continue;
2052 :
2053 47517 : const auto side_dim = elem->dim() - 1;
2054 47517 : side_dims.insert(side_dim);
2055 :
2056 : // Check for unit normals on each boundary side
2057 265346 : for (const auto s : elem->side_index_range())
2058 217829 : if (!elem->neighbor_ptr(s))
2059 : {
2060 : // Reinit to get the normal
2061 52932 : fe_faces[side_dim]->reinit(elem, s);
2062 52932 : const auto & normal = fe_faces[side_dim]->get_normals()[0];
2063 :
2064 : // Get the boundary ID(s) for this side. If there is more
2065 : // than 1 boundary id, then we already can't determine a
2066 : // unique pairing of sides in this direction, but we'll just
2067 : // keep going to keep the logic simple.
2068 52932 : boundary_info.boundary_ids(elem, s, face_ids);
2069 :
2070 52932 : bool found = false;
2071 211728 : for (const auto unit_dim : unit_dims)
2072 : {
2073 158796 : if (libMesh::absolute_fuzzy_equals(normal(unit_dim), 0.0))
2074 103964 : continue;
2075 54832 : nonzero_dims[unit_dim] = true;
2076 54832 : if (!found)
2077 88072 : for (const auto plus : {false, true})
2078 : {
2079 84272 : if (libMesh::absolute_fuzzy_equals(normal(unit_dim), plus ? 1.0 : -1.0))
2080 : {
2081 51032 : ids[side_dim][unit_dim][plus].insert(face_ids.begin(), face_ids.end());
2082 51032 : found = true;
2083 51032 : break;
2084 : }
2085 : }
2086 : }
2087 : }
2088 538 : }
2089 :
2090 : // For a distributed mesh, boundaries may be distributed as well. We therefore collect information
2091 : // from everyone. If the mesh is already serial, then there is no need to do an allgather. Note
2092 : // that this is just going to gather information about what the periodic bc ids are. We are not
2093 : // gathering any remote elements or anything like that. It's just that the GhostPointNeighbors
2094 : // ghosting functor currently relies on the fact that every process agrees on whether we have
2095 : // periodic boundaries; every process that thinks there are periodic boundaries will call
2096 : // MeshBase::sub_point_locator which makes a parallel_object_only() assertion (right or wrong). So
2097 : // we all need to go there (or not go there)
2098 538 : if (_use_distributed_mesh && !_mesh->is_serial())
2099 : {
2100 : // Communicate id data by packing as [side dim, unit dim, plus (as a char), boundary id]
2101 122 : std::vector<std::tuple<unsigned int, unsigned int, unsigned char, boundary_id_type>> id_data;
2102 244 : for (const auto side_dim : side_dims)
2103 488 : for (const auto unit_dim : unit_dims)
2104 1098 : for (const auto plus : {false, true})
2105 1192 : for (const auto bd : ids[side_dim][unit_dim][plus])
2106 460 : id_data.emplace_back(side_dim, unit_dim, plus, bd);
2107 122 : _communicator.allgather(id_data, false);
2108 1170 : for (const auto & [side_dim, unit_dim, plus_char, bd] : id_data)
2109 1048 : ids[side_dim][unit_dim][bool(plus_char)].insert(bd);
2110 :
2111 : // Gather true-ness of nonzero_dims
2112 488 : for (auto & entry : nonzero_dims)
2113 366 : _communicator.max(entry);
2114 :
2115 : // Gather found side dimensions
2116 122 : _communicator.set_union(side_dims);
2117 122 : } // end if (_use_distributed_mesh && !_need_ghost_ghosted_boundaries)
2118 :
2119 : // Find pairings that have exactly one boundary on each side
2120 538 : std::ostringstream oss_found, oss_missing;
2121 1076 : for (const auto side_dim : side_dims)
2122 : {
2123 2152 : for (const auto unit_dim : unit_dims)
2124 1614 : if (nonzero_dims[unit_dim])
2125 : {
2126 1061 : const auto & unit_name = unit_dim_names[unit_dim];
2127 1061 : const auto & minus = ids[side_dim][unit_dim][false];
2128 1061 : const auto & plus = ids[side_dim][unit_dim][true];
2129 :
2130 1061 : if (minus.size() == 1 && plus.size() == 1)
2131 : {
2132 1904 : const auto get_boundary_name = [this](const auto id)
2133 : {
2134 1904 : const auto & name = getBoundaryName(id);
2135 1904 : return name.size() ? name : std::to_string(id);
2136 952 : };
2137 :
2138 952 : oss_found << "\n " << side_dim + 1 << "D " << unit_name
2139 952 : << "-direction: " << get_boundary_name(*minus.begin()) << " <-> "
2140 1904 : << get_boundary_name(*plus.begin());
2141 952 : _paired_boundary->emplace_back(std::make_pair(*minus.begin(), *plus.begin()));
2142 : }
2143 : else
2144 109 : oss_missing << "\n " << side_dim + 1 << "D -" << unit_name << "/+" << unit_name
2145 109 : << ": Found " << minus.size() << " -" << unit_name << " boundaries and "
2146 109 : << plus.size() << " +" << unit_name << " boundaries";
2147 : }
2148 : }
2149 :
2150 538 : std::ostringstream oss;
2151 538 : const auto found = oss_found.str();
2152 538 : const auto missing = oss_missing.str();
2153 538 : if (found.size())
2154 : oss << "The following paired boundaries were automatically detected for periodicity:\n"
2155 504 : << found << "\n";
2156 538 : if (missing.size())
2157 : {
2158 75 : if (found.size())
2159 41 : oss << "\n";
2160 : oss << "Paired boundaries were not automatically detected for the following:\n"
2161 : << missing
2162 : << "\n\nAutomatic detection requires that exactly one boundary is found in each unit "
2163 75 : "direction.\n";
2164 : }
2165 :
2166 538 : mooseInfoRepeated(oss.str());
2167 538 : }
2168 :
2169 : Real
2170 72942 : MooseMesh::dimensionWidth(unsigned int component) const
2171 : {
2172 72942 : return getMaxInDimension(component) - getMinInDimension(component);
2173 : }
2174 :
2175 : Real
2176 33427 : MooseMesh::getMinInDimension(unsigned int component) const
2177 : {
2178 : mooseAssert(_mesh, "The MeshBase has not been constructed");
2179 : mooseAssert(component < _bounds.size(), "Requested dimension out of bounds");
2180 :
2181 33427 : return _bounds[component][MIN];
2182 : }
2183 :
2184 : Real
2185 33427 : MooseMesh::getMaxInDimension(unsigned int component) const
2186 : {
2187 : mooseAssert(_mesh, "The MeshBase has not been constructed");
2188 : mooseAssert(component < _bounds.size(), "Requested dimension out of bounds");
2189 :
2190 33427 : return _bounds[component][MAX];
2191 : }
2192 :
2193 : void
2194 873 : MooseMesh::addPeriodicVariable(const unsigned int sys_num,
2195 : const unsigned int var_num,
2196 : const BoundaryID primary,
2197 : const BoundaryID secondary)
2198 : {
2199 873 : if (!_regular_orthogonal_mesh)
2200 0 : return;
2201 :
2202 873 : const auto key = std::make_pair(sys_num, var_num);
2203 873 : auto & entry = _periodic_dim.try_emplace(key, periodic_dim_default).first->second;
2204 :
2205 873 : _half_range = Point(dimensionWidth(0) / 2.0, dimensionWidth(1) / 2.0, dimensionWidth(2) / 2.0);
2206 :
2207 873 : bool component_found = false;
2208 2707 : for (const auto component : make_range(dimension()))
2209 : {
2210 1834 : const std::pair<BoundaryID, BoundaryID> * boundary_ids = getPairedBoundaryMapping(component);
2211 :
2212 1834 : if (boundary_ids && ((boundary_ids->first == primary && boundary_ids->second == secondary) ||
2213 976 : (boundary_ids->first == secondary && boundary_ids->second == primary)))
2214 : {
2215 858 : entry[component] = true;
2216 858 : component_found = true;
2217 : }
2218 : }
2219 :
2220 873 : if (!component_found)
2221 30 : mooseWarning("Could not find a match between boundary '",
2222 15 : getBoundaryName(primary),
2223 : "' and '",
2224 15 : getBoundaryName(secondary),
2225 : "' to set periodic boundary conditions for variable (index:",
2226 : var_num,
2227 : ") in either the X, Y or Z direction. The periodic dimension of the mesh for this "
2228 : "variable will not be stored.");
2229 : }
2230 :
2231 : const std::array<bool, 3> &
2232 5074751 : MooseMesh::queryPeriodicDimensions(const unsigned int sys_num, const unsigned int var_num) const
2233 : {
2234 5074751 : const auto key = std::make_pair(sys_num, var_num);
2235 5074751 : if (const auto it = _periodic_dim.find(key); it != _periodic_dim.end())
2236 4971067 : return it->second;
2237 103684 : return periodic_dim_default;
2238 : }
2239 :
2240 : const std::array<bool, 3> &
2241 27 : MooseMesh::queryPeriodicDimensions(const MooseVariableBase & var) const
2242 : {
2243 27 : return queryPeriodicDimensions(var.sys().number(), var.number());
2244 : }
2245 :
2246 : bool
2247 0 : MooseMesh::isTranslatedPeriodic(const unsigned int sys_num,
2248 : const unsigned int var_num,
2249 : const unsigned int component) const
2250 : {
2251 : mooseAssert(component < dimension(), "Requested dimension out of bounds");
2252 0 : return queryPeriodicDimensions(sys_num, var_num)[component];
2253 : }
2254 :
2255 : bool
2256 0 : MooseMesh::isTranslatedPeriodic(const MooseVariableBase & var, const unsigned int component) const
2257 : {
2258 0 : return isTranslatedPeriodic(var.sys().number(), var.number(), component);
2259 : }
2260 :
2261 : bool
2262 0 : MooseMesh::isTranslatedPeriodic(const unsigned int var_num, const unsigned int component) const
2263 : {
2264 0 : mooseDoOnce(mooseDeprecated(
2265 : "MooseMesh::isTranslatedPeriodic(const unsigned int, const unsigned int) is deprecated. Use "
2266 : "the method that additionally takes the system number or the MooseVariableBase instead."));
2267 0 : return isTranslatedPeriodic(0, var_num, component);
2268 : }
2269 :
2270 : RealVectorValue
2271 5074724 : MooseMesh::minPeriodicVector(const unsigned int sys_num,
2272 : const unsigned int var_num,
2273 : Point p,
2274 : Point q) const
2275 : {
2276 5074724 : const auto & periodic_dims = queryPeriodicDimensions(sys_num, var_num);
2277 :
2278 15159128 : for (const auto i : make_range(dimension()))
2279 : {
2280 : // check to see if we're closer in real or periodic space in x, y, and z
2281 10084404 : if (periodic_dims[i])
2282 : {
2283 : // Need to test order before differencing
2284 9877036 : if (p(i) > q(i))
2285 : {
2286 6390026 : if (p(i) - q(i) > _half_range(i))
2287 2344164 : p(i) -= _half_range(i) * 2;
2288 : }
2289 : else
2290 : {
2291 3487010 : if (q(i) - p(i) > _half_range(i))
2292 926262 : p(i) += _half_range(i) * 2;
2293 : }
2294 : }
2295 : }
2296 :
2297 5074724 : return q - p;
2298 : }
2299 :
2300 : RealVectorValue
2301 0 : MooseMesh::minPeriodicVector(const MooseVariableBase & var, const Point & p, const Point & q) const
2302 : {
2303 0 : return minPeriodicVector(var.sys().number(), var.number(), p, q);
2304 : }
2305 :
2306 : RealVectorValue
2307 0 : MooseMesh::minPeriodicVector(const unsigned int var_num, const Point & p, const Point & q) const
2308 : {
2309 0 : mooseDoOnce(mooseDeprecated("MooseMesh::minPeriodicVector(const unsigned int, const Point &, "
2310 : "const Point &) is deprecated. Use the method that additionally "
2311 : "takes the system number or the MooseVariableBase instead."));
2312 0 : return minPeriodicVector(0, var_num, p, q);
2313 : }
2314 :
2315 : Real
2316 5074724 : MooseMesh::minPeriodicDistance(const unsigned int sys_num,
2317 : const unsigned int var_num,
2318 : const Point & p,
2319 : const Point & q) const
2320 : {
2321 5074724 : return minPeriodicVector(sys_num, var_num, p, q).norm();
2322 : }
2323 :
2324 : Real
2325 25600 : MooseMesh::minPeriodicDistance(const MooseVariableBase & var,
2326 : const Point & p,
2327 : const Point & q) const
2328 : {
2329 25600 : return minPeriodicDistance(var.sys().number(), var.number(), p, q);
2330 : }
2331 :
2332 : Real
2333 0 : MooseMesh::minPeriodicDistance(const unsigned int var_num, const Point & p, const Point & q) const
2334 : {
2335 0 : mooseDoOnce(mooseDeprecated("MooseMesh::minPeriodicDistance(const unsigned int, const Point &, "
2336 : "const Point &) is deprecated. Use the method that additionally "
2337 : "takes the system number or the MooseVariableBase instead."));
2338 0 : return minPeriodicDistance(0, var_num, p, q);
2339 : }
2340 :
2341 : const std::pair<BoundaryID, BoundaryID> *
2342 2383 : MooseMesh::getPairedBoundaryMapping(unsigned int component) const
2343 : {
2344 2383 : if (!_regular_orthogonal_mesh)
2345 0 : mooseError("Trying to retrieve automatic paired mapping for a mesh that is not regular and "
2346 : "orthogonal");
2347 :
2348 : mooseAssert(component < dimension(), "Requested dimension out of bounds");
2349 :
2350 2383 : if (!hasDetectedPairedSidesets())
2351 0 : mooseError("MooseMesh::getPairedBoundaryMapping(): Paired boundaries not built; must call "
2352 : "detectPairedSidesets() first");
2353 :
2354 2383 : if (component < _paired_boundary->size())
2355 2380 : return &(*_paired_boundary)[component];
2356 : else
2357 3 : return nullptr;
2358 : }
2359 :
2360 : void
2361 33 : MooseMesh::buildHRefinementAndCoarseningMaps(Assembly * const assembly)
2362 : {
2363 33 : std::map<ElemType, Elem *> canonical_elems;
2364 :
2365 : // First, loop over all elements and find a canonical element for each type
2366 : // Doing it this way guarantees that this is going to work in parallel
2367 19937 : for (const auto & elem : getMesh().element_ptr_range()) // TODO: Thread this
2368 : {
2369 9952 : ElemType type = elem->type();
2370 :
2371 9952 : if (canonical_elems.find(type) ==
2372 19904 : canonical_elems.end()) // If we haven't seen this type of elem before save it
2373 42 : canonical_elems[type] = elem;
2374 : else
2375 : {
2376 9910 : Elem * stored = canonical_elems[type];
2377 9910 : if (elem->id() < stored->id()) // Arbitrarily keep the one with a lower id
2378 0 : canonical_elems[type] = elem;
2379 : }
2380 33 : }
2381 : // Now build the maps using these templates
2382 : // Note: This MUST be done NOT threaded!
2383 75 : for (const auto & can_it : canonical_elems)
2384 : {
2385 42 : Elem * elem = can_it.second;
2386 :
2387 : // Need to do this just once to get the right qrules put in place
2388 42 : assembly->setCurrentSubdomainID(elem->subdomain_id());
2389 42 : assembly->reinit(elem);
2390 42 : assembly->reinit(elem, 0);
2391 42 : auto && qrule = assembly->writeableQRule();
2392 42 : auto && qrule_face = assembly->writeableQRuleFace();
2393 :
2394 : // Volume to volume projection for refinement
2395 42 : buildRefinementMap(*elem, *qrule, *qrule_face, -1, -1, -1);
2396 :
2397 : // Volume to volume projection for coarsening
2398 42 : buildCoarseningMap(*elem, *qrule, *qrule_face, -1);
2399 :
2400 : // Map the sides of children
2401 216 : for (unsigned int side = 0; side < elem->n_sides(); side++)
2402 : {
2403 : // Side to side for sides that match parent's sides
2404 174 : buildRefinementMap(*elem, *qrule, *qrule_face, side, -1, side);
2405 174 : buildCoarseningMap(*elem, *qrule, *qrule_face, side);
2406 : }
2407 :
2408 : // Child side to parent volume mapping for "internal" child sides
2409 240 : for (unsigned int child = 0; child < elem->n_children(); ++child)
2410 1146 : for (unsigned int side = 0; side < elem->n_sides();
2411 : ++side) // Assume children have the same number of sides!
2412 948 : if (!elem->is_child_on_side(child, side)) // Otherwise we already computed that map
2413 474 : buildRefinementMap(*elem, *qrule, *qrule_face, -1, child, side);
2414 : }
2415 33 : }
2416 :
2417 : void
2418 90 : MooseMesh::buildPRefinementAndCoarseningMaps(Assembly * const assembly)
2419 : {
2420 90 : _elem_type_to_p_refinement_map.clear();
2421 90 : _elem_type_to_p_refinement_side_map.clear();
2422 90 : _elem_type_to_p_coarsening_map.clear();
2423 90 : _elem_type_to_p_coarsening_side_map.clear();
2424 :
2425 90 : std::map<ElemType, std::pair<Elem *, unsigned int>> elems_and_max_p_level;
2426 :
2427 32218 : for (const auto & elem : getMesh().active_element_ptr_range())
2428 : {
2429 32128 : const auto type = elem->type();
2430 32128 : auto & [picked_elem, max_p_level] = elems_and_max_p_level[type];
2431 32128 : if (!picked_elem)
2432 90 : picked_elem = elem;
2433 32128 : max_p_level = std::max(max_p_level, elem->p_level());
2434 90 : }
2435 :
2436 : // The only requirement on the FEType is that it can be arbitrarily p-refined
2437 90 : const FEType p_refinable_fe_type(CONSTANT, libMesh::MONOMIAL);
2438 90 : std::vector<Point> volume_ref_points_coarse, volume_ref_points_fine, face_ref_points_coarse,
2439 90 : face_ref_points_fine;
2440 90 : std::vector<unsigned int> p_levels;
2441 :
2442 180 : for (auto & [elem_type, elem_p_level_pair] : elems_and_max_p_level)
2443 : {
2444 90 : auto & [moose_elem, max_p_level] = elem_p_level_pair;
2445 90 : const auto dim = moose_elem->dim();
2446 : // Need to do this just once to get the right qrules put in place
2447 90 : assembly->setCurrentSubdomainID(moose_elem->subdomain_id());
2448 90 : assembly->reinit(moose_elem);
2449 90 : assembly->reinit(moose_elem, 0);
2450 90 : auto & qrule = assembly->writeableQRule();
2451 90 : auto & qrule_face = assembly->writeableQRuleFace();
2452 :
2453 90 : libMesh::Parallel::Communicator self_comm{};
2454 90 : ReplicatedMesh mesh(self_comm);
2455 90 : mesh.set_mesh_dimension(dim);
2456 630 : for (const auto & nd : moose_elem->node_ref_range())
2457 540 : mesh.add_point(nd);
2458 :
2459 90 : Elem * const elem = mesh.add_elem(Elem::build(elem_type).release());
2460 630 : for (const auto i : elem->node_index_range())
2461 540 : elem->set_node(i, mesh.node_ptr(i));
2462 :
2463 90 : std::unique_ptr<FEBase> fe_face(FEBase::build(dim, p_refinable_fe_type));
2464 90 : fe_face->get_phi();
2465 90 : const auto & face_phys_points = fe_face->get_xyz();
2466 90 : fe_face->attach_quadrature_rule(qrule_face);
2467 :
2468 90 : qrule->init(*elem);
2469 90 : volume_ref_points_coarse = qrule->get_points();
2470 90 : fe_face->reinit(elem, (unsigned int)0);
2471 90 : libMesh::FEMap::inverse_map(dim, elem, face_phys_points, face_ref_points_coarse);
2472 :
2473 90 : p_levels.resize(max_p_level + 1);
2474 90 : std::iota(p_levels.begin(), p_levels.end(), 0);
2475 90 : libMesh::MeshRefinement mesh_refinement(mesh);
2476 :
2477 306 : for (const auto p_level : p_levels)
2478 : {
2479 216 : mesh_refinement.uniformly_p_refine(1);
2480 216 : qrule->init(*elem);
2481 216 : volume_ref_points_fine = qrule->get_points();
2482 216 : fe_face->reinit(elem, (unsigned int)0);
2483 216 : libMesh::FEMap::inverse_map(dim, elem, face_phys_points, face_ref_points_fine);
2484 :
2485 216 : const auto map_key = std::make_pair(elem_type, p_level);
2486 216 : auto & volume_refine_map = _elem_type_to_p_refinement_map[map_key];
2487 216 : auto & face_refine_map = _elem_type_to_p_refinement_side_map[map_key];
2488 216 : auto & volume_coarsen_map = _elem_type_to_p_coarsening_map[map_key];
2489 216 : auto & face_coarsen_map = _elem_type_to_p_coarsening_side_map[map_key];
2490 :
2491 432 : auto fill_maps = [this](const auto & coarse_ref_points,
2492 : const auto & fine_ref_points,
2493 : auto & coarsen_map,
2494 : auto & refine_map)
2495 : {
2496 432 : mapPoints(fine_ref_points, coarse_ref_points, refine_map);
2497 432 : mapPoints(coarse_ref_points, fine_ref_points, coarsen_map);
2498 648 : };
2499 :
2500 216 : fill_maps(
2501 : volume_ref_points_coarse, volume_ref_points_fine, volume_coarsen_map, volume_refine_map);
2502 216 : fill_maps(face_ref_points_coarse, face_ref_points_fine, face_coarsen_map, face_refine_map);
2503 :
2504 : // With this level's maps filled our fine points now become our coarse points
2505 216 : volume_ref_points_fine.swap(volume_ref_points_coarse);
2506 216 : face_ref_points_fine.swap(face_ref_points_coarse);
2507 : }
2508 90 : }
2509 90 : }
2510 :
2511 : void
2512 57 : MooseMesh::buildRefinementAndCoarseningMaps(Assembly * const assembly)
2513 : {
2514 285 : TIME_SECTION("buildRefinementAndCoarseningMaps", 5, "Building Refinement And Coarsening Maps");
2515 57 : if (doingPRefinement())
2516 24 : buildPRefinementAndCoarseningMaps(assembly);
2517 : else
2518 33 : buildHRefinementAndCoarseningMaps(assembly);
2519 57 : }
2520 :
2521 : void
2522 690 : MooseMesh::buildRefinementMap(const Elem & elem,
2523 : QBase & qrule,
2524 : QBase & qrule_face,
2525 : int parent_side,
2526 : int child,
2527 : int child_side)
2528 : {
2529 3450 : TIME_SECTION("buildRefinementMap", 5, "Building Refinement Map");
2530 :
2531 690 : if (child == -1) // Doing volume mapping or parent side mapping
2532 : {
2533 : mooseAssert(parent_side == child_side,
2534 : "Parent side must match child_side if not passing a specific child!");
2535 :
2536 216 : std::pair<int, ElemType> the_pair(parent_side, elem.type());
2537 :
2538 216 : if (_elem_type_to_refinement_map.find(the_pair) != _elem_type_to_refinement_map.end())
2539 0 : mooseError("Already built a qp refinement map!");
2540 :
2541 216 : std::vector<std::pair<unsigned int, QpMap>> coarsen_map;
2542 216 : std::vector<std::vector<QpMap>> & refinement_map = _elem_type_to_refinement_map[the_pair];
2543 216 : findAdaptivityQpMaps(
2544 : &elem, qrule, qrule_face, refinement_map, coarsen_map, parent_side, child, child_side);
2545 216 : }
2546 : else // Need to map a child side to parent volume qps
2547 : {
2548 474 : std::pair<int, int> child_pair(child, child_side);
2549 :
2550 474 : if (_elem_type_to_child_side_refinement_map.find(elem.type()) !=
2551 1380 : _elem_type_to_child_side_refinement_map.end() &&
2552 432 : _elem_type_to_child_side_refinement_map[elem.type()].find(child_pair) !=
2553 906 : _elem_type_to_child_side_refinement_map[elem.type()].end())
2554 0 : mooseError("Already built a qp refinement map!");
2555 :
2556 474 : std::vector<std::pair<unsigned int, QpMap>> coarsen_map;
2557 : std::vector<std::vector<QpMap>> & refinement_map =
2558 474 : _elem_type_to_child_side_refinement_map[elem.type()][child_pair];
2559 474 : findAdaptivityQpMaps(
2560 : &elem, qrule, qrule_face, refinement_map, coarsen_map, parent_side, child, child_side);
2561 474 : }
2562 690 : }
2563 :
2564 : const std::vector<std::vector<QpMap>> &
2565 3422 : MooseMesh::getRefinementMap(const Elem & elem, int parent_side, int child, int child_side)
2566 : {
2567 3422 : if (child == -1) // Doing volume mapping or parent side mapping
2568 : {
2569 : mooseAssert(parent_side == child_side,
2570 : "Parent side must match child_side if not passing a specific child!");
2571 :
2572 3422 : std::pair<int, ElemType> the_pair(parent_side, elem.type());
2573 :
2574 3422 : if (_elem_type_to_refinement_map.find(the_pair) == _elem_type_to_refinement_map.end())
2575 0 : mooseError("Could not find a suitable qp refinement map!");
2576 :
2577 3422 : return _elem_type_to_refinement_map[the_pair];
2578 : }
2579 : else // Need to map a child side to parent volume qps
2580 : {
2581 0 : std::pair<int, int> child_pair(child, child_side);
2582 :
2583 0 : if (_elem_type_to_child_side_refinement_map.find(elem.type()) ==
2584 0 : _elem_type_to_child_side_refinement_map.end() ||
2585 0 : _elem_type_to_child_side_refinement_map[elem.type()].find(child_pair) ==
2586 0 : _elem_type_to_child_side_refinement_map[elem.type()].end())
2587 0 : mooseError("Could not find a suitable qp refinement map!");
2588 :
2589 0 : return _elem_type_to_child_side_refinement_map[elem.type()][child_pair];
2590 : }
2591 :
2592 : /**
2593 : * TODO: When running with parallel mesh + stateful adaptivty we will need to make sure that each
2594 : * processor has a complete map. This may require parallel communication. This is likely to
2595 : * happen
2596 : * when running on a mixed element mesh.
2597 : */
2598 : }
2599 :
2600 : void
2601 216 : MooseMesh::buildCoarseningMap(const Elem & elem, QBase & qrule, QBase & qrule_face, int input_side)
2602 : {
2603 1080 : TIME_SECTION("buildCoarseningMap", 5, "Building Coarsening Map");
2604 :
2605 216 : std::pair<int, ElemType> the_pair(input_side, elem.type());
2606 :
2607 216 : if (_elem_type_to_coarsening_map.find(the_pair) != _elem_type_to_coarsening_map.end())
2608 0 : mooseError("Already built a qp coarsening map!");
2609 :
2610 216 : std::vector<std::vector<QpMap>> refinement_map;
2611 : std::vector<std::pair<unsigned int, QpMap>> & coarsen_map =
2612 216 : _elem_type_to_coarsening_map[the_pair];
2613 :
2614 : // The -1 here is for a specific child. We don't do that for coarsening maps
2615 : // Also note that we're always mapping the same side to the same side (which is guaranteed by
2616 : // libMesh).
2617 216 : findAdaptivityQpMaps(
2618 : &elem, qrule, qrule_face, refinement_map, coarsen_map, input_side, -1, input_side);
2619 :
2620 : /**
2621 : * TODO: When running with parallel mesh + stateful adaptivty we will need to make sure that each
2622 : * processor has a complete map. This may require parallel communication. This is likely to
2623 : * happen
2624 : * when running on a mixed element mesh.
2625 : */
2626 216 : }
2627 :
2628 : const std::vector<std::pair<unsigned int, QpMap>> &
2629 1288 : MooseMesh::getCoarseningMap(const Elem & elem, int input_side)
2630 : {
2631 1288 : std::pair<int, ElemType> the_pair(input_side, elem.type());
2632 :
2633 1288 : if (_elem_type_to_coarsening_map.find(the_pair) == _elem_type_to_coarsening_map.end())
2634 0 : mooseError("Could not find a suitable qp refinement map!");
2635 :
2636 2576 : return _elem_type_to_coarsening_map[the_pair];
2637 : }
2638 :
2639 : void
2640 7038 : MooseMesh::mapPoints(const std::vector<Point> & from,
2641 : const std::vector<Point> & to,
2642 : std::vector<QpMap> & qp_map)
2643 : {
2644 7038 : unsigned int n_from = from.size();
2645 7038 : unsigned int n_to = to.size();
2646 :
2647 7038 : qp_map.resize(n_from);
2648 :
2649 61340 : for (unsigned int i = 0; i < n_from; ++i)
2650 : {
2651 54302 : const Point & from_point = from[i];
2652 :
2653 54302 : QpMap & current_map = qp_map[i];
2654 :
2655 1247054 : for (unsigned int j = 0; j < n_to; ++j)
2656 : {
2657 1192752 : const Point & to_point = to[j];
2658 1192752 : Real distance = (from_point - to_point).norm();
2659 :
2660 1192752 : if (distance < current_map._distance)
2661 : {
2662 167558 : current_map._distance = distance;
2663 167558 : current_map._from = i;
2664 167558 : current_map._to = j;
2665 : }
2666 : }
2667 : }
2668 7038 : }
2669 :
2670 : void
2671 906 : MooseMesh::findAdaptivityQpMaps(const Elem * template_elem,
2672 : QBase & qrule,
2673 : QBase & qrule_face,
2674 : std::vector<std::vector<QpMap>> & refinement_map,
2675 : std::vector<std::pair<unsigned int, QpMap>> & coarsen_map,
2676 : int parent_side,
2677 : int child,
2678 : int child_side)
2679 : {
2680 2718 : TIME_SECTION("findAdaptivityQpMaps", 5);
2681 :
2682 906 : ReplicatedMesh mesh(_communicator);
2683 906 : mesh.skip_partitioning(true);
2684 :
2685 906 : unsigned int dim = template_elem->dim();
2686 906 : mesh.set_mesh_dimension(dim);
2687 :
2688 7092 : for (unsigned int i = 0; i < template_elem->n_nodes(); ++i)
2689 6186 : mesh.add_point(template_elem->point(i));
2690 :
2691 906 : Elem * elem = mesh.add_elem(Elem::build(template_elem->type()).release());
2692 :
2693 7092 : for (unsigned int i = 0; i < template_elem->n_nodes(); ++i)
2694 6186 : elem->set_node(i, mesh.node_ptr(i));
2695 :
2696 906 : std::unique_ptr<FEBase> fe(FEBase::build(dim, FEType()));
2697 906 : fe->get_phi();
2698 906 : const std::vector<Point> & q_points_volume = fe->get_xyz();
2699 :
2700 906 : std::unique_ptr<FEBase> fe_face(FEBase::build(dim, FEType()));
2701 906 : fe_face->get_phi();
2702 906 : const std::vector<Point> & q_points_face = fe_face->get_xyz();
2703 :
2704 906 : fe->attach_quadrature_rule(&qrule);
2705 906 : fe_face->attach_quadrature_rule(&qrule_face);
2706 :
2707 : // The current q_points (locations in *physical* space)
2708 : const std::vector<Point> * q_points;
2709 :
2710 906 : if (parent_side != -1)
2711 : {
2712 348 : fe_face->reinit(elem, parent_side);
2713 348 : q_points = &q_points_face;
2714 : }
2715 : else
2716 : {
2717 558 : fe->reinit(elem);
2718 558 : q_points = &q_points_volume;
2719 : }
2720 :
2721 906 : std::vector<Point> parent_ref_points;
2722 :
2723 906 : libMesh::FEMap::inverse_map(elem->dim(), elem, *q_points, parent_ref_points);
2724 906 : libMesh::MeshRefinement mesh_refinement(mesh);
2725 906 : mesh_refinement.uniformly_refine(1);
2726 :
2727 : // A map from the child element index to the locations of all the child's quadrature points in
2728 : // *reference* space. Note that we use a map here instead of a vector because the caller can
2729 : // pass an explicit child index. We are not guaranteed to have a sequence from [0, n_children)
2730 906 : std::map<unsigned int, std::vector<Point>> child_to_ref_points;
2731 :
2732 906 : unsigned int n_children = elem->n_children();
2733 :
2734 906 : refinement_map.resize(n_children);
2735 :
2736 906 : std::vector<unsigned int> children;
2737 :
2738 906 : if (child != -1) // Passed in a child explicitly
2739 474 : children.push_back(child);
2740 : else
2741 : {
2742 432 : children.resize(n_children);
2743 2724 : for (unsigned int child = 0; child < n_children; ++child)
2744 2292 : children[child] = child;
2745 : }
2746 :
2747 3672 : for (unsigned int i = 0; i < children.size(); ++i)
2748 : {
2749 2766 : unsigned int child = children[i];
2750 :
2751 2766 : if ((parent_side != -1 && !elem->is_child_on_side(child, parent_side)))
2752 948 : continue;
2753 :
2754 1818 : const Elem * child_elem = elem->child_ptr(child);
2755 :
2756 1818 : if (child_side != -1)
2757 : {
2758 1422 : fe_face->reinit(child_elem, child_side);
2759 1422 : q_points = &q_points_face;
2760 : }
2761 : else
2762 : {
2763 396 : fe->reinit(child_elem);
2764 396 : q_points = &q_points_volume;
2765 : }
2766 :
2767 1818 : std::vector<Point> child_ref_points;
2768 :
2769 1818 : libMesh::FEMap::inverse_map(elem->dim(), elem, *q_points, child_ref_points);
2770 1818 : child_to_ref_points[child] = child_ref_points;
2771 :
2772 1818 : std::vector<QpMap> & qp_map = refinement_map[child];
2773 :
2774 : // Find the closest parent_qp to each child_qp
2775 1818 : mapPoints(child_ref_points, parent_ref_points, qp_map);
2776 1818 : }
2777 :
2778 906 : coarsen_map.resize(parent_ref_points.size());
2779 :
2780 : // For each parent qp find the closest child qp
2781 6210 : for (unsigned int child = 0; child < n_children; child++)
2782 : {
2783 5304 : if (parent_side != -1 && !elem->is_child_on_side(child, child_side))
2784 948 : continue;
2785 :
2786 4356 : std::vector<Point> & child_ref_points = child_to_ref_points[child];
2787 :
2788 4356 : std::vector<QpMap> qp_map;
2789 :
2790 : // Find all of the closest points from parent_qp to _THIS_ child's qp
2791 4356 : mapPoints(parent_ref_points, child_ref_points, qp_map);
2792 :
2793 : // Check those to see if they are closer than what we currently have for each point
2794 32856 : for (unsigned int parent_qp = 0; parent_qp < parent_ref_points.size(); ++parent_qp)
2795 : {
2796 28500 : std::pair<unsigned int, QpMap> & child_and_map = coarsen_map[parent_qp];
2797 28500 : unsigned int & closest_child = child_and_map.first;
2798 28500 : QpMap & closest_map = child_and_map.second;
2799 :
2800 28500 : QpMap & current_map = qp_map[parent_qp];
2801 :
2802 28500 : if (current_map._distance < closest_map._distance)
2803 : {
2804 6300 : closest_child = child;
2805 6300 : closest_map = current_map;
2806 : }
2807 : }
2808 4356 : }
2809 906 : }
2810 :
2811 : void
2812 0 : MooseMesh::changeBoundaryId(const boundary_id_type old_id,
2813 : const boundary_id_type new_id,
2814 : bool delete_prev)
2815 : {
2816 0 : TIME_SECTION("changeBoundaryId", 6);
2817 0 : changeBoundaryId(getMesh(), old_id, new_id, delete_prev);
2818 0 : }
2819 :
2820 : void
2821 0 : MooseMesh::changeBoundaryId(MeshBase & mesh,
2822 : const boundary_id_type old_id,
2823 : const boundary_id_type new_id,
2824 : bool delete_prev)
2825 : {
2826 : // Get a reference to our BoundaryInfo object, we will use it several times below...
2827 0 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
2828 :
2829 : // Container to catch ids passed back from BoundaryInfo
2830 0 : std::vector<boundary_id_type> old_ids;
2831 :
2832 : // Only level-0 elements store BCs. Loop over them.
2833 0 : for (auto & elem : as_range(mesh.level_elements_begin(0), mesh.level_elements_end(0)))
2834 : {
2835 0 : unsigned int n_sides = elem->n_sides();
2836 0 : for (unsigned int s = 0; s != n_sides; ++s)
2837 : {
2838 0 : boundary_info.boundary_ids(elem, s, old_ids);
2839 0 : if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
2840 : {
2841 0 : std::vector<boundary_id_type> new_ids(old_ids);
2842 0 : std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
2843 0 : if (delete_prev)
2844 : {
2845 0 : boundary_info.remove_side(elem, s);
2846 0 : boundary_info.add_side(elem, s, new_ids);
2847 : }
2848 : else
2849 0 : boundary_info.add_side(elem, s, new_ids);
2850 0 : }
2851 : }
2852 0 : }
2853 :
2854 : // Remove any remaining references to the old ID from the
2855 : // BoundaryInfo object. This prevents things like empty sidesets
2856 : // from showing up when printing information, etc.
2857 0 : if (delete_prev)
2858 0 : boundary_info.remove_id(old_id);
2859 :
2860 : // The cached boundary id sets will need re-preparation
2861 0 : mesh.unset_has_boundary_id_sets();
2862 0 : }
2863 :
2864 : const RealVectorValue &
2865 0 : MooseMesh::getNormalByBoundaryID(BoundaryID id) const
2866 : {
2867 : mooseAssert(_boundary_to_normal_map.get() != nullptr, "Boundary To Normal Map not built!");
2868 :
2869 : // Note: Boundaries that are not in the map (existing boundaries) will default
2870 : // construct a new RealVectorValue - (x,y,z)=(0, 0, 0)
2871 0 : return (*_boundary_to_normal_map)[id];
2872 : }
2873 :
2874 : MooseMesh &
2875 0 : MooseMesh::clone() const
2876 : {
2877 0 : mooseError("MooseMesh::clone() is no longer supported, use MooseMesh::safeClone() instead.");
2878 : }
2879 :
2880 : void
2881 74277 : MooseMesh::determineUseDistributedMesh()
2882 : {
2883 74277 : switch (_parallel_type)
2884 : {
2885 69415 : case ParallelType::DEFAULT:
2886 : // The user did not specify 'parallel_type = XYZ' in the input file,
2887 : // so we allow the --distributed-mesh command line arg to possibly turn
2888 : // on DistributedMesh. If the command line arg is not present, we pick ReplicatedMesh.
2889 69415 : if (_app.getDistributedMeshOnCommandLine())
2890 10876 : _use_distributed_mesh = true;
2891 69415 : break;
2892 3515 : case ParallelType::REPLICATED:
2893 3515 : if (_app.getDistributedMeshOnCommandLine() || _is_nemesis || _is_split)
2894 741 : _parallel_type_overridden = true;
2895 3515 : _use_distributed_mesh = false;
2896 3515 : break;
2897 1347 : case ParallelType::DISTRIBUTED:
2898 1347 : _use_distributed_mesh = true;
2899 1347 : break;
2900 : }
2901 :
2902 : // If the user specifies 'nemesis = true' in the Mesh block, or they are using --use-split,
2903 : // we must use DistributedMesh.
2904 74277 : if (_is_nemesis || _is_split)
2905 562 : _use_distributed_mesh = true;
2906 74277 : }
2907 :
2908 : std::unique_ptr<MeshBase>
2909 69440 : MooseMesh::buildMeshBaseObject(unsigned int dim)
2910 : {
2911 69440 : std::unique_ptr<MeshBase> mesh;
2912 69440 : if (_use_distributed_mesh)
2913 10902 : mesh = buildTypedMesh<DistributedMesh>(dim);
2914 : else
2915 58538 : mesh = buildTypedMesh<ReplicatedMesh>(dim);
2916 :
2917 69440 : return mesh;
2918 0 : }
2919 :
2920 : void
2921 66366 : MooseMesh::setMeshBase(std::unique_ptr<MeshBase> mesh_base)
2922 : {
2923 66366 : _mesh = std::move(mesh_base);
2924 66366 : _mesh->allow_remote_element_removal(_allow_remote_element_removal);
2925 66366 : }
2926 :
2927 : void
2928 65878 : MooseMesh::init()
2929 : {
2930 : /**
2931 : * If the mesh base hasn't been constructed by the time init is called, just do it here.
2932 : * This can happen if somebody builds a mesh outside of the normal Action system. Forcing
2933 : * developers to create, construct the MeshBase, and then init separately is a bit much for casual
2934 : * use but it gives us the ability to run MeshGenerators in-between.
2935 : */
2936 65878 : if (!_mesh)
2937 10 : _mesh = buildMeshBaseObject();
2938 :
2939 65878 : if (_app.isSplitMesh() && _use_distributed_mesh)
2940 0 : mooseError("You cannot use the mesh splitter capability with DistributedMesh!");
2941 :
2942 197634 : TIME_SECTION("init", 2);
2943 :
2944 65878 : if (_app.isRecovering() && _allow_recovery && _app.isUltimateMaster())
2945 : {
2946 : // Some partitioners are not idempotent. Some recovery data
2947 : // files require partitioning to match mesh partitioning. This
2948 : // means that, when recovering, we can't safely repartition.
2949 3349 : const bool skip_partitioning_later = getMesh().skip_partitioning();
2950 3349 : getMesh().skip_partitioning(true);
2951 3349 : const bool allow_renumbering_later = getMesh().allow_renumbering();
2952 3349 : getMesh().allow_renumbering(false);
2953 :
2954 : // For now, only read the recovery mesh on the Ultimate Master..
2955 : // sub-apps need to just build their mesh like normal
2956 : {
2957 10047 : TIME_SECTION("readRecoveredMesh", 2);
2958 3349 : getMesh().read(_app.getRestartRecoverFileBase() + MooseApp::checkpointSuffix());
2959 3349 : }
2960 :
2961 3349 : getMesh().allow_renumbering(allow_renumbering_later);
2962 3349 : getMesh().skip_partitioning(skip_partitioning_later);
2963 : }
2964 : else // Normally just build the mesh
2965 : {
2966 : // Don't allow partitioning during building
2967 62529 : if (_app.isSplitMesh())
2968 89 : getMesh().skip_partitioning(true);
2969 62529 : buildMesh();
2970 :
2971 187569 : if (getParam<bool>("build_all_side_lowerd_mesh"))
2972 205 : buildLowerDMesh();
2973 : }
2974 65872 : }
2975 :
2976 : std::vector<std::filesystem::path>
2977 13111 : MooseMesh::writeRecoveryFiles(const std::filesystem::path & file_base)
2978 : {
2979 13111 : CheckpointIO io(getMesh(), false);
2980 13111 : io.write(file_base);
2981 26222 : return {};
2982 13111 : }
2983 :
2984 : unsigned int
2985 91569640 : MooseMesh::dimension() const
2986 : {
2987 91569640 : return getMesh().mesh_dimension();
2988 : }
2989 :
2990 : unsigned int
2991 37845 : MooseMesh::effectiveSpatialDimension() const
2992 : {
2993 37845 : const Real abs_zero = 1e-12;
2994 :
2995 : // See if the mesh is completely containd in the z and y planes to calculate effective spatial
2996 : // dim
2997 69777 : for (unsigned int dim = LIBMESH_DIM; dim >= 1; --dim)
2998 69777 : if (dimensionWidth(dim - 1) >= abs_zero)
2999 37845 : return dim;
3000 :
3001 : // If we get here, we have a 1D mesh on the x-axis.
3002 0 : return 1;
3003 : }
3004 :
3005 : unsigned int
3006 89991 : MooseMesh::getBlocksMaxDimension(const std::vector<SubdomainName> & blocks) const
3007 : {
3008 89991 : const auto & mesh = getMesh();
3009 :
3010 : // Take a shortcut if possible
3011 89991 : if (const auto & elem_dims = mesh.elem_dimensions(); mesh.is_prepared() && elem_dims.size() == 1)
3012 78847 : return *elem_dims.begin();
3013 :
3014 11144 : unsigned short dim = 0;
3015 11144 : const auto subdomain_ids = getSubdomainIDs(blocks);
3016 11144 : const std::set<SubdomainID> subdomain_ids_set(subdomain_ids.begin(), subdomain_ids.end());
3017 2328536 : for (const auto & elem : mesh.active_subdomain_set_elements_ptr_range(subdomain_ids_set))
3018 2328536 : dim = std::max(dim, elem->dim());
3019 :
3020 : // Get the maximumal globally
3021 11144 : _communicator.max(dim);
3022 11144 : return dim;
3023 11144 : }
3024 :
3025 : std::vector<BoundaryID>
3026 106488560 : MooseMesh::getBoundaryIDs(const Elem * const elem, const unsigned short int side) const
3027 : {
3028 106488560 : std::vector<BoundaryID> ids;
3029 106488560 : getMesh().get_boundary_info().boundary_ids(elem, side, ids);
3030 106488560 : return ids;
3031 0 : }
3032 :
3033 : std::vector<std::vector<BoundaryID>>
3034 407162082 : MooseMesh::getBoundaryIDs(const Elem * const elem) const
3035 : {
3036 407162082 : std::vector<std::vector<BoundaryID>> ids;
3037 407162082 : getMesh().get_boundary_info().side_boundary_ids(elem, ids);
3038 407162082 : return ids;
3039 0 : }
3040 :
3041 : const std::set<BoundaryID> &
3042 506989 : MooseMesh::getBoundaryIDs() const
3043 : {
3044 506989 : return getMesh().get_boundary_info().get_boundary_ids();
3045 : }
3046 :
3047 : void
3048 224181 : MooseMesh::buildNodeListFromSideList()
3049 : {
3050 224181 : auto & boundary_info = getMesh().get_boundary_info();
3051 :
3052 224181 : if (_construct_node_list_from_side_list)
3053 : {
3054 224155 : const std::set<boundary_id_type> & side_bcids = boundary_info.get_side_boundary_ids();
3055 :
3056 224155 : if (_displace_node_list_by_side_list)
3057 : {
3058 : // Don't want to use auto here - the rbegin trick relies on a
3059 : // sorted set and we want the compiler to scream if libMesh ever
3060 : // switches type
3061 224155 : const std::set<boundary_id_type> & node_bcids = boundary_info.get_node_boundary_ids();
3062 :
3063 : // If we've got a reasonable largest BC id, we can just use the
3064 : // subsequent unused ones
3065 224155 : boundary_id_type next_bcid = 0;
3066 224155 : if (!node_bcids.empty())
3067 215729 : next_bcid = std::max(next_bcid, cast_int<boundary_id_type>(*node_bcids.rbegin() + 1));
3068 224155 : if (!side_bcids.empty())
3069 219194 : next_bcid = std::max(next_bcid, cast_int<boundary_id_type>(*side_bcids.rbegin() + 1));
3070 :
3071 : // We need all processors to agree on the id to use, even when
3072 : // each only sees the bcids on their own portions of a
3073 : // distributed mesh.
3074 224155 : _communicator.max(next_bcid);
3075 :
3076 : // If we've got an unreasonably high largest BC id, we should
3077 : // probably just search for unused ones with moderate values, so we
3078 : // don't risk wrapping.
3079 224155 : if (next_bcid > 1000 || next_bcid <= 0)
3080 3047 : next_bcid = 1000;
3081 :
3082 : // If any side bcid is already a node bcid with a different name,
3083 : // that's a different boundary condition that we need to reassign
3084 : // rather than overwrite or merge to.
3085 1080063 : for (auto bcid : side_bcids)
3086 1679714 : if (node_bcids.count(bcid) &&
3087 823806 : (boundary_info.get_sideset_name(bcid) != boundary_info.get_nodeset_name(bcid)))
3088 : {
3089 2002 : boundary_info.renumber_node_id(bcid, next_bcid);
3090 : do
3091 : {
3092 2012 : ++next_bcid;
3093 2012 : } while (node_bcids.count(next_bcid) || side_bcids.count(next_bcid));
3094 : }
3095 : }
3096 :
3097 : // For any side bcid that has a name, make sure that our new node
3098 : // bcid is given the same name. We need to iterate over the
3099 : // actual name map (which is global) here, not over side_bcids
3100 : // (which only includes local ids on a distributed mesh).
3101 1067423 : for (auto & [id, name] : boundary_info.get_sideset_name_map())
3102 843268 : boundary_info.nodeset_name(id) = name;
3103 :
3104 224155 : boundary_info.build_node_list_from_side_list();
3105 : }
3106 224181 : }
3107 :
3108 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
3109 218 : MooseMesh::buildSideList()
3110 : {
3111 218 : return getMesh().get_boundary_info().build_side_list();
3112 : }
3113 :
3114 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
3115 4742 : MooseMesh::buildActiveSideList() const
3116 : {
3117 4742 : return getMesh().get_boundary_info().build_active_side_list();
3118 : }
3119 :
3120 : unsigned int
3121 26176 : MooseMesh::sideWithBoundaryID(const Elem * const elem, const BoundaryID boundary_id) const
3122 : {
3123 26176 : return getMesh().get_boundary_info().side_with_boundary_id(elem, boundary_id);
3124 : }
3125 :
3126 : MeshBase::node_iterator
3127 3402 : MooseMesh::localNodesBegin()
3128 : {
3129 3402 : return getMesh().local_nodes_begin();
3130 : }
3131 :
3132 : MeshBase::node_iterator
3133 3402 : MooseMesh::localNodesEnd()
3134 : {
3135 3402 : return getMesh().local_nodes_end();
3136 : }
3137 :
3138 : MeshBase::const_node_iterator
3139 0 : MooseMesh::localNodesBegin() const
3140 : {
3141 0 : return getMesh().local_nodes_begin();
3142 : }
3143 :
3144 : MeshBase::const_node_iterator
3145 0 : MooseMesh::localNodesEnd() const
3146 : {
3147 0 : return getMesh().local_nodes_end();
3148 : }
3149 :
3150 : MeshBase::element_iterator
3151 148750 : MooseMesh::activeLocalElementsBegin()
3152 : {
3153 148750 : return getMesh().active_local_elements_begin();
3154 : }
3155 :
3156 : const MeshBase::element_iterator
3157 148750 : MooseMesh::activeLocalElementsEnd()
3158 : {
3159 148750 : return getMesh().active_local_elements_end();
3160 : }
3161 :
3162 : MeshBase::const_element_iterator
3163 0 : MooseMesh::activeLocalElementsBegin() const
3164 : {
3165 0 : return getMesh().active_local_elements_begin();
3166 : }
3167 :
3168 : const MeshBase::const_element_iterator
3169 0 : MooseMesh::activeLocalElementsEnd() const
3170 : {
3171 0 : return getMesh().active_local_elements_end();
3172 : }
3173 :
3174 : dof_id_type
3175 55393 : MooseMesh::nNodes() const
3176 : {
3177 55393 : return getMesh().n_nodes();
3178 : }
3179 :
3180 : dof_id_type
3181 1562 : MooseMesh::nElem() const
3182 : {
3183 1562 : return getMesh().n_elem();
3184 : }
3185 :
3186 : dof_id_type
3187 0 : MooseMesh::maxNodeId() const
3188 : {
3189 0 : return getMesh().max_node_id();
3190 : }
3191 :
3192 : dof_id_type
3193 0 : MooseMesh::maxElemId() const
3194 : {
3195 0 : return getMesh().max_elem_id();
3196 : }
3197 :
3198 : Elem *
3199 0 : MooseMesh::elem(const dof_id_type i)
3200 : {
3201 0 : mooseDeprecated("MooseMesh::elem() is deprecated, please use MooseMesh::elemPtr() instead");
3202 0 : return elemPtr(i);
3203 : }
3204 :
3205 : const Elem *
3206 0 : MooseMesh::elem(const dof_id_type i) const
3207 : {
3208 0 : mooseDeprecated("MooseMesh::elem() is deprecated, please use MooseMesh::elemPtr() instead");
3209 0 : return elemPtr(i);
3210 : }
3211 :
3212 : Elem *
3213 15691943 : MooseMesh::elemPtr(const dof_id_type i)
3214 : {
3215 15691943 : return getMesh().elem_ptr(i);
3216 : }
3217 :
3218 : const Elem *
3219 1245362 : MooseMesh::elemPtr(const dof_id_type i) const
3220 : {
3221 1245362 : return getMesh().elem_ptr(i);
3222 : }
3223 :
3224 : Elem *
3225 21675 : MooseMesh::queryElemPtr(const dof_id_type i)
3226 : {
3227 21675 : return getMesh().query_elem_ptr(i);
3228 : }
3229 :
3230 : const Elem *
3231 36392 : MooseMesh::queryElemPtr(const dof_id_type i) const
3232 : {
3233 36392 : return getMesh().query_elem_ptr(i);
3234 : }
3235 :
3236 : bool
3237 0 : MooseMesh::prepared() const
3238 : {
3239 0 : return _mesh->is_prepared() && _moose_mesh_prepared;
3240 : }
3241 :
3242 : void
3243 0 : MooseMesh::prepared(bool state)
3244 : {
3245 0 : if (state)
3246 0 : mooseError("We don't have any right to tell the libmesh mesh that it *is* prepared. Only a "
3247 : "call to prepare_for_use should tell us that");
3248 :
3249 : // Some people may call this even before we have a MeshBase object. This isn't dangerous really
3250 : // because when the MeshBase object is born, it knows it's in an unprepared state
3251 0 : if (_mesh)
3252 0 : _mesh->unset_is_prepared();
3253 :
3254 : // If the libMesh mesh isn't preparead, then our MooseMesh wrapper is also no longer prepared
3255 0 : _moose_mesh_prepared = false;
3256 :
3257 : /**
3258 : * If we are explicitly setting the mesh to not prepared, then we've likely modified the mesh
3259 : * and can no longer make assumptions about orthogonality. We really should recheck.
3260 : */
3261 0 : _regular_orthogonal_mesh = false;
3262 0 : }
3263 :
3264 : void
3265 0 : MooseMesh::needsPrepareForUse()
3266 : {
3267 0 : prepared(false);
3268 0 : }
3269 :
3270 : const std::set<SubdomainID> &
3271 8073881 : MooseMesh::meshSubdomains() const
3272 : {
3273 8073881 : return _mesh_subdomains;
3274 : }
3275 :
3276 : const std::set<BoundaryID> &
3277 12498 : MooseMesh::meshBoundaryIds() const
3278 : {
3279 12498 : return _mesh_boundary_ids;
3280 : }
3281 :
3282 : const std::set<BoundaryID> &
3283 29636 : MooseMesh::meshSidesetIds() const
3284 : {
3285 29636 : return _mesh_sideset_ids;
3286 : }
3287 :
3288 : const std::set<BoundaryID> &
3289 149568 : MooseMesh::meshNodesetIds() const
3290 : {
3291 149568 : return _mesh_nodeset_ids;
3292 : }
3293 :
3294 : void
3295 0 : MooseMesh::setMeshBoundaryIDs(std::set<BoundaryID> boundary_IDs)
3296 : {
3297 0 : _mesh_boundary_ids = boundary_IDs;
3298 0 : }
3299 :
3300 : void
3301 0 : MooseMesh::setBoundaryToNormalMap(
3302 : std::unique_ptr<std::map<BoundaryID, RealVectorValue>> boundary_map)
3303 : {
3304 0 : _boundary_to_normal_map = std::move(boundary_map);
3305 0 : }
3306 :
3307 : void
3308 0 : MooseMesh::setBoundaryToNormalMap(std::map<BoundaryID, RealVectorValue> * boundary_map)
3309 : {
3310 0 : mooseDeprecated("setBoundaryToNormalMap(std::map<BoundaryID, RealVectorValue> * boundary_map) is "
3311 : "deprecated, use the unique_ptr version instead");
3312 0 : _boundary_to_normal_map.reset(boundary_map);
3313 0 : }
3314 :
3315 : unsigned int
3316 125591 : MooseMesh::uniformRefineLevel() const
3317 : {
3318 125591 : return _uniform_refine_level;
3319 : }
3320 :
3321 : void
3322 67894 : MooseMesh::setUniformRefineLevel(unsigned int level, bool deletion)
3323 : {
3324 67894 : _uniform_refine_level = level;
3325 67894 : _skip_deletion_repartition_after_refine = deletion;
3326 67894 : }
3327 :
3328 : void
3329 56144 : MooseMesh::addGhostedBoundary(BoundaryID boundary_id)
3330 : {
3331 56144 : _ghosted_boundaries.insert(boundary_id);
3332 56144 : }
3333 :
3334 : void
3335 0 : MooseMesh::setGhostedBoundaryInflation(const std::vector<Real> & inflation)
3336 : {
3337 0 : _ghosted_boundaries_inflation = inflation;
3338 0 : }
3339 :
3340 : const std::set<unsigned int> &
3341 0 : MooseMesh::getGhostedBoundaries() const
3342 : {
3343 0 : return _ghosted_boundaries;
3344 : }
3345 :
3346 : const std::vector<Real> &
3347 11484 : MooseMesh::getGhostedBoundaryInflation() const
3348 : {
3349 11484 : return _ghosted_boundaries_inflation;
3350 : }
3351 :
3352 : namespace // Anonymous namespace for helpers
3353 : {
3354 : // A class for templated methods that expect output iterator
3355 : // arguments, which adds objects to the Mesh.
3356 : // Although extra_ghost_elem_inserter can add any object, we
3357 : // template it around object type so that type inference and
3358 : // iterator_traits will work.
3359 : // This object specifically is used to insert extra ghost elems into the mesh
3360 : template <typename T>
3361 : struct extra_ghost_elem_inserter
3362 : {
3363 : using iterator_category = std::output_iterator_tag;
3364 : using value_type = T;
3365 :
3366 44944 : extra_ghost_elem_inserter(DistributedMesh & m) : mesh(m) {}
3367 :
3368 19101 : void operator=(const Elem * e) { mesh.add_extra_ghost_elem(const_cast<Elem *>(e)); }
3369 :
3370 35698 : void operator=(Node * n) { mesh.add_node(n); }
3371 :
3372 : void operator=(Point * p) { mesh.add_point(*p); }
3373 :
3374 : extra_ghost_elem_inserter & operator++() { return *this; }
3375 :
3376 54799 : extra_ghost_elem_inserter operator++(int) { return extra_ghost_elem_inserter(*this); }
3377 :
3378 : // We don't return a reference-to-T here because we don't want to
3379 : // construct one or have any of its methods called. We just want
3380 : // to allow the returned object to be able to do mesh insertions
3381 : // with operator=().
3382 54799 : extra_ghost_elem_inserter & operator*() { return *this; }
3383 :
3384 : private:
3385 : DistributedMesh & mesh;
3386 : };
3387 :
3388 : /**
3389 : * Specific weak ordering for Elem *'s to be used in a set.
3390 : * We use the id, but first sort by level. This guarantees
3391 : * when traversing the set from beginning to end the lower
3392 : * level (parent) elements are encountered first.
3393 : *
3394 : * This was swiped from libMesh mesh_communication.C, and ought to be
3395 : * replaced with libMesh::CompareElemIdsByLevel just as soon as I refactor to
3396 : * create that - @roystgnr
3397 : */
3398 : struct CompareElemsByLevel
3399 : {
3400 104881 : bool operator()(const Elem * a, const Elem * b) const
3401 : {
3402 : libmesh_assert(a);
3403 : libmesh_assert(b);
3404 104881 : const unsigned int al = a->level(), bl = b->level();
3405 104881 : const dof_id_type aid = a->id(), bid = b->id();
3406 :
3407 104881 : return (al == bl) ? aid < bid : al < bl;
3408 : }
3409 : };
3410 :
3411 : } // anonymous namespace
3412 :
3413 : void
3414 138917 : MooseMesh::ghostGhostedBoundaries()
3415 : {
3416 : // No need to do this if using a serial mesh
3417 : // We do not need to ghost boundary elements when _need_ghost_ghosted_boundaries
3418 : // is not true. _need_ghost_ghosted_boundaries can be set by a mesh generator
3419 : // where boundaries are already ghosted accordingly
3420 138917 : if (!_use_distributed_mesh || !_need_ghost_ghosted_boundaries)
3421 116445 : return;
3422 :
3423 67416 : TIME_SECTION("GhostGhostedBoundaries", 3);
3424 :
3425 : parallel_object_only();
3426 :
3427 22472 : DistributedMesh & mesh = dynamic_cast<DistributedMesh &>(getMesh());
3428 :
3429 : // We clear ghosted elements that were added by previous invocations of this
3430 : // method but leave ghosted elements that were added by other code, e.g.
3431 : // OversampleOutput, untouched
3432 22472 : mesh.clear_extra_ghost_elems(_ghost_elems_from_ghost_boundaries);
3433 22472 : _ghost_elems_from_ghost_boundaries.clear();
3434 :
3435 22472 : std::set<const Elem *, CompareElemsByLevel> boundary_elems_to_ghost;
3436 22472 : std::set<Node *> connected_nodes_to_ghost;
3437 :
3438 22472 : std::vector<const Elem *> family_tree;
3439 :
3440 904489 : for (const auto & t : mesh.get_boundary_info().build_side_list())
3441 : {
3442 882017 : auto elem_id = std::get<0>(t);
3443 882017 : auto bc_id = std::get<2>(t);
3444 :
3445 882017 : if (_ghosted_boundaries.find(bc_id) != _ghosted_boundaries.end())
3446 : {
3447 5695 : Elem * elem = mesh.elem_ptr(elem_id);
3448 :
3449 : #ifdef LIBMESH_ENABLE_AMR
3450 5695 : elem->family_tree(family_tree);
3451 5695 : Elem * parent = elem->parent();
3452 5695 : while (parent)
3453 : {
3454 0 : family_tree.push_back(parent);
3455 0 : parent = parent->parent();
3456 : }
3457 : #else
3458 : family_tree.clear();
3459 : family_tree.push_back(elem);
3460 : #endif
3461 15230 : for (const auto & felem : family_tree)
3462 : {
3463 9535 : boundary_elems_to_ghost.insert(felem);
3464 :
3465 : // The entries of connected_nodes_to_ghost need to be
3466 : // non-constant, so that they will work in things like
3467 : // UpdateDisplacedMeshThread. The container returned by
3468 : // family_tree contains const Elems even when the Elem
3469 : // it is called on is non-const, so once that interface
3470 : // gets fixed we can remove this const_cast.
3471 57258 : for (unsigned int n = 0; n < felem->n_nodes(); ++n)
3472 47723 : connected_nodes_to_ghost.insert(const_cast<Node *>(felem->node_ptr(n)));
3473 : }
3474 : }
3475 22472 : }
3476 :
3477 : // We really do want to store this by value instead of by reference
3478 22472 : const auto prior_ghost_elems = mesh.extra_ghost_elems();
3479 :
3480 22472 : mesh.comm().allgather_packed_range(&mesh,
3481 : connected_nodes_to_ghost.begin(),
3482 : connected_nodes_to_ghost.end(),
3483 : extra_ghost_elem_inserter<Node>(mesh));
3484 :
3485 22472 : mesh.comm().allgather_packed_range(&mesh,
3486 : boundary_elems_to_ghost.begin(),
3487 : boundary_elems_to_ghost.end(),
3488 : extra_ghost_elem_inserter<Elem>(mesh));
3489 :
3490 22472 : const auto & current_ghost_elems = mesh.extra_ghost_elems();
3491 :
3492 44944 : std::set_difference(current_ghost_elems.begin(),
3493 : current_ghost_elems.end(),
3494 : prior_ghost_elems.begin(),
3495 : prior_ghost_elems.end(),
3496 22472 : std::inserter(_ghost_elems_from_ghost_boundaries,
3497 : _ghost_elems_from_ghost_boundaries.begin()));
3498 22472 : }
3499 :
3500 : unsigned int
3501 11484 : MooseMesh::getPatchSize() const
3502 : {
3503 11484 : return _patch_size;
3504 : }
3505 :
3506 : void
3507 0 : MooseMesh::setPatchUpdateStrategy(Moose::PatchUpdateType patch_update_strategy)
3508 : {
3509 0 : _patch_update_strategy = patch_update_strategy;
3510 0 : }
3511 :
3512 : const Moose::PatchUpdateType &
3513 36409 : MooseMesh::getPatchUpdateStrategy() const
3514 : {
3515 36409 : return _patch_update_strategy;
3516 : }
3517 :
3518 : BoundingBox
3519 114404 : MooseMesh::getInflatedProcessorBoundingBox(Real inflation_multiplier) const
3520 : {
3521 : // Grab a bounding box to speed things up. Note that
3522 : // local_bounding_box is *not* equivalent to processor_bounding_box
3523 : // with processor_id() except in serial.
3524 114404 : BoundingBox bbox = MeshTools::create_local_bounding_box(getMesh());
3525 :
3526 : // Inflate the bbox just a bit to deal with roundoff
3527 : // Adding 1% of the diagonal size in each direction on each end
3528 114404 : Real inflation_amount = inflation_multiplier * (bbox.max() - bbox.min()).norm();
3529 114404 : Point inflation(inflation_amount, inflation_amount, inflation_amount);
3530 :
3531 114404 : bbox.first -= inflation; // min
3532 114404 : bbox.second += inflation; // max
3533 :
3534 228808 : return bbox;
3535 : }
3536 :
3537 160588 : MooseMesh::operator libMesh::MeshBase &() { return getMesh(); }
3538 :
3539 2990 : MooseMesh::operator const libMesh::MeshBase &() const { return getMesh(); }
3540 :
3541 : const MeshBase *
3542 440171 : MooseMesh::getMeshPtr() const
3543 : {
3544 440171 : return _mesh.get();
3545 : }
3546 :
3547 : MeshBase &
3548 59423200 : MooseMesh::getMesh()
3549 : {
3550 : mooseAssert(_mesh, "Mesh hasn't been created");
3551 59423200 : return *_mesh;
3552 : }
3553 :
3554 : const MeshBase &
3555 707885737 : MooseMesh::getMesh() const
3556 : {
3557 : mooseAssert(_mesh, "Mesh hasn't been created");
3558 707885737 : return *_mesh;
3559 : }
3560 :
3561 : void
3562 0 : MooseMesh::printInfo(std::ostream & os, const unsigned int verbosity /* = 0 */) const
3563 : {
3564 0 : os << '\n';
3565 0 : getMesh().print_info(os, verbosity);
3566 0 : os << std::flush;
3567 0 : }
3568 :
3569 : const std::vector<dof_id_type> &
3570 229 : MooseMesh::getNodeList(boundary_id_type nodeset_id) const
3571 : {
3572 : std::map<boundary_id_type, std::vector<dof_id_type>>::const_iterator it =
3573 229 : _node_set_nodes.find(nodeset_id);
3574 :
3575 229 : if (it == _node_set_nodes.end())
3576 : {
3577 : // On a distributed mesh we might not know about a remote nodeset,
3578 : // so we'll return an empty vector and hope the nodeset exists
3579 : // elsewhere.
3580 0 : if (!getMesh().is_serial())
3581 : {
3582 0 : static const std::vector<dof_id_type> empty_vec;
3583 0 : return empty_vec;
3584 : }
3585 : // On a replicated mesh we should know about every nodeset and if
3586 : // we're asked for one that doesn't exist then it must be a bug.
3587 : else
3588 : {
3589 0 : mooseError("Unable to nodeset ID: ", nodeset_id, '.');
3590 : }
3591 : }
3592 :
3593 229 : return it->second;
3594 : }
3595 :
3596 : const std::set<BoundaryID> &
3597 4683485 : MooseMesh::getSubdomainBoundaryIds(const SubdomainID subdomain_id) const
3598 : {
3599 4683485 : const auto it = _sub_to_data.find(subdomain_id);
3600 :
3601 4683485 : if (it == _sub_to_data.end())
3602 0 : mooseError("Unable to find subdomain ID: ", subdomain_id, '.');
3603 :
3604 9366970 : return it->second.boundary_ids;
3605 : }
3606 :
3607 : std::set<BoundaryID>
3608 22 : MooseMesh::getSubdomainInterfaceBoundaryIds(const SubdomainID subdomain_id) const
3609 : {
3610 22 : const auto & bnd_ids = getSubdomainBoundaryIds(subdomain_id);
3611 22 : std::set<BoundaryID> boundary_ids(bnd_ids.begin(), bnd_ids.end());
3612 : std::unordered_map<SubdomainID, std::set<BoundaryID>>::const_iterator it =
3613 22 : _neighbor_subdomain_boundary_ids.find(subdomain_id);
3614 :
3615 22 : boundary_ids.insert(it->second.begin(), it->second.end());
3616 :
3617 44 : return boundary_ids;
3618 0 : }
3619 :
3620 : std::set<SubdomainID>
3621 203 : MooseMesh::getBoundaryConnectedBlocks(const BoundaryID bid) const
3622 : {
3623 203 : std::set<SubdomainID> subdomain_ids;
3624 763 : for (const auto & [sub_id, data] : _sub_to_data)
3625 560 : if (data.boundary_ids.find(bid) != data.boundary_ids.end())
3626 203 : subdomain_ids.insert(sub_id);
3627 :
3628 203 : return subdomain_ids;
3629 0 : }
3630 :
3631 : std::set<SubdomainID>
3632 169 : MooseMesh::getBoundaryConnectedSecondaryBlocks(const BoundaryID bid) const
3633 : {
3634 169 : std::set<SubdomainID> subdomain_ids;
3635 507 : for (const auto & it : _neighbor_subdomain_boundary_ids)
3636 338 : if (it.second.find(bid) != it.second.end())
3637 169 : subdomain_ids.insert(it.first);
3638 :
3639 169 : return subdomain_ids;
3640 0 : }
3641 :
3642 : std::set<SubdomainID>
3643 11 : MooseMesh::getInterfaceConnectedBlocks(const BoundaryID bid) const
3644 : {
3645 11 : std::set<SubdomainID> subdomain_ids = getBoundaryConnectedBlocks(bid);
3646 110 : for (const auto & it : _neighbor_subdomain_boundary_ids)
3647 99 : if (it.second.find(bid) != it.second.end())
3648 44 : subdomain_ids.insert(it.first);
3649 :
3650 11 : return subdomain_ids;
3651 0 : }
3652 :
3653 : const std::set<SubdomainID> &
3654 0 : MooseMesh::getBlockConnectedBlocks(const SubdomainID subdomain_id) const
3655 : {
3656 0 : const auto it = _sub_to_data.find(subdomain_id);
3657 :
3658 0 : if (it == _sub_to_data.end())
3659 0 : mooseError("Unable to find subdomain ID: ", subdomain_id, '.');
3660 :
3661 0 : return it->second.neighbor_subs;
3662 : }
3663 :
3664 : bool
3665 1216264 : MooseMesh::isBoundaryNode(dof_id_type node_id) const
3666 : {
3667 1216264 : bool found_node = false;
3668 4992112 : for (const auto & it : _bnd_node_ids)
3669 : {
3670 4053776 : if (it.second.find(node_id) != it.second.end())
3671 : {
3672 277928 : found_node = true;
3673 277928 : break;
3674 : }
3675 : }
3676 1216264 : return found_node;
3677 : }
3678 :
3679 : bool
3680 995742 : MooseMesh::isBoundaryNode(dof_id_type node_id, BoundaryID bnd_id) const
3681 : {
3682 995742 : bool found_node = false;
3683 995742 : std::map<boundary_id_type, std::set<dof_id_type>>::const_iterator it = _bnd_node_ids.find(bnd_id);
3684 995742 : if (it != _bnd_node_ids.end())
3685 935442 : if (it->second.find(node_id) != it->second.end())
3686 11620 : found_node = true;
3687 995742 : return found_node;
3688 : }
3689 :
3690 : bool
3691 0 : MooseMesh::isBoundaryElem(dof_id_type elem_id) const
3692 : {
3693 0 : bool found_elem = false;
3694 0 : for (const auto & it : _bnd_elem_ids)
3695 : {
3696 0 : if (it.second.find(elem_id) != it.second.end())
3697 : {
3698 0 : found_elem = true;
3699 0 : break;
3700 : }
3701 : }
3702 0 : return found_elem;
3703 : }
3704 :
3705 : bool
3706 425114 : MooseMesh::isBoundaryElem(dof_id_type elem_id, BoundaryID bnd_id) const
3707 : {
3708 425114 : bool found_elem = false;
3709 425114 : auto it = _bnd_elem_ids.find(bnd_id);
3710 425114 : if (it != _bnd_elem_ids.end())
3711 393181 : if (it->second.find(elem_id) != it->second.end())
3712 22342 : found_elem = true;
3713 425114 : return found_elem;
3714 : }
3715 :
3716 : void
3717 1276 : MooseMesh::errorIfDistributedMesh(std::string name) const
3718 : {
3719 1276 : if (_use_distributed_mesh)
3720 0 : mooseError("Cannot use ",
3721 : name,
3722 : " with DistributedMesh!\n",
3723 : "Consider specifying parallel_type = 'replicated' in your input file\n",
3724 : "to prevent it from being run with DistributedMesh.");
3725 1276 : }
3726 :
3727 : void
3728 70165 : MooseMesh::setPartitionerHelper(MeshBase * const mesh)
3729 : {
3730 70165 : if (_use_distributed_mesh && (_partitioner_name != "default" && _partitioner_name != "parmetis"))
3731 : {
3732 16 : _partitioner_name = "parmetis";
3733 16 : _partitioner_overridden = true;
3734 : }
3735 :
3736 70165 : setPartitioner(mesh ? *mesh : getMesh(), _partitioner_name, _use_distributed_mesh, _pars, *this);
3737 70165 : }
3738 :
3739 : void
3740 70165 : MooseMesh::setPartitioner(MeshBase & mesh_base,
3741 : MooseEnum & partitioner,
3742 : bool use_distributed_mesh,
3743 : const InputParameters & params,
3744 : MooseObject & context_obj)
3745 : {
3746 : // Set the partitioner based on partitioner name
3747 70165 : switch (partitioner)
3748 : {
3749 65165 : case -3: // default
3750 : // We'll use the default partitioner, but notify the user of which one is being used...
3751 65165 : if (use_distributed_mesh)
3752 21414 : partitioner = "parmetis";
3753 : else
3754 108916 : partitioner = "metis";
3755 65165 : break;
3756 :
3757 : // No need to explicitily create the metis or parmetis partitioners,
3758 : // They are the default for serial and parallel mesh respectively
3759 4888 : case -2: // metis
3760 : case -1: // parmetis
3761 4888 : break;
3762 :
3763 60 : case 0: // linear
3764 60 : mesh_base.partitioner().reset(new libMesh::LinearPartitioner);
3765 60 : break;
3766 52 : case 1: // centroid
3767 : {
3768 104 : if (!params.isParamValid("centroid_partitioner_direction"))
3769 0 : context_obj.paramError(
3770 : "centroid_partitioner_direction",
3771 : "If using the centroid partitioner you _must_ specify centroid_partitioner_direction!");
3772 :
3773 52 : MooseEnum direction = params.get<MooseEnum>("centroid_partitioner_direction");
3774 :
3775 52 : if (direction == "x")
3776 32 : mesh_base.partitioner().reset(
3777 16 : new libMesh::CentroidPartitioner(libMesh::CentroidPartitioner::X));
3778 36 : else if (direction == "y")
3779 72 : mesh_base.partitioner().reset(
3780 36 : new libMesh::CentroidPartitioner(libMesh::CentroidPartitioner::Y));
3781 0 : else if (direction == "z")
3782 0 : mesh_base.partitioner().reset(
3783 0 : new libMesh::CentroidPartitioner(libMesh::CentroidPartitioner::Z));
3784 0 : else if (direction == "radial")
3785 0 : mesh_base.partitioner().reset(
3786 0 : new libMesh::CentroidPartitioner(libMesh::CentroidPartitioner::RADIAL));
3787 52 : break;
3788 52 : }
3789 0 : case 2: // hilbert_sfc
3790 0 : mesh_base.partitioner().reset(new libMesh::HilbertSFCPartitioner);
3791 0 : break;
3792 0 : case 3: // morton_sfc
3793 0 : mesh_base.partitioner().reset(new libMesh::MortonSFCPartitioner);
3794 0 : break;
3795 : }
3796 70165 : }
3797 :
3798 : void
3799 1497 : MooseMesh::setCustomPartitioner(Partitioner * partitioner)
3800 : {
3801 1497 : _custom_partitioner = partitioner->clone();
3802 1497 : setIsCustomPartitionerRequested(true);
3803 1497 : if (_mesh)
3804 12 : _mesh->partitioner() = _custom_partitioner->clone();
3805 1497 : _partitioner_name = "custom";
3806 1497 : }
3807 :
3808 : bool
3809 0 : MooseMesh::isCustomPartitionerRequested() const
3810 : {
3811 0 : return _custom_partitioner_requested;
3812 : }
3813 :
3814 : bool
3815 146491 : MooseMesh::hasSecondOrderElements()
3816 : {
3817 146491 : bool mesh_has_second_order_elements = false;
3818 46339193 : for (auto it = activeLocalElementsBegin(), end = activeLocalElementsEnd(); it != end; ++it)
3819 23112963 : if ((*it)->default_order() == SECOND)
3820 : {
3821 16612 : mesh_has_second_order_elements = true;
3822 16612 : break;
3823 146491 : }
3824 :
3825 : // We checked our local elements, so take the max over all processors.
3826 146491 : comm().max(mesh_has_second_order_elements);
3827 146491 : return mesh_has_second_order_elements;
3828 : }
3829 :
3830 : void
3831 3001 : MooseMesh::setIsCustomPartitionerRequested(bool cpr)
3832 : {
3833 3001 : _custom_partitioner_requested = cpr;
3834 3001 : }
3835 :
3836 : std::unique_ptr<libMesh::PointLocatorBase>
3837 7000 : MooseMesh::getPointLocator() const
3838 : {
3839 7000 : return getMesh().sub_point_locator();
3840 : }
3841 :
3842 : void
3843 4742 : MooseMesh::buildFiniteVolumeInfo() const
3844 : {
3845 : mooseAssert(!Threads::in_threads,
3846 : "This routine has not been implemented for threads. Please query this routine before "
3847 : "a threaded region or contact a MOOSE developer to discuss.");
3848 4742 : _finite_volume_info_dirty = false;
3849 :
3850 : using Keytype = std::pair<const Elem *, unsigned short int>;
3851 :
3852 : // create a map from elem/side --> boundary ids
3853 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> side_list =
3854 4742 : buildActiveSideList();
3855 4742 : std::map<Keytype, std::set<boundary_id_type>> side_map;
3856 164212 : for (auto & [elem_id, side, bc_id] : side_list)
3857 : {
3858 159470 : const Elem * elem = _mesh->elem_ptr(elem_id);
3859 159470 : Keytype key(elem, side);
3860 159470 : auto & bc_set = side_map[key];
3861 159470 : bc_set.insert(bc_id);
3862 : }
3863 :
3864 4742 : _face_info.clear();
3865 4742 : _all_face_info.clear();
3866 4742 : _elem_side_to_face_info.clear();
3867 :
3868 4742 : _elem_to_elem_info.clear();
3869 4742 : _elem_info.clear();
3870 :
3871 : // by performing the element ID comparison check in the below loop, we are ensuring that we never
3872 : // double count face contributions. If a face lies along a process boundary, the only process that
3873 : // will contribute to both sides of the face residuals/Jacobians will be the process that owns the
3874 : // element with the lower ID.
3875 4742 : auto begin = getMesh().active_elements_begin();
3876 4742 : auto end = getMesh().active_elements_end();
3877 :
3878 : // We prepare a map connecting the Elem* and the corresponding ElemInfo
3879 : // for the active elements.
3880 4742 : _elem_to_elem_info.reserve(nActiveLocalElem());
3881 4742 : unsigned int num_sides = 0;
3882 1137318 : for (const Elem * elem : as_range(begin, end))
3883 : {
3884 1132576 : _elem_to_elem_info.emplace(elem->id(), elem);
3885 1132576 : num_sides += elem->n_sides();
3886 4742 : }
3887 :
3888 : // Used to speed up FaceInfo creation:
3889 : // - element side builder that caches per type of element
3890 4742 : libMesh::ElemSideBuilder side_builder;
3891 :
3892 4742 : _all_face_info.reserve(num_sides / 2);
3893 4742 : dof_id_type face_index = 0;
3894 2269894 : for (const Elem * elem : as_range(begin, end))
3895 : {
3896 5035168 : for (unsigned int side = 0; side < elem->n_sides(); ++side)
3897 : {
3898 : // get the neighbor element
3899 3902592 : const Elem * neighbor = elem->neighbor_ptr(side);
3900 :
3901 : // Check if the FaceInfo shall belong to the element. If yes,
3902 : // create and initialize the FaceInfo. We need this to ensure that
3903 : // we do not duplicate FaceInfo-s.
3904 3902592 : if (Moose::FV::elemHasFaceInfo(*elem, neighbor))
3905 : {
3906 : mooseAssert(!neighbor || (neighbor->level() < elem->level() ? neighbor->active() : true),
3907 : "If the neighbor is coarser than the element, we expect that the neighbor must "
3908 : "be active.");
3909 :
3910 : // We construct the faceInfo using the elementinfo and side index
3911 : mooseAssert(elem->default_order() < 4, "Did not expect such high element orders in FV");
3912 4057954 : _all_face_info.emplace_back(
3913 2028977 : &_elem_to_elem_info[elem->id()], side, face_index++, side_builder);
3914 :
3915 2028977 : auto & fi = _all_face_info.back();
3916 :
3917 : // get all the sidesets that this face is contained in and cache them
3918 : // in the face info.
3919 2028977 : std::set<boundary_id_type> & boundary_ids = fi.boundaryIDs();
3920 2028977 : boundary_ids.clear();
3921 :
3922 : // We initialize the weights/other information in faceInfo. If the neighbor does not exist
3923 : // or is remote (so when we are on some sort of mesh boundary), we initialize the ghost
3924 : // cell and use it to compute the weights corresponding to the faceInfo.
3925 2028977 : if (!neighbor || neighbor == libMesh::remote_elem)
3926 152539 : fi.computeBoundaryCoefficients();
3927 : else
3928 1876438 : fi.computeInternalCoefficients(&_elem_to_elem_info[neighbor->id()]);
3929 :
3930 2028977 : auto lit = side_map.find(Keytype(&fi.elem(), fi.elemSideID()));
3931 2028977 : if (lit != side_map.end())
3932 152633 : boundary_ids.insert(lit->second.begin(), lit->second.end());
3933 :
3934 2028977 : if (fi.neighborPtr())
3935 : {
3936 1876438 : auto rit = side_map.find(Keytype(fi.neighborPtr(), fi.neighborSideID()));
3937 1876438 : if (rit != side_map.end())
3938 4220 : boundary_ids.insert(rit->second.begin(), rit->second.end());
3939 : }
3940 : }
3941 : }
3942 4742 : }
3943 :
3944 : // Build the local face info and elem_side to face info maps. We need to do this after
3945 : // _all_face_info is finished being constructed because emplace_back invalidates all iterators and
3946 : // references if ever the new size exceeds capacity
3947 4742 : _elem_side_to_face_info.reserve(_all_face_info.size());
3948 : // heuristic to avoid resizing too much
3949 4742 : _face_info.reserve(_all_face_info.size());
3950 2033719 : for (auto & fi : _all_face_info)
3951 : {
3952 2028977 : const Elem * const elem = &fi.elem();
3953 2028977 : const auto side = fi.elemSideID();
3954 :
3955 : #ifndef NDEBUG
3956 : auto pair_it =
3957 : #endif
3958 2028977 : _elem_side_to_face_info.emplace(std::make_pair(elem, side), &fi);
3959 : mooseAssert(pair_it.second, "We should be adding unique FaceInfo objects.");
3960 :
3961 : // We will add the faces on processor boundaries to the list of face infos on each
3962 : // associated processor.
3963 2583648 : if (fi.elem().processor_id() == this->processor_id() ||
3964 554671 : (fi.neighborPtr() && (fi.neighborPtr()->processor_id() == this->processor_id())))
3965 1745141 : _face_info.push_back(&fi);
3966 : }
3967 :
3968 4742 : _elem_info.reserve(nActiveLocalElem());
3969 1137318 : for (auto & ei : _elem_to_elem_info)
3970 1132576 : if (ei.second.elem()->processor_id() == this->processor_id())
3971 978244 : _elem_info.push_back(&ei.second);
3972 4742 : }
3973 :
3974 : const FaceInfo *
3975 122662892 : MooseMesh::faceInfo(const Elem * elem, unsigned int side) const
3976 : {
3977 122662892 : auto it = _elem_side_to_face_info.find(std::make_pair(elem, side));
3978 :
3979 122662892 : if (it == _elem_side_to_face_info.end())
3980 792 : return nullptr;
3981 : else
3982 : {
3983 : mooseAssert(it->second,
3984 : "For some reason, the FaceInfo object is NULL! Try calling "
3985 : "`buildFiniteVolumeInfo()` before using this accessor!");
3986 122662100 : return it->second;
3987 : }
3988 : }
3989 :
3990 : const ElemInfo &
3991 108299835 : MooseMesh::elemInfo(const dof_id_type id) const
3992 : {
3993 108299835 : return libmesh_map_find(_elem_to_elem_info, id);
3994 : }
3995 :
3996 : void
3997 4722 : MooseMesh::computeFiniteVolumeCoords() const
3998 : {
3999 4722 : if (_finite_volume_info_dirty)
4000 0 : mooseError("Trying to compute face- and elem-info coords when the information is dirty");
4001 :
4002 2032979 : for (auto & fi : _all_face_info)
4003 : {
4004 : // get elem & neighbor elements, and set subdomain ids
4005 2028257 : const SubdomainID elem_subdomain_id = fi.elemSubdomainID();
4006 2028257 : const SubdomainID neighbor_subdomain_id = fi.neighborSubdomainID();
4007 :
4008 2028257 : coordTransformFactor(
4009 2028257 : *this, elem_subdomain_id, fi.faceCentroid(), fi.faceCoord(), neighbor_subdomain_id);
4010 : }
4011 :
4012 1137138 : for (auto & ei : _elem_to_elem_info)
4013 1132416 : coordTransformFactor(
4014 2264832 : *this, ei.second.subdomain_id(), ei.second.centroid(), ei.second.coordFactor());
4015 4722 : }
4016 :
4017 : MooseEnum
4018 204250 : MooseMesh::partitioning()
4019 : {
4020 : MooseEnum partitioning(
4021 612750 : "default=-3 metis=-2 parmetis=-1 linear=0 centroid hilbert_sfc morton_sfc custom", "default");
4022 204250 : return partitioning;
4023 : }
4024 :
4025 : MooseEnum
4026 3425 : MooseMesh::elemTypes()
4027 : {
4028 : MooseEnum elemTypes(
4029 : "EDGE EDGE2 EDGE3 EDGE4 QUAD QUAD4 QUAD8 QUAD9 TRI3 TRI6 HEX HEX8 HEX20 HEX27 TET4 TET10 "
4030 10275 : "PRISM6 PRISM15 PRISM18 PYRAMID5 PYRAMID13 PYRAMID14");
4031 3425 : return elemTypes;
4032 : }
4033 :
4034 : void
4035 35755 : MooseMesh::allowRemoteElementRemoval(const bool allow_remote_element_removal)
4036 : {
4037 35755 : _allow_remote_element_removal = allow_remote_element_removal;
4038 35755 : if (_mesh)
4039 16324 : _mesh->allow_remote_element_removal(allow_remote_element_removal);
4040 :
4041 35755 : if (!allow_remote_element_removal)
4042 : // If we're not allowing remote element removal now, then we will need deletion later after
4043 : // late geoemetric ghosting functors have been added (late geometric ghosting functor addition
4044 : // happens when algebraic ghosting functors are added)
4045 35755 : _need_delete = true;
4046 35755 : }
4047 :
4048 : void
4049 17623 : MooseMesh::deleteRemoteElements()
4050 : {
4051 17623 : _allow_remote_element_removal = true;
4052 17623 : if (!_mesh)
4053 0 : mooseError("Cannot delete remote elements because we have not yet attached a MeshBase");
4054 :
4055 17623 : _mesh->allow_remote_element_removal(true);
4056 :
4057 17623 : _mesh->delete_remote_elements();
4058 17623 : }
4059 :
4060 : void
4061 4718 : MooseMesh::cacheFaceInfoVariableOwnership() const
4062 : {
4063 : mooseAssert(
4064 : !Threads::in_threads,
4065 : "Performing writes to faceInfo variable association maps. This must be done unthreaded!");
4066 :
4067 4718 : const unsigned int num_eqs = _app.feProblem().es().n_systems();
4068 :
4069 4060395 : auto face_lambda = [this](const SubdomainID elem_subdomain_id,
4070 : const SubdomainID neighbor_subdomain_id,
4071 : SystemBase & sys,
4072 : std::vector<std::vector<FaceInfo::VarFaceNeighbors>> & face_type_vector)
4073 : {
4074 4060395 : face_type_vector[sys.number()].resize(sys.nVariables(), FaceInfo::VarFaceNeighbors::NEITHER);
4075 4060395 : const auto & variables = sys.getVariables(0);
4076 :
4077 8732670 : for (const auto & var : variables)
4078 : {
4079 4672275 : const unsigned int var_num = var->number();
4080 4672275 : const unsigned int sys_num = var->sys().number();
4081 4672275 : std::set<SubdomainID> var_subdomains = var->blockIDs();
4082 : /**
4083 : * The following paragraph of code assigns the VarFaceNeighbors
4084 : * 1. The face is an internal face of this variable if it is defined on
4085 : * the elem and neighbor subdomains
4086 : * 2. The face is an invalid face of this variable if it is neither defined
4087 : * on the elem nor the neighbor subdomains
4088 : * 3. If not 1. or 2. then this is a boundary for this variable and the else clause
4089 : * applies
4090 : */
4091 4672275 : bool var_defined_elem = var_subdomains.find(elem_subdomain_id) != var_subdomains.end();
4092 : bool var_defined_neighbor =
4093 4672275 : var_subdomains.find(neighbor_subdomain_id) != var_subdomains.end();
4094 4672275 : if (var_defined_elem && var_defined_neighbor)
4095 3885061 : face_type_vector[sys_num][var_num] = FaceInfo::VarFaceNeighbors::BOTH;
4096 787214 : else if (!var_defined_elem && !var_defined_neighbor)
4097 349377 : face_type_vector[sys_num][var_num] = FaceInfo::VarFaceNeighbors::NEITHER;
4098 : else
4099 : {
4100 : // this is a boundary face for this variable, set elem or neighbor
4101 437837 : if (var_defined_elem)
4102 432405 : face_type_vector[sys_num][var_num] = FaceInfo::VarFaceNeighbors::ELEM;
4103 5432 : else if (var_defined_neighbor)
4104 5432 : face_type_vector[sys_num][var_num] = FaceInfo::VarFaceNeighbors::NEIGHBOR;
4105 : else
4106 0 : mooseError("Should never get here");
4107 : }
4108 4672275 : }
4109 4060395 : };
4110 :
4111 : // We loop through the faces and check if they are internal, boundary or external to
4112 : // the variables in the problem
4113 2032927 : for (FaceInfo & face : _all_face_info)
4114 : {
4115 2028209 : const SubdomainID elem_subdomain_id = face.elemSubdomainID();
4116 2028209 : const SubdomainID neighbor_subdomain_id = face.neighborSubdomainID();
4117 :
4118 2028209 : auto & face_type_vector = face.faceType();
4119 :
4120 2028209 : face_type_vector.clear();
4121 2028209 : face_type_vector.resize(num_eqs);
4122 :
4123 : // First, we check the variables in the solver systems (linear/nonlinear)
4124 4060395 : for (const auto i : make_range(_app.feProblem().numSolverSystems()))
4125 2032186 : face_lambda(elem_subdomain_id,
4126 : neighbor_subdomain_id,
4127 2032186 : _app.feProblem().getSolverSystem(i),
4128 : face_type_vector);
4129 :
4130 : // Then we check the variables in the auxiliary system
4131 2028209 : face_lambda(elem_subdomain_id,
4132 : neighbor_subdomain_id,
4133 2028209 : _app.feProblem().getAuxiliarySystem(),
4134 : face_type_vector);
4135 : }
4136 4718 : }
4137 :
4138 : void
4139 4718 : MooseMesh::cacheFVElementalDoFs() const
4140 : {
4141 : mooseAssert(!Threads::in_threads,
4142 : "Performing writes to elemInfo dof indices. This must be done unthreaded!");
4143 :
4144 2267443 : auto elem_lambda = [](const ElemInfo & elem_info,
4145 : SystemBase & sys,
4146 : std::vector<std::vector<dof_id_type>> & dof_vector)
4147 : {
4148 2267443 : if (sys.nFVVariables())
4149 : {
4150 1199180 : dof_vector[sys.number()].resize(sys.nVariables(), libMesh::DofObject::invalid_id);
4151 1199180 : const auto & variables = sys.getVariables(0);
4152 :
4153 3685546 : for (const auto & var : variables)
4154 2486366 : if (var->isFV())
4155 : {
4156 1440695 : const auto & var_subdomains = var->blockIDs();
4157 :
4158 : // We will only cache for FV variables and if they live on the current subdomain
4159 1440695 : if (var_subdomains.find(elem_info.subdomain_id()) != var_subdomains.end())
4160 : {
4161 1333516 : std::vector<dof_id_type> indices;
4162 1333516 : var->dofMap().dof_indices(elem_info.elem(), indices, var->number());
4163 : mooseAssert(indices.size() == 1, "We expect to have only one dof per element!");
4164 1333516 : dof_vector[sys.number()][var->number()] = indices[0];
4165 1333516 : }
4166 : }
4167 : }
4168 2267443 : };
4169 :
4170 4718 : const unsigned int num_eqs = _app.feProblem().es().n_systems();
4171 :
4172 : // We loop through the elements in the mesh and cache the dof indices
4173 : // for the corresponding variables.
4174 1137118 : for (auto & ei_pair : _elem_to_elem_info)
4175 : {
4176 1132400 : auto & elem_info = ei_pair.second;
4177 1132400 : auto & dof_vector = elem_info.dofIndices();
4178 :
4179 1132400 : dof_vector.clear();
4180 1132400 : dof_vector.resize(num_eqs);
4181 :
4182 : // First, we cache the dof indices for the variables in the solver systems (linear, nonlinear)
4183 2267443 : for (const auto i : make_range(_app.feProblem().numSolverSystems()))
4184 1135043 : elem_lambda(elem_info, _app.feProblem().getSolverSystem(i), dof_vector);
4185 :
4186 : // Then we cache the dof indices for the auxvariables
4187 1132400 : elem_lambda(elem_info, _app.feProblem().getAuxiliarySystem(), dof_vector);
4188 : }
4189 4718 : }
4190 :
4191 : void
4192 4718 : MooseMesh::setupFiniteVolumeMeshData() const
4193 : {
4194 4718 : buildFiniteVolumeInfo();
4195 4718 : computeFiniteVolumeCoords();
4196 4718 : cacheFaceInfoVariableOwnership();
4197 4718 : cacheFVElementalDoFs();
4198 4718 : }
4199 :
4200 : void
4201 65884 : MooseMesh::setCoordSystem(const std::vector<SubdomainName> & blocks,
4202 : const MultiMooseEnum & coord_sys)
4203 : {
4204 329420 : TIME_SECTION("setCoordSystem", 5, "Setting Coordinate System");
4205 65884 : if (!_provided_coord_blocks.empty() && (_provided_coord_blocks != blocks))
4206 : {
4207 0 : const std::string param_name = isParamValid("coord_block") ? "coord_block" : "block";
4208 0 : mooseWarning("Supplied blocks in the 'setCoordSystem' method do not match the value of the "
4209 : "'Mesh/",
4210 : param_name,
4211 : "' parameter. Did you provide different parameter values for 'Mesh/",
4212 : param_name,
4213 : "' and 'Problem/block'?. We will honor the parameter value from 'Mesh/",
4214 : param_name,
4215 : "'");
4216 : mooseAssert(_coord_system_set,
4217 : "If we are arriving here due to a bad specification in the Problem block, then we "
4218 : "should have already set our coordinate system subdomains from the Mesh block");
4219 0 : return;
4220 0 : }
4221 198810 : if (_pars.isParamSetByUser("coord_type") && getParam<MultiMooseEnum>("coord_type") != coord_sys)
4222 0 : mooseError("Supplied coordinate systems in the 'setCoordSystem' method do not match the value "
4223 : "of the 'Mesh/coord_type' parameter. Did you provide different parameter values for "
4224 : "'coord_type' to 'Mesh' and 'Problem'?");
4225 :
4226 : // If blocks contain ANY_BLOCK_ID, it should be the only block specified, and coord_sys should
4227 : // have one and only one entry. In that case, the same coordinate system will be set for all
4228 : // subdomains.
4229 65884 : if (blocks.size() == 1 && blocks[0] == "ANY_BLOCK_ID")
4230 : {
4231 0 : if (coord_sys.size() > 1)
4232 0 : mooseError("If you specify ANY_BLOCK_ID as the only block, you must also specify a single "
4233 : "coordinate system for it.");
4234 0 : if (!_mesh->is_prepared())
4235 0 : mooseError(
4236 : "You cannot set the coordinate system for ANY_BLOCK_ID before the mesh is prepared. "
4237 : "Please call this method after the mesh is prepared.");
4238 0 : const auto coord_type = coord_sys.size() == 0
4239 0 : ? Moose::COORD_XYZ
4240 0 : : Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[0]);
4241 0 : for (const auto sid : meshSubdomains())
4242 0 : _coord_sys[sid] = coord_type;
4243 0 : return;
4244 : }
4245 :
4246 : // If multiple blocks are specified, but one of them is ANY_BLOCK_ID, let's emit a helpful error
4247 65884 : if (std::find(blocks.begin(), blocks.end(), "ANY_BLOCK_ID") != blocks.end())
4248 0 : mooseError("You cannot specify ANY_BLOCK_ID together with other blocks in the "
4249 : "setCoordSystem() method. If you want to set the same coordinate system for all "
4250 : "blocks, use ANY_BLOCK_ID as the only block.");
4251 :
4252 65884 : auto subdomains = meshSubdomains();
4253 : // It's possible that a user has called this API before the mesh is prepared and consequently we
4254 : // don't yet have the subdomains in meshSubdomains()
4255 66283 : for (const auto & sub_name : blocks)
4256 : {
4257 399 : const auto sub_id = getSubdomainID(sub_name);
4258 399 : subdomains.insert(sub_id);
4259 : }
4260 :
4261 65884 : if (coord_sys.size() <= 1)
4262 : {
4263 : // We will specify the same coordinate system for all blocks
4264 65860 : const auto coord_type = coord_sys.size() == 0
4265 65860 : ? Moose::COORD_XYZ
4266 65860 : : Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[0]);
4267 158221 : for (const auto sid : subdomains)
4268 92361 : _coord_sys[sid] = coord_type;
4269 : }
4270 : else
4271 : {
4272 24 : if (blocks.size() != coord_sys.size())
4273 0 : mooseError("Number of blocks and coordinate systems does not match.");
4274 :
4275 96 : for (const auto i : index_range(blocks))
4276 : {
4277 72 : SubdomainID sid = getSubdomainID(blocks[i]);
4278 : Moose::CoordinateSystemType coord_type =
4279 72 : Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[i]);
4280 72 : _coord_sys[sid] = coord_type;
4281 : }
4282 :
4283 96 : for (const auto & sid : subdomains)
4284 72 : if (_coord_sys.find(sid) == _coord_sys.end())
4285 0 : mooseError("Subdomain '" + Moose::stringify(sid) +
4286 : "' does not have a coordinate system specified.");
4287 : }
4288 :
4289 65884 : _coord_system_set = true;
4290 :
4291 65884 : updateCoordTransform();
4292 65884 : }
4293 :
4294 : Moose::CoordinateSystemType
4295 2291096004 : MooseMesh::getCoordSystem(SubdomainID sid) const
4296 : {
4297 2291096004 : auto it = _coord_sys.find(sid);
4298 2291096004 : if (it != _coord_sys.end())
4299 4582192008 : return (*it).second;
4300 : else
4301 0 : mooseError("Requested subdomain ", sid, " does not exist.");
4302 : }
4303 :
4304 : Moose::CoordinateSystemType
4305 55126 : MooseMesh::getUniqueCoordSystem() const
4306 : {
4307 55126 : const auto unique_system = _coord_sys.find(*meshSubdomains().begin())->second;
4308 : // Check that it is actually unique
4309 55126 : bool result = std::all_of(
4310 55126 : std::next(_coord_sys.begin()),
4311 55126 : _coord_sys.end(),
4312 4346 : [unique_system](
4313 : typename std::unordered_map<SubdomainID, Moose::CoordinateSystemType>::const_reference
4314 4346 : item) { return (item.second == unique_system); });
4315 55126 : if (!result)
4316 0 : mooseError("The unique coordinate system of the mesh was requested by the mesh contains "
4317 : "multiple blocks with different coordinate systems");
4318 :
4319 55126 : if (usingGeneralAxisymmetricCoordAxes())
4320 0 : mooseError("General axisymmetric coordinate axes are being used, and it is currently "
4321 : "conservatively assumed that in this case there is no unique coordinate system.");
4322 :
4323 55126 : return unique_system;
4324 : }
4325 :
4326 : const std::map<SubdomainID, Moose::CoordinateSystemType> &
4327 68888 : MooseMesh::getCoordSystem() const
4328 : {
4329 68888 : return _coord_sys;
4330 : }
4331 :
4332 : void
4333 0 : MooseMesh::setAxisymmetricCoordAxis(const MooseEnum & rz_coord_axis)
4334 : {
4335 0 : _rz_coord_axis = rz_coord_axis;
4336 :
4337 0 : updateCoordTransform();
4338 0 : }
4339 :
4340 : void
4341 17 : MooseMesh::setGeneralAxisymmetricCoordAxes(
4342 : const std::vector<SubdomainName> & blocks,
4343 : const std::vector<std::pair<Point, RealVectorValue>> & axes)
4344 : {
4345 : // Set the axes for the given blocks
4346 : mooseAssert(blocks.size() == axes.size(), "Blocks and axes vectors must be the same length.");
4347 58 : for (const auto i : index_range(blocks))
4348 : {
4349 41 : const auto subdomain_id = getSubdomainID(blocks[i]);
4350 41 : const auto it = _coord_sys.find(subdomain_id);
4351 41 : if (it == _coord_sys.end())
4352 0 : mooseError("The block '",
4353 0 : blocks[i],
4354 : "' has not set a coordinate system. Make sure to call setCoordSystem() before "
4355 : "setGeneralAxisymmetricCoordAxes().");
4356 : else
4357 : {
4358 41 : if (it->second == Moose::COORD_RZ)
4359 : {
4360 41 : const auto direction = axes[i].second;
4361 41 : if (direction.is_zero())
4362 0 : mooseError("Only nonzero vectors may be supplied for RZ directions.");
4363 :
4364 41 : _subdomain_id_to_rz_coord_axis[subdomain_id] =
4365 82 : std::make_pair(axes[i].first, direction.unit());
4366 : }
4367 : else
4368 0 : mooseError("The block '",
4369 0 : blocks[i],
4370 : "' was provided in setGeneralAxisymmetricCoordAxes(), but the coordinate system "
4371 : "for this block is not 'RZ'.");
4372 : }
4373 : }
4374 :
4375 : // Make sure there are no RZ blocks that still do not have axes
4376 17 : const auto all_subdomain_ids = meshSubdomains();
4377 70 : for (const auto subdomain_id : all_subdomain_ids)
4378 94 : if (getCoordSystem(subdomain_id) == Moose::COORD_RZ &&
4379 41 : !_subdomain_id_to_rz_coord_axis.count(subdomain_id))
4380 0 : mooseError("The block '",
4381 0 : getSubdomainName(subdomain_id),
4382 : "' was specified to use the 'RZ' coordinate system but was not given in "
4383 : "setGeneralAxisymmetricCoordAxes().");
4384 :
4385 17 : updateCoordTransform();
4386 17 : }
4387 :
4388 : const std::pair<Point, RealVectorValue> &
4389 1041955 : MooseMesh::getGeneralAxisymmetricCoordAxis(SubdomainID subdomain_id) const
4390 : {
4391 1041955 : auto it = _subdomain_id_to_rz_coord_axis.find(subdomain_id);
4392 1041955 : if (it != _subdomain_id_to_rz_coord_axis.end())
4393 2083910 : return (*it).second;
4394 : else
4395 0 : mooseError("Requested subdomain ", subdomain_id, " does not exist.");
4396 : }
4397 :
4398 : bool
4399 24298712 : MooseMesh::usingGeneralAxisymmetricCoordAxes() const
4400 : {
4401 24298712 : return _subdomain_id_to_rz_coord_axis.size() > 0;
4402 : }
4403 :
4404 : void
4405 68888 : MooseMesh::updateCoordTransform()
4406 : {
4407 68888 : if (!_coord_transform)
4408 68867 : _coord_transform = std::make_unique<MooseAppCoordTransform>(*this);
4409 : else
4410 21 : _coord_transform->setCoordinateSystem(*this);
4411 68888 : }
4412 :
4413 : unsigned int
4414 20684596 : MooseMesh::getAxisymmetricRadialCoord() const
4415 : {
4416 20684596 : if (usingGeneralAxisymmetricCoordAxes())
4417 0 : mooseError("getAxisymmetricRadialCoord() should not be called if "
4418 : "setGeneralAxisymmetricCoordAxes() has been called.");
4419 :
4420 20684596 : if (_rz_coord_axis == 0)
4421 133200 : return 1; // if the rotation axis is x (0), then the radial direction is y (1)
4422 : else
4423 20551396 : return 0; // otherwise the radial direction is assumed to be x, i.e., the rotation axis is y
4424 : }
4425 :
4426 : void
4427 61282 : MooseMesh::checkCoordinateSystems()
4428 : {
4429 27928712 : for (const auto & elem : getMesh().element_ptr_range())
4430 : {
4431 13933718 : SubdomainID sid = elem->subdomain_id();
4432 13933718 : if (_coord_sys[sid] == Moose::COORD_RZ && elem->dim() == 3)
4433 3 : mooseError("An RZ coordinate system was requested for subdomain " + Moose::stringify(sid) +
4434 : " which contains 3D elements.");
4435 13933715 : if (_coord_sys[sid] == Moose::COORD_RSPHERICAL && elem->dim() > 1)
4436 0 : mooseError("An RSPHERICAL coordinate system was requested for subdomain " +
4437 0 : Moose::stringify(sid) + " which contains 2D or 3D elements.");
4438 61279 : }
4439 61279 : }
4440 :
4441 : void
4442 2022 : MooseMesh::setCoordData(const MooseMesh & other_mesh)
4443 : {
4444 2022 : _coord_sys = other_mesh._coord_sys;
4445 2022 : _rz_coord_axis = other_mesh._rz_coord_axis;
4446 2022 : _subdomain_id_to_rz_coord_axis = other_mesh._subdomain_id_to_rz_coord_axis;
4447 2022 : }
4448 :
4449 : const MooseUnits &
4450 2 : MooseMesh::lengthUnit() const
4451 : {
4452 : mooseAssert(_coord_transform, "This must be non-null");
4453 2 : return _coord_transform->lengthUnit();
4454 : }
4455 :
4456 : void
4457 68769 : MooseMesh::checkDuplicateSubdomainNames()
4458 : {
4459 68769 : std::map<SubdomainName, SubdomainID> subdomain;
4460 165270 : for (const auto & sbd_id : _mesh_subdomains)
4461 : {
4462 96504 : std::string sub_name = getSubdomainName(sbd_id);
4463 96504 : if (!sub_name.empty() && subdomain.count(sub_name) > 0)
4464 6 : mooseError("The subdomain name ",
4465 : sub_name,
4466 : " is used for both subdomain with ID=",
4467 3 : subdomain[sub_name],
4468 : " and ID=",
4469 : sbd_id,
4470 : ", Please rename one of them!");
4471 : else
4472 96501 : subdomain[sub_name] = sbd_id;
4473 96501 : }
4474 68766 : }
4475 :
4476 : const std::vector<QpMap> &
4477 800 : MooseMesh::getPRefinementMapHelper(
4478 : const Elem & elem,
4479 : const std::map<std::pair<ElemType, unsigned int>, std::vector<QpMap>> & map) const
4480 : {
4481 : // We are actually seeking the map stored with the p_level - 1 key, e.g. the refinement map that
4482 : // maps from the previous p_level to this element's p_level
4483 800 : return libmesh_map_find(map,
4484 : std::make_pair(elem.type(), cast_int<unsigned int>(elem.p_level() - 1)));
4485 : }
4486 :
4487 : const std::vector<QpMap> &
4488 0 : MooseMesh::getPCoarseningMapHelper(
4489 : const Elem & elem,
4490 : const std::map<std::pair<ElemType, unsigned int>, std::vector<QpMap>> & map) const
4491 : {
4492 : mooseAssert(elem.active() && elem.p_refinement_flag() == Elem::JUST_COARSENED,
4493 : "These are the conditions that should be met for requesting a coarsening map");
4494 0 : return libmesh_map_find(map, std::make_pair(elem.type(), elem.p_level()));
4495 : }
4496 :
4497 : const std::vector<QpMap> &
4498 800 : MooseMesh::getPRefinementMap(const Elem & elem) const
4499 : {
4500 800 : return getPRefinementMapHelper(elem, _elem_type_to_p_refinement_map);
4501 : }
4502 :
4503 : const std::vector<QpMap> &
4504 0 : MooseMesh::getPRefinementSideMap(const Elem & elem) const
4505 : {
4506 0 : return getPRefinementMapHelper(elem, _elem_type_to_p_refinement_side_map);
4507 : }
4508 :
4509 : const std::vector<QpMap> &
4510 0 : MooseMesh::getPCoarseningMap(const Elem & elem) const
4511 : {
4512 0 : return getPCoarseningMapHelper(elem, _elem_type_to_p_coarsening_map);
4513 : }
4514 :
4515 : const std::vector<QpMap> &
4516 0 : MooseMesh::getPCoarseningSideMap(const Elem & elem) const
4517 : {
4518 0 : return getPCoarseningMapHelper(elem, _elem_type_to_p_coarsening_side_map);
4519 : }
4520 :
4521 : bool
4522 27055 : MooseMesh::skipNoncriticalPartitioning() const
4523 : {
4524 27055 : return _mesh->skip_noncritical_partitioning();
4525 : }
|
