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 254 : 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 254 : const auto & side_bd_name_map = mesh.get_boundary_info().get_sideset_name_map();
40 254 : const auto & node_bd_name_map = mesh.get_boundary_info().get_nodeset_name_map();
41 254 : std::map<boundary_id_type, boundary_id_type> same_name_ids;
42 :
43 508 : 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 2515 : for (const auto & pair_outer : map)
47 16200 : for (const auto & pair_inner : map)
48 : // The last condition is needed to make sure we only store one combination
49 15589 : if (pair_outer.second == pair_inner.second && pair_outer.first != pair_inner.first &&
50 15589 : same_ids.find(pair_inner.first) == same_ids.end())
51 698 : same_ids[pair_outer.first] = pair_inner.first;
52 508 : };
53 :
54 254 : populate_map(side_bd_name_map, same_name_ids);
55 254 : populate_map(node_bd_name_map, same_name_ids);
56 :
57 639 : for (const auto & [id1, id2] : same_name_ids)
58 385 : mesh.get_boundary_info().renumber_id(id2, id1);
59 254 : }
60 :
61 : void
62 541 : 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 541 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
69 :
70 : // Container to catch ids passed back from BoundaryInfo
71 541 : std::vector<boundary_id_type> old_ids;
72 :
73 : // Only level-0 elements store BCs. Loop over them.
74 666721 : for (auto & elem : as_range(mesh.level_elements_begin(0), mesh.level_elements_end(0)))
75 : {
76 333090 : unsigned int n_sides = elem->n_sides();
77 2296104 : for (const auto s : make_range(n_sides))
78 : {
79 1963014 : boundary_info.boundary_ids(elem, s, old_ids);
80 1963014 : if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
81 : {
82 18752 : std::vector<boundary_id_type> new_ids(old_ids);
83 18752 : std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
84 18752 : if (delete_prev)
85 : {
86 18299 : boundary_info.remove_side(elem, s);
87 18299 : boundary_info.add_side(elem, s, new_ids);
88 : }
89 : else
90 453 : boundary_info.add_side(elem, s, new_ids);
91 18752 : }
92 : }
93 541 : }
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 541 : if (delete_prev)
99 439 : boundary_info.remove_id(old_id);
100 :
101 : // global information may now be out of sync
102 541 : mesh.set_isnt_prepared();
103 541 : }
104 :
105 : std::vector<boundary_id_type>
106 10438 : getBoundaryIDs(const MeshBase & mesh,
107 : const std::vector<BoundaryName> & boundary_name,
108 : bool generate_unknown)
109 : {
110 : return getBoundaryIDs(
111 10438 : mesh, boundary_name, generate_unknown, mesh.get_boundary_info().get_boundary_ids());
112 : }
113 :
114 : std::vector<boundary_id_type>
115 133075 : 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 133075 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
121 133075 : const std::map<BoundaryID, std::string> & sideset_map = boundary_info.get_sideset_name_map();
122 133075 : const std::map<BoundaryID, std::string> & nodeset_map = boundary_info.get_nodeset_name_map();
123 :
124 133075 : 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 133075 : if (generate_unknown)
131 : {
132 128152 : const auto & bids = mesh.is_prepared() ? mesh.get_boundary_info().get_global_boundary_ids()
133 857 : : mesh.get_boundary_info().get_boundary_ids();
134 128152 : max_boundary_local_id = bids.empty() ? 0 : *(bids.rbegin());
135 : /* We should not hit this often */
136 128152 : if (!mesh.is_prepared() && !mesh.is_serial())
137 91 : mesh.comm().max(max_boundary_local_id);
138 : }
139 :
140 133075 : BoundaryID max_boundary_id = mesh_boundary_ids.empty() ? 0 : *(mesh_boundary_ids.rbegin());
141 :
142 133075 : max_boundary_id =
143 133075 : max_boundary_id > max_boundary_local_id ? max_boundary_id : max_boundary_local_id;
144 :
145 133075 : std::vector<BoundaryID> ids(boundary_name.size());
146 300736 : for (const auto i : index_range(boundary_name))
147 : {
148 167788 : if (boundary_name[i] == "ANY_BOUNDARY_ID")
149 : {
150 127 : ids.assign(mesh_boundary_ids.begin(), mesh_boundary_ids.end());
151 127 : if (i)
152 0 : mooseWarning("You passed \"ANY_BOUNDARY_ID\" in addition to other boundary_names. This "
153 : "may be a logic error.");
154 127 : break;
155 : }
156 :
157 167661 : 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 167661 : 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 333893 : if (generate_unknown &&
171 238695 : !MooseUtils::doesMapContainValue(sideset_map, std::string(boundary_name[i])) &&
172 123908 : !MooseUtils::doesMapContainValue(nodeset_map, std::string(boundary_name[i])))
173 8886 : id = ++max_boundary_id;
174 : else
175 105901 : id = boundary_info.get_id_by_name(boundary_name[i]);
176 : }
177 : else
178 52874 : id = getIDFromName<BoundaryName, BoundaryID>(boundary_name[i]);
179 :
180 167661 : ids[i] = id;
181 : }
182 :
183 266150 : return ids;
184 0 : }
185 :
186 : std::set<BoundaryID>
187 188 : getBoundaryIDSet(const MeshBase & mesh,
188 : const std::vector<BoundaryName> & boundary_name,
189 : bool generate_unknown)
190 : {
191 188 : auto boundaries = getBoundaryIDs(mesh, boundary_name, generate_unknown);
192 376 : return std::set<BoundaryID>(boundaries.begin(), boundaries.end());
193 188 : }
194 :
195 : std::vector<subdomain_id_type>
196 316844 : getSubdomainIDs(const MeshBase & mesh, const std::vector<SubdomainName> & subdomain_names)
197 : {
198 316844 : std::vector<subdomain_id_type> ids;
199 :
200 : // shortcut for "ANY_BLOCK_ID"
201 316844 : 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 61021 : ids.assign(mesh.get_mesh_subdomains().begin(), mesh.get_mesh_subdomains().end());
208 61021 : return ids;
209 : }
210 :
211 : // loop through subdomain names and get IDs (this preserves the order of subdomain_names)
212 255823 : ids.resize(subdomain_names.size());
213 381349 : for (auto i : index_range(subdomain_names))
214 : {
215 125526 : 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 125526 : ids[i] = MooseMeshUtils::getSubdomainID(subdomain_names[i], mesh);
219 : }
220 :
221 255823 : 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 559311 : getBoundaryID(const BoundaryName & boundary_name, const MeshBase & mesh)
234 : {
235 559311 : BoundaryID id = Moose::INVALID_BOUNDARY_ID;
236 559311 : if (boundary_name.empty())
237 0 : return id;
238 :
239 559311 : if (!MooseUtils::isDigits(boundary_name))
240 362627 : id = mesh.get_boundary_info().get_id_by_name(boundary_name);
241 : else
242 196684 : id = getIDFromName<BoundaryName, BoundaryID>(boundary_name);
243 :
244 559307 : return id;
245 : }
246 :
247 : SubdomainID
248 673059 : getSubdomainID(const SubdomainName & subdomain_name, const MeshBase & mesh)
249 : {
250 673059 : if (subdomain_name == "ANY_BLOCK_ID")
251 0 : mooseError("getSubdomainID() does not work with \"ANY_BLOCK_ID\"");
252 :
253 673059 : SubdomainID id = Moose::INVALID_BLOCK_ID;
254 673059 : if (subdomain_name.empty())
255 0 : return id;
256 :
257 673059 : if (!MooseUtils::isDigits(subdomain_name))
258 395044 : id = mesh.get_id_by_name(subdomain_name);
259 : else
260 278015 : id = getIDFromName<SubdomainName, SubdomainID>(subdomain_name);
261 :
262 673059 : 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.set_isnt_prepared();
274 0 : }
275 :
276 : Point
277 391 : meshCentroidCalculator(const MeshBase & mesh)
278 : {
279 391 : Point centroid_pt = Point(0.0, 0.0, 0.0);
280 391 : Real vol_tmp = 0.0;
281 391 : for (const auto & elem :
282 7770 : as_range(mesh.active_local_elements_begin(), mesh.active_local_elements_end()))
283 : {
284 6988 : Real elem_vol = elem->volume();
285 6988 : centroid_pt += (elem->true_centroid()) * elem_vol;
286 6988 : vol_tmp += elem_vol;
287 391 : }
288 391 : mesh.comm().sum(centroid_pt);
289 391 : mesh.comm().sum(vol_tmp);
290 391 : centroid_pt /= vol_tmp;
291 782 : return centroid_pt;
292 : }
293 :
294 : std::unordered_map<dof_id_type, dof_id_type>
295 280 : getExtraIDUniqueCombinationMap(const MeshBase & mesh,
296 : const std::set<SubdomainID> & block_ids,
297 : std::vector<ExtraElementIDName> extra_ids)
298 : {
299 : // check block restriction
300 280 : const bool block_restricted = !block_ids.empty();
301 : // get element id name of interest in recursive parsing algorithm
302 280 : ExtraElementIDName id_name = extra_ids.back();
303 280 : extra_ids.pop_back();
304 280 : const auto id_index = mesh.get_elem_integer_index(id_name);
305 :
306 : // create base parsed id set
307 280 : if (extra_ids.empty())
308 : {
309 : // get set of extra id values;
310 164 : std::vector<dof_id_type> ids;
311 : {
312 164 : std::set<dof_id_type> ids_set;
313 564284 : for (const auto & elem : mesh.active_element_ptr_range())
314 : {
315 564120 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
316 578 : continue;
317 563542 : const auto id = elem->get_extra_integer(id_index);
318 563542 : ids_set.insert(id);
319 164 : }
320 164 : mesh.comm().set_union(ids_set);
321 164 : ids.assign(ids_set.begin(), ids_set.end());
322 164 : }
323 :
324 : // determine new extra id values;
325 164 : std::unordered_map<dof_id_type, dof_id_type> parsed_ids;
326 564284 : for (auto & elem : mesh.active_element_ptr_range())
327 : {
328 564120 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
329 578 : continue;
330 1127084 : parsed_ids[elem->id()] = std::distance(
331 1127084 : ids.begin(), std::lower_bound(ids.begin(), ids.end(), elem->get_extra_integer(id_index)));
332 164 : }
333 164 : return parsed_ids;
334 164 : }
335 :
336 : // if extra_ids is not empty, recursively call getExtraIDUniqueCombinationMap
337 : const auto base_parsed_ids =
338 116 : MooseMeshUtils::getExtraIDUniqueCombinationMap(mesh, block_ids, extra_ids);
339 : // parsing extra ids based on ref_parsed_ids
340 116 : std::vector<std::pair<dof_id_type, dof_id_type>> unique_ids;
341 : {
342 116 : std::set<std::pair<dof_id_type, dof_id_type>> unique_ids_set;
343 428212 : for (const auto & elem : mesh.active_element_ptr_range())
344 : {
345 428096 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
346 528 : continue;
347 427568 : const dof_id_type id1 = libmesh_map_find(base_parsed_ids, elem->id());
348 427568 : const dof_id_type id2 = elem->get_extra_integer(id_index);
349 427568 : const std::pair<dof_id_type, dof_id_type> ids = std::make_pair(id1, id2);
350 427568 : unique_ids_set.insert(ids);
351 116 : }
352 116 : mesh.comm().set_union(unique_ids_set);
353 116 : unique_ids.assign(unique_ids_set.begin(), unique_ids_set.end());
354 116 : }
355 :
356 116 : std::unordered_map<dof_id_type, dof_id_type> parsed_ids;
357 :
358 428212 : for (const auto & elem : mesh.active_element_ptr_range())
359 : {
360 428096 : if (block_restricted && block_ids.find(elem->subdomain_id()) == block_ids.end())
361 528 : continue;
362 427568 : const dof_id_type id1 = libmesh_map_find(base_parsed_ids, elem->id());
363 427568 : const dof_id_type id2 = elem->get_extra_integer(id_index);
364 427568 : const dof_id_type new_id = std::distance(
365 : unique_ids.begin(),
366 427568 : std::lower_bound(unique_ids.begin(), unique_ids.end(), std::make_pair(id1, id2)));
367 427568 : parsed_ids[elem->id()] = new_id;
368 116 : }
369 :
370 116 : return parsed_ids;
371 280 : }
372 :
373 : bool
374 306 : isCoPlanar(const std::vector<Point> vec_pts, const Point plane_nvec, const Point fixed_pt)
375 : {
376 2661 : for (const auto & pt : vec_pts)
377 2359 : if (!MooseUtils::absoluteFuzzyEqual((pt - fixed_pt) * plane_nvec, 0.0))
378 4 : return false;
379 302 : return true;
380 : }
381 :
382 : bool
383 0 : isCoPlanar(const std::vector<Point> vec_pts, const Point plane_nvec)
384 : {
385 0 : return isCoPlanar(vec_pts, plane_nvec, vec_pts.front());
386 : }
387 :
388 : bool
389 0 : isCoPlanar(const std::vector<Point> vec_pts)
390 : {
391 : // Assuming that overlapped Points are allowed, the Points that are overlapped with vec_pts[0] are
392 : // removed before further calculation.
393 0 : std::vector<Point> vec_pts_nonzero{vec_pts[0]};
394 0 : for (const auto i : index_range(vec_pts))
395 0 : if (!MooseUtils::absoluteFuzzyEqual((vec_pts[i] - vec_pts[0]).norm(), 0.0))
396 0 : vec_pts_nonzero.push_back(vec_pts[i]);
397 : // 3 or fewer points are always coplanar
398 0 : if (vec_pts_nonzero.size() <= 3)
399 0 : return true;
400 : else
401 : {
402 0 : for (const auto i : make_range(vec_pts_nonzero.size() - 1))
403 : {
404 0 : const Point tmp_pt = (vec_pts_nonzero[i] - vec_pts_nonzero[0])
405 0 : .cross(vec_pts_nonzero[i + 1] - vec_pts_nonzero[0]);
406 : // if the three points are not collinear, use cross product as the normal vector of the plane
407 0 : if (!MooseUtils::absoluteFuzzyEqual(tmp_pt.norm(), 0.0))
408 0 : return isCoPlanar(vec_pts_nonzero, tmp_pt.unit());
409 : }
410 : }
411 : // If all the points are collinear, they are also coplanar
412 0 : return true;
413 0 : }
414 :
415 : SubdomainID
416 871 : getNextFreeSubdomainID(MeshBase & input_mesh)
417 : {
418 : // Call this to get most up to date block id information
419 871 : input_mesh.cache_elem_data();
420 :
421 871 : std::set<SubdomainID> preexisting_subdomain_ids;
422 871 : input_mesh.subdomain_ids(preexisting_subdomain_ids);
423 871 : if (preexisting_subdomain_ids.empty())
424 0 : return 0;
425 : else
426 : {
427 : const auto highest_subdomain_id =
428 871 : *std::max_element(preexisting_subdomain_ids.begin(), preexisting_subdomain_ids.end());
429 : mooseAssert(highest_subdomain_id < std::numeric_limits<SubdomainID>::max(),
430 : "A SubdomainID with max possible value was found");
431 871 : return highest_subdomain_id + 1;
432 : }
433 871 : }
434 :
435 : BoundaryID
436 775 : getNextFreeBoundaryID(MeshBase & input_mesh)
437 : {
438 775 : auto boundary_ids = input_mesh.get_boundary_info().get_boundary_ids();
439 775 : if (boundary_ids.empty())
440 116 : return 0;
441 659 : return (*boundary_ids.rbegin() + 1);
442 775 : }
443 :
444 : bool
445 13246 : hasSubdomainID(const MeshBase & input_mesh, const SubdomainID & id)
446 : {
447 13246 : std::set<SubdomainID> mesh_blocks;
448 13246 : input_mesh.subdomain_ids(mesh_blocks);
449 :
450 : // On a distributed mesh we may have sideset IDs that only exist on
451 : // other processors
452 13246 : if (!input_mesh.is_replicated())
453 1238 : input_mesh.comm().set_union(mesh_blocks);
454 :
455 26492 : return mesh_blocks.count(id) && (id != Moose::INVALID_BLOCK_ID);
456 13246 : }
457 :
458 : bool
459 10734 : hasSubdomainName(const MeshBase & input_mesh, const SubdomainName & name)
460 : {
461 10734 : const auto id = getSubdomainID(name, input_mesh);
462 21468 : return hasSubdomainID(input_mesh, id);
463 : }
464 :
465 : bool
466 4378 : hasBoundaryID(const MeshBase & input_mesh, const BoundaryID id)
467 : {
468 4378 : const BoundaryInfo & boundary_info = input_mesh.get_boundary_info();
469 4378 : std::set<boundary_id_type> boundary_ids = boundary_info.get_boundary_ids();
470 :
471 : // On a distributed mesh we may have boundary IDs that only exist on
472 : // other processors
473 4378 : if (!input_mesh.is_replicated())
474 726 : input_mesh.comm().set_union(boundary_ids);
475 :
476 8756 : return boundary_ids.count(id) && (id != Moose::INVALID_BOUNDARY_ID);
477 4378 : }
478 :
479 : bool
480 4204 : hasBoundaryName(const MeshBase & input_mesh, const BoundaryName & name)
481 : {
482 4204 : const auto id = getBoundaryID(name, input_mesh);
483 4204 : return hasBoundaryID(input_mesh, id);
484 : }
485 :
486 : void
487 496 : makeOrderedNodeList(std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
488 : std::vector<dof_id_type> & elem_id_list,
489 : std::vector<dof_id_type> & midpoint_node_list,
490 : std::vector<dof_id_type> & ordered_node_list,
491 : std::vector<dof_id_type> & ordered_elem_id_list)
492 : {
493 : // a flag to indicate if the ordered_node_list has been reversed
494 496 : bool is_flipped = false;
495 : // Start from the first element, try to find a chain of nodes
496 : mooseAssert(node_assm.size(), "Node list must not be empty");
497 496 : ordered_node_list.push_back(node_assm.front().first);
498 496 : if (midpoint_node_list.front() != DofObject::invalid_id)
499 0 : ordered_node_list.push_back(midpoint_node_list.front());
500 496 : ordered_node_list.push_back(node_assm.front().second);
501 496 : ordered_elem_id_list.push_back(elem_id_list.front());
502 : // Remove the element that has just been added to ordered_node_list
503 496 : node_assm.erase(node_assm.begin());
504 496 : midpoint_node_list.erase(midpoint_node_list.begin());
505 496 : elem_id_list.erase(elem_id_list.begin());
506 496 : const unsigned int node_assm_size_0 = node_assm.size();
507 2889 : for (unsigned int i = 0; i < node_assm_size_0; i++)
508 : {
509 : // Find nodes to expand the chain
510 2397 : dof_id_type end_node_id = ordered_node_list.back();
511 8058 : auto isMatch1 = [end_node_id](std::pair<dof_id_type, dof_id_type> old_id_pair)
512 8058 : { return old_id_pair.first == end_node_id; };
513 2146 : auto isMatch2 = [end_node_id](std::pair<dof_id_type, dof_id_type> old_id_pair)
514 2146 : { return old_id_pair.second == end_node_id; };
515 2397 : auto result = std::find_if(node_assm.begin(), node_assm.end(), isMatch1);
516 : bool match_first;
517 2397 : if (result == node_assm.end())
518 : {
519 1214 : match_first = false;
520 1214 : result = std::find_if(node_assm.begin(), node_assm.end(), isMatch2);
521 : }
522 : else
523 : {
524 1183 : match_first = true;
525 : }
526 : // If found, add the node to boundary_ordered_node_list
527 2397 : if (result != node_assm.end())
528 : {
529 2139 : const auto elem_index = std::distance(node_assm.begin(), result);
530 2139 : if (midpoint_node_list[elem_index] != DofObject::invalid_id)
531 0 : ordered_node_list.push_back(midpoint_node_list[elem_index]);
532 2139 : ordered_node_list.push_back(match_first ? (*result).second : (*result).first);
533 2139 : node_assm.erase(result);
534 2139 : midpoint_node_list.erase(midpoint_node_list.begin() + elem_index);
535 2139 : ordered_elem_id_list.push_back(elem_id_list[elem_index]);
536 2139 : elem_id_list.erase(elem_id_list.begin() + elem_index);
537 : }
538 : // If there are still elements in node_assm and result ==
539 : // node_assm.end(), this means the curve is not a loop, the
540 : // ordered_node_list is flipped and try the other direction that has not
541 : // been examined yet.
542 : else
543 : {
544 258 : if (is_flipped)
545 : // Flipped twice; this means the node list has at least two segments.
546 4 : throw MooseException("The node list provided has more than one segments.");
547 :
548 : // mark the first flip event.
549 254 : is_flipped = true;
550 254 : std::reverse(ordered_node_list.begin(), ordered_node_list.end());
551 254 : std::reverse(midpoint_node_list.begin(), midpoint_node_list.end());
552 254 : std::reverse(ordered_elem_id_list.begin(), ordered_elem_id_list.end());
553 : // As this iteration is wasted, set the iterator backward
554 254 : i--;
555 : }
556 : }
557 492 : }
558 :
559 : void
560 234 : makeOrderedNodeList(std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
561 : std::vector<dof_id_type> & elem_id_list,
562 : std::vector<dof_id_type> & ordered_node_list,
563 : std::vector<dof_id_type> & ordered_elem_id_list)
564 : {
565 234 : std::vector<dof_id_type> dummy_midpoint_node_list(node_assm.size(), DofObject::invalid_id);
566 234 : makeOrderedNodeList(
567 : node_assm, elem_id_list, dummy_midpoint_node_list, ordered_node_list, ordered_elem_id_list);
568 234 : }
569 :
570 : void
571 9036 : swapNodesInElem(Elem & elem, const unsigned int nd1, const unsigned int nd2)
572 : {
573 9036 : Node * n_temp = elem.node_ptr(nd1);
574 9036 : elem.set_node(nd1, elem.node_ptr(nd2));
575 9036 : elem.set_node(nd2, n_temp);
576 9036 : }
577 :
578 : void
579 241 : extraElemIntegerSwapParametersProcessor(
580 : const std::string & class_name,
581 : const unsigned int num_sections,
582 : const unsigned int num_integers,
583 : const std::vector<std::vector<std::vector<dof_id_type>>> & elem_integers_swaps,
584 : std::vector<std::unordered_map<dof_id_type, dof_id_type>> & elem_integers_swap_pairs)
585 : {
586 241 : elem_integers_swap_pairs.reserve(num_sections * num_integers);
587 267 : for (const auto i : make_range(num_integers))
588 : {
589 26 : const auto & elem_integer_swaps = elem_integers_swaps[i];
590 26 : std::vector<std::unordered_map<dof_id_type, dof_id_type>> elem_integer_swap_pairs;
591 : try
592 : {
593 52 : MooseMeshUtils::idSwapParametersProcessor(class_name,
594 : "elem_integers_swaps",
595 : elem_integer_swaps,
596 : elem_integer_swap_pairs,
597 : i * num_sections);
598 : }
599 0 : catch (const MooseException & e)
600 : {
601 0 : throw MooseException(e.what());
602 0 : }
603 :
604 26 : elem_integers_swap_pairs.insert(elem_integers_swap_pairs.end(),
605 : elem_integer_swap_pairs.begin(),
606 : elem_integer_swap_pairs.end());
607 26 : }
608 241 : }
609 :
610 : std::unique_ptr<ReplicatedMesh>
611 0 : buildBoundaryMesh(const ReplicatedMesh & input_mesh, const boundary_id_type boundary_id)
612 : {
613 0 : auto poly_mesh = std::make_unique<ReplicatedMesh>(input_mesh.comm());
614 :
615 0 : auto side_list = input_mesh.get_boundary_info().build_side_list();
616 :
617 0 : std::unordered_map<dof_id_type, dof_id_type> old_new_node_map;
618 0 : for (const auto & bside : side_list)
619 : {
620 0 : if (std::get<2>(bside) != boundary_id)
621 0 : continue;
622 :
623 0 : const Elem * elem = input_mesh.elem_ptr(std::get<0>(bside));
624 0 : const auto side = std::get<1>(bside);
625 0 : auto side_elem = elem->build_side_ptr(side);
626 0 : auto copy = side_elem->build(side_elem->type());
627 :
628 0 : for (const auto i : side_elem->node_index_range())
629 : {
630 0 : auto & n = side_elem->node_ref(i);
631 :
632 0 : if (old_new_node_map.count(n.id()))
633 0 : copy->set_node(i, poly_mesh->node_ptr(old_new_node_map[n.id()]));
634 : else
635 : {
636 0 : Node * node = poly_mesh->add_point(side_elem->point(i));
637 0 : copy->set_node(i, node);
638 0 : old_new_node_map[n.id()] = node->id();
639 : }
640 : }
641 0 : poly_mesh->add_elem(copy.release());
642 0 : }
643 0 : poly_mesh->skip_partitioning(true);
644 0 : poly_mesh->prepare_for_use();
645 0 : if (poly_mesh->n_elem() == 0)
646 0 : mooseError("The input mesh does not have a boundary with id ", boundary_id);
647 :
648 0 : return poly_mesh;
649 0 : }
650 :
651 : void
652 2435 : createSubdomainFromSidesets(std::unique_ptr<MeshBase> & mesh,
653 : std::vector<BoundaryName> boundary_names,
654 : const SubdomainID new_subdomain_id,
655 : const SubdomainName new_subdomain_name,
656 : const std::string type_name)
657 : {
658 : // Generate a new block id if one isn't supplied.
659 2435 : SubdomainID new_block_id = new_subdomain_id;
660 :
661 : // Make sure our boundary info and parallel counts are setup
662 2435 : if (!mesh->is_prepared())
663 : {
664 539 : const bool allow_remote_element_removal = mesh->allow_remote_element_removal();
665 : // We want all of our boundary elements available, so avoid removing them if they haven't
666 : // already been so
667 539 : mesh->allow_remote_element_removal(false);
668 539 : mesh->prepare_for_use();
669 539 : mesh->allow_remote_element_removal(allow_remote_element_removal);
670 : }
671 :
672 : // Check that the sidesets are present in the mesh
673 5039 : for (const auto & sideset : boundary_names)
674 2616 : if (!MooseMeshUtils::hasBoundaryName(*mesh, sideset))
675 12 : mooseException("The sideset '", sideset, "' was not found within the mesh");
676 :
677 2423 : auto sideset_ids = MooseMeshUtils::getBoundaryIDs(*mesh, boundary_names, true);
678 2423 : std::set<boundary_id_type> sidesets(sideset_ids.begin(), sideset_ids.end());
679 2423 : auto side_list = mesh->get_boundary_info().build_side_list();
680 2423 : if (!mesh->is_serial() && mesh->comm().size() > 1)
681 : {
682 10 : std::vector<Elem *> elements_to_send;
683 10 : unsigned short i_need_boundary_elems = 0;
684 172 : for (const auto & [elem_id, side, bc_id] : side_list)
685 : {
686 162 : libmesh_ignore(side);
687 162 : if (sidesets.count(bc_id))
688 : {
689 : // Whether we have this boundary information through our locally owned element or a ghosted
690 : // element, we'll need the boundary elements for parallel consistent addition
691 92 : i_need_boundary_elems = 1;
692 92 : auto * elem = mesh->elem_ptr(elem_id);
693 92 : if (elem->processor_id() == mesh->processor_id())
694 80 : elements_to_send.push_back(elem);
695 : }
696 : }
697 :
698 : std::set<const Elem *, libMesh::CompareElemIdsByLevel> connected_elements(
699 10 : elements_to_send.begin(), elements_to_send.end());
700 10 : std::set<const Node *> connected_nodes;
701 10 : reconnect_nodes(connected_elements, connected_nodes);
702 10 : std::set<dof_id_type> connected_node_ids;
703 178 : for (auto * nd : connected_nodes)
704 168 : connected_node_ids.insert(nd->id());
705 :
706 10 : std::vector<unsigned short> need_boundary_elems(mesh->comm().size());
707 10 : mesh->comm().allgather(i_need_boundary_elems, need_boundary_elems);
708 10 : std::unordered_map<processor_id_type, decltype(elements_to_send)> push_element_data;
709 10 : std::unordered_map<processor_id_type, decltype(connected_nodes)> push_node_data;
710 :
711 30 : for (const auto pid : index_range(mesh->comm()))
712 : // Don't need to send to self
713 20 : if (pid != mesh->processor_id() && need_boundary_elems[pid])
714 : {
715 10 : if (elements_to_send.size())
716 10 : push_element_data[pid] = elements_to_send;
717 10 : if (connected_nodes.size())
718 10 : push_node_data[pid] = connected_nodes;
719 : }
720 :
721 10 : auto node_action_functor = [](processor_id_type, const auto &)
722 : {
723 : // Node packing specialization already has unpacked node into mesh, so nothing to do
724 10 : };
725 20 : Parallel::push_parallel_packed_range(
726 10 : mesh->comm(), push_node_data, mesh.get(), node_action_functor);
727 10 : auto elem_action_functor = [](processor_id_type, const auto &)
728 : {
729 : // Elem packing specialization already has unpacked elem into mesh, so nothing to do
730 10 : };
731 20 : TIMPI::push_parallel_packed_range(
732 10 : mesh->comm(), push_element_data, mesh.get(), elem_action_functor);
733 :
734 : // now that we've gathered everything, we need to rebuild the side list
735 10 : side_list = mesh->get_boundary_info().build_side_list();
736 10 : }
737 :
738 2423 : std::vector<std::pair<dof_id_type, ElemSidePair>> element_sides_on_boundary;
739 2423 : dof_id_type counter = 0;
740 155305 : for (const auto & triple : side_list)
741 152886 : if (sidesets.count(std::get<2>(triple)))
742 : {
743 24511 : if (auto elem = mesh->query_elem_ptr(std::get<0>(triple)))
744 : {
745 24511 : if (!elem->active())
746 4 : mooseError(
747 : "Only active, level 0 elements can be made interior parents of new level 0 lower-d "
748 : "elements. Make sure that ",
749 : type_name,
750 : "s are run before any refinement generators");
751 24507 : element_sides_on_boundary.push_back(
752 49014 : std::make_pair(counter, ElemSidePair(elem, std::get<1>(triple))));
753 : }
754 24507 : ++counter;
755 : }
756 :
757 2419 : dof_id_type max_elem_id = mesh->max_elem_id();
758 2419 : unique_id_type max_unique_id = mesh->parallel_max_unique_id();
759 :
760 : // Making an important assumption that at least our boundary elements are the same on all
761 : // processes even in distributed mesh mode (this is reliant on the correct ghosting functors
762 : // existing on the mesh)
763 26926 : for (auto & [i, elem_side] : element_sides_on_boundary)
764 : {
765 24507 : Elem * elem = elem_side.elem;
766 :
767 24507 : const auto side = elem_side.side;
768 :
769 : // Build a non-proxy element from this side.
770 24507 : std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
771 :
772 : // The side will be added with the same processor id as the parent.
773 24507 : side_elem->processor_id() = elem->processor_id();
774 :
775 : // Add subdomain ID
776 24507 : side_elem->subdomain_id() = new_block_id;
777 :
778 : // Also assign the side's interior parent, so it is always
779 : // easy to figure out the Elem we came from.
780 24507 : side_elem->set_interior_parent(elem);
781 :
782 : // Add id
783 24507 : side_elem->set_id(max_elem_id + i);
784 24507 : side_elem->set_unique_id(max_unique_id + i);
785 :
786 : // Finally, add the lower-dimensional element to the mesh->
787 24507 : mesh->add_elem(side_elem.release());
788 24507 : };
789 :
790 : // Assign block name, if provided
791 2419 : if (new_subdomain_name.size())
792 1485 : mesh->subdomain_name(new_block_id) = new_subdomain_name;
793 :
794 2419 : const bool skip_partitioning_old = mesh->skip_partitioning();
795 2419 : mesh->skip_partitioning(true);
796 2419 : mesh->prepare_for_use();
797 2419 : mesh->skip_partitioning(skip_partitioning_old);
798 2419 : }
799 :
800 : void
801 316 : convertBlockToMesh(std::unique_ptr<MeshBase> & source_mesh,
802 : std::unique_ptr<MeshBase> & target_mesh,
803 : const std::vector<SubdomainName> & target_blocks)
804 : {
805 316 : if (!source_mesh->is_replicated())
806 0 : mooseError("This generator does not support distributed meshes.");
807 :
808 316 : const auto target_block_ids = MooseMeshUtils::getSubdomainIDs(*source_mesh, target_blocks);
809 :
810 : // Check that the block ids/names exist in the mesh
811 316 : std::set<SubdomainID> mesh_blocks;
812 316 : source_mesh->subdomain_ids(mesh_blocks);
813 :
814 706 : for (const auto i : index_range(target_block_ids))
815 390 : if (target_block_ids[i] == Moose::INVALID_BLOCK_ID || !mesh_blocks.count(target_block_ids[i]))
816 : {
817 0 : mooseException("The target_block '", target_blocks[i], "' was not found within the mesh.");
818 : }
819 :
820 : // know which nodes have already been inserted, by tracking the old mesh's node's ids'
821 316 : std::unordered_map<dof_id_type, dof_id_type> old_new_node_map;
822 :
823 686 : for (const auto target_block_id : target_block_ids)
824 : {
825 :
826 8848 : for (auto elem : source_mesh->active_subdomain_elements_ptr_range(target_block_id))
827 : {
828 4241 : if (elem->level() != 0)
829 4 : mooseError("Refined blocks are not supported by this generator. "
830 : "Can you re-organize mesh generators to refine after converting the block?");
831 :
832 : // make a deep copy so that mutiple meshes' destructors don't segfault at program termination
833 4237 : auto copy = elem->build(elem->type());
834 :
835 : // index of node in the copy element must be managed manually as there is no intelligent
836 : // insert method
837 4237 : dof_id_type copy_n_index = 0;
838 :
839 : // correctly assign new copies of nodes, loop over nodes
840 19819 : for (dof_id_type i : elem->node_index_range())
841 : {
842 15582 : auto & n = elem->node_ref(i);
843 :
844 15582 : if (old_new_node_map.count(n.id()))
845 : {
846 : // case where we have already inserted this particular point before
847 : // then we need to find the already-inserted one and hook it up right
848 : // to it's respective element
849 8921 : copy->set_node(copy_n_index++, target_mesh->node_ptr(old_new_node_map[n.id()]));
850 : }
851 : else
852 : {
853 : // case where we've NEVER inserted this particular point before
854 : // add them both to the element and the mesh
855 :
856 : // Nodes' IDs are their indexes in the nodes' respective mesh
857 : // If we set them as invalid they are automatically assigned
858 : // Add to mesh, auto-assigning a new id.
859 6661 : Node * node = target_mesh->add_point(elem->point(i));
860 :
861 : // Add to element copy (manually)
862 6661 : copy->set_node(copy_n_index++, node);
863 :
864 : // remember the (old) ID
865 6661 : old_new_node_map[n.id()] = node->id();
866 : }
867 : }
868 :
869 : // it is ok to release the copy element into the mesh because derived meshes class
870 : // (ReplicatedMesh, DistributedMesh) manage their own elements, will delete them
871 4237 : target_mesh->add_elem(copy.release());
872 4607 : }
873 : }
874 312 : }
875 : }
|
