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