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