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 "BoundaryLayerUtils.h"
11 : #include "MeshTriangulationUtils.h"
12 : #include "MooseMeshUtils.h"
13 : #include "CastUniquePointer.h"
14 : #include "GeometryUtils.h"
15 :
16 : #include "libmesh/elem.h"
17 : #include "libmesh/boundary_info.h"
18 : #include "libmesh/edge_edge3.h"
19 : #include "libmesh/fe_base.h"
20 : #include "libmesh/int_range.h"
21 : #include "libmesh/mesh_serializer.h"
22 : #include "libmesh/mesh_triangle_holes.h"
23 : #include "libmesh/poly2tri_triangulator.h"
24 : #include "libmesh/unstructured_mesh.h"
25 :
26 : using namespace libMesh;
27 :
28 : namespace BoundaryLayerUtils
29 : {
30 :
31 : std::unique_ptr<MeshBase>
32 130 : buildBoundaryLayerRing(MeshGenerator & mg,
33 : MeshBase & input_mesh,
34 : const std::vector<BoundaryName> & boundary_names,
35 : unsigned int num_layers,
36 : Real thickness,
37 : Real layer_bias,
38 : bool outward,
39 : const MooseEnum & tri_elem_type,
40 : SubdomainID output_subdomain_id,
41 : const SubdomainName & output_subdomain_name)
42 : {
43 : // Extract seed boundary polyline points (vertices + optional midpoints) from input_mesh.
44 130 : std::set<std::size_t> bdry_id_set;
45 130 : if (!boundary_names.empty())
46 : {
47 : auto ids =
48 10 : MooseMeshUtils::getBoundaryIDs(input_mesh, boundary_names, /*generate unknown*/ false);
49 10 : bdry_id_set.insert(ids.begin(), ids.end());
50 10 : }
51 130 : TriangulatorInterface::MeshedHole bdry_mh(input_mesh, bdry_id_set);
52 130 : std::vector<Point> cur_pts;
53 130 : std::vector<Point> cur_mids;
54 : // Preserve every input boundary node (including colinear interior nodes on straight edges) so
55 : // the ring's innermost polyline can be stitched back to the input mesh's boundary exactly when
56 : // keep_input is requested by a downstream caller.
57 130 : collectExteriorVertexPointsFromMesh(bdry_mh,
58 : cur_pts,
59 : cur_mids,
60 : /*skip_node_reduction=*/true);
61 :
62 : // Geometric progression of incremental thicknesses. layer_thicknesses[i] is the offset
63 : // distance from polyline i-1 to polyline i (for i >= 1); layer_thicknesses[0] is unused.
64 130 : std::vector<Real> layer_thicknesses(num_layers + 1);
65 : const Real unit_thickness =
66 : (layer_bias == 1.0)
67 130 : ? (thickness / num_layers)
68 60 : : (thickness / (std::pow(layer_bias, num_layers) - 1.0) * (layer_bias - 1.0));
69 130 : layer_thicknesses[0] = 0.0;
70 500 : for (auto i : make_range(std::vector<Real>::size_type(1), layer_thicknesses.size()))
71 370 : layer_thicknesses[i] = unit_thickness * std::pow(layer_bias, i - 1);
72 :
73 : // Generate the N+1 polyline meshes. Polylines are indexed so that polyline 0 is geometrically
74 : // innermost (smallest) and polyline N is outermost (largest). Iteration order depends on
75 : // direction: for outward, build polyline 0 first (= input boundary); for inward, build polyline
76 : // N first.
77 130 : std::vector<std::unique_ptr<MeshBase>> polylines(num_layers + 1);
78 630 : for (auto layer_i : make_range(num_layers + 1))
79 : {
80 500 : const unsigned int layer_index = outward ? layer_i : (num_layers - layer_i);
81 :
82 : // Build the polyline mesh for storage as polylines[layer_index].
83 500 : auto ply = std::make_unique<ReplicatedMesh>(mg.comm());
84 500 : MooseMeshUtils::buildPolyLineMesh(*ply,
85 : cur_pts,
86 : cur_mids,
87 : /*loop=*/true,
88 1000 : BoundaryName(),
89 1000 : BoundaryName(),
90 1000 : std::vector<unsigned int>({1}));
91 500 : polylines[layer_index] = std::move(ply);
92 :
93 500 : if (layer_i + 1 < num_layers + 1)
94 : {
95 : // Build a sacrificial polyline mesh to feed generateOffsetPolyline (it triangulates the
96 : // input mesh in place to compute side normals; we discard it afterwards).
97 370 : auto ply_for_offset = std::make_unique<ReplicatedMesh>(mg.comm());
98 370 : MooseMeshUtils::buildPolyLineMesh(*ply_for_offset,
99 : cur_pts,
100 : cur_mids,
101 : /*loop=*/true,
102 740 : BoundaryName(),
103 740 : BoundaryName(),
104 740 : std::vector<unsigned int>({1}));
105 : std::unique_ptr<UnstructuredMesh> ply_for_offset_u =
106 370 : dynamic_pointer_cast<UnstructuredMesh>(std::move(ply_for_offset));
107 :
108 : std::vector<Point> next_combined = generateOffsetPolyline(
109 370 : &mg, ply_for_offset_u, cur_pts, cur_mids, outward, layer_thicknesses[layer_i + 1]);
110 370 : if (cur_mids.empty())
111 340 : cur_pts = std::move(next_combined);
112 : else
113 : {
114 30 : const auto n_vert = cur_pts.size();
115 30 : cur_pts.assign(next_combined.begin(), next_combined.begin() + n_vert);
116 30 : cur_mids.assign(next_combined.begin() + n_vert, next_combined.end());
117 : }
118 370 : }
119 500 : }
120 :
121 : // Triangulate each annulus (between polylines[i] and polylines[i+1]). Stitch each annulus to
122 : // the accumulating ring (except the first one).
123 130 : std::unique_ptr<MeshBase> ring;
124 500 : for (auto i : make_range(num_layers))
125 : {
126 370 : MeshTriangulationUtils::XYDelaunayOptions xyd_opts;
127 370 : xyd_opts.refine_bdy = false;
128 370 : xyd_opts.verify_holes = false;
129 370 : xyd_opts.stitch_holes = {i > 0};
130 370 : xyd_opts.refine_holes = {false};
131 370 : xyd_opts.tri_elem_type = std::string(tri_elem_type);
132 370 : if (output_subdomain_id != 0)
133 : {
134 310 : xyd_opts.has_output_subdomain_id = true;
135 310 : xyd_opts.output_subdomain_id = output_subdomain_id;
136 : }
137 370 : if (output_subdomain_name.size())
138 : {
139 310 : xyd_opts.has_output_subdomain_name = true;
140 310 : xyd_opts.output_subdomain_name = output_subdomain_name;
141 : }
142 :
143 370 : std::vector<std::unique_ptr<MeshBase>> holes;
144 370 : holes.reserve(1);
145 370 : if (i == 0)
146 130 : holes.push_back(std::move(polylines[0]));
147 : else
148 240 : holes.push_back(std::move(ring));
149 :
150 1110 : ring = MeshTriangulationUtils::triangulateWithDelaunay(
151 1110 : mg, std::move(polylines[i + 1]), std::move(holes), xyd_opts);
152 370 : }
153 : // We now have 2 * num_layers boundaries
154 : // Let's only keep the innermost (1) and outermost (2 * num_layers) boundaries, and remove all
155 : // intermediate ring bcids
156 130 : std::vector<BoundaryID> bids_to_delete;
157 870 : for (const auto b : make_range(2 * num_layers))
158 : {
159 740 : if (b == 1 || b == (num_layers - 1) * 2)
160 260 : continue;
161 480 : bids_to_delete.push_back(b);
162 : }
163 130 : auto & bi = ring->get_boundary_info();
164 610 : for (auto b : bids_to_delete)
165 480 : bi.remove_id(b);
166 :
167 260 : return ring;
168 130 : }
169 :
170 : std::vector<Point>
171 370 : generateOffsetPolyline(MeshGenerator * mg,
172 : std::unique_ptr<libMesh::UnstructuredMesh> & ply_mesh_u,
173 : std::vector<Point> & points,
174 : std::vector<Point> & mid_points,
175 : const bool outward,
176 : const Real thickness)
177 : {
178 : // If the input points are empty, we will extract them from the input 1D mesh
179 370 : if (points.empty())
180 : {
181 : mooseAssert(mid_points.empty(),
182 : "If the input points are empty, the input mid_points must be also empty.");
183 :
184 0 : TriangulatorInterface::MeshedHole bdry_mh(*ply_mesh_u);
185 0 : collectExteriorVertexPointsFromMesh(bdry_mh, points, mid_points);
186 0 : }
187 : // Generate a very simple triangulation mesh so that we can get the outward normal vectors
188 370 : libMesh::Poly2TriTriangulator poly2tri(*ply_mesh_u);
189 370 : poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
190 :
191 370 : poly2tri.set_interpolate_boundary_points(0);
192 370 : poly2tri.set_refine_boundary_allowed(false);
193 370 : poly2tri.set_verify_hole_boundaries(false);
194 370 : poly2tri.desired_area() = 0;
195 370 : poly2tri.minimum_angle() = 0; // Not yet supported
196 370 : poly2tri.smooth_after_generating() = false;
197 370 : if (mid_points.size())
198 30 : poly2tri.elem_type() = libMesh::ElemType::TRI6;
199 370 : poly2tri.triangulate();
200 :
201 : // We need to serialize the mesh for next steps
202 370 : libMesh::MeshSerializer serial(*ply_mesh_u);
203 : // The mesh now only contains one side set that corresponds to the outer boundary with an ID of 0
204 370 : auto bdry_list(ply_mesh_u->get_boundary_info().build_side_list());
205 :
206 : // For each vertex, the shifting direction to form the offset is defined by the normal vectors of
207 : // the two sides that contain the vertex.
208 : // We gather the normals for each node on the boundary here, which are all vertices because of
209 : // the pre-selection of nodes.
210 370 : std::map<dof_id_type, std::vector<Point>> node_normal_map;
211 370 : std::map<dof_id_type, Point> mid_node_normal_map;
212 9210 : for (const auto & bside : bdry_list)
213 : {
214 8840 : const auto & side = ply_mesh_u->elem_ptr(std::get<0>(bside))->side_ptr(std::get<1>(bside));
215 : // For linear elements, the side normal is constant and can be obtained straightforwardly;
216 : // For quadratic elements, we need the normal at the shared vertices
217 : const Point side_normal_0 =
218 8840 : mid_points.size()
219 8840 : ? getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 0)
220 8000 : : ply_mesh_u->elem_ptr(std::get<0>(bside))
221 8000 : ->side_vertex_average_normal(std::get<1>(bside));
222 : const Point side_normal_1 =
223 8840 : mid_points.size()
224 8840 : ? getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 1)
225 8000 : : side_normal_0;
226 :
227 8840 : if (node_normal_map.count(side->node_ptr(0)->id()))
228 4260 : node_normal_map[side->node_ptr(0)->id()].push_back(side_normal_0);
229 : else
230 4580 : node_normal_map[side->node_ptr(0)->id()] = {side_normal_0};
231 8840 : if (node_normal_map.count(side->node_ptr(1)->id()))
232 4580 : node_normal_map[side->node_ptr(1)->id()].push_back(side_normal_1);
233 : else
234 4260 : node_normal_map[side->node_ptr(1)->id()] = {side_normal_1};
235 :
236 8840 : if (mid_points.size())
237 : {
238 : const Point mid_node_normal =
239 840 : getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 2);
240 : mid_node_normal_map
241 840 : [ply_mesh_u->elem_ptr(std::get<0>(bside))->node_ptr(std::get<1>(bside) + 3)->id()] =
242 : mid_node_normal;
243 : }
244 8840 : }
245 :
246 370 : std::vector<Point> mod_reduced_pts_list(points);
247 9210 : for (const auto & [node_id, normal_vecs] : node_normal_map)
248 : {
249 : mooseAssert(normal_vecs.size() == 2,
250 : "Each vertex should be connected to exactly two sides in a polygon.");
251 :
252 8840 : const Point original_pt = *(ply_mesh_u->node_ptr(node_id));
253 : // Form an average normal at the vertex from the two connected sides' normals
254 : const Point move_dir =
255 8840 : (normal_vecs.front() + normal_vecs.back()).unit() * (outward ? 1.0 : -1.0);
256 : // Consider four points of interest to determine the moving distance
257 : // 1. the vertex point
258 : // 2. point along normal_vecs.front() from the vertex with a distance thickness
259 : // 3. point along normal_vecs.back() from the vertex with a distance thickness
260 : // 4. point along move_dir from the vertex with a distance mov_dist
261 : // The four points form a kite shape with its symmetry axis along 1-4 direction
262 : // Angle 1-2-4 and angle 1-3-4 are right angles
263 : // Angle 2-1-3 (i.e., theta) can be calculated using the interior product of
264 : // v1 = normal_vecs.front() and v2 = normal_vecs.back():
265 : // cos(theta) = (v1 . v2) / (|v1| * |v2|)
266 : // Because of the right angles, the distance between 1 and 4 can be calculated as:
267 : // mov_dist = thickness / cos(theta/2)
268 : // and cos(theta/2) = sqrt((1 + cos(theta)) / 2)
269 : const Real mov_dist =
270 8840 : thickness / std::sqrt((1.0 + (normal_vecs.front() * normal_vecs.back()) /
271 8840 : (normal_vecs.front().norm() * normal_vecs.back().norm())) /
272 8840 : 2.0);
273 : mooseAssert(std::count(points.begin(), points.end(), original_pt) == 1,
274 : "The original point should be found exactly once in the reduced points list.");
275 8840 : mod_reduced_pts_list[std::distance(points.begin(),
276 17680 : std::find(points.begin(), points.end(), original_pt))] =
277 17680 : original_pt + move_dir * mov_dist;
278 : }
279 :
280 : // To ensure no overlapping, we need to check set of four points
281 : // p1 and p2 should be a pair of points before and after shifting
282 : // p3 and p4 should be a pair of adjacent shifted points that are neither not p2
283 9210 : for (const auto & i_node_1 : make_range(mod_reduced_pts_list.size()))
284 : {
285 8840 : const Point & p1 = points[i_node_1];
286 8840 : const Point & p2 = mod_reduced_pts_list[i_node_1];
287 239160 : for (const auto & i_node_2 : make_range(mod_reduced_pts_list.size()))
288 : {
289 230320 : if (i_node_2 == i_node_1 || (i_node_2 + 1) % mod_reduced_pts_list.size() == i_node_1)
290 17680 : continue;
291 212640 : const Point & p3 = mod_reduced_pts_list[i_node_2];
292 212640 : const Point & p4 = mod_reduced_pts_list[(i_node_2 + 1) % mod_reduced_pts_list.size()];
293 212640 : if (thickness > 0)
294 212640 : if (geom_utils::segmentsIntersect(p1, p2, p3, p4))
295 0 : mg->paramError(
296 : "thickness",
297 : "The thickness is so large that the mesh is tangled because the offset nodes "
298 : "are no longer in the same order when following the original boundary. Please "
299 : "reduce the thickness value.");
300 : }
301 : }
302 :
303 370 : std::vector<Point> mid_mod_reduced_pts_list(mid_points.size());
304 1210 : for (const auto & [node_id, normal_vec] : mid_node_normal_map)
305 : {
306 840 : const Point original_pt = *(ply_mesh_u->node_ptr(node_id));
307 840 : const Point move_dir = normal_vec.unit() * (outward ? 1.0 : -1.0);
308 840 : mid_mod_reduced_pts_list[std::distance(
309 1680 : mid_points.begin(), std::find(mid_points.begin(), mid_points.end(), original_pt))] =
310 1680 : original_pt + move_dir * thickness;
311 : }
312 :
313 : // combine mod_reduced_pts_list and mid_mod_reduced_pts_list to get the final list of points for
314 : // the layer mesh
315 370 : std::vector<Point> layer_pts_list(mod_reduced_pts_list);
316 740 : layer_pts_list.insert(
317 370 : layer_pts_list.end(), mid_mod_reduced_pts_list.begin(), mid_mod_reduced_pts_list.end());
318 :
319 740 : return layer_pts_list;
320 370 : }
321 :
322 : void
323 130 : collectExteriorVertexPointsFromMesh(libMesh::TriangulatorInterface::MeshedHole & bdry_mh,
324 : std::vector<Point> & points,
325 : std::vector<Point> & mid_points,
326 : const bool skip_node_reduction)
327 : {
328 3170 : for (const auto i : make_range(bdry_mh.n_points()))
329 : {
330 3040 : if (skip_node_reduction || !geom_utils::arePointsColinear(
331 0 : bdry_mh.point((i - 1 + bdry_mh.n_points()) % bdry_mh.n_points()),
332 0 : bdry_mh.point(i),
333 3040 : bdry_mh.point((i + 1) % bdry_mh.n_points())))
334 : {
335 3040 : points.push_back(bdry_mh.point(i));
336 3040 : if (bdry_mh.n_midpoints() == 1 && skip_node_reduction)
337 280 : mid_points.push_back(bdry_mh.midpoint(0, i));
338 : }
339 : }
340 130 : }
341 :
342 : Point
343 2520 : getKeyNormal(const Elem * elem, const unsigned int s, const unsigned int node_index)
344 : {
345 2520 : const std::unique_ptr<const Elem> face = elem->build_side_ptr(s);
346 : mooseAssert(face->type() == ElemType::EDGE3,
347 : "Only elements with EDGE3 sides are supported in this function.");
348 : mooseAssert(node_index < 3,
349 : "The node index for an EDGE3 side should be 0, 1, or 2 (for the two vertices and the "
350 : "midpoint).");
351 : std::unique_ptr<libMesh::FEBase> fe(
352 2520 : libMesh::FEBase::build(2, libMesh::FEType(elem->default_order())));
353 2520 : const std::vector<Point> & normals = fe->get_normals();
354 5040 : std::vector<Point> ref_pts = {face->reference_elem()->point(node_index)};
355 2520 : fe->reinit(elem, s, TOLERANCE, &ref_pts);
356 5040 : return normals[0];
357 2520 : }
358 : }
|
