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