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