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 "XYDelaunayGenerator.h"
11 :
12 : #include "CastUniquePointer.h"
13 : #include "MooseMeshUtils.h"
14 : #include "MeshTriangulationUtils.h"
15 : #include "BoundaryLayerUtils.h"
16 : #include "MooseUtils.h"
17 :
18 : #include "libmesh/int_range.h"
19 : #include "libmesh/mesh_modification.h"
20 : #include "libmesh/mesh_serializer.h"
21 : #include "libmesh/unstructured_mesh.h"
22 : #include "DelimitedFileReader.h"
23 :
24 : registerMooseObject("MooseApp", XYDelaunayGenerator);
25 :
26 : InputParameters
27 3798 : XYDelaunayGenerator::validParams()
28 : {
29 3798 : InputParameters params = SurfaceDelaunayGeneratorBase::validParams();
30 :
31 15192 : MooseEnum algorithm("BINARY EXHAUSTIVE", "BINARY");
32 15192 : MooseEnum tri_elem_type("TRI3 TRI6 TRI7 DEFAULT", "DEFAULT");
33 :
34 15192 : params.addRequiredParam<MeshGeneratorName>(
35 : "boundary",
36 : "The input MeshGenerator defining the output outer boundary and required Steiner points.");
37 15192 : params.addParam<std::vector<BoundaryName>>(
38 : "input_boundary_names", "2D-input-mesh boundaries defining the output mesh outer boundary");
39 15192 : params.addParam<std::vector<SubdomainName>>(
40 : "input_subdomain_names", "1D-input-mesh subdomains defining the output mesh outer boundary");
41 11394 : params.addParam<unsigned int>("add_nodes_per_boundary_segment",
42 7596 : 0,
43 : "How many more nodes to add in each outer boundary segment.");
44 11394 : params.addParam<bool>(
45 7596 : "refine_boundary", true, "Whether to allow automatically refining the outer boundary.");
46 :
47 15192 : params.addParam<SubdomainName>("output_subdomain_name",
48 : "Subdomain name to set on new triangles.");
49 15192 : params.addParam<SubdomainID>("output_subdomain_id", "Subdomain id to set on new triangles.");
50 :
51 15192 : params.addParam<BoundaryName>(
52 : "output_boundary",
53 : "Boundary name to set on new outer boundary. Default ID: 0 if no hole meshes are stitched; "
54 : "or maximum boundary ID of all the stitched hole meshes + 1.");
55 15192 : params.addParam<std::vector<BoundaryName>>(
56 : "hole_boundaries",
57 : "Boundary names to set on holes. Default IDs are numbered up from 1 if no hole meshes are "
58 : "stitched; or from maximum boundary ID of all the stitched hole meshes + 2.");
59 :
60 11394 : params.addParam<bool>(
61 : "verify_holes",
62 7596 : true,
63 : "Verify holes do not intersect boundary or each other. Asymptotically costly.");
64 :
65 11394 : params.addParam<bool>("smooth_triangulation",
66 7596 : false,
67 : "Whether to do Laplacian mesh smoothing on the generated triangles.");
68 11394 : params.addParam<std::vector<MeshGeneratorName>>(
69 7596 : "holes", std::vector<MeshGeneratorName>(), "The MeshGenerators that define mesh holes.");
70 11394 : params.addParam<std::vector<bool>>(
71 7596 : "stitch_holes", std::vector<bool>(), "Whether to stitch to the mesh defining each hole.");
72 11394 : params.addParam<std::vector<bool>>("refine_holes",
73 7596 : std::vector<bool>(),
74 : "Whether to allow automatically refining each hole boundary.");
75 22788 : params.addRangeCheckedParam<Real>(
76 : "desired_area",
77 : 0,
78 : "desired_area>=0",
79 : "Desired (maximum) triangle area, or 0 to skip uniform refinement");
80 11394 : params.addParam<std::string>(
81 : "desired_area_func",
82 7596 : std::string(),
83 : "Desired area as a function of x,y; omit to skip non-uniform refinement");
84 :
85 15192 : params.addParam<MooseEnum>(
86 : "algorithm",
87 : algorithm,
88 : "Control the use of binary search for the nodes of the stitched surfaces.");
89 15192 : params.addParam<MooseEnum>(
90 : "tri_element_type", tri_elem_type, "Type of the triangular elements to be generated.");
91 11394 : params.addParam<bool>(
92 7596 : "verbose_stitching", false, "Whether mesh stitching should have verbose output.");
93 15192 : params.addParam<std::vector<Point>>("interior_points",
94 : {},
95 : "Interior node locations, if no smoothing is used. Any point "
96 : "outside the surface will not be meshed.");
97 15192 : params.addParam<std::vector<FileName>>(
98 : "interior_point_files", {}, "Text file(s) with the interior points, one per line");
99 7596 : params.addClassDescription("Triangulates meshes within boundaries defined by input meshes.");
100 :
101 22788 : params.addRangeCheckedParam<Real>("outer_boundary_layer_thickness",
102 : 0,
103 : "outer_boundary_layer_thickness>=0",
104 : "Thickness of the outer boundary layer to be added.");
105 11394 : params.addParam<unsigned int>(
106 7596 : "outer_boundary_layer_num", 0, "Number of layers for the outer boundary layer.");
107 18990 : params.addRangeCheckedParam<Real>(
108 : "outer_boundary_layer_bias",
109 7596 : 1.0,
110 : "outer_boundary_layer_bias>0",
111 : "Bias factor for the thickness of each layer in the outer boundary layer.");
112 :
113 22788 : params.addRangeCheckedParam<std::vector<Real>>(
114 : "holes_boundary_layer_thickness",
115 : "holes_boundary_layer_thickness>=0",
116 : "Thickness of the hole boundary layers to be added.");
117 15192 : params.addParam<std::vector<unsigned int>>("holes_boundary_layer_num",
118 : "Numbers of layers for the hole boundary layers.");
119 22788 : params.addRangeCheckedParam<std::vector<Real>>(
120 : "holes_boundary_layer_bias",
121 : "holes_boundary_layer_bias>0",
122 : "Bias factors for the thickness of each layer in the hole boundary layers.");
123 :
124 15192 : params.addParamNamesToGroup("interior_points interior_point_files",
125 : "Mandatory mesh interior nodes");
126 :
127 11394 : params.addParamNamesToGroup("outer_boundary_layer_thickness outer_boundary_layer_num "
128 : "outer_boundary_layer_bias holes_boundary_layer_thickness "
129 : "holes_boundary_layer_num holes_boundary_layer_bias",
130 : "Boundary layer");
131 :
132 7596 : return params;
133 3798 : }
134 :
135 373 : XYDelaunayGenerator::XYDelaunayGenerator(const InputParameters & parameters)
136 : : SurfaceDelaunayGeneratorBase(parameters),
137 373 : _bdy_ptr(getMesh("boundary")),
138 746 : _add_nodes_per_boundary_segment(getParam<unsigned int>("add_nodes_per_boundary_segment")),
139 746 : _refine_bdy(getParam<bool>("refine_boundary")),
140 373 : _output_subdomain_id(0),
141 746 : _smooth_tri(getParam<bool>("smooth_triangulation")),
142 746 : _verify_holes(getParam<bool>("verify_holes")),
143 746 : _hole_ptrs(getMeshes("holes")),
144 746 : _stitch_holes(getParam<std::vector<bool>>("stitch_holes")),
145 746 : _refine_holes(getParam<std::vector<bool>>("refine_holes")),
146 746 : _desired_area(getParam<Real>("desired_area")),
147 746 : _desired_area_func(getParam<std::string>("desired_area_func")),
148 373 : _algorithm(parameters.get<MooseEnum>("algorithm")),
149 373 : _tri_elem_type(parameters.get<MooseEnum>("tri_element_type")),
150 373 : _verbose_stitching(parameters.get<bool>("verbose_stitching")),
151 746 : _interior_points(getParam<std::vector<Point>>("interior_points")),
152 746 : _outer_boundary_layer_thickness(getParam<Real>("outer_boundary_layer_thickness")),
153 746 : _outer_boundary_layer_num(getParam<unsigned int>("outer_boundary_layer_num")),
154 746 : _outer_boundary_layer_bias(getParam<Real>("outer_boundary_layer_bias")),
155 383 : _holes_boundary_layer_thickness(
156 1119 : isParamValid("holes_boundary_layer_thickness")
157 373 : ? getParam<std::vector<Real>>("holes_boundary_layer_thickness")
158 : : std::vector<Real>()),
159 1139 : _holes_boundary_layer_num(isParamValid("holes_boundary_layer_num")
160 373 : ? getParam<std::vector<unsigned int>>("holes_boundary_layer_num")
161 : : std::vector<unsigned int>()),
162 1139 : _holes_boundary_layer_bias(isParamValid("holes_boundary_layer_bias")
163 373 : ? getParam<std::vector<Real>>("holes_boundary_layer_bias")
164 373 : : std::vector<Real>())
165 : {
166 310 : if ((_desired_area > 0.0 && !_desired_area_func.empty()) ||
167 1053 : (_desired_area > 0.0 && _use_auto_area_func) ||
168 370 : (!_desired_area_func.empty() && _use_auto_area_func))
169 6 : paramError("desired_area_func",
170 : "Only one of the three methods ('desired_area', 'desired_area_func', and "
171 : "'_use_auto_area_func') to set element area limit should be used.");
172 :
173 370 : if (!_use_auto_area_func)
174 1080 : if (isParamSetByUser("auto_area_func_default_size") ||
175 1440 : isParamSetByUser("auto_area_func_default_size_dist") ||
176 2520 : isParamSetByUser("auto_area_function_num_points") ||
177 1440 : isParamSetByUser("auto_area_function_power"))
178 6 : paramError("use_auto_area_func",
179 : "If this parameter is set to false, the following parameters should not be set: "
180 : "'auto_area_func_default_size', 'auto_area_func_default_size_dist', "
181 : "'auto_area_function_num_points', 'auto_area_function_power'.");
182 :
183 367 : if (!_stitch_holes.empty() && _stitch_holes.size() != _hole_ptrs.size())
184 0 : paramError("stitch_holes", "Need one stitch_holes entry per hole, if specified.");
185 :
186 546 : for (auto hole_i : index_range(_stitch_holes))
187 179 : if (_stitch_holes[hole_i] && (hole_i >= _refine_holes.size() || _refine_holes[hole_i]))
188 0 : paramError("refine_holes", "Disable auto refine of any hole boundary to be stitched.");
189 :
190 1101 : if (isParamValid("hole_boundaries"))
191 : {
192 76 : auto & hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
193 38 : if (hole_boundaries.size() != _hole_ptrs.size())
194 0 : paramError("hole_boundaries", "Need one hole_boundaries entry per hole, if specified.");
195 : }
196 : // Copied from MultiApp.C
197 734 : const auto & positions_files = getParam<std::vector<FileName>>("interior_point_files");
198 380 : for (const auto p_file_it : index_range(positions_files))
199 : {
200 13 : const std::string positions_file = positions_files[p_file_it];
201 13 : MooseUtils::DelimitedFileReader file(positions_file, &_communicator);
202 13 : file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
203 13 : file.read();
204 :
205 13 : const std::vector<Point> & data = file.getDataAsPoints();
206 65 : for (const auto & d : data)
207 52 : _interior_points.push_back(d);
208 13 : }
209 : bool has_duplicates =
210 367 : std::any_of(_interior_points.begin(),
211 : _interior_points.end(),
212 113 : [&](const Point & p)
213 113 : { return std::count(_interior_points.begin(), _interior_points.end(), p) > 1; });
214 367 : if (has_duplicates)
215 6 : paramError("interior_points", "Duplicate points were found in the provided interior points.");
216 :
217 364 : if ((_outer_boundary_layer_thickness > 0 && _outer_boundary_layer_num == 0) ||
218 364 : (_outer_boundary_layer_thickness == 0 && _outer_boundary_layer_num > 0))
219 0 : paramError(
220 : "outer_boundary_layer_thickness",
221 : "this parameter must be set as non-zero along with a non-zero outer_boundary_layer_num.");
222 :
223 364 : if (_outer_boundary_layer_num > 0)
224 : {
225 10 : if (_add_nodes_per_boundary_segment)
226 0 : paramError("add_nodes_per_boundary_segment",
227 : "Cannot add nodes per boundary segment when using an outer boundary layer.");
228 10 : if (_refine_bdy)
229 0 : paramError("refine_boundary", "Cannot refine boundary when using an outer boundary layer.");
230 : }
231 :
232 374 : if ((_holes_boundary_layer_thickness.size() && _holes_boundary_layer_num.size() &&
233 10 : _holes_boundary_layer_bias.size()) !=
234 718 : (_holes_boundary_layer_thickness.size() || _holes_boundary_layer_num.size() ||
235 354 : _holes_boundary_layer_bias.size()))
236 0 : paramError("holes_boundary_layer_num",
237 : "holes_boundary_layer_thickness, holes_boundary_layer_bias and this parameter must "
238 : "be specified or not specified together.");
239 374 : if (_holes_boundary_layer_thickness.size() &&
240 10 : _holes_boundary_layer_thickness.size() != _hole_ptrs.size())
241 0 : paramError("holes_boundary_layer_thickness",
242 : "If specified, this parameter must have the same length as 'holes'.");
243 364 : if (_holes_boundary_layer_num.size() && _holes_boundary_layer_num.size() != _hole_ptrs.size())
244 0 : paramError("holes_boundary_layer_num",
245 : "If specified, this parameter must have the same length as 'holes'.");
246 364 : if (_holes_boundary_layer_bias.size() && _holes_boundary_layer_bias.size() != _hole_ptrs.size())
247 0 : paramError("holes_boundary_layer_bias",
248 : "If specified, this parameter must have the same length as 'holes'.");
249 384 : for (const auto & i : index_range(_holes_boundary_layer_thickness))
250 : {
251 40 : if ((_holes_boundary_layer_thickness[i] > 0 && _holes_boundary_layer_num[i] == 0) ||
252 20 : (_holes_boundary_layer_thickness[i] == 0 && _holes_boundary_layer_num[i] > 0))
253 0 : paramError(
254 : "holes_boundary_layer_thickness",
255 0 : "entry " + std::to_string(i) +
256 : " must be set as non-zero along with a non-zero holes_boundary_layer_num entry.");
257 20 : if (_holes_boundary_layer_thickness[i] > 0)
258 : {
259 20 : if ((_refine_holes.size() > i && _refine_holes[i]) || _refine_holes.empty())
260 0 : paramError("refine_holes", "Cannot refine hole boundary when using a hole boundary layer.");
261 : }
262 : }
263 364 : }
264 :
265 : std::unique_ptr<MeshBase>
266 352 : XYDelaunayGenerator::generate()
267 : {
268 352 : MeshTriangulationUtils::XYDelaunayOptions xyd_opts;
269 1056 : if (isParamValid("input_boundary_names"))
270 30 : xyd_opts.input_boundary_names = getParam<std::vector<BoundaryName>>("input_boundary_names");
271 1056 : if (isParamValid("input_subdomain_names"))
272 30 : xyd_opts.input_subdomain_names = getParam<std::vector<SubdomainName>>("input_subdomain_names");
273 352 : xyd_opts.add_nodes_per_boundary_segment = _add_nodes_per_boundary_segment;
274 352 : xyd_opts.refine_bdy = _refine_bdy;
275 352 : xyd_opts.verify_holes = _verify_holes;
276 352 : xyd_opts.smooth_tri = _smooth_tri;
277 352 : xyd_opts.desired_area = _desired_area;
278 352 : xyd_opts.desired_area_func = _desired_area_func;
279 352 : xyd_opts.use_auto_area_func = _use_auto_area_func;
280 352 : xyd_opts.auto_area_func_default_size = _auto_area_func_default_size;
281 352 : xyd_opts.auto_area_func_default_size_dist = _auto_area_func_default_size_dist;
282 352 : xyd_opts.auto_area_function_num_points = _auto_area_function_num_points;
283 352 : xyd_opts.auto_area_function_power = _auto_area_function_power;
284 352 : xyd_opts.interior_points = _interior_points;
285 352 : xyd_opts.tri_elem_type = std::string(_tri_elem_type);
286 352 : xyd_opts.stitch_holes = _stitch_holes;
287 352 : xyd_opts.refine_holes = _refine_holes;
288 352 : xyd_opts.use_binary_search = (_algorithm == "BINARY");
289 352 : xyd_opts.verbose_stitching = _verbose_stitching;
290 1056 : if (isParamValid("output_subdomain_id"))
291 : {
292 20 : xyd_opts.has_output_subdomain_id = true;
293 60 : xyd_opts.output_subdomain_id = getParam<SubdomainID>("output_subdomain_id");
294 : }
295 1056 : if (isParamValid("output_subdomain_name"))
296 : {
297 97 : xyd_opts.has_output_subdomain_name = true;
298 291 : xyd_opts.output_subdomain_name = getParam<SubdomainName>("output_subdomain_name");
299 : }
300 1056 : if (isParamValid("output_boundary"))
301 : {
302 10 : xyd_opts.has_output_boundary = true;
303 30 : xyd_opts.output_boundary = getParam<BoundaryName>("output_boundary");
304 : }
305 1056 : if (isParamValid("hole_boundaries"))
306 96 : xyd_opts.hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
307 :
308 352 : std::vector<std::unique_ptr<MeshBase>> hole_meshes(_hole_ptrs.size());
309 673 : for (auto hole_i : index_range(_hole_ptrs))
310 321 : hole_meshes[hole_i] = std::move(*_hole_ptrs[hole_i]);
311 :
312 352 : std::unique_ptr<MeshBase> boundary_mesh = std::move(_bdy_ptr);
313 :
314 : // Preserve the user-facing output_boundary so we can restore it after the outer-ring stitch-back
315 : // (we'll temporarily replace it with a sentinel name to locate the seam).
316 352 : const bool user_has_output_boundary = xyd_opts.has_output_boundary;
317 352 : const BoundaryName user_output_boundary = xyd_opts.output_boundary;
318 352 : const BoundaryName tmp_outer_name("__xyd_bdry_layer_tmp_outer__");
319 :
320 352 : const bool using_outer_layer = (_outer_boundary_layer_num > 0);
321 352 : const SubdomainID layer_sd_id =
322 352 : xyd_opts.has_output_subdomain_id ? xyd_opts.output_subdomain_id : SubdomainID(0);
323 : const SubdomainName layer_sd_name =
324 352 : xyd_opts.has_output_subdomain_name ? xyd_opts.output_subdomain_name : SubdomainName();
325 :
326 : // Build outer boundary-layer ring if requested. The ring's innermost side (bcid 1) becomes the
327 : // outer constraint for the interior triangulation; we then stitch a clone of the ring back onto
328 : // the result at that seam.
329 352 : std::unique_ptr<MeshBase> outer_ring_clone;
330 352 : if (using_outer_layer)
331 : {
332 : auto outer_ring = BoundaryLayerUtils::buildBoundaryLayerRing(*this,
333 10 : *boundary_mesh,
334 : xyd_opts.input_boundary_names,
335 10 : _outer_boundary_layer_num,
336 10 : _outer_boundary_layer_thickness,
337 10 : _outer_boundary_layer_bias,
338 : /*outward=*/false,
339 10 : _tri_elem_type,
340 : layer_sd_id,
341 10 : layer_sd_name);
342 10 : outer_ring_clone = outer_ring->clone();
343 10 : boundary_mesh = std::move(outer_ring);
344 20 : xyd_opts.input_boundary_names = {BoundaryName("1")};
345 10 : xyd_opts.has_output_boundary = true;
346 10 : xyd_opts.output_boundary = tmp_outer_name;
347 10 : }
348 :
349 : // Build hole boundary-layer rings if requested. Each ring is stitched into the result by the
350 : // standard hole-stitching path in triangulateWithDelaunay (with stitch_holes forced true).
351 673 : for (auto hole_i : index_range(_hole_ptrs))
352 : {
353 321 : if (hole_i < _holes_boundary_layer_num.size() && _holes_boundary_layer_num[hole_i] > 0)
354 : {
355 20 : const bool keep_input = (hole_i < _stitch_holes.size() && _stitch_holes[hole_i]);
356 : auto hole_ring =
357 : BoundaryLayerUtils::buildBoundaryLayerRing(*this,
358 20 : *hole_meshes[hole_i],
359 20 : std::vector<BoundaryName>(),
360 20 : _holes_boundary_layer_num[hole_i],
361 20 : _holes_boundary_layer_thickness[hole_i],
362 20 : _holes_boundary_layer_bias[hole_i],
363 : /*outward=*/true,
364 20 : _tri_elem_type,
365 : layer_sd_id,
366 40 : layer_sd_name);
367 :
368 20 : if (keep_input)
369 : {
370 : // Stitch the original hole mesh into the ring at ring's innermost side (bcid 1).
371 10 : auto & ring_u = dynamic_cast<UnstructuredMesh &>(*hole_ring);
372 10 : auto & inp_u = dynamic_cast<UnstructuredMesh &>(*hole_meshes[hole_i]);
373 10 : libMesh::MeshSerializer s1(ring_u), s2(inp_u);
374 10 : if (!ring_u.is_prepared())
375 10 : ring_u.prepare_for_use();
376 10 : if (!inp_u.is_prepared())
377 10 : inp_u.prepare_for_use();
378 : // Renumber input mesh boundary ids so they don't overlap with the ring's, mirroring the
379 : // approach in StitchMeshGenerator. This avoids degenerate stitching when both meshes
380 : // contain the same bcids on different sides.
381 10 : const auto & ring_bids = ring_u.get_boundary_info().get_global_boundary_ids();
382 10 : const auto inp_bids = inp_u.get_boundary_info().get_global_boundary_ids();
383 10 : const auto max_bid = std::max(*ring_bids.rbegin(),
384 10 : inp_bids.empty() ? boundary_id_type(0) : *inp_bids.rbegin());
385 10 : BoundaryID ext_id = 1;
386 10 : bool overlap = false;
387 50 : for (auto b : inp_bids)
388 40 : if (ring_bids.count(b))
389 10 : overlap = true;
390 10 : if (overlap)
391 : {
392 10 : BoundaryID idx = 1;
393 50 : for (auto b : inp_bids)
394 : {
395 40 : const auto new_b = max_bid + (idx++);
396 40 : inp_u.get_boundary_info().renumber_id(b, new_b);
397 : }
398 10 : ext_id = max_bid + idx;
399 : }
400 : else
401 0 : ext_id = max_bid + 1;
402 10 : inp_u.comm().max(ext_id);
403 10 : bool has_ext = false;
404 10 : MooseMeshUtils::addExternalBoundary(inp_u, ext_id, has_ext);
405 : mooseAssert(has_ext, "A 2D-XY mesh should have an external boundary.");
406 : const auto ring_u_ext_id =
407 10 : std::max(MooseMeshUtils::getNextFreeBoundaryID(ring_u), BoundaryID(ext_id + 1));
408 10 : MooseMeshUtils::changeBoundaryId(ring_u, 1, ring_u_ext_id, false);
409 10 : if (xyd_opts.hole_boundary_inner_id_defaults.size() <= hole_i)
410 10 : xyd_opts.hole_boundary_inner_id_defaults.resize(_hole_ptrs.size());
411 10 : xyd_opts.hole_boundary_inner_id_defaults[hole_i] = {ring_u_ext_id};
412 : // we want to keep the ring's original inner bcid (1) for later use.
413 10 : ring_u.stitch_meshes(inp_u,
414 : 1,
415 : ext_id,
416 : TOLERANCE,
417 : /*clear_stitched_bcids=*/true,
418 10 : _verbose_stitching,
419 10 : _algorithm == "BINARY",
420 : /*enforce_all_nodes_match_on_boundaries=*/false,
421 : /*merge_boundary_nodes_all_or_nothing=*/false,
422 : /*remap_subdomain_ids=*/false);
423 10 : }
424 :
425 20 : hole_meshes[hole_i] = std::move(hole_ring);
426 20 : if (xyd_opts.stitch_holes.size() <= hole_i)
427 0 : xyd_opts.stitch_holes.resize(_hole_ptrs.size(), false);
428 20 : xyd_opts.stitch_holes[hole_i] = true;
429 : // The ring presents multiple external manifolds; tell triangulateWithDelaunay to use the
430 : // outermost (bcid (num_layers - 1) * 2) as the hole's outer boundary.
431 20 : if (xyd_opts.hole_boundary_id_filters.size() <= hole_i)
432 10 : xyd_opts.hole_boundary_id_filters.resize(_hole_ptrs.size());
433 20 : xyd_opts.hole_boundary_id_filters[hole_i] = {
434 20 : std::size_t((_holes_boundary_layer_num[hole_i] - 1) * 2)};
435 20 : }
436 : }
437 :
438 : auto result = MeshTriangulationUtils::triangulateWithDelaunay(
439 364 : *this, std::move(boundary_mesh), std::move(hole_meshes), xyd_opts);
440 :
441 : // Stitch the outer ring clone back onto the interior triangulation at the recorded seam.
442 340 : if (outer_ring_clone)
443 : {
444 30 : auto sentinel_ids = MooseMeshUtils::getBoundaryIDs(*result, {tmp_outer_name}, false);
445 10 : const boundary_id_type sentinel_id = sentinel_ids[0];
446 :
447 : // Preserve the ring's outermost bcid (= (num_layers - 1) * 2) by renaming it to a high temp
448 : // value so the post-stitch rename can recover it as the final outer bcid.
449 10 : const boundary_id_type ring_outermost_orig =
450 10 : boundary_id_type((_outer_boundary_layer_num - 1) * 2);
451 : // The maximum boundary ID in the outer ring is _outer_boundary_layer_num * 2 - 1, so
452 : // _outer_boundary_layer_num * 2 is safe for itself. We need the maximum boundary ID of the
453 : // result mesh too.
454 : const boundary_id_type ring_outermost_temp =
455 20 : std::max(boundary_id_type(_outer_boundary_layer_num * 2),
456 10 : MooseMeshUtils::getNextFreeBoundaryID(*result));
457 10 : libMesh::MeshTools::Modification::change_boundary_id(
458 10 : *outer_ring_clone, ring_outermost_orig, ring_outermost_temp);
459 :
460 10 : auto & result_u = dynamic_cast<UnstructuredMesh &>(*result);
461 10 : auto & clone_u = dynamic_cast<UnstructuredMesh &>(*outer_ring_clone);
462 10 : libMesh::MeshSerializer s1(result_u), s2(clone_u);
463 10 : result_u.stitch_meshes(clone_u,
464 : sentinel_id,
465 : 1,
466 : TOLERANCE,
467 : /*clear_stitched_bcids=*/true,
468 10 : _verbose_stitching,
469 10 : _algorithm == "BINARY");
470 :
471 10 : libMesh::MeshTools::Modification::change_boundary_id(*result, ring_outermost_temp, sentinel_id);
472 :
473 10 : if (user_has_output_boundary)
474 : {
475 0 : auto user_id = MooseMeshUtils::getBoundaryIDs(*result, {user_output_boundary}, true).front();
476 0 : if (user_id != sentinel_id)
477 0 : libMesh::MeshTools::Modification::change_boundary_id(*result, sentinel_id, user_id);
478 0 : result->get_boundary_info().sideset_name(user_id) = user_output_boundary;
479 : }
480 :
481 10 : result->unset_is_prepared();
482 10 : }
483 :
484 680 : return result;
485 402 : }
|
