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