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 : // MOOSE includes
11 : #include "MooseMeshUtils.h"
12 :
13 : #include "libmesh/elem.h"
14 : #include "libmesh/boundary_info.h"
15 : #include "libmesh/id_types.h"
16 : #include "libmesh/int_range.h"
17 : #include "libmesh/parallel.h"
18 : #include "libmesh/parallel_algebra.h"
19 : #include "libmesh/utility.h"
20 :
21 : #include "libmesh/distributed_mesh.h"
22 : #include "libmesh/parallel_elem.h"
23 : #include "libmesh/parallel_node.h"
24 : #include "libmesh/compare_elems_by_level.h"
25 : #include "libmesh/mesh_communication.h"
26 : #include "libmesh/edge_edge3.h"
27 : #include "libmesh/enum_to_string.h"
28 : #include "libmesh/unstructured_mesh.h"
29 :
30 : #include "timpi/parallel_sync.h"
31 :
32 : using namespace libMesh;
33 :
34 : namespace MooseMeshUtils
35 : {
36 :
37 : void
38 410 : mergeBoundaryIDsWithSameName(MeshBase & mesh)
39 : {
40 : // We check if we have the same boundary name with different IDs. If we do, we assign the
41 : // first ID to every occurrence.
42 410 : const auto & side_bd_name_map = mesh.get_boundary_info().get_sideset_name_map();
43 410 : const auto & node_bd_name_map = mesh.get_boundary_info().get_nodeset_name_map();
44 410 : std::map<boundary_id_type, boundary_id_type> same_name_ids;
45 :
46 820 : auto populate_map = [](const std::map<boundary_id_type, std::string> & map,
47 : std::map<boundary_id_type, boundary_id_type> & same_ids)
48 : {
49 4036 : for (const auto & pair_outer : map)
50 26256 : for (const auto & pair_inner : map)
51 : // The last condition is needed to make sure we only store one combination
52 24240 : if (pair_outer.second == pair_inner.second && pair_outer.first != pair_inner.first &&
53 24240 : same_ids.find(pair_inner.first) == same_ids.end())
54 600 : same_ids[pair_outer.first] = pair_inner.first;
55 820 : };
56 :
57 410 : populate_map(side_bd_name_map, same_name_ids);
58 410 : populate_map(node_bd_name_map, same_name_ids);
59 :
60 742 : for (const auto & [id1, id2] : same_name_ids)
61 332 : mesh.get_boundary_info().renumber_id(id2, id1);
62 410 : }
63 :
64 : void
65 508 : changeBoundaryId(MeshBase & mesh,
66 : const boundary_id_type old_id,
67 : const boundary_id_type new_id,
68 : bool delete_prev)
69 : {
70 : // Get a reference to our BoundaryInfo object, we will use it several times below...
71 508 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
72 :
73 : // Container to catch ids passed back from BoundaryInfo
74 508 : std::vector<boundary_id_type> old_ids;
75 :
76 : // Only level-0 elements store BCs. Loop over them.
77 562116 : for (auto & elem : as_range(mesh.level_elements_begin(0), mesh.level_elements_end(0)))
78 : {
79 280804 : unsigned int n_sides = elem->n_sides();
80 1930316 : for (const auto s : make_range(n_sides))
81 : {
82 1649512 : boundary_info.boundary_ids(elem, s, old_ids);
83 1649512 : if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
84 : {
85 16581 : std::vector<boundary_id_type> new_ids(old_ids);
86 16581 : std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
87 16581 : if (delete_prev)
88 : {
89 15585 : boundary_info.remove_side(elem, s);
90 15585 : boundary_info.add_side(elem, s, new_ids);
91 : }
92 : else
93 996 : boundary_info.add_side(elem, s, new_ids);
94 16581 : }
95 : }
96 508 : }
97 :
98 : // Remove any remaining references to the old ID from the
99 : // BoundaryInfo object. This prevents things like empty sidesets
100 : // from showing up when printing information, etc.
101 508 : if (delete_prev)
102 399 : boundary_info.remove_id(old_id);
103 :
104 : // global information may now be out of sync
105 508 : mesh.unset_is_prepared();
106 508 : }
107 :
108 : std::vector<boundary_id_type>
109 9968 : getBoundaryIDs(const MeshBase & mesh,
110 : const std::vector<BoundaryName> & boundary_name,
111 : bool generate_unknown)
112 : {
113 : return getBoundaryIDs(
114 9968 : mesh, boundary_name, generate_unknown, mesh.get_boundary_info().get_boundary_ids());
115 : }
116 :
117 : std::vector<boundary_id_type>
118 125908 : getBoundaryIDs(const MeshBase & mesh,
119 : const std::vector<BoundaryName> & boundary_name,
120 : bool generate_unknown,
121 : const std::set<BoundaryID> & mesh_boundary_ids)
122 : {
123 125908 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
124 125908 : const std::map<BoundaryID, std::string> & sideset_map = boundary_info.get_sideset_name_map();
125 125908 : const std::map<BoundaryID, std::string> & nodeset_map = boundary_info.get_nodeset_name_map();
126 :
127 125908 : BoundaryID max_boundary_local_id = 0;
128 : /* It is required to generate a new ID for a given name. It is used often in mesh modifiers such
129 : * as SideSetsBetweenSubdomains. Then we need to check the current boundary ids since they are
130 : * changing during "mesh modify()", and figure out the right max boundary ID. Most of mesh
131 : * modifiers are running in serial, and we won't involve a global communication.
132 : */
133 125908 : if (generate_unknown)
134 : {
135 120838 : const auto & bids = mesh.is_prepared() ? mesh.get_boundary_info().get_global_boundary_ids()
136 1002 : : mesh.get_boundary_info().get_boundary_ids();
137 120838 : max_boundary_local_id = bids.empty() ? 0 : *(bids.rbegin());
138 : /* We should not hit this often */
139 120838 : if (!mesh.is_prepared() && !mesh.is_serial())
140 93 : mesh.comm().max(max_boundary_local_id);
141 : }
142 :
143 125908 : BoundaryID max_boundary_id = mesh_boundary_ids.empty() ? 0 : *(mesh_boundary_ids.rbegin());
144 :
145 125908 : max_boundary_id =
146 125908 : max_boundary_id > max_boundary_local_id ? max_boundary_id : max_boundary_local_id;
147 :
148 125908 : std::vector<BoundaryID> ids(boundary_name.size());
149 283389 : for (const auto i : index_range(boundary_name))
150 : {
151 157597 : if (boundary_name[i] == "ANY_BOUNDARY_ID")
152 : {
153 116 : ids.assign(mesh_boundary_ids.begin(), mesh_boundary_ids.end());
154 116 : if (i)
155 0 : mooseWarning("You passed \"ANY_BOUNDARY_ID\" in addition to other boundary_names. This "
156 : "may be a logic error.");
157 116 : break;
158 : }
159 :
160 157481 : if (boundary_name[i].empty() && !generate_unknown)
161 0 : mooseError("Incoming boundary name is empty and we are not generating unknown boundary IDs. "
162 : "This is invalid.");
163 :
164 : BoundaryID id;
165 :
166 157481 : if (boundary_name[i].empty() || !MooseUtils::isDigits(boundary_name[i]))
167 : {
168 : /**
169 : * If the conversion from a name to a number fails, that means that this must be a named
170 : * boundary. We will look in the complete map for this sideset and create a new name/ID pair
171 : * if requested.
172 : */
173 307584 : if (generate_unknown &&
174 221280 : !MooseUtils::doesMapContainValue(sideset_map, std::string(boundary_name[i])) &&
175 114822 : !MooseUtils::doesMapContainValue(nodeset_map, std::string(boundary_name[i])))
176 8146 : id = ++max_boundary_id;
177 : else
178 98312 : id = boundary_info.get_id_by_name(boundary_name[i]);
179 : }
180 : else
181 51023 : id = getIDFromName<BoundaryName, BoundaryID>(boundary_name[i]);
182 :
183 157481 : ids[i] = id;
184 : }
185 :
186 251816 : return ids;
187 0 : }
188 :
189 : std::set<BoundaryID>
190 169 : getBoundaryIDSet(const MeshBase & mesh,
191 : const std::vector<BoundaryName> & boundary_name,
192 : bool generate_unknown)
193 : {
194 169 : auto boundaries = getBoundaryIDs(mesh, boundary_name, generate_unknown);
195 338 : return std::set<BoundaryID>(boundaries.begin(), boundaries.end());
196 169 : }
197 :
198 : std::vector<subdomain_id_type>
199 306239 : getSubdomainIDs(const MeshBase & mesh, const std::vector<SubdomainName> & subdomain_names)
200 : {
201 306239 : std::vector<subdomain_id_type> ids;
202 :
203 : // shortcut for "ANY_BLOCK_ID"
204 306239 : if (subdomain_names.size() == 1 && subdomain_names[0] == "ANY_BLOCK_ID")
205 : {
206 : // since get_mesh_subdomains() requires a prepared mesh, we need to check that here
207 : mooseAssert(mesh.is_prepared(),
208 : "getSubdomainIDs() should only be called on a prepared mesh if ANY_BLOCK_ID is "
209 : "used to query all block IDs");
210 60117 : ids.assign(mesh.get_mesh_subdomains().begin(), mesh.get_mesh_subdomains().end());
211 60117 : return ids;
212 : }
213 :
214 : // loop through subdomain names and get IDs (this preserves the order of subdomain_names)
215 246122 : ids.resize(subdomain_names.size());
216 365396 : for (auto i : index_range(subdomain_names))
217 : {
218 119274 : if (subdomain_names[i] == "ANY_BLOCK_ID")
219 0 : mooseError("getSubdomainIDs() accepts \"ANY_BLOCK_ID\" if and only if it is the only "
220 : "subdomain name being queried.");
221 119274 : ids[i] = MooseMeshUtils::getSubdomainID(subdomain_names[i], mesh);
222 : }
223 :
224 246122 : return ids;
225 0 : }
226 :
227 : std::set<subdomain_id_type>
228 0 : getSubdomainIDs(const MeshBase & mesh, const std::set<SubdomainName> & subdomain_names)
229 : {
230 : const auto blk_ids = getSubdomainIDs(
231 0 : mesh, std::vector<SubdomainName>(subdomain_names.begin(), subdomain_names.end()));
232 0 : return {blk_ids.begin(), blk_ids.end()};
233 0 : }
234 :
235 : BoundaryID
236 373782 : getBoundaryID(const BoundaryName & boundary_name, const MeshBase & mesh)
237 : {
238 373782 : BoundaryID id = Moose::INVALID_BOUNDARY_ID;
239 373782 : if (boundary_name.empty())
240 0 : return id;
241 :
242 373782 : if (!MooseUtils::isDigits(boundary_name))
243 195349 : id = mesh.get_boundary_info().get_id_by_name(boundary_name);
244 : else
245 178433 : id = getIDFromName<BoundaryName, BoundaryID>(boundary_name);
246 :
247 373779 : return id;
248 : }
249 :
250 : SubdomainID
251 482715 : getSubdomainID(const SubdomainName & subdomain_name, const MeshBase & mesh)
252 : {
253 482715 : if (subdomain_name == "ANY_BLOCK_ID")
254 0 : mooseError("getSubdomainID() does not work with \"ANY_BLOCK_ID\"");
255 :
256 482715 : SubdomainID id = Moose::INVALID_BLOCK_ID;
257 482715 : if (subdomain_name.empty())
258 0 : return id;
259 :
260 482715 : if (!MooseUtils::isDigits(subdomain_name))
261 208078 : id = mesh.get_id_by_name(subdomain_name);
262 : else
263 274637 : id = getIDFromName<SubdomainName, SubdomainID>(subdomain_name);
264 :
265 482715 : return id;
266 : }
267 :
268 : void
269 0 : changeSubdomainId(MeshBase & mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
270 : {
271 0 : for (const auto & elem : mesh.element_ptr_range())
272 0 : if (elem->subdomain_id() == old_id)
273 0 : elem->subdomain_id() = new_id;
274 :
275 : // global cached information may now be out of sync
276 0 : mesh.unset_is_prepared();
277 0 : }
278 :
279 : Point
280 922 : meshCentroidCalculator(const MeshBase & mesh)
281 : {
282 922 : Point centroid_pt = Point(0.0, 0.0, 0.0);
283 922 : Real vol_tmp = 0.0;
284 19798 : for (const auto & elem : mesh.active_local_element_ptr_range())
285 : {
286 18876 : Real elem_vol = elem->volume();
287 18876 : centroid_pt += (elem->true_centroid()) * elem_vol;
288 18876 : vol_tmp += elem_vol;
289 922 : }
290 922 : mesh.comm().sum(centroid_pt);
291 922 : mesh.comm().sum(vol_tmp);
292 922 : centroid_pt /= vol_tmp;
293 1844 : return centroid_pt;
294 : }
295 :
296 : Point
297 152 : boundaryCentroidCalculator(const BoundaryName & boundary, MeshBase & mesh)
298 : {
299 : // Need boundaries to be synchronized
300 152 : if (!mesh.preparation().has_boundary_id_sets)
301 152 : mesh.get_boundary_info().synchronize_global_id_set();
302 152 : BoundaryInfo & mesh_boundary_info = mesh.get_boundary_info();
303 152 : boundary_id_type boundary_id = mesh_boundary_info.get_id_by_name(boundary);
304 152 : const auto side_list = mesh_boundary_info.build_side_list();
305 :
306 : // Initialize sums
307 152 : Real volume_sum = 0;
308 152 : Point volume_weighted_centroid_sum(0, 0, 0);
309 :
310 1208 : for (const auto & [eid, side_i, bid] : side_list)
311 : {
312 1056 : if (bid != boundary_id)
313 672 : continue;
314 :
315 : // Get the side
316 384 : const auto elem = mesh.elem_ptr(eid);
317 384 : const auto side = elem->side_ptr(side_i);
318 :
319 384 : volume_sum += side->volume();
320 384 : volume_weighted_centroid_sum += side->volume() * side->true_centroid();
321 384 : }
322 : // Sum across processes
323 152 : mesh.comm().sum(volume_weighted_centroid_sum);
324 152 : mesh.comm().sum(volume_sum);
325 :
326 304 : return volume_weighted_centroid_sum / volume_sum;
327 152 : }
328 :
329 : RealVectorValue
330 20 : boundaryWeightedNormal(const BoundaryName & boundary, MeshBase & mesh)
331 : {
332 : // Need boundaries to be synchronized
333 20 : if (!mesh.preparation().has_boundary_id_sets)
334 0 : mesh.get_boundary_info().synchronize_global_id_set();
335 20 : BoundaryInfo & mesh_boundary_info = mesh.get_boundary_info();
336 20 : boundary_id_type boundary_id = mesh_boundary_info.get_id_by_name(boundary);
337 20 : const auto side_list = mesh_boundary_info.build_side_list();
338 :
339 : // Initialize sums
340 20 : Real volume_sum = 0;
341 20 : RealVectorValue volume_weighted_normal_sum(0, 0, 0);
342 :
343 180 : for (const auto & [eid, side_i, bid] : side_list)
344 : {
345 160 : if (bid != boundary_id)
346 120 : continue;
347 :
348 : // Get the side
349 40 : const auto elem = mesh.elem_ptr(eid);
350 40 : const auto side = elem->side_ptr(side_i);
351 :
352 40 : volume_sum += side->volume();
353 40 : volume_weighted_normal_sum += side->volume() * elem->side_vertex_average_normal(side_i);
354 40 : }
355 : // Sum across processes
356 20 : mesh.comm().sum(volume_weighted_normal_sum);
357 20 : mesh.comm().sum(volume_sum);
358 :
359 40 : return volume_weighted_normal_sum / volume_sum;
360 20 : }
361 :
362 : Real
363 40 : computeMaxDistanceToAxis(const MeshBase & mesh,
364 : const Point & origin,
365 : const RealVectorValue & direction)
366 : {
367 40 : Real distance = 0;
368 : mooseAssert(MooseUtils::absoluteFuzzyEqual(direction.norm_sq(), 1),
369 : "Direction should be normalized");
370 1384 : for (const auto & node : mesh.node_ptr_range())
371 672 : if (const auto dist_node = (*node - origin).cross(direction).norm(); dist_node > distance)
372 94 : distance = dist_node;
373 40 : mesh.comm().max(distance);
374 40 : return distance;
375 : }
376 :
377 : std::unordered_map<dof_id_type, dof_id_type>
378 476 : getExtraIDUniqueCombinationMap(const MeshBase & mesh,
379 : const std::set<SubdomainID> & block_ids,
380 : std::vector<ExtraElementIDName> extra_ids)
381 : {
382 : // check block restriction
383 476 : const bool block_restricted = !block_ids.empty();
384 : // get element id name of interest in recursive parsing algorithm
385 476 : ExtraElementIDName id_name = extra_ids.back();
386 476 : extra_ids.pop_back();
387 476 : const auto id_index = mesh.get_elem_integer_index(id_name);
388 :
389 : // create base parsed id set
390 476 : if (extra_ids.empty())
391 : {
392 : // get set of extra id values;
393 278 : std::vector<dof_id_type> ids;
394 : {
395 278 : std::set<dof_id_type> ids_set;
396 1129933 : for (const auto & elem : mesh.active_element_ptr_range())
397 : {
398 1129655 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
399 520 : continue;
400 1129135 : const auto id = elem->get_extra_integer(id_index);
401 1129135 : ids_set.insert(id);
402 278 : }
403 278 : mesh.comm().set_union(ids_set);
404 278 : ids.assign(ids_set.begin(), ids_set.end());
405 278 : }
406 :
407 : // determine new extra id values;
408 278 : std::unordered_map<dof_id_type, dof_id_type> parsed_ids;
409 1129933 : for (auto & elem : mesh.active_element_ptr_range())
410 : {
411 1129655 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
412 520 : continue;
413 2258270 : parsed_ids[elem->id()] = std::distance(
414 2258270 : ids.begin(), std::lower_bound(ids.begin(), ids.end(), elem->get_extra_integer(id_index)));
415 278 : }
416 278 : return parsed_ids;
417 278 : }
418 :
419 : // if extra_ids is not empty, recursively call getExtraIDUniqueCombinationMap
420 : const auto base_parsed_ids =
421 198 : MooseMeshUtils::getExtraIDUniqueCombinationMap(mesh, block_ids, extra_ids);
422 : // parsing extra ids based on ref_parsed_ids
423 198 : std::vector<std::pair<dof_id_type, dof_id_type>> unique_ids;
424 : {
425 198 : std::set<std::pair<dof_id_type, dof_id_type>> unique_ids_set;
426 870470 : for (const auto & elem : mesh.active_element_ptr_range())
427 : {
428 870272 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
429 480 : continue;
430 869792 : const dof_id_type id1 = libmesh_map_find(base_parsed_ids, elem->id());
431 869792 : const dof_id_type id2 = elem->get_extra_integer(id_index);
432 869792 : const std::pair<dof_id_type, dof_id_type> ids = std::make_pair(id1, id2);
433 869792 : unique_ids_set.insert(ids);
434 198 : }
435 198 : mesh.comm().set_union(unique_ids_set);
436 198 : unique_ids.assign(unique_ids_set.begin(), unique_ids_set.end());
437 198 : }
438 :
439 198 : std::unordered_map<dof_id_type, dof_id_type> parsed_ids;
440 :
441 870470 : for (const auto & elem : mesh.active_element_ptr_range())
442 : {
443 870272 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
444 480 : continue;
445 869792 : const dof_id_type id1 = libmesh_map_find(base_parsed_ids, elem->id());
446 869792 : const dof_id_type id2 = elem->get_extra_integer(id_index);
447 869792 : const dof_id_type new_id = std::distance(
448 : unique_ids.begin(),
449 869792 : std::lower_bound(unique_ids.begin(), unique_ids.end(), std::make_pair(id1, id2)));
450 869792 : parsed_ids[elem->id()] = new_id;
451 198 : }
452 :
453 198 : return parsed_ids;
454 476 : }
455 :
456 : bool
457 9372 : isCoPlanar(const std::vector<Point> & vec_pts, const Point plane_nvec, const Point fixed_pt)
458 : {
459 38888 : for (const auto & pt : vec_pts)
460 38543 : if (!MooseUtils::absoluteFuzzyEqual((pt - fixed_pt) * plane_nvec, 0.0))
461 9027 : return false;
462 345 : return true;
463 : }
464 :
465 : bool
466 9104 : isCoPlanar(const std::vector<Point> & vec_pts, const Point plane_nvec)
467 : {
468 9104 : return isCoPlanar(vec_pts, plane_nvec, vec_pts.front());
469 : }
470 :
471 : bool
472 9104 : isCoPlanar(const std::vector<Point> & vec_pts)
473 : {
474 : // Assuming that overlapped Points are allowed, the Points that are overlapped with vec_pts[0] are
475 : // removed before further calculation.
476 18208 : std::vector<Point> vec_pts_nonzero{vec_pts[0]};
477 81936 : for (const auto i : index_range(vec_pts))
478 72832 : if (!MooseUtils::absoluteFuzzyEqual((vec_pts[i] - vec_pts[0]).norm(), 0.0))
479 36332 : vec_pts_nonzero.push_back(vec_pts[i]);
480 : // 3 or fewer points are always coplanar
481 9104 : if (vec_pts_nonzero.size() <= 3)
482 0 : return true;
483 : else
484 : {
485 18228 : for (const auto i : make_range(vec_pts_nonzero.size() - 1))
486 : {
487 18228 : const Point tmp_pt = (vec_pts_nonzero[i] - vec_pts_nonzero[0])
488 18228 : .cross(vec_pts_nonzero[i + 1] - vec_pts_nonzero[0]);
489 : // if the three points are not collinear, use cross product as the normal vector of the plane
490 18228 : if (!MooseUtils::absoluteFuzzyEqual(tmp_pt.norm(), 0.0))
491 9104 : return isCoPlanar(vec_pts_nonzero, tmp_pt.unit());
492 : }
493 : }
494 : // If all the points are collinear, they are also coplanar
495 0 : return true;
496 9104 : }
497 :
498 : SubdomainID
499 2187 : getNextFreeSubdomainID(MeshBase & input_mesh)
500 : {
501 : // Call this to get most up to date block id information
502 2187 : input_mesh.cache_elem_data();
503 :
504 2187 : std::set<SubdomainID> preexisting_subdomain_ids;
505 2187 : input_mesh.subdomain_ids(preexisting_subdomain_ids);
506 2187 : if (preexisting_subdomain_ids.empty())
507 0 : return 0;
508 : else
509 : {
510 : const auto highest_subdomain_id =
511 2187 : *std::max_element(preexisting_subdomain_ids.begin(), preexisting_subdomain_ids.end());
512 : mooseAssert(highest_subdomain_id < std::numeric_limits<SubdomainID>::max(),
513 : "A SubdomainID with max possible value was found");
514 2187 : return highest_subdomain_id + 1;
515 : }
516 2187 : }
517 :
518 : BoundaryID
519 1342 : getNextFreeBoundaryID(MeshBase & input_mesh)
520 : {
521 1342 : if (!input_mesh.preparation().has_boundary_id_sets)
522 345 : input_mesh.get_boundary_info().regenerate_id_sets();
523 :
524 1342 : auto boundary_ids = input_mesh.get_boundary_info().get_boundary_ids();
525 1342 : if (boundary_ids.empty())
526 406 : return 0;
527 936 : return (*boundary_ids.rbegin() + 1);
528 1342 : }
529 :
530 : bool
531 16101 : hasSubdomainID(const MeshBase & input_mesh, const SubdomainID & id)
532 : {
533 16101 : std::set<SubdomainID> mesh_blocks;
534 16101 : input_mesh.subdomain_ids(mesh_blocks);
535 :
536 : // On a distributed mesh we may have sideset IDs that only exist on
537 : // other processors
538 16101 : if (!input_mesh.is_replicated())
539 2164 : input_mesh.comm().set_union(mesh_blocks);
540 :
541 32202 : return mesh_blocks.count(id) && (id != Moose::INVALID_BLOCK_ID);
542 16101 : }
543 :
544 : bool
545 13683 : hasSubdomainName(const MeshBase & input_mesh, const SubdomainName & name)
546 : {
547 13683 : const auto id = getSubdomainID(name, input_mesh);
548 27366 : return hasSubdomainID(input_mesh, id);
549 : }
550 :
551 : bool
552 7776 : hasBoundaryID(const MeshBase & input_mesh, const BoundaryID id)
553 : {
554 7776 : const BoundaryInfo & boundary_info = input_mesh.get_boundary_info();
555 7776 : std::set<boundary_id_type> boundary_ids = boundary_info.get_boundary_ids();
556 :
557 : // On a distributed mesh we may have boundary IDs that only exist on
558 : // other processors
559 7776 : if (!input_mesh.is_replicated())
560 1054 : input_mesh.comm().set_union(boundary_ids);
561 :
562 15552 : return boundary_ids.count(id) && (id != Moose::INVALID_BOUNDARY_ID);
563 7776 : }
564 :
565 : bool
566 0 : hasBoundaryName(const MeshBase & mesh, const BoundaryName & name)
567 : {
568 0 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
569 0 : const BoundaryID id = boundary_info.get_id_by_name(name);
570 0 : return id != Moose::INVALID_BOUNDARY_ID;
571 : }
572 :
573 : bool
574 7633 : hasBoundaryNameOrID(const MeshBase & mesh, const BoundaryName & name_or_id)
575 : {
576 7633 : const auto id = getBoundaryID(name_or_id, mesh);
577 7633 : return hasBoundaryID(mesh, id);
578 : }
579 :
580 : void
581 436 : makeOrderedNodeList(std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
582 : std::vector<dof_id_type> & elem_id_list,
583 : std::vector<dof_id_type> & midpoint_node_list,
584 : std::vector<dof_id_type> & ordered_node_list,
585 : std::vector<dof_id_type> & ordered_elem_id_list)
586 : {
587 : // a flag to indicate if the ordered_node_list has been reversed
588 436 : bool is_flipped = false;
589 : // Start from the first element, try to find a chain of nodes
590 : mooseAssert(node_assm.size(), "Node list must not be empty");
591 436 : ordered_node_list.push_back(node_assm.front().first);
592 436 : if (midpoint_node_list.front() != DofObject::invalid_id)
593 0 : ordered_node_list.push_back(midpoint_node_list.front());
594 436 : ordered_node_list.push_back(node_assm.front().second);
595 436 : ordered_elem_id_list.push_back(elem_id_list.front());
596 : // Remove the element that has just been added to ordered_node_list
597 436 : node_assm.erase(node_assm.begin());
598 436 : midpoint_node_list.erase(midpoint_node_list.begin());
599 436 : elem_id_list.erase(elem_id_list.begin());
600 436 : const unsigned int node_assm_size_0 = node_assm.size();
601 2521 : for (unsigned int i = 0; i < node_assm_size_0; i++)
602 : {
603 : // Find nodes to expand the chain
604 2088 : dof_id_type end_node_id = ordered_node_list.back();
605 6994 : auto isMatch1 = [end_node_id](std::pair<dof_id_type, dof_id_type> old_id_pair)
606 6994 : { return old_id_pair.first == end_node_id; };
607 1894 : auto isMatch2 = [end_node_id](std::pair<dof_id_type, dof_id_type> old_id_pair)
608 1894 : { return old_id_pair.second == end_node_id; };
609 2088 : auto result = std::find_if(node_assm.begin(), node_assm.end(), isMatch1);
610 : bool match_first;
611 2088 : if (result == node_assm.end())
612 : {
613 1070 : match_first = false;
614 1070 : result = std::find_if(node_assm.begin(), node_assm.end(), isMatch2);
615 : }
616 : else
617 : {
618 1018 : match_first = true;
619 : }
620 : // If found, add the node to boundary_ordered_node_list
621 2088 : if (result != node_assm.end())
622 : {
623 1861 : const auto elem_index = std::distance(node_assm.begin(), result);
624 1861 : if (midpoint_node_list[elem_index] != DofObject::invalid_id)
625 0 : ordered_node_list.push_back(midpoint_node_list[elem_index]);
626 1861 : ordered_node_list.push_back(match_first ? (*result).second : (*result).first);
627 1861 : node_assm.erase(result);
628 1861 : midpoint_node_list.erase(midpoint_node_list.begin() + elem_index);
629 1861 : ordered_elem_id_list.push_back(elem_id_list[elem_index]);
630 1861 : elem_id_list.erase(elem_id_list.begin() + elem_index);
631 : }
632 : // If there are still elements in node_assm and result ==
633 : // node_assm.end(), this means the curve is not a loop, the
634 : // ordered_node_list is flipped and try the other direction that has not
635 : // been examined yet.
636 : else
637 : {
638 227 : if (is_flipped)
639 : // Flipped twice; this means the node list has at least two segments.
640 3 : throw MooseException("The node list provided has more than one segments.");
641 :
642 : // mark the first flip event.
643 224 : is_flipped = true;
644 224 : std::reverse(ordered_node_list.begin(), ordered_node_list.end());
645 224 : std::reverse(midpoint_node_list.begin(), midpoint_node_list.end());
646 224 : std::reverse(ordered_elem_id_list.begin(), ordered_elem_id_list.end());
647 : // As this iteration is wasted, set the iterator backward
648 224 : i--;
649 : }
650 : }
651 433 : }
652 :
653 : void
654 208 : makeOrderedNodeList(std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
655 : std::vector<dof_id_type> & elem_id_list,
656 : std::vector<dof_id_type> & ordered_node_list,
657 : std::vector<dof_id_type> & ordered_elem_id_list)
658 : {
659 208 : std::vector<dof_id_type> dummy_midpoint_node_list(node_assm.size(), DofObject::invalid_id);
660 208 : makeOrderedNodeList(
661 : node_assm, elem_id_list, dummy_midpoint_node_list, ordered_node_list, ordered_elem_id_list);
662 208 : }
663 :
664 : void
665 8160 : swapNodesInElem(Elem & elem, const unsigned int nd1, const unsigned int nd2)
666 : {
667 8160 : Node * n_temp = elem.node_ptr(nd1);
668 8160 : elem.set_node(nd1, elem.node_ptr(nd2));
669 8160 : elem.set_node(nd2, n_temp);
670 8160 : }
671 :
672 : void
673 394 : extraElemIntegerSwapParametersProcessor(
674 : const std::string & class_name,
675 : const unsigned int num_sections,
676 : const unsigned int num_integers,
677 : const std::vector<std::vector<std::vector<dof_id_type>>> & elem_integers_swaps,
678 : std::vector<std::unordered_map<dof_id_type, dof_id_type>> & elem_integers_swap_pairs)
679 : {
680 394 : elem_integers_swap_pairs.reserve(num_sections * num_integers);
681 418 : for (const auto i : make_range(num_integers))
682 : {
683 24 : const auto & elem_integer_swaps = elem_integers_swaps[i];
684 24 : std::vector<std::unordered_map<dof_id_type, dof_id_type>> elem_integer_swap_pairs;
685 : try
686 : {
687 48 : MooseMeshUtils::idSwapParametersProcessor(class_name,
688 : "elem_integers_swaps",
689 : elem_integer_swaps,
690 : elem_integer_swap_pairs,
691 : i * num_sections);
692 : }
693 0 : catch (const MooseException & e)
694 : {
695 0 : throw MooseException(e.what());
696 0 : }
697 :
698 24 : elem_integers_swap_pairs.insert(elem_integers_swap_pairs.end(),
699 : elem_integer_swap_pairs.begin(),
700 : elem_integer_swap_pairs.end());
701 24 : }
702 394 : }
703 :
704 : std::unique_ptr<ReplicatedMesh>
705 0 : buildBoundaryMesh(const MeshBase & input_mesh, const boundary_id_type boundary_id)
706 : {
707 0 : if (!input_mesh.is_serial())
708 0 : ::mooseError("Input mesh should be serialized for extracting the boundary mesh.\nInput mesh:" +
709 0 : input_mesh.get_info());
710 0 : auto poly_mesh = std::make_unique<ReplicatedMesh>(input_mesh.comm());
711 :
712 0 : auto side_list = input_mesh.get_boundary_info().build_side_list();
713 :
714 0 : std::unordered_map<dof_id_type, dof_id_type> old_new_node_map;
715 0 : for (const auto & [eid, side_i, bid] : side_list)
716 : {
717 0 : if (bid != boundary_id)
718 0 : continue;
719 :
720 : // Get the side
721 0 : const auto elem = input_mesh.elem_ptr(eid);
722 0 : const auto side = elem->side_ptr(side_i);
723 0 : auto side_elem = elem->build_side_ptr(side_i);
724 0 : auto copy = side_elem->build(side_elem->type());
725 :
726 0 : for (const auto i : side_elem->node_index_range())
727 : {
728 0 : auto & n = side_elem->node_ref(i);
729 :
730 0 : if (old_new_node_map.count(n.id()))
731 0 : copy->set_node(i, poly_mesh->node_ptr(old_new_node_map[n.id()]));
732 : else
733 : {
734 0 : Node * node = poly_mesh->add_point(side_elem->point(i));
735 0 : copy->set_node(i, node);
736 0 : old_new_node_map[n.id()] = node->id();
737 : }
738 : }
739 0 : poly_mesh->add_elem(copy.release());
740 0 : }
741 0 : poly_mesh->skip_partitioning(true);
742 0 : poly_mesh->prepare_for_use();
743 0 : if (poly_mesh->n_elem() == 0)
744 0 : mooseError("The input mesh to extract the boundary from does not have a boundary with id ",
745 : boundary_id,
746 : ".\n",
747 : input_mesh);
748 :
749 0 : return poly_mesh;
750 0 : }
751 :
752 : std::unique_ptr<ReplicatedMesh>
753 88 : buildLoopBoundaryOf2DMesh(const MeshBase & input_mesh, const boundary_id_type boundary_id)
754 : {
755 88 : if (!input_mesh.is_serial())
756 0 : ::mooseError(
757 0 : "Input 2D mesh should be serialized for extracting the loop boundary mesh.\nInput mesh:" +
758 0 : input_mesh.get_info());
759 88 : auto edge_mesh = std::make_unique<ReplicatedMesh>(input_mesh.comm());
760 88 : auto side_list = input_mesh.get_boundary_info().build_side_list();
761 88 : std::set<BoundaryInfo::BCTuple> visited;
762 88 : bool already_seen_this_side_tuple = false;
763 88 : BoundaryInfo::BCTuple first_side_visited = {libMesh::invalid_uint, 0, 0};
764 :
765 : // Helps move elem to elem at a given node
766 88 : const auto node_to_elem_map = buildBoundaryNodeToElemMap(input_mesh, boundary_id);
767 : // Helps check if a node is part of a boundary
768 88 : const auto & node_to_bids = input_mesh.get_boundary_info().get_nodeset_map();
769 :
770 : // Traverse from the first side (edge) in the side_list that matches the boundary_id
771 648 : for (const auto & bside : side_list)
772 : {
773 560 : if (std::get<2>(bside) != boundary_id)
774 472 : continue;
775 :
776 : // Check that we are not starting 'another' loop
777 560 : if (bside != first_side_visited)
778 : {
779 560 : if (visited.size() && !visited.count(bside))
780 0 : mooseWarning(
781 0 : "Boundary " + std::to_string(boundary_id) +
782 0 : " is not a (contiguous) loop. Boundary side: (" + Moose::stringify(bside) +
783 0 : ") was not visited after a single pass around the boundary. Boundary sides visited: " +
784 0 : Moose::stringify(visited));
785 560 : else if (visited.empty())
786 88 : first_side_visited = bside;
787 : else
788 472 : continue;
789 : }
790 :
791 : // Form the element to be able to find the side
792 : // These three variables will be updated while traversing the loop boundary
793 88 : const Elem * elem = input_mesh.elem_ptr(std::get<0>(bside));
794 88 : auto current_side = std::get<1>(bside);
795 88 : auto side_elem = elem->build_side_ptr(current_side);
796 :
797 : // 3D elements should not be part of this boundary
798 88 : if (elem->dim() != 2)
799 0 : mooseError(
800 : "Finding the loop boundary of a 2D mesh cannot be done with non-2D elements such as ",
801 : *elem);
802 :
803 : // Start from node 0 of the side (on the boundary), set the next node as the other node
804 : // one that side, and keep going from tht next node
805 88 : bool looped_back = false;
806 88 : const Node * starting_node = side_elem->node_ptr(0);
807 88 : const auto new_mesh_starting_node = edge_mesh->add_point(side_elem->point(0));
808 88 : Node * new_first_node = new_mesh_starting_node;
809 88 : [[maybe_unused]] dof_id_type first_node_index = starting_node->id();
810 88 : dof_id_type second_node_index = input_mesh.node_ptr(side_elem->node_id(1))->id();
811 :
812 648 : while (!looped_back && !already_seen_this_side_tuple)
813 : {
814 560 : if (MooseUtils::absoluteFuzzyEqual(input_mesh.point(second_node_index),
815 1120 : Point(*starting_node)))
816 88 : looped_back = true;
817 :
818 : // Get the opposite node (the next node) and add it to the edge mesh
819 : Node * new_second_node = looped_back
820 560 : ? new_mesh_starting_node
821 472 : : edge_mesh->add_point(input_mesh.point(second_node_index));
822 :
823 : // Add a copy of the edge side element to the mesh
824 560 : side_elem = elem->build_side_ptr(current_side);
825 560 : auto copy = side_elem->build(side_elem->type());
826 560 : copy->set_node(0, new_first_node);
827 560 : copy->set_node(1, new_second_node);
828 560 : edge_mesh->add_elem(copy.release());
829 :
830 : // Make this side as 'visited'
831 : std::tuple<dof_id_type, unsigned short int, boundary_id_type> bc_tuple = {
832 560 : elem->id(), current_side, boundary_id};
833 560 : const auto & visit_iter = visited.insert(bc_tuple);
834 560 : if (!looped_back && !visit_iter.second)
835 0 : already_seen_this_side_tuple = true;
836 :
837 : // Find the next element and side_elem
838 560 : auto & connected_elems = libmesh_map_find(node_to_elem_map, second_node_index);
839 560 : bool found_match = false;
840 560 : const auto current_eid = elem->id();
841 :
842 1032 : for (const auto eid : connected_elems)
843 : {
844 : mooseAssert(!found_match,
845 : "We should only find one node on a connected element on this boundary");
846 856 : if (eid != current_eid)
847 : {
848 : // Update the element (on the input mesh)
849 496 : elem = input_mesh.elem_ptr(eid);
850 :
851 : // Find the side and the opposite node index in that side
852 1032 : for (const auto si : elem->side_index_range())
853 : {
854 : // Check that second node is on the side
855 : const auto local_second_node_index =
856 920 : elem->get_node_index(input_mesh.node_ptr(second_node_index));
857 : // 2 sides should match this
858 920 : if (elem->is_node_on_side(local_second_node_index, si))
859 : {
860 : // Only one side should be on the same boundary (node is connected to two elements)
861 : // Form a bc_tuple and check the list of boundary sides
862 : std::tuple<dof_id_type, unsigned short int, boundary_id_type> side_bc_tuple = {
863 760 : elem->id(), si, boundary_id};
864 :
865 760 : if (std::find(side_list.begin(), side_list.end(), side_bc_tuple) == side_list.end())
866 376 : continue;
867 :
868 : // We are on the right boundary, just need to get the other node
869 768 : for (const auto local_side_node_id : elem->nodes_on_side(si))
870 : {
871 768 : const auto side_node_id = elem->node_id(local_side_node_id);
872 :
873 : // Skip current node (use global index to compare)
874 768 : if (side_node_id == second_node_index)
875 384 : continue;
876 : mooseAssert(side_node_id != first_node_index,
877 : "Somehow looped back in a single element");
878 :
879 384 : current_side = si;
880 384 : second_node_index = side_node_id;
881 384 : found_match = true;
882 384 : break;
883 384 : }
884 : }
885 : // No need to examine more sides
886 544 : if (found_match)
887 384 : break;
888 : }
889 :
890 : // No need to examine more elements
891 496 : if (found_match)
892 384 : break;
893 : }
894 : // next node could be on the same element, just moving on to the next side
895 360 : else if (connected_elems.size() == 1)
896 : {
897 176 : elem = input_mesh.elem_ptr(eid);
898 : const auto local_second_node_index =
899 176 : elem->get_node_index(input_mesh.node_ptr(second_node_index));
900 :
901 : // Move on to the next side
902 416 : for (const auto si : elem->side_index_range())
903 416 : if (si != current_side && elem->is_node_on_side(local_second_node_index, si))
904 : {
905 : // Check all nodes on that next side
906 528 : for (const auto local_side_node_id : elem->nodes_on_side(si))
907 : {
908 352 : const auto side_node_id = elem->node_id(local_side_node_id);
909 : // Skip current node
910 352 : if (side_node_id == second_node_index)
911 176 : continue;
912 : mooseAssert((side_node_id != first_node_index) ||
913 : (side_list.size() == elem->n_sides()),
914 : "Somehow looped back in a single element");
915 :
916 : // Check all the boundaries the other node (on the edge side) is part of
917 176 : const auto bids_range = node_to_bids.equal_range(input_mesh.node_ptr(side_node_id));
918 :
919 352 : for (auto iter = bids_range.first; iter != bids_range.second; iter++)
920 176 : if (iter->second == boundary_id)
921 : {
922 176 : current_side = si;
923 176 : second_node_index = side_node_id;
924 176 : found_match = true;
925 : }
926 176 : }
927 :
928 : // no need to examine other sides
929 176 : if (found_match)
930 176 : break;
931 : }
932 : }
933 : }
934 :
935 : // Set current node to opposite node of new element
936 : // NOTE: do not use new_first_node or new_second_node to search in the input mesh!
937 560 : new_first_node = new_second_node;
938 560 : first_node_index = second_node_index;
939 :
940 : // Handle loop ending criterion
941 560 : if (!found_match)
942 : {
943 0 : mooseWarning("Search for next element in loop boundary failed. Is boundary '" +
944 0 : std::to_string(boundary_id) + "' of mesh ",
945 : input_mesh,
946 : " a loop boundary?");
947 0 : break;
948 : }
949 560 : }
950 88 : }
951 :
952 88 : if (already_seen_this_side_tuple)
953 0 : mooseWarning("Boundary " + std::to_string(boundary_id) +
954 : " seems to have cycles. A single-cycle loop should be used");
955 :
956 88 : edge_mesh->skip_partitioning(true);
957 88 : edge_mesh->prepare_for_use();
958 88 : if (edge_mesh->n_elem() == 0)
959 0 : mooseError("The input mesh to extract the boundary from does not have a boundary with id ",
960 : boundary_id,
961 : "\n",
962 : input_mesh);
963 :
964 176 : return edge_mesh;
965 88 : }
966 :
967 : std::unordered_map<dof_id_type, std::unordered_set<dof_id_type>>
968 88 : buildBoundaryNodeToElemMap(const MeshBase & input_mesh, const boundary_id_type boundary_id)
969 : {
970 88 : if (!input_mesh.is_serial())
971 0 : ::mooseError(
972 0 : "Input 2D mesh should be serialized for extracting the loop boundary mesh.\nInput mesh:" +
973 0 : input_mesh.get_info());
974 :
975 : // Get all nodes on that boundary
976 : // Boundary ID might be a sideset or a nodeset, get nodes regardless
977 88 : const auto particular_node_ids = getBoundaryNodes(input_mesh, boundary_id);
978 :
979 88 : std::unordered_map<dof_id_type, std::unordered_set<dof_id_type>> nid_to_eids_map;
980 : // Fill the map from looping over elements
981 88 : for (const auto & elem :
982 944 : as_range(input_mesh.active_elements_begin(), input_mesh.active_elements_end()))
983 : {
984 1536 : for (const auto & nd : elem->node_ref_range())
985 : {
986 : // Only add the element id if the node is on the boundary
987 1152 : if (!particular_node_ids.count(nd.id()))
988 0 : continue;
989 :
990 1152 : auto & elem_ids = nid_to_eids_map[nd.id()];
991 1152 : elem_ids.insert(elem->id());
992 : }
993 88 : }
994 176 : return nid_to_eids_map;
995 88 : }
996 :
997 : std::set<dof_id_type>
998 88 : getBoundaryNodes(const MeshBase & mesh, const BoundaryID boundary_id)
999 : {
1000 88 : std::set<dof_id_type> boundary_node_ids;
1001 88 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1002 :
1003 : // Get all nodes from the sideset with ID of boundary_id
1004 : const auto & bc_sides =
1005 88 : boundary_info.build_side_list(libMesh::BoundaryInfo::BCTupleSortBy::BOUNDARY_ID);
1006 648 : for (const auto & [elem_id, side, bc_id] : bc_sides)
1007 : {
1008 560 : if (bc_id == boundary_id)
1009 : {
1010 560 : const auto elem = mesh.elem_ptr(elem_id);
1011 1680 : for (const auto ni : elem->nodes_on_side(side))
1012 1680 : boundary_node_ids.insert(elem->node_id(ni));
1013 : }
1014 : }
1015 :
1016 : // Get all nodes from nodeset with ID of boundary_id
1017 88 : const auto & bc_nodes = boundary_info.build_node_list();
1018 648 : for (const auto & [n_id, bc_id] : bc_nodes)
1019 560 : if (bc_id == boundary_id)
1020 560 : boundary_node_ids.insert(n_id);
1021 :
1022 176 : return boundary_node_ids;
1023 88 : }
1024 :
1025 : void
1026 2553 : createSubdomainFromSidesets(MeshBase & mesh,
1027 : std::vector<BoundaryName> boundary_names,
1028 : const SubdomainID new_subdomain_id,
1029 : const SubdomainName new_subdomain_name,
1030 : const std::string type_name)
1031 : {
1032 : // Generate a new block id if one isn't supplied.
1033 2553 : SubdomainID new_block_id = new_subdomain_id;
1034 :
1035 : // Make sure our boundary info and parallel counts are setup
1036 2553 : if (!mesh.is_prepared())
1037 : {
1038 475 : const bool allow_remote_element_removal = mesh.allow_remote_element_removal();
1039 : // We want all of our boundary elements available, so avoid removing them if they haven't
1040 : // already been so
1041 475 : mesh.allow_remote_element_removal(false);
1042 475 : mesh.prepare_for_use();
1043 475 : mesh.allow_remote_element_removal(allow_remote_element_removal);
1044 : }
1045 :
1046 : // Check that the sidesets are present in the mesh
1047 5258 : for (const auto & sideset : boundary_names)
1048 2714 : if (!MooseMeshUtils::hasBoundaryNameOrID(mesh, sideset))
1049 9 : mooseException("The sideset '", sideset, "' was not found within the mesh");
1050 :
1051 2544 : auto sideset_ids = MooseMeshUtils::getBoundaryIDs(mesh, boundary_names, true);
1052 2544 : std::set<boundary_id_type> sidesets(sideset_ids.begin(), sideset_ids.end());
1053 2544 : auto side_list = mesh.get_boundary_info().build_side_list();
1054 2544 : if (!mesh.is_serial() && mesh.comm().size() > 1)
1055 : {
1056 10 : std::vector<Elem *> elements_to_send;
1057 10 : unsigned short i_need_boundary_elems = 0;
1058 172 : for (const auto & [elem_id, side, bc_id] : side_list)
1059 : {
1060 162 : libmesh_ignore(side);
1061 162 : if (sidesets.count(bc_id))
1062 : {
1063 : // Whether we have this boundary information through our locally owned element or a ghosted
1064 : // element, we'll need the boundary elements for parallel consistent addition
1065 92 : i_need_boundary_elems = 1;
1066 92 : auto * elem = mesh.elem_ptr(elem_id);
1067 92 : if (elem->processor_id() == mesh.processor_id())
1068 80 : elements_to_send.push_back(elem);
1069 : }
1070 : }
1071 :
1072 : std::set<const Elem *, libMesh::CompareElemIdsByLevel> connected_elements(
1073 10 : elements_to_send.begin(), elements_to_send.end());
1074 10 : std::set<const Node *> connected_nodes;
1075 10 : reconnect_nodes(connected_elements, connected_nodes);
1076 10 : std::set<dof_id_type> connected_node_ids;
1077 178 : for (auto * nd : connected_nodes)
1078 168 : connected_node_ids.insert(nd->id());
1079 :
1080 10 : std::vector<unsigned short> need_boundary_elems(mesh.comm().size());
1081 10 : mesh.comm().allgather(i_need_boundary_elems, need_boundary_elems);
1082 10 : std::unordered_map<processor_id_type, decltype(elements_to_send)> push_element_data;
1083 10 : std::unordered_map<processor_id_type, decltype(connected_nodes)> push_node_data;
1084 :
1085 30 : for (const auto pid : index_range(mesh.comm()))
1086 : // Don't need to send to self
1087 20 : if (pid != mesh.processor_id() && need_boundary_elems[pid])
1088 : {
1089 10 : if (elements_to_send.size())
1090 10 : push_element_data[pid] = elements_to_send;
1091 10 : if (connected_nodes.size())
1092 10 : push_node_data[pid] = connected_nodes;
1093 : }
1094 :
1095 10 : auto node_action_functor = [](processor_id_type, const auto &)
1096 : {
1097 : // Node packing specialization already has unpacked node into mesh, so nothing to do
1098 10 : };
1099 10 : Parallel::push_parallel_packed_range(mesh.comm(), push_node_data, &mesh, node_action_functor);
1100 10 : auto elem_action_functor = [](processor_id_type, const auto &)
1101 : {
1102 : // Elem packing specialization already has unpacked elem into mesh, so nothing to do
1103 10 : };
1104 10 : TIMPI::push_parallel_packed_range(mesh.comm(), push_element_data, &mesh, elem_action_functor);
1105 :
1106 : // now that we've gathered everything, we need to rebuild the side list
1107 10 : side_list = mesh.get_boundary_info().build_side_list();
1108 10 : }
1109 :
1110 2544 : std::vector<std::pair<dof_id_type, ElemSidePair>> element_sides_on_boundary;
1111 2544 : dof_id_type counter = 0;
1112 150766 : for (const auto & [eid, side, bid] : side_list)
1113 148225 : if (sidesets.count(bid))
1114 : {
1115 26609 : if (auto elem = mesh.query_elem_ptr(eid))
1116 : {
1117 26609 : if (!elem->active())
1118 3 : mooseError(
1119 : "Only active, level 0 elements can be made interior parents of new level 0 lower-d "
1120 : "elements. Make sure that ",
1121 : type_name,
1122 : "s are run before any refinement generators");
1123 26606 : element_sides_on_boundary.push_back(std::make_pair(counter, ElemSidePair(elem, side)));
1124 : }
1125 26606 : ++counter;
1126 : }
1127 :
1128 2541 : dof_id_type max_elem_id = mesh.max_elem_id();
1129 2541 : unique_id_type max_unique_id = mesh.parallel_max_unique_id();
1130 :
1131 : // Making an important assumption that at least our boundary elements are the same on all
1132 : // processes even in distributed mesh mode (this is reliant on the correct ghosting functors
1133 : // existing on the mesh)
1134 29147 : for (auto & [i, elem_side] : element_sides_on_boundary)
1135 : {
1136 26606 : Elem * elem = elem_side.elem;
1137 :
1138 26606 : const auto side = elem_side.side;
1139 :
1140 : // Build a non-proxy element from this side.
1141 26606 : std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
1142 :
1143 : // The side will be added with the same processor id as the parent.
1144 26606 : side_elem->processor_id() = elem->processor_id();
1145 :
1146 : // Add subdomain ID
1147 26606 : side_elem->subdomain_id() = new_block_id;
1148 :
1149 : // Also assign the side's interior parent, so it is always
1150 : // easy to figure out the Elem we came from.
1151 26606 : side_elem->set_interior_parent(elem);
1152 :
1153 : // Add id
1154 26606 : side_elem->set_id(max_elem_id + i);
1155 26606 : side_elem->set_unique_id(max_unique_id + i);
1156 :
1157 : // Finally, add the lower-dimensional element to the mesh.
1158 26606 : mesh.add_elem(side_elem.release());
1159 26606 : };
1160 :
1161 : // Assign block name, if provided
1162 2541 : if (new_subdomain_name.size())
1163 1407 : mesh.subdomain_name(new_block_id) = new_subdomain_name;
1164 :
1165 2541 : const bool skip_partitioning_old = mesh.skip_partitioning();
1166 2541 : mesh.skip_partitioning(true);
1167 2541 : mesh.prepare_for_use();
1168 2541 : mesh.skip_partitioning(skip_partitioning_old);
1169 2541 : }
1170 :
1171 : void
1172 929 : convertBlockToMesh(MeshBase & source_mesh,
1173 : MeshBase & target_mesh,
1174 : const std::vector<SubdomainName> & target_blocks)
1175 : {
1176 929 : if (!source_mesh.is_replicated())
1177 0 : mooseError("This generator does not support distributed meshes.");
1178 :
1179 929 : const auto target_block_ids = MooseMeshUtils::getSubdomainIDs(source_mesh, target_blocks);
1180 :
1181 : // Check that the block ids/names exist in the mesh
1182 929 : std::set<SubdomainID> mesh_blocks;
1183 929 : source_mesh.subdomain_ids(mesh_blocks);
1184 :
1185 1921 : for (const auto i : index_range(target_block_ids))
1186 992 : if (target_block_ids[i] == Moose::INVALID_BLOCK_ID || !mesh_blocks.count(target_block_ids[i]))
1187 : {
1188 0 : mooseException("The target_block '", target_blocks[i], "' was not found within the mesh.");
1189 : }
1190 :
1191 : // know which nodes have already been inserted, by tracking the old mesh's node's ids'
1192 929 : std::unordered_map<dof_id_type, dof_id_type> old_new_node_map;
1193 :
1194 1906 : for (const auto target_block_id : target_block_ids)
1195 : {
1196 :
1197 17232 : for (auto elem : source_mesh.active_subdomain_elements_ptr_range(target_block_id))
1198 : {
1199 8129 : if (elem->level() != 0)
1200 3 : mooseError("Refined blocks are not supported by this generator. "
1201 : "Can you re-organize mesh generators to refine after converting the block?");
1202 :
1203 : // make a deep copy so that mutiple meshes' destructors don't segfault at program termination
1204 8126 : auto copy = elem->build(elem->type());
1205 :
1206 : // Keep the subdomain id
1207 8126 : copy->subdomain_id() = elem->subdomain_id();
1208 :
1209 : // index of node in the copy element must be managed manually as there is no intelligent
1210 : // insert method
1211 8126 : dof_id_type copy_n_index = 0;
1212 :
1213 : // correctly assign new copies of nodes, loop over nodes
1214 35950 : for (dof_id_type i : elem->node_index_range())
1215 : {
1216 27824 : auto & n = elem->node_ref(i);
1217 :
1218 27824 : if (old_new_node_map.count(n.id()))
1219 : {
1220 : // case where we have already inserted this particular point before
1221 : // then we need to find the already-inserted one and hook it up right
1222 : // to it's respective element
1223 16164 : copy->set_node(copy_n_index++, target_mesh.node_ptr(old_new_node_map[n.id()]));
1224 : }
1225 : else
1226 : {
1227 : // case where we've NEVER inserted this particular point before
1228 : // add them both to the element and the mesh
1229 :
1230 : // Nodes' IDs are their indexes in the nodes' respective mesh
1231 : // If we set them as invalid they are automatically assigned
1232 : // Add to mesh, auto-assigning a new id.
1233 11660 : Node * node = target_mesh.add_point(elem->point(i));
1234 :
1235 : // Add to element copy (manually)
1236 11660 : copy->set_node(copy_n_index++, node);
1237 :
1238 : // remember the (old) ID
1239 11660 : old_new_node_map[n.id()] = node->id();
1240 : }
1241 : }
1242 :
1243 : // it is ok to release the copy element into the mesh because derived meshes class
1244 : // (ReplicatedMesh, DistributedMesh) manage their own elements, will delete them
1245 8126 : target_mesh.add_elem(copy.release());
1246 9103 : }
1247 : }
1248 :
1249 : // Move subdomain names
1250 1900 : for (const auto sbd_id : target_block_ids)
1251 974 : target_mesh.subdomain_name(sbd_id) = source_mesh.subdomain_name(sbd_id);
1252 926 : }
1253 :
1254 : void
1255 1967 : copyIntoMesh(MeshGenerator & mg,
1256 : UnstructuredMesh & destination,
1257 : const UnstructuredMesh & source,
1258 : const bool avoid_merging_subdomains,
1259 : const bool avoid_merging_boundaries,
1260 : const Parallel::Communicator & communicator)
1261 : {
1262 1967 : dof_id_type node_delta = destination.max_node_id();
1263 1967 : dof_id_type elem_delta = destination.max_elem_id();
1264 :
1265 : unique_id_type unique_delta =
1266 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1267 1967 : destination.parallel_max_unique_id();
1268 : #else
1269 : 0;
1270 : #endif
1271 :
1272 : // Prevent overlaps by offsetting the subdomains in
1273 1967 : std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
1274 1967 : unsigned int block_offset = 0;
1275 1967 : if (avoid_merging_subdomains)
1276 : {
1277 : // Note: if performance becomes an issue, this is overkill for just getting the max node id
1278 132 : std::set<subdomain_id_type> source_ids;
1279 132 : std::set<subdomain_id_type> dest_ids;
1280 :
1281 : // We need source subdomain ids already cached; libMesh will
1282 : // scream otherwise
1283 132 : source.subdomain_ids(source_ids, true);
1284 :
1285 : // Our destination is non-const, so we can fix any missing caches
1286 132 : if (!destination.preparation().has_cached_elem_data)
1287 132 : destination.cache_elem_data();
1288 :
1289 132 : destination.subdomain_ids(dest_ids, true);
1290 :
1291 : mooseAssert(source_ids.size(), "Should have a subdomain");
1292 : mooseAssert(dest_ids.size(), "Should have a subdomain");
1293 132 : unsigned int max_dest_bid = *dest_ids.rbegin();
1294 132 : unsigned int min_source_bid = *source_ids.begin();
1295 132 : communicator.max(max_dest_bid);
1296 132 : communicator.min(min_source_bid);
1297 132 : block_offset = 1 + max_dest_bid - min_source_bid;
1298 264 : for (const auto bid : source_ids)
1299 132 : id_remapping[bid] = block_offset + bid;
1300 132 : }
1301 :
1302 : // Copy mesh data over from the other mesh
1303 1967 : destination.copy_nodes_and_elements(source,
1304 : // Skipping this should cause the neighbors
1305 : // to simply be copied from the other mesh
1306 : // (which makes sense and is way faster)
1307 : /*skip_find_neighbors = */ true,
1308 : elem_delta,
1309 : node_delta,
1310 : unique_delta,
1311 : avoid_merging_subdomains ? &id_remapping : nullptr);
1312 :
1313 : // Get an offset to prevent overlaps / wild merging between boundaries
1314 1967 : BoundaryInfo & boundary = destination.get_boundary_info();
1315 1967 : const BoundaryInfo & other_boundary = source.get_boundary_info();
1316 :
1317 1967 : unsigned int bid_offset = 0;
1318 1967 : if (avoid_merging_boundaries)
1319 : {
1320 162 : const auto boundary_ids = boundary.get_boundary_ids();
1321 162 : const auto other_boundary_ids = other_boundary.get_boundary_ids();
1322 162 : unsigned int max_dest_bid = boundary_ids.size() ? *boundary_ids.rbegin() : 0;
1323 162 : unsigned int min_source_bid = other_boundary_ids.size() ? *other_boundary_ids.begin() : 0;
1324 162 : communicator.max(max_dest_bid);
1325 162 : communicator.min(min_source_bid);
1326 162 : bid_offset = 1 + max_dest_bid - min_source_bid;
1327 162 : }
1328 :
1329 : // Note: the code below originally came from ReplicatedMesh::stitch_mesh_helper()
1330 : // in libMesh replicated_mesh.C around line 1203
1331 :
1332 : // Copy BoundaryInfo from other_mesh too. We do this via the
1333 : // list APIs rather than element-by-element for speed.
1334 55013 : for (const auto & t : other_boundary.build_node_list())
1335 55013 : boundary.add_node(std::get<0>(t) + node_delta, bid_offset + std::get<1>(t));
1336 :
1337 37877 : for (const auto & t : other_boundary.build_side_list())
1338 37877 : boundary.add_side(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
1339 :
1340 1967 : for (const auto & t : other_boundary.build_edge_list())
1341 1967 : boundary.add_edge(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
1342 :
1343 1967 : for (const auto & t : other_boundary.build_shellface_list())
1344 0 : boundary.add_shellface(
1345 1967 : std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
1346 :
1347 : // Check for the case with two block ids sharing the same name
1348 1967 : if (avoid_merging_subdomains)
1349 : {
1350 : mooseAssert(mg.parameters().isParamDefined("avoid_merging_subdomains"),
1351 : "Missing parameter in the mesh generator calling this function: "
1352 : "avoid_merging_subdomains. Considering setting avoid_merging_subdomains to true.");
1353 304 : for (const auto & [block_id, block_name] : destination.get_subdomain_name_map())
1354 233 : for (const auto & [source_id, source_name] : source.get_subdomain_name_map())
1355 61 : if (block_name == source_name)
1356 0 : mg.paramWarning(
1357 : "avoid_merging_subdomains",
1358 0 : "Not merging subdomains is creating two subdomains with the same name '" +
1359 0 : block_name + "' but different ids: " + std::to_string(source_id) + " & " +
1360 0 : std::to_string(block_id + block_offset) +
1361 : ".\n We recommend using a RenameBlockGenerator to prevent this as you "
1362 : "will get errors reading the Exodus output later.");
1363 : }
1364 :
1365 2137 : for (const auto & [block_id, block_name] : source.get_subdomain_name_map())
1366 340 : destination.set_subdomain_name_map().insert(
1367 340 : std::make_pair<SubdomainID, SubdomainName>(block_id + block_offset, block_name));
1368 :
1369 : // Check for the case with two boundary ids sharing the same name
1370 1967 : if (avoid_merging_boundaries)
1371 : {
1372 : mooseAssert(mg.parameters().isParamDefined("avoid_merging_boundaries"),
1373 : "Missing parameter in the mesh generator calling this function: "
1374 : "avoid_merging_boundaries. Considering setting avoid_merging_boundaries to true.");
1375 830 : for (const auto & [b_id, b_name] : other_boundary.get_sideset_name_map())
1376 4740 : for (const auto & [source_id, source_name] : boundary.get_sideset_name_map())
1377 4072 : if (b_name == source_name)
1378 0 : mg.paramWarning(
1379 : "avoid_merging_boundaries",
1380 0 : "Not merging boundaries is creating two sidesets with the same name '" + b_name +
1381 0 : "' but different ids: " + std::to_string(source_id) + " & " +
1382 0 : std::to_string(b_id + bid_offset) +
1383 : ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
1384 : "will get errors reading the Exodus output later.");
1385 830 : for (const auto & [b_id, b_name] : other_boundary.get_nodeset_name_map())
1386 4740 : for (const auto & [source_id, source_name] : boundary.get_nodeset_name_map())
1387 4072 : if (b_name == source_name)
1388 0 : mg.paramWarning(
1389 : "avoid_merging_boundaries",
1390 0 : "Not merging boundaries is creating two nodesets with the same name '" + b_name +
1391 0 : "' but different ids: " + std::to_string(source_id) + " & " +
1392 0 : std::to_string(b_id + bid_offset) +
1393 : ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
1394 : "will get errors reading the Exodus output later.");
1395 : }
1396 :
1397 9043 : for (const auto & [nodeset_id, nodeset_name] : other_boundary.get_nodeset_name_map())
1398 14152 : boundary.set_nodeset_name_map().insert(
1399 14152 : std::make_pair<BoundaryID, BoundaryName>(nodeset_id + bid_offset, nodeset_name));
1400 :
1401 9303 : for (const auto & [sideset_id, sideset_name] : other_boundary.get_sideset_name_map())
1402 14672 : boundary.set_sideset_name_map().insert(
1403 14672 : std::make_pair<BoundaryID, BoundaryName>(sideset_id + bid_offset, sideset_name));
1404 :
1405 1967 : for (const auto & [edgeset_id, edgeset_name] : other_boundary.get_edgeset_name_map())
1406 0 : boundary.set_edgeset_name_map().insert(
1407 0 : std::make_pair<BoundaryID, BoundaryName>(edgeset_id + bid_offset, edgeset_name));
1408 1967 : }
1409 :
1410 : void
1411 1364 : buildPolyLineMesh(MeshBase & mesh,
1412 : const std::vector<Point> & points,
1413 : const std::vector<Point> & mid_points,
1414 : const bool loop,
1415 : const BoundaryName & start_boundary,
1416 : const BoundaryName & end_boundary,
1417 : const std::vector<unsigned int> & nums_edges_between_points)
1418 : {
1419 : mooseAssert(nums_edges_between_points.size() == 1 ||
1420 : nums_edges_between_points.size() == points.size() - 1 + loop,
1421 : "nums_edges_between_points must be either a single value or have the same number of "
1422 : "entries as segments defined by the points.");
1423 : mooseAssert(
1424 : mid_points.size() == 0 || mid_points.size() == points.size() - (loop ? 0 : 1),
1425 : "mid_points must be either empty or have the consistent number of entries as points.");
1426 : mooseAssert(
1427 : mid_points.size() == 0 ||
1428 : (nums_edges_between_points.size() == 1 && nums_edges_between_points.front() == 1) ||
1429 : (nums_edges_between_points.size() == points.size() - 1 + loop &&
1430 : std::all_of(nums_edges_between_points.begin(),
1431 : nums_edges_between_points.end(),
1432 : [](unsigned int n) { return n == 1; })),
1433 : "mid_points can only be provided if each segment has exactly one edge.");
1434 :
1435 1364 : const auto n_points = points.size();
1436 24256 : for (auto i : make_range(n_points))
1437 : {
1438 : const auto & num_edges_between_points =
1439 22902 : (nums_edges_between_points.size() == 1)
1440 23522 : ? nums_edges_between_points[0]
1441 620 : : (i == nums_edges_between_points.size() ? 0 : nums_edges_between_points[i]);
1442 :
1443 22902 : Point p = points[i];
1444 22902 : const auto pt_counter = (nums_edges_between_points.size() == 1)
1445 23522 : ? i
1446 620 : : std::accumulate(nums_edges_between_points.begin(),
1447 620 : nums_edges_between_points.begin() + i,
1448 22902 : 0);
1449 45184 : mesh.add_point(
1450 45184 : p, nums_edges_between_points.size() == 1 ? (i * num_edges_between_points) : pt_counter);
1451 :
1452 22902 : if (num_edges_between_points > 1)
1453 : {
1454 750 : if (!loop && (i + 1) == n_points)
1455 10 : break;
1456 :
1457 740 : const auto ip1 = (i + 1) % n_points;
1458 740 : const Point pvec = (points[ip1] - p) / num_edges_between_points;
1459 :
1460 1650 : for (auto j : make_range(1u, num_edges_between_points))
1461 : {
1462 910 : p += pvec;
1463 1820 : mesh.add_point(
1464 : p,
1465 910 : (nums_edges_between_points.size() == 1 ? (i * num_edges_between_points) : pt_counter) +
1466 910 : j);
1467 : }
1468 : }
1469 : }
1470 : // Add mid points if applicable. When mid points are provided, each segment has exactly one edge,
1471 : // so the midpoint node ids follow the vertex node ids.
1472 3324 : for (const auto & i : make_range(mid_points.size()))
1473 1960 : mesh.add_point(mid_points[i], n_points + i);
1474 :
1475 1364 : const auto n_segments = loop ? n_points : (n_points - 1);
1476 : const auto n_elem =
1477 1364 : nums_edges_between_points.size() == 1
1478 1472 : ? n_segments * nums_edges_between_points[0]
1479 108 : : std::accumulate(nums_edges_between_points.begin(), nums_edges_between_points.end(), 0);
1480 : const auto max_nodes =
1481 1472 : (nums_edges_between_points.size() == 1 ? n_segments * nums_edges_between_points[0]
1482 108 : : std::accumulate(nums_edges_between_points.begin(),
1483 : nums_edges_between_points.end(),
1484 : 0)) +
1485 1364 : (loop ? 0 : 1);
1486 25127 : for (auto i : make_range(n_elem))
1487 : {
1488 23763 : std::unique_ptr<Elem> elem;
1489 23763 : if (mid_points.size())
1490 : {
1491 1960 : elem = std::make_unique<Edge3>();
1492 1960 : elem->set_node(2, mesh.node_ptr(n_points + i));
1493 : }
1494 : else
1495 21803 : elem = Elem::build(EDGE2);
1496 23763 : const auto ip1 = (i + 1) % max_nodes;
1497 23763 : elem->set_node(0, mesh.node_ptr(i));
1498 23763 : elem->set_node(1, mesh.node_ptr(ip1));
1499 23763 : elem->set_id() = i;
1500 23763 : mesh.add_elem(std::move(elem));
1501 23763 : }
1502 :
1503 1364 : if (!loop)
1504 : {
1505 49 : BoundaryInfo & bi = mesh.get_boundary_info();
1506 196 : std::vector<BoundaryName> bdy_names{start_boundary, end_boundary};
1507 49 : std::vector<boundary_id_type> ids = MooseMeshUtils::getBoundaryIDs(mesh, bdy_names, true);
1508 49 : bi.add_side(mesh.elem_ptr(0), 0, ids[0]);
1509 49 : bi.add_side(mesh.elem_ptr(n_elem - 1), 1, ids[1]);
1510 49 : }
1511 : else
1512 : mooseAssert(start_boundary.empty() && end_boundary.empty(),
1513 : "Cannot assign start/end boundaries on a looped polyline.");
1514 :
1515 1364 : mesh.prepare_for_use();
1516 1413 : }
1517 :
1518 : void
1519 406 : buildPolyLineMesh(MeshBase & mesh,
1520 : const std::vector<Point> & points,
1521 : const bool loop,
1522 : const BoundaryName & start_boundary,
1523 : const BoundaryName & end_boundary,
1524 : const std::vector<unsigned int> & nums_edges_between_points)
1525 : {
1526 406 : buildPolyLineMesh(
1527 : mesh, points, {}, loop, start_boundary, end_boundary, nums_edges_between_points);
1528 406 : }
1529 :
1530 : void
1531 88 : buildPolyLineMesh(MeshBase & mesh,
1532 : const std::vector<Point> & points,
1533 : const bool loop,
1534 : const BoundaryName & start_boundary,
1535 : const BoundaryName & end_boundary,
1536 : const Real max_elem_size)
1537 : {
1538 88 : std::vector<unsigned int> nums_edges_between_points;
1539 88 : const auto n_points = points.size();
1540 648 : for (auto i : make_range(n_points))
1541 : {
1542 560 : if (!loop && (i + 1) == n_points)
1543 0 : break;
1544 :
1545 560 : const auto ip1 = (i + 1) % n_points;
1546 560 : const Real length = (points[ip1] - points[i]).norm();
1547 1120 : const unsigned int n_elems = std::max(
1548 560 : static_cast<unsigned int>(std::ceil(length / max_elem_size)), static_cast<unsigned int>(1));
1549 560 : nums_edges_between_points.push_back(n_elems);
1550 : }
1551 :
1552 88 : buildPolyLineMesh(
1553 : mesh, points, {}, loop, start_boundary, end_boundary, nums_edges_between_points);
1554 88 : }
1555 :
1556 : void
1557 510 : addExternalBoundary(MeshBase & mesh, const BoundaryID extern_bid, bool & has_external_bid)
1558 : {
1559 510 : auto & binfo = mesh.get_boundary_info();
1560 66734 : for (const auto & elem : mesh.active_element_ptr_range())
1561 133368 : for (const auto & i_side : elem->side_index_range())
1562 100256 : if (elem->neighbor_ptr(i_side) == nullptr)
1563 : {
1564 13144 : has_external_bid = true;
1565 13144 : binfo.add_side(elem, i_side, extern_bid);
1566 510 : }
1567 510 : }
1568 : }
|
