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