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 "MeshTriangulationUtils.h"
11 : #include "MooseMeshUtils.h"
12 : #include "CastUniquePointer.h"
13 :
14 : #include "libmesh/elem.h"
15 : #include "libmesh/boundary_info.h"
16 : #include "libmesh/int_range.h"
17 : #include "libmesh/enum_to_string.h"
18 : #include "libmesh/mesh_modification.h"
19 : #include "libmesh/mesh_serializer.h"
20 : #include "libmesh/mesh_triangle_holes.h"
21 : #include "libmesh/parsed_function.h"
22 : #include "libmesh/poly2tri_triangulator.h"
23 : #include "libmesh/unstructured_mesh.h"
24 :
25 : using namespace libMesh;
26 :
27 : namespace MeshTriangulationUtils
28 : {
29 :
30 : std::unique_ptr<MeshBase>
31 722 : triangulateWithDelaunay(MeshGenerator & mg,
32 : std::unique_ptr<MeshBase> boundary_mesh,
33 : std::vector<std::unique_ptr<MeshBase>> hole_meshes,
34 : const XYDelaunayOptions & xyd_opts)
35 : {
36 : // Put the boundary mesh in a local pointer
37 : std::unique_ptr<UnstructuredMesh> mesh =
38 722 : dynamic_pointer_cast<UnstructuredMesh>(std::move(boundary_mesh));
39 :
40 : // Get ready to triangulate the line segments we extract from it
41 722 : libMesh::Poly2TriTriangulator poly2tri(*mesh);
42 722 : poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
43 :
44 : // If we're using a user-requested subset of boundaries on that
45 : // mesh, get their ids.
46 722 : std::set<std::size_t> bdy_ids;
47 :
48 722 : if (!xyd_opts.input_boundary_names.empty())
49 : {
50 20 : if (!xyd_opts.input_subdomain_names.empty())
51 0 : mg.paramError(
52 : "input_subdomain_names",
53 : "input_boundary_names and input_subdomain_names cannot both specify an outer boundary.");
54 :
55 40 : for (const auto & name : xyd_opts.input_boundary_names)
56 : {
57 20 : auto bcid = MooseMeshUtils::getBoundaryID(name, *mesh);
58 20 : if (bcid == BoundaryInfo::invalid_id)
59 0 : mg.paramError("input_boundary_names", name, " is not a boundary name in the input mesh");
60 :
61 20 : bdy_ids.insert(bcid);
62 : }
63 : }
64 :
65 722 : if (!xyd_opts.input_subdomain_names.empty())
66 : {
67 : // Make sure subdomain info caches are up to date
68 10 : if (!mesh->preparation().has_cached_elem_data)
69 10 : mesh->cache_elem_data();
70 :
71 : const auto subdomain_ids =
72 10 : MooseMeshUtils::getSubdomainIDs(*mesh, xyd_opts.input_subdomain_names);
73 :
74 : // Check that the requested subdomains exist in the mesh
75 10 : std::set<SubdomainID> subdomains;
76 10 : mesh->subdomain_ids(subdomains);
77 :
78 20 : for (auto i : index_range(subdomain_ids))
79 : {
80 10 : if (subdomain_ids[i] == Moose::INVALID_BLOCK_ID || !subdomains.count(subdomain_ids[i]))
81 0 : mg.paramError("input_subdomain_names",
82 0 : xyd_opts.input_subdomain_names[i],
83 : " was not found in the boundary mesh");
84 :
85 10 : bdy_ids.insert(subdomain_ids[i]);
86 : }
87 10 : }
88 :
89 722 : if (!bdy_ids.empty())
90 30 : poly2tri.set_outer_boundary_ids(bdy_ids);
91 :
92 722 : poly2tri.set_interpolate_boundary_points(xyd_opts.add_nodes_per_boundary_segment);
93 722 : poly2tri.set_refine_boundary_allowed(xyd_opts.refine_bdy);
94 722 : poly2tri.set_verify_hole_boundaries(xyd_opts.verify_holes);
95 :
96 722 : poly2tri.desired_area() = xyd_opts.desired_area;
97 722 : poly2tri.minimum_angle() = 0; // Not yet supported
98 722 : poly2tri.smooth_after_generating() = xyd_opts.smooth_tri;
99 :
100 722 : std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
101 1444 : std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(hole_meshes.size());
102 : // This tells us the element orders of the hole meshes
103 : // For the boundary meshes, it can be access through poly2tri.segment_midpoints.
104 722 : std::vector<bool> holes_with_midpoints(hole_meshes.size());
105 722 : bool stitch_second_order_holes(false);
106 :
107 : // Make sure pointers here aren't invalidated by a resize
108 722 : meshed_holes.reserve(hole_meshes.size());
109 1413 : for (auto hole_i : index_range(hole_meshes))
110 : {
111 691 : if (!hole_meshes[hole_i]->is_prepared())
112 438 : hole_meshes[hole_i]->prepare_for_use();
113 711 : if (hole_i < xyd_opts.hole_boundary_id_filters.size() &&
114 20 : !xyd_opts.hole_boundary_id_filters[hole_i].empty())
115 20 : meshed_holes.emplace_back(*hole_meshes[hole_i], xyd_opts.hole_boundary_id_filters[hole_i]);
116 : else
117 671 : meshed_holes.emplace_back(*hole_meshes[hole_i]);
118 691 : holes_with_midpoints[hole_i] = meshed_holes.back().n_midpoints();
119 691 : stitch_second_order_holes =
120 691 : xyd_opts.stitch_holes.empty()
121 1240 : ? false
122 549 : : ((holes_with_midpoints[hole_i] && xyd_opts.stitch_holes[hole_i]) ||
123 : stitch_second_order_holes);
124 691 : if (hole_i < xyd_opts.refine_holes.size())
125 591 : meshed_holes.back().set_refine_boundary_allowed(xyd_opts.refine_holes[hole_i]);
126 :
127 691 : triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
128 : }
129 765 : if (stitch_second_order_holes &&
130 43 : (xyd_opts.tri_elem_type == "TRI3" || xyd_opts.tri_elem_type == "DEFAULT"))
131 6 : mg.paramError(
132 : "tri_element_type",
133 : "Cannot use first order elements with stitched quadratic element holes. Please try "
134 : "to specify a higher-order tri_element_type or reduce the order of the hole inputs.");
135 :
136 719 : if (!triangulator_hole_ptrs.empty())
137 585 : poly2tri.attach_hole_list(&triangulator_hole_ptrs);
138 :
139 719 : if (xyd_opts.desired_area_func != "")
140 : {
141 : // poly2tri will clone this so it's fine going out of scope
142 10 : libMesh::ParsedFunction<Real> area_func{xyd_opts.desired_area_func};
143 10 : poly2tri.set_desired_area_function(&area_func);
144 10 : }
145 709 : else if (xyd_opts.use_auto_area_func)
146 : {
147 30 : poly2tri.set_auto_area_function(
148 : mg.comm(),
149 10 : xyd_opts.auto_area_function_num_points,
150 10 : xyd_opts.auto_area_function_power,
151 10 : xyd_opts.auto_area_func_default_size > 0.0 ? xyd_opts.auto_area_func_default_size : 0.0,
152 10 : xyd_opts.auto_area_func_default_size_dist > 0.0 ? xyd_opts.auto_area_func_default_size_dist
153 : : -1.0);
154 : }
155 :
156 719 : if (xyd_opts.tri_elem_type == "TRI6")
157 60 : poly2tri.elem_type() = libMesh::ElemType::TRI6;
158 659 : else if (xyd_opts.tri_elem_type == "TRI7")
159 20 : poly2tri.elem_type() = libMesh::ElemType::TRI7;
160 : // Add interior points before triangulating. Only points inside the boundaries
161 : // will be meshed.
162 829 : for (const auto & point : xyd_opts.interior_points)
163 110 : mesh->add_point(point);
164 :
165 719 : poly2tri.triangulate();
166 :
167 713 : SubdomainID output_subdomain_id =
168 713 : xyd_opts.has_output_subdomain_id ? xyd_opts.output_subdomain_id : 0;
169 :
170 713 : if (xyd_opts.has_output_subdomain_name)
171 : {
172 407 : auto id = MooseMeshUtils::getSubdomainID(xyd_opts.output_subdomain_name, *mesh);
173 :
174 407 : if (id == Elem::invalid_subdomain_id)
175 : {
176 353 : if (!xyd_opts.has_output_subdomain_id)
177 : {
178 : // We'll probably need to make a new ID, then
179 33 : output_subdomain_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
180 :
181 : // But check the hole meshes for our output subdomain name too
182 99 : for (auto & hole_ptr : hole_meshes)
183 : {
184 : auto possible_sbdid =
185 66 : MooseMeshUtils::getSubdomainID(xyd_opts.output_subdomain_name, *hole_ptr);
186 : // Huh, it was in one of them
187 66 : if (possible_sbdid != Elem::invalid_subdomain_id)
188 : {
189 0 : output_subdomain_id = possible_sbdid;
190 0 : break;
191 : }
192 66 : output_subdomain_id =
193 66 : std::max(output_subdomain_id, MooseMeshUtils::getNextFreeSubdomainID(*hole_ptr));
194 : }
195 : }
196 : }
197 : else
198 : {
199 54 : if (xyd_opts.has_output_subdomain_id)
200 : {
201 10 : if (id != output_subdomain_id)
202 0 : mg.paramError("output_subdomain_name",
203 : "name has been used by the input meshes and the corresponding id is not "
204 : "equal to 'output_subdomain_id'");
205 : }
206 : else
207 44 : output_subdomain_id = id;
208 : }
209 : // We do not want to set an empty subdomain name
210 407 : if (xyd_opts.output_subdomain_name.size())
211 407 : mesh->subdomain_name(output_subdomain_id) = xyd_opts.output_subdomain_name;
212 : }
213 :
214 713 : if (xyd_opts.smooth_tri || output_subdomain_id)
215 22990 : for (auto elem : mesh->element_ptr_range())
216 : {
217 : mooseAssert(elem->type() == (xyd_opts.tri_elem_type == "TRI6"
218 : ? TRI6
219 : : (xyd_opts.tri_elem_type == "TRI7" ? TRI7 : TRI3)),
220 : "Unexpected element type " << libMesh::Utility::enum_to_string(elem->type())
221 : << " found in triangulation");
222 :
223 22563 : elem->subdomain_id() = output_subdomain_id;
224 :
225 : // I do not trust Laplacian mesh smoothing not to invert
226 : // elements near reentrant corners. Eventually we'll add better
227 : // smoothing options, but even those might have failure cases.
228 : // Better to always do extra tests here than to ever let users
229 : // try to run on a degenerate mesh.
230 22563 : if (xyd_opts.smooth_tri)
231 : {
232 2176 : auto cross_prod = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
233 :
234 2176 : if (cross_prod(2) <= 0)
235 0 : mooseError("Inverted element found in triangulation.\n"
236 : "Laplacian smoothing can create these at reentrant corners; disable it?");
237 : }
238 427 : }
239 :
240 : // The hole meshes are specified by the user, so they could have any
241 : // BCID or no BCID or any combination of BCIDs on their outer
242 : // boundary, so we'll have to set our own BCID to use for stitching
243 : // there. We'll need to check all the holes for used BCIDs, if we
244 : // want to pick a new ID on hole N that doesn't conflict with any
245 : // IDs on hole M < N (or with the IDs on the new triangulation)
246 :
247 : // The new triangulation by default assigns BCID i+1 to hole i ...
248 : // but we can't even use this for mesh stitching, because we can't
249 : // be sure it isn't also already in use on the hole's mesh and so we
250 : // won't be able to safely clear it afterwards.
251 713 : const boundary_id_type end_bcid = hole_meshes.size() + 1;
252 :
253 : // If a hole has its boundary layer mesh, we need to move the hole bcid to the "real" hole
254 : // boundary in the boundary layer mesh. So we need to record them here.
255 713 : std::vector<BoundaryID> hole_boundary_rec(hole_meshes.size());
256 713 : std::iota(hole_boundary_rec.begin(), hole_boundary_rec.end(), 1);
257 :
258 : // For the hole meshes that need to be stitched, we would like to make sure the hole boundary ids
259 : // and output boundary id are not conflicting with the existing boundary ids of the hole meshes to
260 : // be stitched.
261 713 : BoundaryID free_boundary_id = 0;
262 713 : if (xyd_opts.stitch_holes.size())
263 : {
264 1049 : for (auto hole_i : index_range(hole_meshes))
265 : {
266 546 : if (xyd_opts.stitch_holes[hole_i])
267 : {
268 376 : free_boundary_id =
269 376 : std::max(free_boundary_id, MooseMeshUtils::getNextFreeBoundaryID(*hole_meshes[hole_i]));
270 376 : hole_meshes[hole_i]->comm().max(free_boundary_id);
271 : }
272 : }
273 1049 : for (auto h : index_range(hole_meshes))
274 : {
275 546 : libMesh::MeshTools::Modification::change_boundary_id(*mesh, h + 1, h + 1 + free_boundary_id);
276 546 : hole_boundary_rec[h] = h + 1 + free_boundary_id;
277 : }
278 : }
279 713 : boundary_id_type new_hole_bcid = end_bcid + free_boundary_id;
280 :
281 : // We might be overriding the default bcid numbers. We have to be
282 : // careful about how we renumber, though. We pick unused temporary
283 : // numbers because e.g. "0->2, 2->0" is impossible to do
284 : // sequentially, but "0->N, 2->N+2, N->2, N+2->0" works.
285 713 : libMesh::MeshTools::Modification::change_boundary_id(
286 713 : *mesh, 0, (xyd_opts.has_output_boundary ? end_bcid : 0) + free_boundary_id);
287 :
288 713 : if (!xyd_opts.hole_boundaries.empty())
289 : {
290 32 : auto hole_boundary_ids = MooseMeshUtils::getBoundaryIDs(*mesh, xyd_opts.hole_boundaries, true);
291 :
292 74 : for (auto h : index_range(hole_meshes))
293 42 : libMesh::MeshTools::Modification::change_boundary_id(
294 42 : *mesh, h + 1 + free_boundary_id, h + 1 + free_boundary_id + end_bcid);
295 :
296 74 : for (auto h : index_range(hole_meshes))
297 : {
298 42 : libMesh::MeshTools::Modification::change_boundary_id(
299 42 : *mesh, h + 1 + free_boundary_id + end_bcid, hole_boundary_ids[h]);
300 42 : hole_boundary_rec[h] = hole_boundary_ids[h];
301 42 : mesh->get_boundary_info().sideset_name(hole_boundary_ids[h]) = xyd_opts.hole_boundaries[h];
302 42 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(hole_boundary_ids[h] + 1));
303 : }
304 32 : }
305 :
306 713 : if (xyd_opts.has_output_boundary)
307 : {
308 : const std::vector<BoundaryID> output_boundary_id =
309 60 : MooseMeshUtils::getBoundaryIDs(*mesh, {xyd_opts.output_boundary}, true);
310 :
311 20 : libMesh::MeshTools::Modification::change_boundary_id(
312 20 : *mesh, end_bcid + free_boundary_id, output_boundary_id[0]);
313 20 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = xyd_opts.output_boundary;
314 :
315 20 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(output_boundary_id[0] + 1));
316 20 : }
317 :
318 713 : bool doing_stitching = false;
319 :
320 1401 : for (auto hole_i : index_range(hole_meshes))
321 : {
322 688 : const MeshBase & hole_mesh = *hole_meshes[hole_i];
323 688 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
324 688 : const std::set<boundary_id_type> & local_hole_bcids = hole_boundary_info.get_boundary_ids();
325 :
326 688 : if (!local_hole_bcids.empty())
327 438 : new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(*local_hole_bcids.rbegin() + 1));
328 688 : hole_mesh.comm().max(new_hole_bcid);
329 :
330 688 : if (hole_i < xyd_opts.stitch_holes.size() && xyd_opts.stitch_holes[hole_i])
331 376 : doing_stitching = true;
332 : }
333 :
334 713 : const boundary_id_type inner_bcid = new_hole_bcid + 1;
335 :
336 : // libMesh mesh stitching still requires a serialized mesh, and it's
337 : // cheaper to do that once than to do it once-per-hole
338 713 : libMesh::MeshSerializer serial(*mesh, doing_stitching);
339 :
340 : // Define a reference map variable for subdomain map
341 713 : auto & main_subdomain_map = mesh->set_subdomain_name_map();
342 1401 : for (auto hole_i : index_range(hole_meshes))
343 : {
344 688 : if (hole_i < xyd_opts.stitch_holes.size() && xyd_opts.stitch_holes[hole_i])
345 : {
346 376 : UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(*hole_meshes[hole_i]);
347 : // increase hole mesh order if the triangulation mesh has higher order
348 376 : if (!holes_with_midpoints[hole_i])
349 : {
350 336 : if (xyd_opts.tri_elem_type == "TRI6")
351 10 : hole_mesh.all_second_order();
352 326 : else if (xyd_opts.tri_elem_type == "TRI7")
353 10 : hole_mesh.all_complete_order();
354 : }
355 376 : auto & hole_boundary_info = hole_mesh.get_boundary_info();
356 :
357 : // Our algorithm here requires a serialized Mesh. To avoid
358 : // redundant serialization and deserialization (libMesh
359 : // MeshedHole and stitch_meshes still also require
360 : // serialization) we'll do the serialization up front.
361 376 : libMesh::MeshSerializer serial_hole(hole_mesh);
362 :
363 : // It would have been nicer for MeshedHole to add the BCID
364 : // itself, but we want MeshedHole to work with a const mesh.
365 : // We'll still use MeshedHole, for its code distinguishing
366 : // outer boundaries from inner boundaries on a
367 : // hole-with-holes.
368 376 : const auto & hole_bdy_id_filter = (hole_i < xyd_opts.hole_boundary_id_filters.size())
369 376 : ? xyd_opts.hole_boundary_id_filters[hole_i]
370 376 : : std::set<std::size_t>();
371 376 : libMesh::TriangulatorInterface::MeshedHole mh{hole_mesh, hole_bdy_id_filter};
372 :
373 : // We have to translate from MeshedHole points to mesh
374 : // sides.
375 376 : std::unordered_map<Point, Point> next_hole_boundary_point;
376 376 : const int np = mh.n_points();
377 7252 : for (auto pi : make_range(1, np))
378 6876 : next_hole_boundary_point[mh.point(pi - 1)] = mh.point(pi);
379 376 : next_hole_boundary_point[mh.point(np - 1)] = mh.point(0);
380 :
381 : #ifndef NDEBUG
382 : int found_hole_sides = 0;
383 : #endif
384 23020 : for (auto elem : hole_mesh.element_ptr_range())
385 : {
386 22644 : if (elem->dim() != 2)
387 0 : mooseError("Non 2-D element found in hole; stitching is not supported.");
388 :
389 22644 : auto ns = elem->n_sides();
390 91460 : for (auto s : make_range(ns))
391 : {
392 68816 : auto it_s = next_hole_boundary_point.find(elem->point(s));
393 68816 : if (it_s != next_hole_boundary_point.end())
394 21320 : if (it_s->second == elem->point((s + 1) % ns))
395 : {
396 7252 : hole_boundary_info.add_side(elem, s, new_hole_bcid);
397 : #ifndef NDEBUG
398 : ++found_hole_sides;
399 : #endif
400 : }
401 : }
402 376 : }
403 : mooseAssert(found_hole_sides == np, "Failed to find full outer boundary of meshed hole");
404 :
405 376 : auto & mesh_boundary_info = mesh->get_boundary_info();
406 : #ifndef NDEBUG
407 : int found_inner_sides = 0;
408 : #endif
409 24319 : for (auto elem : mesh->element_ptr_range())
410 : {
411 23943 : auto ns = elem->n_sides();
412 95799 : for (auto s : make_range(ns))
413 : {
414 71856 : auto it_s = next_hole_boundary_point.find(elem->point((s + 1) % ns));
415 71856 : if (it_s != next_hole_boundary_point.end())
416 22511 : if (it_s->second == elem->point(s))
417 : {
418 7252 : mesh_boundary_info.add_side(elem, s, inner_bcid);
419 : #ifndef NDEBUG
420 : ++found_inner_sides;
421 : #endif
422 : }
423 : }
424 376 : }
425 : mooseAssert(found_inner_sides == np, "Failed to find full boundary around meshed hole");
426 :
427 : // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
428 : // subdomain map
429 376 : const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
430 376 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
431 :
432 : // We do not need the hole_bdy_id_filter anymore
433 396 : for (const auto & bcid : hole_bdy_id_filter)
434 20 : hole_boundary_info.remove_id(bcid);
435 : // If we are stitching a hole boundary layer mesh, we need to reassign the bcid
436 376 : if (hole_bdy_id_filter.size())
437 : {
438 30 : MooseMeshUtils::changeBoundaryId(
439 : hole_mesh,
440 20 : xyd_opts.hole_boundary_inner_id_defaults[hole_i].empty()
441 : ? 1
442 10 : : *xyd_opts.hole_boundary_inner_id_defaults[hole_i].begin(),
443 20 : hole_boundary_rec[hole_i],
444 : true);
445 20 : hole_mesh.get_boundary_info().sideset_name(hole_boundary_rec[hole_i]) =
446 40 : mesh->get_boundary_info().sideset_name(hole_boundary_rec[hole_i]);
447 20 : mesh->get_boundary_info().remove_id(hole_boundary_rec[hole_i]);
448 : }
449 :
450 376 : mesh->stitch_meshes(hole_mesh,
451 : inner_bcid,
452 : new_hole_bcid,
453 : TOLERANCE,
454 : /*clear_stitched_bcids*/ true,
455 376 : xyd_opts.verbose_stitching,
456 376 : xyd_opts.use_binary_search);
457 376 : }
458 : }
459 : // Check if one SubdomainName is shared by more than one subdomain ids
460 713 : std::set<SubdomainName> main_subdomain_map_name_list;
461 1176 : for (auto const & id_name_pair : main_subdomain_map)
462 463 : main_subdomain_map_name_list.emplace(id_name_pair.second);
463 713 : if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
464 6 : mg.paramError("holes", "The hole meshes contain subdomain name maps with conflicts.");
465 :
466 710 : mesh->unset_is_prepared();
467 1420 : return mesh;
468 766 : }
469 : }
|
