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