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 17640 : std::size_t operator()(const std::tuple<libMesh::Point, libMesh::Point, libMesh::Point> & p) const
33 : {
34 17640 : std::size_t seed = 0;
35 17640 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<0>(p)));
36 17640 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<1>(p)));
37 17640 : libMesh::boostcopy::hash_combine(seed, std::hash<libMesh::Point>()(std::get<2>(p)));
38 17640 : return seed;
39 : }
40 : };
41 :
42 : } // namespace std
43 :
44 : InputParameters
45 14761 : XYZDelaunayGenerator::validParams()
46 : {
47 14761 : InputParameters params = MeshGenerator::validParams();
48 :
49 59044 : MooseEnum algorithm("BINARY EXHAUSTIVE", "BINARY");
50 :
51 59044 : 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 59044 : params.addParam<SubdomainName>("output_subdomain_name",
58 : "Subdomain name to set on new triangles.");
59 :
60 59044 : 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 59044 : 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 44283 : params.addParam<bool>("smooth_triangulation",
70 29522 : false,
71 : "Whether to do Laplacian mesh smoothing on the generated triangles.");
72 59044 : 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 44283 : params.addParam<std::vector<bool>>(
80 29522 : "stitch_holes", std::vector<bool>(), "Whether to stitch to the mesh defining each hole.");
81 44283 : params.addParam<bool>("convert_holes_for_stitching",
82 29522 : false,
83 : "Whether to convert the 3D hole meshes with non-TRI3 surface sides into "
84 : "a compatible form.");
85 :
86 59044 : MooseEnum conversion_method("ALL SURFACE", "ALL");
87 59044 : 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 44283 : params.addParam<bool>(
95 : "combined_stitching",
96 29522 : false,
97 : "Whether to stitch all holes in one combined stitching step. This is efficient if a great "
98 : "number of holes are to be stitched. But it is hard to debug if problems arise.");
99 :
100 88566 : params.addRangeCheckedParam<Real>(
101 : "desired_volume",
102 : 0,
103 : "desired_volume>=0",
104 : "Desired (maximum) tetrahedral volume, or 0 to skip uniform refinement");
105 :
106 59044 : params.addParam<MooseEnum>(
107 : "algorithm",
108 : algorithm,
109 : "Control the use of binary search for the nodes of the stitched surfaces.");
110 44283 : params.addParam<bool>(
111 29522 : "verbose_stitching", false, "Whether mesh hole stitching should have verbose output.");
112 :
113 14761 : params.addClassDescription(
114 : "Creates tetrahedral 3D meshes within boundaries defined by input meshes.");
115 :
116 29522 : return params;
117 14761 : }
118 :
119 251 : XYZDelaunayGenerator::XYZDelaunayGenerator(const InputParameters & parameters)
120 : : MeshGenerator(parameters),
121 251 : _bdy_ptr(getMesh("boundary")),
122 502 : _desired_volume(getParam<Real>("desired_volume")),
123 251 : _output_subdomain_id(0),
124 502 : _smooth_tri(getParam<bool>("smooth_triangulation")),
125 502 : _hole_ptrs(getMeshes("holes")),
126 502 : _stitch_holes(getParam<std::vector<bool>>("stitch_holes")),
127 502 : _convert_holes_for_stitching(getParam<bool>("convert_holes_for_stitching")),
128 251 : _conversion_method(parameters.get<MooseEnum>("conversion_method")),
129 251 : _combined_stitching(parameters.get<bool>("combined_stitching")),
130 251 : _algorithm(parameters.get<MooseEnum>("algorithm")),
131 502 : _verbose_stitching(parameters.get<bool>("verbose_stitching"))
132 : {
133 251 : if (!_stitch_holes.empty() && _stitch_holes.size() != _hole_ptrs.size())
134 8 : paramError("stitch_holes", "Need one stitch_holes entry per hole, if specified.");
135 :
136 741 : if (isParamValid("hole_boundaries"))
137 : {
138 74 : auto & hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
139 37 : if (hole_boundaries.size() != _hole_ptrs.size())
140 8 : paramError("hole_boundaries", "Need one hole_boundaries entry per hole, if specified.");
141 : }
142 :
143 729 : 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 243 : }
148 :
149 : std::unique_ptr<MeshBase>
150 225 : XYZDelaunayGenerator::generate()
151 : {
152 : #ifdef LIBMESH_HAVE_NETGEN
153 : // Put the boundary mesh in a local pointer
154 : std::unique_ptr<UnstructuredMesh> mesh =
155 225 : 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 225 : mesh->get_boundary_info().clear_boundary_node_ids();
161 :
162 : // Get ready to triangulate its boundary
163 225 : libMesh::NetGenMeshInterface ngint(*mesh);
164 :
165 225 : ngint.smooth_after_generating() = _smooth_tri;
166 :
167 225 : 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 497 : for (const auto hole_i : index_range(_hole_ptrs))
174 : {
175 288 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
176 288 : 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 288 : std::set<ElemType> hole_elem_types;
180 288 : std::set<unsigned short> hole_elem_dims;
181 288 : std::vector<std::pair<dof_id_type, unsigned int>> hole_elem_external_sides;
182 9088 : for (auto elem : hole_mesh.element_ptr_range())
183 : {
184 8800 : 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 8800 : if (elem->dim() == 3)
189 44184 : 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 35604 : if (!elem->neighbor_ptr(s))
193 : {
194 5380 : hole_elem_types.emplace(elem->side_ptr(s)->type());
195 5380 : 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 220 : hole_elem_types.emplace(elem->type());
201 288 : }
202 288 : if (hole_elem_dims.size() != 1 || *hole_elem_dims.begin() < 2)
203 16 : paramError(
204 : "holes",
205 : "All elements in a hole mesh must have the same dimension that is either 2D or 3D.");
206 280 : 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 254 : if (*hole_elem_types.begin() != ElemType::TRI3 || hole_elem_types.size() > 1)
213 : {
214 70 : if (_stitch_holes.size() && _stitch_holes[hole_i] && !_convert_holes_for_stitching)
215 8 : 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 66 : 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 11 : BoundaryID temp_ext_bid = MooseMeshUtils::getNextFreeBoundaryID(hole_mesh);
222 605 : for (const auto & hees : hole_elem_external_sides)
223 594 : hole_mesh.get_boundary_info().add_side(hees.first, hees.second, temp_ext_bid);
224 11 : MooseMeshElementConversionUtils::transitionLayerGenerator(
225 22 : hole_mesh, std::vector<BoundaryName>({std::to_string(temp_ext_bid)}), 1, false);
226 11 : 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 11 : if (!hole_mesh.is_replicated())
230 2 : hole_mesh.prepare_for_use();
231 : else
232 9 : hole_mesh.find_neighbors();
233 : }
234 : else
235 55 : MeshTools::Modification::all_tri(**_hole_ptrs[hole_i]);
236 : }
237 : }
238 : else // if (*hole_elem_dims.begin() == 2)
239 : {
240 : // There is no point to stitch a 2D hole mesh, so we throw an error if this is attempted
241 26 : if (_stitch_holes.size() && _stitch_holes[hole_i])
242 4 : paramError("holes",
243 4 : "the hole mesh with index " + std::to_string(hole_i) +
244 : " is a 2D mesh, for which stitching onto a 3D mesh does not make sense.");
245 : }
246 272 : }
247 :
248 : std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh>>> ngholes =
249 209 : std::make_unique<std::vector<std::unique_ptr<UnstructuredMesh>>>();
250 :
251 : // The libMesh Netgen interface will modify hole meshes in-place, so
252 : // we make copies to pass in.
253 481 : for (std::unique_ptr<MeshBase> * hole_ptr : _hole_ptrs)
254 : {
255 : // How did we never add a ReplicatedMesh(MeshBase&) constructor in
256 : // libMesh?
257 272 : const UnstructuredMesh & hole = dynamic_cast<UnstructuredMesh &>(**hole_ptr);
258 272 : ngholes->push_back(std::make_unique<ReplicatedMesh>(hole));
259 : }
260 :
261 209 : if (!_hole_ptrs.empty())
262 169 : ngint.attach_hole_list(std::move(ngholes));
263 :
264 209 : ngint.triangulate();
265 :
266 627 : if (isParamValid("output_subdomain_name"))
267 : {
268 114 : auto output_subdomain_name = getParam<SubdomainName>("output_subdomain_name");
269 57 : _output_subdomain_id = MooseMeshUtils::getSubdomainID(output_subdomain_name, *mesh);
270 :
271 57 : if (_output_subdomain_id == Elem::invalid_subdomain_id)
272 : {
273 : // We'll probably need to make a new ID, then
274 57 : _output_subdomain_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
275 :
276 : // But check the hole meshes for our output subdomain name too
277 153 : for (auto & hole_ptr : _hole_ptrs)
278 : {
279 96 : auto possible_sbdid = MooseMeshUtils::getSubdomainID(output_subdomain_name, **hole_ptr);
280 : // Huh, it was in one of them
281 96 : if (possible_sbdid != Elem::invalid_subdomain_id)
282 : {
283 0 : _output_subdomain_id = possible_sbdid;
284 0 : break;
285 : }
286 96 : _output_subdomain_id =
287 96 : std::max(_output_subdomain_id, MooseMeshUtils::getNextFreeSubdomainID(**hole_ptr));
288 : }
289 :
290 57 : mesh->subdomain_name(_output_subdomain_id) = output_subdomain_name;
291 : }
292 57 : }
293 :
294 209 : if (_smooth_tri || _output_subdomain_id)
295 30721 : for (auto elem : mesh->element_ptr_range())
296 : {
297 30631 : elem->subdomain_id() = _output_subdomain_id;
298 :
299 : // I do not trust Laplacian mesh smoothing not to invert
300 : // elements near reentrant corners. Eventually we'll add better
301 : // smoothing options, but even those might have failure cases.
302 : // Better to always do extra tests here than to ever let users
303 : // try to run on a degenerate mesh.
304 30631 : if (elem->is_flipped())
305 : {
306 0 : if (_smooth_tri)
307 0 : mooseError("Inverted element found in triangulation.\n"
308 : "Laplacian smoothing can create these at reentrant corners; disable it?");
309 : else
310 0 : mooseError("Unexplained inverted element found in triangulation.\n");
311 : }
312 90 : }
313 :
314 209 : const bool use_binary_search = (_algorithm == "BINARY");
315 :
316 : // The hole meshes are specified by the user, so they could have any
317 : // BCID or no BCID or any combination of BCIDs on their outer
318 : // boundary, so we'll have to set our own BCID to use for stitching
319 : // there. We'll need to check all the holes for used BCIDs, if we
320 : // want to pick a new ID on hole N that doesn't conflict with any
321 : // IDs on hole M < N (or with the IDs on the new triangulation)
322 :
323 : // The Netgen generated mesh does not have hole and output boundary ids,
324 : // which is different from its 2D counterpart.
325 : // So, we need to assign boundary ids to the hole boundaries and output boundary
326 : // By default, we assign the hole boundary ids to be 1,2,..., N
327 : // and the output boundary id to be 0
328 : // We will need a different set of boundary ids if we are stitching holes
329 :
330 : // end_bcid is one more than the maximum of default output boundary id and hole boundary ids
331 209 : const boundary_id_type end_bcid = _hole_ptrs.size() + 1;
332 :
333 : // For the hole meshes that need to be stitched, we would like to make sure the hole boundary ids
334 : // and output boundary id are not conflicting with the existing boundary ids of the hole meshes to
335 : // be stitched.
336 209 : BoundaryID free_boundary_id = 0;
337 : // We need to get a free boundary id that is larger than the existing boundary ids in any hole
338 : // meshes
339 209 : if (_stitch_holes.size())
340 : {
341 309 : for (auto hole_i : index_range(_hole_ptrs))
342 : {
343 184 : if (_stitch_holes[hole_i])
344 : {
345 151 : free_boundary_id =
346 151 : std::max(free_boundary_id, MooseMeshUtils::getNextFreeBoundaryID(**_hole_ptrs[hole_i]));
347 151 : (*_hole_ptrs[hole_i])->comm().max(free_boundary_id);
348 : }
349 : }
350 : }
351 : // new_hole_bcid is used to ensure the boundary ids used for stitching have no conflicts
352 : // If we shift the default hole boundary ids and output boundary id by free_boundary_id, we need
353 : // to have new_hole_bcid that is larger than any existing boundary ids in the hole meshes and the
354 : // Netgen generated mesh
355 209 : boundary_id_type new_hole_bcid = end_bcid + free_boundary_id;
356 :
357 418 : auto hole_boundaries = isParamValid("hole_boundaries")
358 209 : ? getParam<std::vector<BoundaryName>>("hole_boundaries")
359 275 : : std::vector<BoundaryName>();
360 : const BoundaryName output_boundary =
361 484 : isParamValid("output_boundary") ? getParam<BoundaryName>("output_boundary") : BoundaryName();
362 :
363 : // if outer boundary id is not specified by a numeric output_boundary, we just assign
364 : // free_boundary_id to the output boundary
365 : const std::vector<BoundaryID> output_boundary_id =
366 418 : isParamValid("output_boundary")
367 33 : ? (MooseUtils::isDigits(output_boundary)
368 33 : ? std::vector<BoundaryID>(
369 0 : 1, MooseMeshUtils::getIDFromName<BoundaryName, BoundaryID>(output_boundary))
370 : : std::vector<BoundaryID>(1, free_boundary_id))
371 275 : : std::vector<BoundaryID>();
372 : // Similarly, set the hole boundary ids one by one
373 209 : std::vector<BoundaryID> hole_boundary_ids;
374 627 : if (isParamValid("hole_boundaries"))
375 99 : for (auto h : index_range(hole_boundaries))
376 66 : hole_boundary_ids.push_back(
377 66 : MooseUtils::isDigits(hole_boundaries[h])
378 22 : ? MooseMeshUtils::getIDFromName<BoundaryName, BoundaryID>(hole_boundaries[h])
379 44 : : h + 1 + free_boundary_id);
380 :
381 : // if large ids are used in the hole boundaries or the output boundary,
382 : // we need to make sure the new_hole_bcid is larger than them
383 209 : if (hole_boundary_ids.size())
384 99 : for (auto h : index_range(_hole_ptrs))
385 66 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(hole_boundary_ids[h] + 1));
386 209 : if (output_boundary_id.size())
387 33 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(output_boundary_id[0] + 1));
388 :
389 209 : bool doing_stitching = false;
390 :
391 481 : for (auto hole_i : index_range(_hole_ptrs))
392 : {
393 272 : const MeshBase & hole_mesh = **_hole_ptrs[hole_i];
394 272 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
395 272 : const std::set<boundary_id_type> & local_hole_bcids = hole_boundary_info.get_boundary_ids();
396 :
397 272 : if (!local_hole_bcids.empty())
398 250 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(*local_hole_bcids.rbegin() + 1));
399 272 : hole_mesh.comm().max(new_hole_bcid);
400 :
401 272 : if (_stitch_holes.size() && _stitch_holes[hole_i])
402 151 : doing_stitching = true;
403 : }
404 : // Now we can define boundaries to assist with stitching
405 209 : const boundary_id_type inner_bcid = new_hole_bcid + 1;
406 :
407 : // libMesh mesh stitching still requires a serialized mesh, and it's
408 : // cheaper to do that once than to do it once-per-hole
409 209 : libMesh::MeshSerializer serial(*mesh, doing_stitching);
410 :
411 : // We'll be looking for any sides that match between hole meshes and
412 : // the newly triangulated mesh, to apply bcids accordingly. We
413 : // can't key on Elem::key() here, because that depends on node ids
414 : // that differ from mesh to mesh. We can't use centroids here,
415 : // because that depends on rounding error in order of operations
416 : // that can differ from mesh to mesh. The node locations themselves
417 : // should always match up exactly, though, so let's use (sorted!)
418 : // tuples of those to map from nodes to elements and side numbers.
419 : // Also added a Boolean to check if the face is already used
420 : std::unordered_map<std::tuple<Point, Point, Point>,
421 : std::pair<std::pair<Elem *, unsigned int>, bool>>
422 209 : mesh_faces;
423 :
424 17640 : auto sorted_point_tuple = [](Elem & elem, unsigned int side)
425 : {
426 17640 : std::vector<unsigned int> nodes_on_side = elem.nodes_on_side(side);
427 : libmesh_assert_equal_to(nodes_on_side.size(), 3);
428 17640 : std::vector<Point> p(3);
429 70560 : for (auto i : index_range(p))
430 52920 : p[i] = elem.point(nodes_on_side[i]);
431 17640 : if (p[0] < p[1])
432 : {
433 8873 : if (p[1] < p[2])
434 2556 : return std::make_tuple(p[0], p[1], p[2]);
435 6317 : else if (p[0] < p[2])
436 3426 : return std::make_tuple(p[0], p[2], p[1]);
437 : else
438 2891 : return std::make_tuple(p[2], p[0], p[1]);
439 : }
440 : else
441 : {
442 8767 : if (p[0] < p[2])
443 2789 : return std::make_tuple(p[1], p[0], p[2]);
444 5978 : else if (p[1] < p[2])
445 3373 : return std::make_tuple(p[1], p[2], p[0]);
446 : else
447 2605 : return std::make_tuple(p[2], p[1], p[0]);
448 : }
449 17640 : };
450 :
451 209 : if (!_hole_ptrs.empty())
452 37887 : for (auto elem : mesh->element_ptr_range())
453 188590 : for (auto s : make_range(elem->n_sides()))
454 150872 : if (!elem->neighbor_ptr(s))
455 : {
456 10848 : auto points = sorted_point_tuple(*elem, s);
457 : libmesh_assert(!mesh_faces.count(points));
458 10848 : mesh_faces.emplace(points, std::make_pair(std::make_pair(elem, s), false));
459 169 : }
460 :
461 209 : auto & mesh_boundary_info = mesh->get_boundary_info();
462 :
463 : // Define a reference map variable for subdomain map
464 209 : auto & main_subdomain_map = mesh->set_subdomain_name_map();
465 481 : for (auto hole_i : index_range(_hole_ptrs))
466 : {
467 272 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
468 272 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
469 :
470 : // Our algorithm here requires a serialized Mesh. To avoid
471 : // redundant serialization and deserialization (libMesh
472 : // MeshedHole and stitch_meshes still also require
473 : // serialization) we'll do the serialization up front.
474 272 : libMesh::MeshSerializer serial_hole(hole_mesh);
475 :
476 : // We'll look for any sides that match between the hole mesh and
477 : // the newly triangulated mesh, and apply bcids accordingly.
478 12705 : for (auto elem : hole_mesh.element_ptr_range())
479 63309 : for (auto s : make_range(elem->n_sides()))
480 50876 : if (!elem->neighbor_ptr(s))
481 : {
482 6792 : auto points = sorted_point_tuple(*elem, s);
483 6792 : auto it = mesh_faces.find(points);
484 :
485 : // I'd love to assert that we don't have any missing
486 : // matches, but our holes might themselves have holes
487 6792 : if (it != mesh_faces.end())
488 : {
489 6792 : auto [main_elem, main_side] = it->second.first;
490 6792 : if (_stitch_holes.size() && _stitch_holes[hole_i])
491 : {
492 3888 : hole_boundary_info.add_side(elem, s, new_hole_bcid);
493 3888 : mesh_boundary_info.add_side(main_elem, main_side, inner_bcid);
494 : }
495 : // We would like to take this opportunity to identify the hole boundary
496 : // We set the boundary id to hole_i + 1 at this stage (the default)
497 6792 : mesh_boundary_info.add_side(main_elem, main_side, hole_i + 1);
498 6792 : it->second.second = true;
499 : }
500 272 : }
501 272 : }
502 :
503 : // Any sides that do not match the hole surfaces belong to the external boundary
504 : // We set the boundary id to 0 at this stage (the default)
505 11057 : for (auto & [points, elem_side] : mesh_faces)
506 10848 : if (!elem_side.second)
507 : {
508 4056 : auto [main_elem, main_side] = elem_side.first;
509 4056 : mesh_boundary_info.add_side(main_elem, main_side, 0);
510 : }
511 :
512 : // Now hole boundary ids are 1,2,..., N
513 : // Now external boundary id is 0
514 : // We will find the upper bound of the boundary ids of the mesh so that we would not overwrite
515 : // things during renumbering (because of the bcids used for stitching, this is different from
516 : // end_bcid)
517 209 : const auto temp_bcid_shift = MooseMeshUtils::getNextFreeBoundaryID(*mesh);
518 : // If we stitch holes, we want to make sure the netgen mesh has no conflicting boundary ids
519 : // with the hole meshes.
520 : // OR, if we assign custom boundary ids, we want to make sure no overwriting occurs.
521 : // So we shift all these boundary ids by temp_bcid_shift + free_boundary_id
522 : // before we assign the new boundary ids
523 209 : if (_stitch_holes.size() || output_boundary_id.size())
524 147 : libMesh::MeshTools::Modification::change_boundary_id(
525 147 : *mesh, 0, temp_bcid_shift + free_boundary_id);
526 209 : if (_stitch_holes.size() || hole_boundary_ids.size())
527 375 : for (auto hole_i : index_range(_hole_ptrs))
528 228 : libMesh::MeshTools::Modification::change_boundary_id(
529 228 : *mesh, hole_i + 1, hole_i + 1 + temp_bcid_shift + free_boundary_id);
530 :
531 : // Now we can reassign the boundary ids
532 : // If these ids are specified, we will use them
533 : // Otherwise, we will use the default ones
534 209 : if (output_boundary_id.size())
535 33 : libMesh::MeshTools::Modification::change_boundary_id(
536 33 : *mesh, temp_bcid_shift + free_boundary_id, output_boundary_id[0]);
537 : else
538 176 : libMesh::MeshTools::Modification::change_boundary_id(
539 176 : *mesh, temp_bcid_shift + free_boundary_id, free_boundary_id);
540 :
541 209 : if (hole_boundary_ids.size())
542 99 : for (auto hole_i : index_range(_hole_ptrs))
543 66 : libMesh::MeshTools::Modification::change_boundary_id(
544 66 : *mesh, hole_i + 1 + temp_bcid_shift + free_boundary_id, hole_boundary_ids[hole_i]);
545 : else
546 382 : for (auto hole_i : index_range(_hole_ptrs))
547 206 : libMesh::MeshTools::Modification::change_boundary_id(
548 206 : *mesh, hole_i + 1 + temp_bcid_shift + free_boundary_id, hole_i + 1 + free_boundary_id);
549 :
550 481 : for (auto hole_i : index_range(_hole_ptrs))
551 : {
552 272 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(**_hole_ptrs[hole_i]);
553 :
554 272 : if (_stitch_holes.size() && _stitch_holes[hole_i])
555 : {
556 : // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
557 : // subdomain map
558 151 : const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
559 151 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
560 :
561 151 : if (_combined_stitching)
562 22 : MooseMeshUtils::copyIntoMesh(*this, *mesh, hole_mesh, false, false, _communicator);
563 : else
564 : {
565 129 : std::size_t n_nodes_stitched = mesh->stitch_meshes(hole_mesh,
566 : inner_bcid,
567 : new_hole_bcid,
568 : TOLERANCE,
569 : /*clear_stitched_bcids*/ true,
570 129 : _verbose_stitching,
571 : use_binary_search);
572 :
573 129 : if (!n_nodes_stitched)
574 0 : mooseError("Failed to stitch hole mesh ", hole_i, " to new tetrahedralization.");
575 : }
576 : }
577 : }
578 209 : if (doing_stitching && _combined_stitching)
579 : {
580 11 : std::size_t n_nodes_stitched = mesh->stitch_surfaces(inner_bcid,
581 : new_hole_bcid,
582 : TOLERANCE,
583 : /*clear_stitched_bcids*/ true,
584 11 : _verbose_stitching,
585 : use_binary_search);
586 11 : if (!n_nodes_stitched)
587 0 : mooseError("Failed to stitch combined hole meshes to new tetrahedralization.");
588 : }
589 :
590 : // Add user-specified sideset names
591 209 : if (hole_boundary_ids.size())
592 99 : for (auto h : index_range(_hole_ptrs))
593 66 : mesh->get_boundary_info().sideset_name(hole_boundary_ids[h]) = hole_boundaries[h];
594 209 : if (output_boundary_id.size())
595 33 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = output_boundary;
596 :
597 : // Check if one SubdomainName is shared by more than one subdomain ids
598 209 : std::set<SubdomainName> main_subdomain_map_name_list;
599 318 : for (auto const & id_name_pair : main_subdomain_map)
600 109 : main_subdomain_map_name_list.emplace(id_name_pair.second);
601 209 : if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
602 8 : paramError("holes", "The hole meshes contain subdomain name maps with conflicts.");
603 :
604 : // We're done with the hole meshes now, and MeshGenerator doesn't
605 : // want them anymore either.
606 469 : for (auto & hole_ptr : _hole_ptrs)
607 264 : hole_ptr->reset();
608 :
609 205 : mesh->set_isnt_prepared();
610 410 : return mesh;
611 : #else
612 : mooseError("Cannot use XYZDelaunayGenerator without NetGen-enabled libMesh.");
613 : return std::unique_ptr<MeshBase>();
614 : #endif
615 216 : }
|
