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