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