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