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 "MooseUtils.h"
15 :
16 : #include "libmesh/elem.h"
17 : #include "libmesh/enum_to_string.h"
18 : #include "libmesh/int_range.h"
19 : #include "libmesh/mesh_modification.h"
20 : #include "libmesh/mesh_serializer.h"
21 : #include "libmesh/mesh_triangle_holes.h"
22 : #include "libmesh/parsed_function.h"
23 : #include "libmesh/poly2tri_triangulator.h"
24 : #include "libmesh/unstructured_mesh.h"
25 : #include "DelimitedFileReader.h"
26 :
27 : registerMooseObject("MooseApp", XYDelaunayGenerator);
28 :
29 : InputParameters
30 15001 : XYDelaunayGenerator::validParams()
31 : {
32 15001 : InputParameters params = MeshGenerator::validParams();
33 :
34 15001 : MooseEnum algorithm("BINARY EXHAUSTIVE", "BINARY");
35 15001 : MooseEnum tri_elem_type("TRI3 TRI6 TRI7 DEFAULT", "DEFAULT");
36 :
37 15001 : params.addRequiredParam<MeshGeneratorName>(
38 : "boundary",
39 : "The input MeshGenerator defining the output outer boundary and required Steiner points.");
40 15001 : params.addParam<std::vector<BoundaryName>>(
41 : "input_boundary_names", "2D-input-mesh boundaries defining the output mesh outer boundary");
42 15001 : params.addParam<std::vector<SubdomainName>>(
43 : "input_subdomain_names", "1D-input-mesh subdomains defining the output mesh outer boundary");
44 45003 : params.addParam<unsigned int>("add_nodes_per_boundary_segment",
45 30002 : 0,
46 : "How many more nodes to add in each outer boundary segment.");
47 45003 : params.addParam<bool>(
48 30002 : "refine_boundary", true, "Whether to allow automatically refining the outer boundary.");
49 :
50 15001 : params.addParam<SubdomainName>("output_subdomain_name",
51 : "Subdomain name to set on new triangles.");
52 15001 : params.addParam<SubdomainID>("output_subdomain_id", "Subdomain id to set on new triangles.");
53 :
54 15001 : params.addParam<BoundaryName>(
55 : "output_boundary",
56 : "Boundary name to set on new outer boundary. Default ID: 0 if no hole meshes are stitched; "
57 : "or maximum boundary ID of all the stitched hole meshes + 1.");
58 15001 : params.addParam<std::vector<BoundaryName>>(
59 : "hole_boundaries",
60 : "Boundary names to set on holes. Default IDs are numbered up from 1 if no hole meshes are "
61 : "stitched; or from maximum boundary ID of all the stitched hole meshes + 2.");
62 :
63 45003 : params.addParam<bool>(
64 : "verify_holes",
65 30002 : true,
66 : "Verify holes do not intersect boundary or each other. Asymptotically costly.");
67 :
68 45003 : params.addParam<bool>("smooth_triangulation",
69 30002 : false,
70 : "Whether to do Laplacian mesh smoothing on the generated triangles.");
71 45003 : params.addParam<std::vector<MeshGeneratorName>>(
72 30002 : "holes", std::vector<MeshGeneratorName>(), "The MeshGenerators that define mesh holes.");
73 45003 : params.addParam<std::vector<bool>>(
74 30002 : "stitch_holes", std::vector<bool>(), "Whether to stitch to the mesh defining each hole.");
75 45003 : params.addParam<std::vector<bool>>("refine_holes",
76 30002 : std::vector<bool>(),
77 : "Whether to allow automatically refining each hole boundary.");
78 15001 : params.addRangeCheckedParam<Real>(
79 : "desired_area",
80 : 0,
81 : "desired_area>=0",
82 : "Desired (maximum) triangle area, or 0 to skip uniform refinement");
83 45003 : params.addParam<std::string>(
84 : "desired_area_func",
85 30002 : std::string(),
86 : "Desired area as a function of x,y; omit to skip non-uniform refinement");
87 :
88 45003 : params.addParam<bool>("use_auto_area_func",
89 30002 : false,
90 : "Use the automatic area function for the triangle meshing region.");
91 45003 : params.addParam<Real>(
92 : "auto_area_func_default_size",
93 30002 : 0,
94 : "Background size for automatic area function, or 0 to use non background size");
95 45003 : params.addParam<Real>("auto_area_func_default_size_dist",
96 30002 : -1.0,
97 : "Effective distance of background size for automatic area "
98 : "function, or negative to use non background size");
99 45003 : params.addParam<unsigned int>("auto_area_function_num_points",
100 30002 : 10,
101 : "Maximum number of nearest points used for the inverse distance "
102 : "interpolation algorithm for automatic area function calculation.");
103 45003 : params.addRangeCheckedParam<Real>(
104 : "auto_area_function_power",
105 30002 : 1.0,
106 : "auto_area_function_power>0",
107 : "Polynomial power of the inverse distance interpolation algorithm for automatic area "
108 : "function calculation.");
109 :
110 15001 : params.addParam<MooseEnum>(
111 : "algorithm",
112 : algorithm,
113 : "Control the use of binary search for the nodes of the stitched surfaces.");
114 15001 : params.addParam<MooseEnum>(
115 : "tri_element_type", tri_elem_type, "Type of the triangular elements to be generated.");
116 45003 : params.addParam<bool>(
117 30002 : "verbose_stitching", false, "Whether mesh stitching should have verbose output.");
118 15001 : params.addParam<std::vector<Point>>("interior_points",
119 : {},
120 : "Interior node locations, if no smoothing is used. Any point "
121 : "outside the surface will not be meshed.");
122 15001 : params.addParam<std::vector<FileName>>(
123 : "interior_point_files", {}, "Text file(s) with the interior points, one per line");
124 15001 : params.addClassDescription("Triangulates meshes within boundaries defined by input meshes.");
125 :
126 15001 : params.addParamNamesToGroup(
127 : "use_auto_area_func auto_area_func_default_size auto_area_func_default_size_dist",
128 : "Automatic triangle meshing area control");
129 15001 : params.addParamNamesToGroup("interior_points interior_point_files",
130 : "Mandatory mesh interior nodes");
131 :
132 30002 : return params;
133 15001 : }
134 :
135 374 : XYDelaunayGenerator::XYDelaunayGenerator(const InputParameters & parameters)
136 : : MeshGenerator(parameters),
137 374 : _bdy_ptr(getMesh("boundary")),
138 374 : _add_nodes_per_boundary_segment(getParam<unsigned int>("add_nodes_per_boundary_segment")),
139 374 : _refine_bdy(getParam<bool>("refine_boundary")),
140 374 : _output_subdomain_id(0),
141 374 : _smooth_tri(getParam<bool>("smooth_triangulation")),
142 374 : _verify_holes(getParam<bool>("verify_holes")),
143 374 : _hole_ptrs(getMeshes("holes")),
144 374 : _stitch_holes(getParam<std::vector<bool>>("stitch_holes")),
145 374 : _refine_holes(getParam<std::vector<bool>>("refine_holes")),
146 374 : _desired_area(getParam<Real>("desired_area")),
147 374 : _desired_area_func(getParam<std::string>("desired_area_func")),
148 374 : _use_auto_area_func(getParam<bool>("use_auto_area_func")),
149 374 : _auto_area_func_default_size(getParam<Real>("auto_area_func_default_size")),
150 374 : _auto_area_func_default_size_dist(getParam<Real>("auto_area_func_default_size_dist")),
151 374 : _auto_area_function_num_points(getParam<unsigned int>("auto_area_function_num_points")),
152 374 : _auto_area_function_power(getParam<Real>("auto_area_function_power")),
153 374 : _algorithm(parameters.get<MooseEnum>("algorithm")),
154 374 : _tri_elem_type(parameters.get<MooseEnum>("tri_element_type")),
155 374 : _verbose_stitching(parameters.get<bool>("verbose_stitching")),
156 748 : _interior_points(getParam<std::vector<Point>>("interior_points"))
157 : {
158 320 : if ((_desired_area > 0.0 && !_desired_area_func.empty()) ||
159 1064 : (_desired_area > 0.0 && _use_auto_area_func) ||
160 370 : (!_desired_area_func.empty() && _use_auto_area_func))
161 4 : paramError("desired_area_func",
162 : "Only one of the three methods ('desired_area', 'desired_area_func', and "
163 : "'_use_auto_area_func') to set element area limit should be used.");
164 :
165 370 : if (!_use_auto_area_func)
166 720 : if (isParamSetByUser("auto_area_func_default_size") ||
167 720 : isParamSetByUser("auto_area_func_default_size_dist") ||
168 1800 : isParamSetByUser("auto_area_function_num_points") ||
169 720 : isParamSetByUser("auto_area_function_power"))
170 4 : paramError("use_auto_area_func",
171 : "If this parameter is set to false, the following parameters should not be set: "
172 : "'auto_area_func_default_size', 'auto_area_func_default_size_dist', "
173 : "'auto_area_function_num_points', 'auto_area_function_power'.");
174 :
175 366 : if (!_stitch_holes.empty() && _stitch_holes.size() != _hole_ptrs.size())
176 0 : paramError("stitch_holes", "Need one stitch_holes entry per hole, if specified.");
177 :
178 528 : for (auto hole_i : index_range(_stitch_holes))
179 162 : if (_stitch_holes[hole_i] && (hole_i >= _refine_holes.size() || _refine_holes[hole_i]))
180 0 : paramError("refine_holes", "Disable auto refine of any hole boundary to be stitched.");
181 :
182 366 : if (isParamValid("hole_boundaries"))
183 : {
184 44 : auto & hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
185 44 : if (hole_boundaries.size() != _hole_ptrs.size())
186 0 : paramError("hole_boundaries", "Need one hole_boundaries entry per hole, if specified.");
187 : }
188 : // Copied from MultiApp.C
189 366 : const auto & positions_files = getParam<std::vector<FileName>>("interior_point_files");
190 380 : for (const auto p_file_it : index_range(positions_files))
191 : {
192 14 : const std::string positions_file = positions_files[p_file_it];
193 14 : MooseUtils::DelimitedFileReader file(positions_file, &_communicator);
194 14 : file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
195 14 : file.read();
196 :
197 14 : const std::vector<Point> & data = file.getDataAsPoints();
198 70 : for (const auto & d : data)
199 56 : _interior_points.push_back(d);
200 14 : }
201 : bool has_duplicates =
202 366 : std::any_of(_interior_points.begin(),
203 : _interior_points.end(),
204 114 : [&](const Point & p)
205 114 : { return std::count(_interior_points.begin(), _interior_points.end(), p) > 1; });
206 366 : if (has_duplicates)
207 4 : paramError("interior_points", "Duplicate points were found in the provided interior points.");
208 362 : }
209 :
210 : std::unique_ptr<MeshBase>
211 346 : XYDelaunayGenerator::generate()
212 : {
213 : // Put the boundary mesh in a local pointer
214 : std::unique_ptr<UnstructuredMesh> mesh =
215 346 : dynamic_pointer_cast<UnstructuredMesh>(std::move(_bdy_ptr));
216 :
217 : // Get ready to triangulate the line segments we extract from it
218 346 : libMesh::Poly2TriTriangulator poly2tri(*mesh);
219 346 : poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
220 :
221 : // If we're using a user-requested subset of boundaries on that
222 : // mesh, get their ids.
223 346 : std::set<std::size_t> bdy_ids;
224 :
225 346 : if (isParamValid("input_boundary_names"))
226 : {
227 10 : if (isParamValid("input_subdomain_names"))
228 0 : paramError(
229 : "input_subdomain_names",
230 : "input_boundary_names and input_subdomain_names cannot both specify an outer boundary.");
231 :
232 10 : const auto & boundary_names = getParam<std::vector<BoundaryName>>("input_boundary_names");
233 20 : for (const auto & name : boundary_names)
234 : {
235 10 : auto bcid = MooseMeshUtils::getBoundaryID(name, *mesh);
236 10 : if (bcid == BoundaryInfo::invalid_id)
237 0 : paramError("input_boundary_names", name, " is not a boundary name in the input mesh");
238 :
239 10 : bdy_ids.insert(bcid);
240 : }
241 : }
242 :
243 346 : if (isParamValid("input_subdomain_names"))
244 : {
245 10 : const auto & subdomain_names = getParam<std::vector<SubdomainName>>("input_subdomain_names");
246 :
247 10 : const auto subdomain_ids = MooseMeshUtils::getSubdomainIDs(*mesh, subdomain_names);
248 :
249 : // Check that the requested subdomains exist in the mesh
250 10 : std::set<SubdomainID> subdomains;
251 10 : mesh->subdomain_ids(subdomains);
252 :
253 20 : for (auto i : index_range(subdomain_ids))
254 : {
255 10 : if (subdomain_ids[i] == Moose::INVALID_BLOCK_ID || !subdomains.count(subdomain_ids[i]))
256 0 : paramError(
257 0 : "input_subdomain_names", subdomain_names[i], " was not found in the boundary mesh");
258 :
259 10 : bdy_ids.insert(subdomain_ids[i]);
260 : }
261 10 : }
262 :
263 346 : if (!bdy_ids.empty())
264 20 : poly2tri.set_outer_boundary_ids(bdy_ids);
265 :
266 346 : poly2tri.set_interpolate_boundary_points(_add_nodes_per_boundary_segment);
267 346 : poly2tri.set_refine_boundary_allowed(_refine_bdy);
268 346 : poly2tri.set_verify_hole_boundaries(_verify_holes);
269 :
270 346 : poly2tri.desired_area() = _desired_area;
271 346 : poly2tri.minimum_angle() = 0; // Not yet supported
272 346 : poly2tri.smooth_after_generating() = _smooth_tri;
273 :
274 346 : std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
275 346 : std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(_hole_ptrs.size());
276 346 : std::vector<std::unique_ptr<MeshBase>> hole_ptrs(_hole_ptrs.size());
277 : // This tells us the element orders of the hole meshes
278 : // For the boundary meshes, it can be access through poly2tri.segment_midpoints.
279 346 : std::vector<bool> holes_with_midpoints(_hole_ptrs.size());
280 346 : bool stitch_second_order_holes(false);
281 :
282 : // Make sure pointers here aren't invalidated by a resize
283 346 : meshed_holes.reserve(_hole_ptrs.size());
284 634 : for (auto hole_i : index_range(_hole_ptrs))
285 : {
286 288 : hole_ptrs[hole_i] = std::move(*_hole_ptrs[hole_i]);
287 288 : if (!hole_ptrs[hole_i]->is_prepared())
288 174 : hole_ptrs[hole_i]->prepare_for_use();
289 288 : meshed_holes.emplace_back(*hole_ptrs[hole_i]);
290 288 : holes_with_midpoints[hole_i] = meshed_holes.back().n_midpoints();
291 288 : stitch_second_order_holes = _stitch_holes.empty()
292 450 : ? false
293 162 : : ((holes_with_midpoints[hole_i] && _stitch_holes[hole_i]) ||
294 : stitch_second_order_holes);
295 288 : if (hole_i < _refine_holes.size())
296 208 : meshed_holes.back().set_refine_boundary_allowed(_refine_holes[hole_i]);
297 :
298 288 : triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
299 : }
300 346 : if (stitch_second_order_holes && (_tri_elem_type == "TRI3" || _tri_elem_type == "DEFAULT"))
301 4 : paramError(
302 : "tri_element_type",
303 : "Cannot use first order elements with stitched quadratic element holes. Please try "
304 : "to specify a higher-order tri_element_type or reduce the order of the hole inputs.");
305 :
306 342 : if (!triangulator_hole_ptrs.empty())
307 200 : poly2tri.attach_hole_list(&triangulator_hole_ptrs);
308 :
309 342 : if (_desired_area_func != "")
310 : {
311 : // poly2tri will clone this so it's fine going out of scope
312 10 : libMesh::ParsedFunction<Real> area_func{_desired_area_func};
313 10 : poly2tri.set_desired_area_function(&area_func);
314 10 : }
315 332 : else if (_use_auto_area_func)
316 : {
317 30 : poly2tri.set_auto_area_function(
318 : this->comm(),
319 10 : _auto_area_function_num_points,
320 10 : _auto_area_function_power,
321 10 : _auto_area_func_default_size > 0.0 ? _auto_area_func_default_size : 0.0,
322 10 : _auto_area_func_default_size_dist > 0.0 ? _auto_area_func_default_size_dist : -1.0);
323 : }
324 :
325 342 : if (_tri_elem_type == "TRI6")
326 20 : poly2tri.elem_type() = libMesh::ElemType::TRI6;
327 322 : else if (_tri_elem_type == "TRI7")
328 20 : poly2tri.elem_type() = libMesh::ElemType::TRI7;
329 : // Add interior points before triangulating. Only points inside the boundaries
330 : // will be meshed.
331 452 : for (auto & point : _interior_points)
332 110 : mesh->add_point(point);
333 :
334 342 : poly2tri.triangulate();
335 :
336 334 : if (isParamValid("output_subdomain_id"))
337 10 : _output_subdomain_id = getParam<SubdomainID>("output_subdomain_id");
338 :
339 334 : if (isParamValid("output_subdomain_name"))
340 : {
341 96 : auto output_subdomain_name = getParam<SubdomainName>("output_subdomain_name");
342 96 : auto id = MooseMeshUtils::getSubdomainID(output_subdomain_name, *mesh);
343 :
344 96 : if (id == Elem::invalid_subdomain_id)
345 : {
346 44 : if (!isParamValid("output_subdomain_id"))
347 : {
348 : // We'll probably need to make a new ID, then
349 34 : _output_subdomain_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
350 :
351 : // But check the hole meshes for our output subdomain name too
352 102 : for (auto & hole_ptr : hole_ptrs)
353 : {
354 68 : auto possible_sbdid = MooseMeshUtils::getSubdomainID(output_subdomain_name, *hole_ptr);
355 : // Huh, it was in one of them
356 68 : if (possible_sbdid != Elem::invalid_subdomain_id)
357 : {
358 0 : _output_subdomain_id = possible_sbdid;
359 0 : break;
360 : }
361 68 : _output_subdomain_id =
362 68 : std::max(_output_subdomain_id, MooseMeshUtils::getNextFreeSubdomainID(*hole_ptr));
363 : }
364 : }
365 : }
366 : else
367 : {
368 52 : if (isParamValid("output_subdomain_id"))
369 : {
370 0 : if (id != _output_subdomain_id)
371 0 : paramError("output_subdomain_name",
372 : "name has been used by the input meshes and the corresponding id is not equal "
373 : "to 'output_subdomain_id'");
374 : }
375 : else
376 52 : _output_subdomain_id = id;
377 : }
378 : // We do not want to set an empty subdomain name
379 96 : if (output_subdomain_name.size())
380 96 : mesh->subdomain_name(_output_subdomain_id) = output_subdomain_name;
381 96 : }
382 :
383 334 : if (_smooth_tri || _output_subdomain_id)
384 7669 : for (auto elem : mesh->element_ptr_range())
385 : {
386 : mooseAssert(elem->type() ==
387 : (_tri_elem_type == "TRI6" ? TRI6 : (_tri_elem_type == "TRI7" ? TRI7 : TRI3)),
388 : "Unexpected element type " << libMesh::Utility::enum_to_string(elem->type())
389 : << " found in triangulation");
390 :
391 7553 : elem->subdomain_id() = _output_subdomain_id;
392 :
393 : // I do not trust Laplacian mesh smoothing not to invert
394 : // elements near reentrant corners. Eventually we'll add better
395 : // smoothing options, but even those might have failure cases.
396 : // Better to always do extra tests here than to ever let users
397 : // try to run on a degenerate mesh.
398 7553 : if (_smooth_tri)
399 : {
400 2176 : auto cross_prod = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
401 :
402 2176 : if (cross_prod(2) <= 0)
403 0 : mooseError("Inverted element found in triangulation.\n"
404 : "Laplacian smoothing can create these at reentrant corners; disable it?");
405 : }
406 116 : }
407 :
408 334 : const bool use_binary_search = (_algorithm == "BINARY");
409 :
410 : // The hole meshes are specified by the user, so they could have any
411 : // BCID or no BCID or any combination of BCIDs on their outer
412 : // boundary, so we'll have to set our own BCID to use for stitching
413 : // there. We'll need to check all the holes for used BCIDs, if we
414 : // want to pick a new ID on hole N that doesn't conflict with any
415 : // IDs on hole M < N (or with the IDs on the new triangulation)
416 :
417 : // The new triangulation by default assigns BCID i+1 to hole i ...
418 : // but we can't even use this for mesh stitching, because we can't
419 : // be sure it isn't also already in use on the hole's mesh and so we
420 : // won't be able to safely clear it afterwards.
421 334 : const boundary_id_type end_bcid = hole_ptrs.size() + 1;
422 :
423 : // For the hole meshes that need to be stitched, we would like to make sure the hole boundary ids
424 : // and output boundary id are not conflicting with the existing boundary ids of the hole meshes to
425 : // be stitched.
426 334 : BoundaryID free_boundary_id = 0;
427 334 : if (_stitch_holes.size())
428 : {
429 282 : for (auto hole_i : index_range(hole_ptrs))
430 : {
431 158 : if (_stitch_holes[hole_i])
432 : {
433 118 : free_boundary_id =
434 118 : std::max(free_boundary_id, MooseMeshUtils::getNextFreeBoundaryID(*hole_ptrs[hole_i]));
435 118 : hole_ptrs[hole_i]->comm().max(free_boundary_id);
436 : }
437 : }
438 282 : for (auto h : index_range(hole_ptrs))
439 158 : libMesh::MeshTools::Modification::change_boundary_id(*mesh, h + 1, h + 1 + free_boundary_id);
440 : }
441 334 : boundary_id_type new_hole_bcid = end_bcid + free_boundary_id;
442 :
443 : // We might be overriding the default bcid numbers. We have to be
444 : // careful about how we renumber, though. We pick unused temporary
445 : // numbers because e.g. "0->2, 2->0" is impossible to do
446 : // sequentially, but "0->N, 2->N+2, N->2, N+2->0" works.
447 334 : libMesh::MeshTools::Modification::change_boundary_id(
448 668 : *mesh, 0, (isParamValid("output_boundary") ? end_bcid : 0) + free_boundary_id);
449 :
450 334 : if (isParamValid("hole_boundaries"))
451 : {
452 36 : auto hole_boundaries = getParam<std::vector<BoundaryName>>("hole_boundaries");
453 36 : auto hole_boundary_ids = MooseMeshUtils::getBoundaryIDs(*mesh, hole_boundaries, true);
454 :
455 82 : for (auto h : index_range(hole_ptrs))
456 46 : libMesh::MeshTools::Modification::change_boundary_id(
457 46 : *mesh, h + 1 + free_boundary_id, h + 1 + free_boundary_id + end_bcid);
458 :
459 82 : for (auto h : index_range(hole_ptrs))
460 : {
461 46 : libMesh::MeshTools::Modification::change_boundary_id(
462 46 : *mesh, h + 1 + free_boundary_id + end_bcid, hole_boundary_ids[h]);
463 46 : mesh->get_boundary_info().sideset_name(hole_boundary_ids[h]) = hole_boundaries[h];
464 46 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(hole_boundary_ids[h] + 1));
465 : }
466 36 : }
467 :
468 334 : if (isParamValid("output_boundary"))
469 : {
470 10 : const BoundaryName output_boundary = getParam<BoundaryName>("output_boundary");
471 : const std::vector<BoundaryID> output_boundary_id =
472 20 : MooseMeshUtils::getBoundaryIDs(*mesh, {output_boundary}, true);
473 :
474 10 : libMesh::MeshTools::Modification::change_boundary_id(
475 10 : *mesh, end_bcid + free_boundary_id, output_boundary_id[0]);
476 10 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = output_boundary;
477 :
478 10 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(output_boundary_id[0] + 1));
479 10 : }
480 :
481 334 : bool doing_stitching = false;
482 :
483 618 : for (auto hole_i : index_range(hole_ptrs))
484 : {
485 284 : const MeshBase & hole_mesh = *hole_ptrs[hole_i];
486 284 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
487 284 : const std::set<boundary_id_type> & local_hole_bcids = hole_boundary_info.get_boundary_ids();
488 :
489 284 : if (!local_hole_bcids.empty())
490 184 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(*local_hole_bcids.rbegin() + 1));
491 284 : hole_mesh.comm().max(new_hole_bcid);
492 :
493 284 : if (hole_i < _stitch_holes.size() && _stitch_holes[hole_i])
494 118 : doing_stitching = true;
495 : }
496 :
497 334 : const boundary_id_type inner_bcid = new_hole_bcid + 1;
498 :
499 : // libMesh mesh stitching still requires a serialized mesh, and it's
500 : // cheaper to do that once than to do it once-per-hole
501 334 : libMesh::MeshSerializer serial(*mesh, doing_stitching);
502 :
503 : // Define a reference map variable for subdomain map
504 334 : auto & main_subdomain_map = mesh->set_subdomain_name_map();
505 618 : for (auto hole_i : index_range(hole_ptrs))
506 : {
507 284 : if (hole_i < _stitch_holes.size() && _stitch_holes[hole_i])
508 : {
509 118 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(*hole_ptrs[hole_i]);
510 : // increase hole mesh order if the triangulation mesh has higher order
511 118 : if (!holes_with_midpoints[hole_i])
512 : {
513 98 : if (_tri_elem_type == "TRI6")
514 10 : hole_mesh.all_second_order();
515 88 : else if (_tri_elem_type == "TRI7")
516 10 : hole_mesh.all_complete_order();
517 : }
518 118 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
519 :
520 : // Our algorithm here requires a serialized Mesh. To avoid
521 : // redundant serialization and deserialization (libMesh
522 : // MeshedHole and stitch_meshes still also require
523 : // serialization) we'll do the serialization up front.
524 118 : libMesh::MeshSerializer serial_hole(hole_mesh);
525 :
526 : // It would have been nicer for MeshedHole to add the BCID
527 : // itself, but we want MeshedHole to work with a const mesh.
528 : // We'll still use MeshedHole, for its code distinguishing
529 : // outer boundaries from inner boundaries on a
530 : // hole-with-holes.
531 118 : libMesh::TriangulatorInterface::MeshedHole mh{hole_mesh};
532 :
533 : // We have to translate from MeshedHole points to mesh
534 : // sides.
535 118 : std::unordered_map<Point, Point> next_hole_boundary_point;
536 118 : const int np = mh.n_points();
537 1336 : for (auto pi : make_range(1, np))
538 1218 : next_hole_boundary_point[mh.point(pi - 1)] = mh.point(pi);
539 118 : next_hole_boundary_point[mh.point(np - 1)] = mh.point(0);
540 :
541 : #ifndef NDEBUG
542 : int found_hole_sides = 0;
543 : #endif
544 4960 : for (auto elem : hole_mesh.element_ptr_range())
545 : {
546 4842 : if (elem->dim() != 2)
547 0 : mooseError("Non 2-D element found in hole; stitching is not supported.");
548 :
549 4842 : auto ns = elem->n_sides();
550 20250 : for (auto s : make_range(ns))
551 : {
552 15408 : auto it_s = next_hole_boundary_point.find(elem->point(s));
553 15408 : if (it_s != next_hole_boundary_point.end())
554 3540 : if (it_s->second == elem->point((s + 1) % ns))
555 : {
556 1336 : hole_boundary_info.add_side(elem, s, new_hole_bcid);
557 : #ifndef NDEBUG
558 : ++found_hole_sides;
559 : #endif
560 : }
561 : }
562 118 : }
563 : mooseAssert(found_hole_sides == np, "Failed to find full outer boundary of meshed hole");
564 :
565 118 : auto & mesh_boundary_info = mesh->get_boundary_info();
566 : #ifndef NDEBUG
567 : int found_inner_sides = 0;
568 : #endif
569 11642 : for (auto elem : mesh->element_ptr_range())
570 : {
571 11524 : auto ns = elem->n_sides();
572 46132 : for (auto s : make_range(ns))
573 : {
574 34608 : auto it_s = next_hole_boundary_point.find(elem->point((s + 1) % ns));
575 34608 : if (it_s != next_hole_boundary_point.end())
576 4618 : if (it_s->second == elem->point(s))
577 : {
578 1336 : mesh_boundary_info.add_side(elem, s, inner_bcid);
579 : #ifndef NDEBUG
580 : ++found_inner_sides;
581 : #endif
582 : }
583 : }
584 118 : }
585 : mooseAssert(found_inner_sides == np, "Failed to find full boundary around meshed hole");
586 :
587 : // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
588 : // subdomain map
589 118 : const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
590 118 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
591 :
592 118 : mesh->stitch_meshes(hole_mesh,
593 : inner_bcid,
594 : new_hole_bcid,
595 : TOLERANCE,
596 : /*clear_stitched_bcids*/ true,
597 118 : _verbose_stitching,
598 : use_binary_search);
599 118 : }
600 : }
601 : // Check if one SubdomainName is shared by more than one subdomain ids
602 334 : std::set<SubdomainName> main_subdomain_map_name_list;
603 488 : for (auto const & id_name_pair : main_subdomain_map)
604 154 : main_subdomain_map_name_list.emplace(id_name_pair.second);
605 334 : if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
606 4 : paramError("holes", "The hole meshes contain subdomain name maps with conflicts.");
607 :
608 330 : mesh->set_isnt_prepared();
609 660 : return mesh;
610 396 : }
|
