Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "XYZDelaunayGenerator.h"
11 :
12 : #include "CastUniquePointer.h"
13 : #include "MooseMeshUtils.h"
14 : #include "MooseUtils.h"
15 : #include "MooseMeshElementConversionUtils.h"
16 :
17 : #include "libmesh/elem.h"
18 : #include "libmesh/int_range.h"
19 : #include "libmesh/mesh_modification.h"
20 : #include "libmesh/mesh_netgen_interface.h"
21 : #include "libmesh/mesh_serializer.h"
22 : #include "libmesh/parsed_function.h"
23 : #include "libmesh/replicated_mesh.h"
24 :
25 : registerMooseObject("MooseApp", XYZDelaunayGenerator);
26 :
27 : namespace std
28 : {
29 : template <>
30 : struct hash<std::tuple<libMesh::Point, libMesh::Point, libMesh::Point>>
31 : {
32 15960 : std::size_t operator()(const std::tuple<libMesh::Point, libMesh::Point, libMesh::Point> & p) const
33 : {
34 15960 : std::size_t seed = 0;
35 15960 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<0>(p)));
36 15960 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<1>(p)));
37 15960 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<2>(p)));
38 15960 : return seed;
39 : }
40 : };
41 :
42 : } // namespace std
43 :
44 : InputParameters
45 3501 : XYZDelaunayGenerator::validParams()
46 : {
47 3501 : InputParameters params = MeshGenerator::validParams();
48 :
49 14004 : MooseEnum algorithm("BINARY EXHAUSTIVE", "BINARY");
50 :
51 14004 : params.addRequiredParam<MeshGeneratorName>(
52 : "boundary",
53 : "The input MeshGenerator defining the output outer boundary. The input mesh (the output mesh "
54 : "of the input mesh generator) can either "
55 : "include 3D volume elements or 2D surface elements.");
56 :
57 14004 : params.addParam<SubdomainName>("output_subdomain_name",
58 : "Subdomain name to set on new triangles.");
59 :
60 14004 : params.addParam<BoundaryName>(
61 : "output_boundary",
62 : "Boundary name to set on new outer boundary. Default ID: 0 if no hole meshes are stitched; "
63 : "or maximum boundary ID of all the stitched hole meshes + 1.");
64 14004 : params.addParam<std::vector<BoundaryName>>(
65 : "hole_boundaries",
66 : "Boundary names to set on holes. Default IDs are numbered up from 1 if no hole meshes are "
67 : "stitched; or from maximum boundary ID of all the stitched hole meshes + 2.");
68 :
69 10503 : params.addParam<bool>("smooth_triangulation",
70 7002 : false,
71 : "Whether to do Laplacian mesh smoothing on the generated triangles.");
72 14004 : params.addParam<std::vector<MeshGeneratorName>>(
73 : "holes",
74 : {},
75 : "The MeshGenerators that create meshes defining the holes. A hole mesh must contain either "
76 : "3D volume "
77 : "elements where the external surface of the mesh works as the closed manifold that defines "
78 : "the hole, or 2D surface elements that form the closed manifold that defines the hole.");
79 10503 : params.addParam<std::vector<bool>>(
80 7002 : "stitch_holes", std::vector<bool>(), "Whether to stitch to the mesh defining each hole.");
81 10503 : params.addParam<bool>("convert_holes_for_stitching",
82 7002 : false,
83 : "Whether to convert the 3D hole meshes with non-TRI3 surface sides into "
84 : "a compatible form.");
85 :
86 14004 : MooseEnum conversion_method("ALL SURFACE", "ALL");
87 14004 : params.addParam<MooseEnum>(
88 : "conversion_method",
89 : conversion_method,
90 : "The method to convert 3D hole meshes into compatible meshes. Options are "
91 : "ALL: convert all elements into TET4; SURFACE: convert only the surface elements "
92 : "that are stitched to the generated Delaunay mesh.");
93 :
94 21006 : params.addRangeCheckedParam<Real>(
95 : "desired_volume",
96 : 0,
97 : "desired_volume>=0",
98 : "Desired (maximum) tetrahedral volume, or 0 to skip uniform refinement");
99 :
100 10503 : params.addParam<bool>(
101 : "combined_stitching",
102 7002 : false,
103 : "Whether to stitch all holes in one combined stitching step. This is efficient if a great "
104 : "number of holes are to be stitched. But it is hard to debug if problems arise.");
105 14004 : params.addParam<MooseEnum>(
106 : "algorithm",
107 : algorithm,
108 : "Control the use of binary search for the nodes of the stitched surfaces.");
109 21006 : params.renameParam("algorithm", "stitching_algorithm", "12/10/2026");
110 10503 : params.addParam<bool>(
111 7002 : "verbose_stitching", false, "Whether mesh hole stitching should have verbose output.");
112 :
113 3501 : params.addClassDescription(
114 : "Creates tetrahedral 3D meshes within boundaries defined by input meshes.");
115 :
116 7002 : return params;
117 3501 : }
118 :
119 223 : XYZDelaunayGenerator::XYZDelaunayGenerator(const InputParameters & parameters)
120 : : MeshGenerator(parameters),
121 223 : _bdy_ptr(getMesh("boundary")),
122 446 : _desired_volume(getParam<Real>("desired_volume")),
123 223 : _output_subdomain_id(0),
124 446 : _smooth_tri(getParam<bool>("smooth_triangulation")),
125 446 : _hole_ptrs(getMeshes("holes")),
126 446 : _stitch_holes(getParam<std::vector<bool>>("stitch_holes")),
127 446 : _convert_holes_for_stitching(getParam<bool>("convert_holes_for_stitching")),
128 223 : _conversion_method(parameters.get<MooseEnum>("conversion_method")),
129 223 : _combined_stitching(parameters.get<bool>("combined_stitching")),
130 223 : _algorithm(parameters.get<MooseEnum>("algorithm")),
131 446 : _verbose_stitching(parameters.get<bool>("verbose_stitching"))
132 : {
133 223 : if (!_stitch_holes.empty() && _stitch_holes.size() != _hole_ptrs.size())
134 6 : paramError("stitch_holes", "Need one stitch_holes entry per hole, if specified.");
135 :
136 660 : if (isParamValid("hole_boundaries"))
137 : {
138 66 : auto & hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
139 33 : if (hole_boundaries.size() != _hole_ptrs.size())
140 6 : paramError("hole_boundaries", "Need one hole_boundaries entry per hole, if specified.");
141 : }
142 :
143 651 : if (isParamSetByUser("conversion_method") && !_convert_holes_for_stitching)
144 0 : paramError(
145 : "conversion_method",
146 : "This parameter is only applicable when convert_holes_for_stitching is set to true.");
147 217 : }
148 :
149 : std::unique_ptr<MeshBase>
150 201 : XYZDelaunayGenerator::generate()
151 : {
152 : #ifdef LIBMESH_HAVE_NETGEN
153 : // Put the boundary mesh in a local pointer
154 : std::unique_ptr<UnstructuredMesh> mesh =
155 201 : dynamic_pointer_cast<UnstructuredMesh>(std::move(_bdy_ptr));
156 :
157 : // The libMesh Netgen removes all the sideset info
158 : // But it keeps some nodeset info for the retained nodes
159 : // We need to clear these nodeset info as they could overlap with upcoming boundary info
160 201 : mesh->get_boundary_info().clear_boundary_node_ids();
161 :
162 : // Get ready to triangulate its boundary
163 201 : libMesh::NetGenMeshInterface ngint(*mesh);
164 :
165 201 : ngint.smooth_after_generating() = _smooth_tri;
166 :
167 201 : ngint.desired_volume() = _desired_volume;
168 :
169 : // The hole meshes will be used for hole boundary identification and optionally for stitching.
170 : // if a hole mesh contains 3D volume elements but has non-TRI3 surface side elements, it cannot be
171 : // used directly for stitching. But it can be converted into an all-TET4 mesh to support hole
172 : // boundary identification
173 447 : for (const auto hole_i : index_range(_hole_ptrs))
174 : {
175 258 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
176 258 : libMesh::MeshSerializer serial_hole(hole_mesh);
177 : // Check the dimension of the hole mesh
178 : // We do not need to worry about element order here as libMesh checks it
179 258 : std::set<ElemType> hole_elem_types;
180 258 : std::set<unsigned short> hole_elem_dims;
181 258 : std::vector<std::pair<dof_id_type, unsigned int>> hole_elem_external_sides;
182 8206 : for (auto elem : hole_mesh.element_ptr_range())
183 : {
184 7948 : hole_elem_dims.emplace(elem->dim());
185 :
186 : // For a 3D element, we need to check the surface side element type instead of the element
187 : // type. It is also a good opportunity to define the external boundary.
188 7948 : if (elem->dim() == 3)
189 39976 : for (auto s : make_range(elem->n_sides()))
190 : {
191 : // Note that the entire external boundary is used for defining the hole at this time
192 32214 : if (!elem->neighbor_ptr(s))
193 : {
194 4854 : hole_elem_types.emplace(elem->side_ptr(s)->type());
195 4854 : hole_elem_external_sides.emplace_back(elem->id(), s);
196 : }
197 : }
198 : // For a non-3D element, we just need to record the element type
199 : else
200 186 : hole_elem_types.emplace(elem->type());
201 258 : }
202 258 : if (hole_elem_dims.size() != 1 || *hole_elem_dims.begin() < 2)
203 12 : paramError(
204 : "holes",
205 : "All elements in a hole mesh must have the same dimension that is either 2D or 3D.");
206 252 : else if (*hole_elem_dims.begin() == 3)
207 : {
208 : // For 3D meshes, if there are non-TRI3 surface side elements
209 : // (1) if no stitching is needed, we can just convert the whole mesh into TET to facilitate
210 : // boundary identification (2) if stitching is needed, we can still convert and stitch, but
211 : // that would modify the input hole mesh
212 229 : if (*hole_elem_types.begin() != ElemType::TRI3 || hole_elem_types.size() > 1)
213 : {
214 63 : if (_stitch_holes.size() && _stitch_holes[hole_i] && !_convert_holes_for_stitching)
215 6 : paramError("holes",
216 : "3D hole meshes with non-TRI3 surface elements cannot be stitched without "
217 : "converting them to TET4. Consider setting convert_holes_for_stitching=true.");
218 60 : else if (_stitch_holes.size() && _stitch_holes[hole_i] && _conversion_method == "SURFACE")
219 : {
220 : // Create a transition layer with triangle sides on the external boundary of the hole
221 10 : BoundaryID temp_ext_bid = MooseMeshUtils::getNextFreeBoundaryID(hole_mesh);
222 550 : for (const auto & hees : hole_elem_external_sides)
223 540 : hole_mesh.get_boundary_info().add_side(hees.first, hees.second, temp_ext_bid);
224 10 : MooseMeshElementConversionUtils::transitionLayerGenerator(
225 20 : hole_mesh, std::vector<BoundaryName>({std::to_string(temp_ext_bid)}), 1, false);
226 10 : hole_mesh.get_boundary_info().remove_id(temp_ext_bid);
227 : // MeshSerializer's destructor calls delete_remote_elements()
228 : // For replicated meshes, we need to refresh the neighbor information
229 : // Also, not doing prepare_for_use() causes an error in DEBUG mode
230 : // in MeshTools::libmesh_assert_topology_consistent_procids()
231 10 : hole_mesh.prepare_for_use();
232 : }
233 : else
234 50 : MeshTools::Modification::all_tri(**_hole_ptrs[hole_i]);
235 : }
236 : }
237 : else // if (*hole_elem_dims.begin() == 2)
238 : {
239 : // There is no point to stitch a 2D hole mesh, so we throw an error if this is attempted
240 23 : if (_stitch_holes.size() && _stitch_holes[hole_i])
241 3 : paramError("holes",
242 3 : "the hole mesh with index " + std::to_string(hole_i) +
243 : " is a 2D mesh, for which stitching onto a 3D mesh does not make sense.");
244 : }
245 246 : }
246 :
247 : std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh>>> ngholes =
248 189 : std::make_unique<std::vector<std::unique_ptr<UnstructuredMesh>>>();
249 :
250 : // The libMesh Netgen interface will modify hole meshes in-place, so
251 : // we make copies to pass in.
252 435 : for (std::unique_ptr<MeshBase> * hole_ptr : _hole_ptrs)
253 : {
254 : // How did we never add a ReplicatedMesh(MeshBase&) constructor in
255 : // libMesh?
256 246 : const UnstructuredMesh & hole = dynamic_cast<UnstructuredMesh &>(**hole_ptr);
257 246 : ngholes->push_back(std::make_unique<ReplicatedMesh>(hole));
258 : }
259 :
260 189 : if (!_hole_ptrs.empty())
261 153 : ngint.attach_hole_list(std::move(ngholes));
262 :
263 189 : ngint.triangulate();
264 :
265 567 : if (isParamValid("output_subdomain_name"))
266 : {
267 102 : auto output_subdomain_name = getParam<SubdomainName>("output_subdomain_name");
268 51 : _output_subdomain_id = MooseMeshUtils::getSubdomainID(output_subdomain_name, *mesh);
269 :
270 51 : if (_output_subdomain_id == Elem::invalid_subdomain_id)
271 : {
272 : // We'll probably need to make a new ID, then
273 51 : _output_subdomain_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
274 :
275 : // But check the hole meshes for our output subdomain name too
276 137 : for (auto & hole_ptr : _hole_ptrs)
277 : {
278 86 : auto possible_sbdid = MooseMeshUtils::getSubdomainID(output_subdomain_name, **hole_ptr);
279 : // Huh, it was in one of them
280 86 : if (possible_sbdid != Elem::invalid_subdomain_id)
281 : {
282 0 : _output_subdomain_id = possible_sbdid;
283 0 : break;
284 : }
285 86 : _output_subdomain_id =
286 86 : std::max(_output_subdomain_id, MooseMeshUtils::getNextFreeSubdomainID(**hole_ptr));
287 : }
288 :
289 51 : mesh->subdomain_name(_output_subdomain_id) = output_subdomain_name;
290 : }
291 51 : }
292 :
293 189 : if (_smooth_tri || _output_subdomain_id)
294 27756 : for (auto elem : mesh->element_ptr_range())
295 : {
296 27675 : elem->subdomain_id() = _output_subdomain_id;
297 :
298 : // I do not trust Laplacian mesh smoothing not to invert
299 : // elements near reentrant corners. Eventually we'll add better
300 : // smoothing options, but even those might have failure cases.
301 : // Better to always do extra tests here than to ever let users
302 : // try to run on a degenerate mesh.
303 27675 : if (elem->is_flipped())
304 : {
305 0 : if (_smooth_tri)
306 0 : mooseError("Inverted element found in triangulation.\n"
307 : "Laplacian smoothing can create these at reentrant corners; disable it?");
308 : else
309 0 : mooseError("Unexplained inverted element found in triangulation.\n");
310 : }
311 81 : }
312 :
313 189 : const bool use_binary_search = (_algorithm == "BINARY");
314 :
315 : // The hole meshes are specified by the user, so they could have any
316 : // BCID or no BCID or any combination of BCIDs on their outer
317 : // boundary, so we'll have to set our own BCID to use for stitching
318 : // there. We'll need to check all the holes for used BCIDs, if we
319 : // want to pick a new ID on hole N that doesn't conflict with any
320 : // IDs on hole M < N (or with the IDs on the new triangulation)
321 :
322 : // The Netgen generated mesh does not have hole and output boundary ids,
323 : // which is different from its 2D counterpart.
324 : // So, we need to assign boundary ids to the hole boundaries and output boundary
325 : // By default, we assign the hole boundary ids to be 1,2,..., N
326 : // and the output boundary id to be 0
327 : // We will need a different set of boundary ids if we are stitching holes
328 :
329 : // end_bcid is one more than the maximum of default output boundary id and hole boundary ids
330 189 : const boundary_id_type end_bcid = _hole_ptrs.size() + 1;
331 :
332 : // For the hole meshes that need to be stitched, we would like to make sure the hole boundary ids
333 : // and output boundary id are not conflicting with the existing boundary ids of the hole meshes to
334 : // be stitched.
335 189 : BoundaryID free_boundary_id = 0;
336 : // We need to get a free boundary id that is larger than the existing boundary ids in any hole
337 : // meshes
338 189 : if (_stitch_holes.size())
339 : {
340 279 : for (auto hole_i : index_range(_hole_ptrs))
341 : {
342 166 : if (_stitch_holes[hole_i])
343 : {
344 136 : free_boundary_id =
345 136 : std::max(free_boundary_id, MooseMeshUtils::getNextFreeBoundaryID(**_hole_ptrs[hole_i]));
346 136 : (*_hole_ptrs[hole_i])->comm().max(free_boundary_id);
347 : }
348 : }
349 : }
350 : // new_hole_bcid is used to ensure the boundary ids used for stitching have no conflicts
351 : // If we shift the default hole boundary ids and output boundary id by free_boundary_id, we need
352 : // to have new_hole_bcid that is larger than any existing boundary ids in the hole meshes and the
353 : // Netgen generated mesh
354 189 : boundary_id_type new_hole_bcid = end_bcid + free_boundary_id;
355 :
356 378 : auto hole_boundaries = isParamValid("hole_boundaries")
357 189 : ? getParam<std::vector<BoundaryName>>("hole_boundaries")
358 249 : : std::vector<BoundaryName>();
359 : const BoundaryName output_boundary =
360 438 : isParamValid("output_boundary") ? getParam<BoundaryName>("output_boundary") : BoundaryName();
361 :
362 : // if outer boundary id is not specified by a numeric output_boundary, we just assign
363 : // free_boundary_id to the output boundary
364 : const std::vector<BoundaryID> output_boundary_id =
365 378 : isParamValid("output_boundary")
366 30 : ? (MooseUtils::isDigits(output_boundary)
367 30 : ? std::vector<BoundaryID>(
368 0 : 1, MooseMeshUtils::getIDFromName<BoundaryName, BoundaryID>(output_boundary))
369 : : std::vector<BoundaryID>(1, free_boundary_id))
370 249 : : std::vector<BoundaryID>();
371 : // Similarly, set the hole boundary ids one by one
372 189 : std::vector<BoundaryID> hole_boundary_ids;
373 567 : if (isParamValid("hole_boundaries"))
374 90 : for (auto h : index_range(hole_boundaries))
375 60 : hole_boundary_ids.push_back(
376 60 : MooseUtils::isDigits(hole_boundaries[h])
377 20 : ? MooseMeshUtils::getIDFromName<BoundaryName, BoundaryID>(hole_boundaries[h])
378 40 : : h + 1 + free_boundary_id);
379 :
380 : // if large ids are used in the hole boundaries or the output boundary,
381 : // we need to make sure the new_hole_bcid is larger than them
382 189 : if (hole_boundary_ids.size())
383 90 : for (auto h : index_range(_hole_ptrs))
384 60 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(hole_boundary_ids[h] + 1));
385 189 : if (output_boundary_id.size())
386 30 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(output_boundary_id[0] + 1));
387 :
388 189 : bool doing_stitching = false;
389 :
390 435 : for (auto hole_i : index_range(_hole_ptrs))
391 : {
392 246 : const MeshBase & hole_mesh = **_hole_ptrs[hole_i];
393 246 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
394 246 : const std::set<boundary_id_type> & local_hole_bcids = hole_boundary_info.get_boundary_ids();
395 :
396 246 : if (!local_hole_bcids.empty())
397 226 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(*local_hole_bcids.rbegin() + 1));
398 246 : hole_mesh.comm().max(new_hole_bcid);
399 :
400 246 : if (_stitch_holes.size() && _stitch_holes[hole_i])
401 136 : doing_stitching = true;
402 : }
403 : // Now we can define boundaries to assist with stitching
404 189 : const boundary_id_type inner_bcid = new_hole_bcid + 1;
405 :
406 : // libMesh mesh stitching still requires a serialized mesh, and it's
407 : // cheaper to do that once than to do it once-per-hole
408 189 : libMesh::MeshSerializer serial(*mesh, doing_stitching);
409 :
410 : // We'll be looking for any sides that match between hole meshes and
411 : // the newly triangulated mesh, to apply bcids accordingly. We
412 : // can't key on Elem::key() here, because that depends on node ids
413 : // that differ from mesh to mesh. We can't use centroids here,
414 : // because that depends on rounding error in order of operations
415 : // that can differ from mesh to mesh. The node locations themselves
416 : // should always match up exactly, though, so let's use (sorted!)
417 : // tuples of those to map from nodes to elements and side numbers.
418 : // Also added a Boolean to check if the face is already used
419 : std::unordered_map<std::tuple<Point, Point, Point>,
420 : std::pair<std::pair<Elem *, unsigned int>, bool>>
421 189 : mesh_faces;
422 :
423 15960 : auto sorted_point_tuple = [](Elem & elem, unsigned int side)
424 : {
425 15960 : std::vector<unsigned int> nodes_on_side = elem.nodes_on_side(side);
426 : libmesh_assert_equal_to(nodes_on_side.size(), 3);
427 15960 : std::vector<Point> p(3);
428 63840 : for (auto i : index_range(p))
429 47880 : p[i] = elem.point(nodes_on_side[i]);
430 15960 : if (p[0] < p[1])
431 : {
432 8024 : if (p[1] < p[2])
433 2312 : return std::make_tuple(p[0], p[1], p[2]);
434 5712 : else if (p[0] < p[2])
435 3095 : return std::make_tuple(p[0], p[2], p[1]);
436 : else
437 2617 : return std::make_tuple(p[2], p[0], p[1]);
438 : }
439 : else
440 : {
441 7936 : if (p[0] < p[2])
442 2524 : return std::make_tuple(p[1], p[0], p[2]);
443 5412 : else if (p[1] < p[2])
444 3051 : return std::make_tuple(p[1], p[2], p[0]);
445 : else
446 2361 : return std::make_tuple(p[2], p[1], p[0]);
447 : }
448 15960 : };
449 :
450 189 : if (!_hole_ptrs.empty())
451 34336 : for (auto elem : mesh->element_ptr_range())
452 170915 : for (auto s : make_range(elem->n_sides()))
453 136732 : if (!elem->neighbor_ptr(s))
454 : {
455 9816 : auto points = sorted_point_tuple(*elem, s);
456 : libmesh_assert(!mesh_faces.count(points));
457 9816 : mesh_faces.emplace(points, std::make_pair(std::make_pair(elem, s), false));
458 153 : }
459 :
460 189 : auto & mesh_boundary_info = mesh->get_boundary_info();
461 :
462 : // Define a reference map variable for subdomain map
463 189 : auto & main_subdomain_map = mesh->set_subdomain_name_map();
464 435 : for (auto hole_i : index_range(_hole_ptrs))
465 : {
466 246 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
467 246 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
468 :
469 : // Our algorithm here requires a serialized Mesh. To avoid
470 : // redundant serialization and deserialization (libMesh
471 : // MeshedHole and stitch_meshes still also require
472 : // serialization) we'll do the serialization up front.
473 246 : libMesh::MeshSerializer serial_hole(hole_mesh);
474 :
475 : // We'll look for any sides that match between the hole mesh and
476 : // the newly triangulated mesh, and apply bcids accordingly.
477 11512 : for (auto elem : hole_mesh.element_ptr_range())
478 57370 : for (auto s : make_range(elem->n_sides()))
479 46104 : if (!elem->neighbor_ptr(s))
480 : {
481 6144 : auto points = sorted_point_tuple(*elem, s);
482 6144 : auto it = mesh_faces.find(points);
483 :
484 : // I'd love to assert that we don't have any missing
485 : // matches, but our holes might themselves have holes
486 6144 : if (it != mesh_faces.end())
487 : {
488 6144 : auto [main_elem, main_side] = it->second.first;
489 6144 : if (_stitch_holes.size() && _stitch_holes[hole_i])
490 : {
491 3504 : hole_boundary_info.add_side(elem, s, new_hole_bcid);
492 3504 : mesh_boundary_info.add_side(main_elem, main_side, inner_bcid);
493 : }
494 : // We would like to take this opportunity to identify the hole boundary
495 : // We set the boundary id to hole_i + 1 at this stage (the default)
496 6144 : mesh_boundary_info.add_side(main_elem, main_side, hole_i + 1);
497 6144 : it->second.second = true;
498 : }
499 246 : }
500 246 : }
501 :
502 : // Any sides that do not match the hole surfaces belong to the external boundary
503 : // We set the boundary id to 0 at this stage (the default)
504 10005 : for (auto & [points, elem_side] : mesh_faces)
505 9816 : if (!elem_side.second)
506 : {
507 3672 : auto [main_elem, main_side] = elem_side.first;
508 3672 : mesh_boundary_info.add_side(main_elem, main_side, 0);
509 : }
510 :
511 : // Now hole boundary ids are 1,2,..., N
512 : // Now external boundary id is 0
513 : // We will find the upper bound of the boundary ids of the mesh so that we would not overwrite
514 : // things during renumbering (because of the bcids used for stitching, this is different from
515 : // end_bcid)
516 189 : const auto temp_bcid_shift = MooseMeshUtils::getNextFreeBoundaryID(*mesh);
517 : // If we stitch holes, we want to make sure the netgen mesh has no conflicting boundary ids
518 : // with the hole meshes.
519 : // OR, if we assign custom boundary ids, we want to make sure no overwriting occurs.
520 : // So we shift all these boundary ids by temp_bcid_shift + free_boundary_id
521 : // before we assign the new boundary ids
522 189 : if (_stitch_holes.size() || output_boundary_id.size())
523 133 : libMesh::MeshTools::Modification::change_boundary_id(
524 133 : *mesh, 0, temp_bcid_shift + free_boundary_id);
525 189 : if (_stitch_holes.size() || hole_boundary_ids.size())
526 339 : for (auto hole_i : index_range(_hole_ptrs))
527 206 : libMesh::MeshTools::Modification::change_boundary_id(
528 206 : *mesh, hole_i + 1, hole_i + 1 + temp_bcid_shift + free_boundary_id);
529 :
530 : // Now we can reassign the boundary ids
531 : // If these ids are specified, we will use them
532 : // Otherwise, we will use the default ones
533 189 : if (output_boundary_id.size())
534 30 : libMesh::MeshTools::Modification::change_boundary_id(
535 30 : *mesh, temp_bcid_shift + free_boundary_id, output_boundary_id[0]);
536 : else
537 159 : libMesh::MeshTools::Modification::change_boundary_id(
538 159 : *mesh, temp_bcid_shift + free_boundary_id, free_boundary_id);
539 :
540 189 : if (hole_boundary_ids.size())
541 90 : for (auto hole_i : index_range(_hole_ptrs))
542 60 : libMesh::MeshTools::Modification::change_boundary_id(
543 60 : *mesh, hole_i + 1 + temp_bcid_shift + free_boundary_id, hole_boundary_ids[hole_i]);
544 : else
545 345 : for (auto hole_i : index_range(_hole_ptrs))
546 186 : libMesh::MeshTools::Modification::change_boundary_id(
547 186 : *mesh, hole_i + 1 + temp_bcid_shift + free_boundary_id, hole_i + 1 + free_boundary_id);
548 :
549 435 : for (auto hole_i : index_range(_hole_ptrs))
550 : {
551 246 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
552 :
553 246 : if (_stitch_holes.size() && _stitch_holes[hole_i])
554 : {
555 : // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
556 : // subdomain map
557 136 : const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
558 136 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
559 :
560 136 : if (_combined_stitching)
561 : {
562 : // The main mesh has been serialized early and will last through the end of this MG
563 : // If the hole mesh to be stitched is not serialized, it will cause parallelization
564 : // issues after combining when the number of processors is large
565 20 : libMesh::MeshSerializer serial_hole(hole_mesh);
566 20 : MooseMeshUtils::copyIntoMesh(*this, *mesh, hole_mesh, false, false, _communicator);
567 20 : }
568 : else
569 : {
570 116 : std::size_t n_nodes_stitched = mesh->stitch_meshes(hole_mesh,
571 : inner_bcid,
572 : new_hole_bcid,
573 : TOLERANCE,
574 : /*clear_stitched_bcids*/ true,
575 116 : _verbose_stitching,
576 : use_binary_search);
577 :
578 116 : if (!n_nodes_stitched)
579 0 : mooseError("Failed to stitch hole mesh ", hole_i, " to new tetrahedralization.");
580 : }
581 : }
582 : }
583 189 : if (doing_stitching && _combined_stitching)
584 : {
585 10 : std::size_t n_nodes_stitched = mesh->stitch_surfaces(inner_bcid,
586 : new_hole_bcid,
587 : TOLERANCE,
588 : /*clear_stitched_bcids*/ true,
589 10 : _verbose_stitching,
590 : use_binary_search);
591 10 : if (!n_nodes_stitched)
592 0 : mooseError("Failed to stitch combined hole meshes to new tetrahedralization.");
593 : }
594 :
595 : // Add user-specified sideset names
596 189 : if (hole_boundary_ids.size())
597 90 : for (auto h : index_range(_hole_ptrs))
598 60 : mesh->get_boundary_info().sideset_name(hole_boundary_ids[h]) = hole_boundaries[h];
599 189 : if (output_boundary_id.size())
600 30 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = output_boundary;
601 :
602 : // Check if one SubdomainName is shared by more than one subdomain ids
603 189 : std::set<SubdomainName> main_subdomain_map_name_list;
604 294 : for (auto const & id_name_pair : main_subdomain_map)
605 105 : main_subdomain_map_name_list.emplace(id_name_pair.second);
606 189 : if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
607 6 : paramError("holes", "The hole meshes contain subdomain name maps with conflicts.");
608 :
609 : // We're done with the hole meshes now, and MeshGenerator doesn't
610 : // want them anymore either.
611 426 : for (auto & hole_ptr : _hole_ptrs)
612 240 : hole_ptr->reset();
613 :
614 186 : mesh->unset_is_prepared();
615 372 : return mesh;
616 : #else
617 : mooseError("Cannot use XYZDelaunayGenerator without NetGen-enabled libMesh.");
618 : return std::unique_ptr<MeshBase>();
619 : #endif
620 196 : }
|
