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