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 : #include "MortarSegmentHelper.h"
10 : #include "MooseError.h"
11 :
12 : #include "libmesh/int_range.h"
13 : #include "libmesh/utility.h"
14 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
15 : #include "libmesh/replicated_mesh.h"
16 : #include "libmesh/mesh_triangle_interface.h"
17 : #include "libmesh/poly2tri_triangulator.h"
18 : #endif
19 :
20 : #include <algorithm>
21 : #include <array>
22 : #include <cmath>
23 : #include <limits>
24 : #include <map>
25 : #include <numeric>
26 : #include <optional>
27 : #include <set>
28 : #include <unordered_map>
29 :
30 : using namespace libMesh;
31 :
32 : namespace
33 : {
34 :
35 : // Signed-area test for the 2D triangle (a, b, c). Returns twice the signed area:
36 : // positive if a->b->c is counter-clockwise, negative if clockwise, zero if
37 : // collinear. Used as the building block for orientation, point-in-triangle, and
38 : // circumcircle predicates.
39 : Real
40 171081 : orient2dHelper(const Point & a, const Point & b, const Point & c)
41 : {
42 171081 : return (b(0) - a(0)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(0) - a(0));
43 : }
44 :
45 : Real
46 68814 : triangleAreaHelper(const Point & a, const Point & b, const Point & c)
47 : {
48 68814 : return 0.5 * std::abs(orient2dHelper(a, b, c));
49 : }
50 :
51 : // Canonical key for an undirected edge: the two endpoint indices sorted so that
52 : // (a, b) and (b, a) hash and compare equal. Used to dedupe / look up edges in
53 : // triangle-adjacency maps.
54 : std::array<unsigned int, 2>
55 55629 : canonicalEdgeHelper(const unsigned int a, const unsigned int b)
56 : {
57 55629 : return {{std::min(a, b), std::max(a, b)}};
58 : }
59 :
60 : // Reorder the three vertex indices (a, b, c) so the resulting triangle is wound
61 : // counter-clockwise (CCW) in the 2D plane spanned by \p nodes. Many of the
62 : // triangulation paths (orientation tests, area accumulation, ear-clipping
63 : // validity checks) assume CCW input, so we normalize before emitting triangles.
64 : std::array<unsigned int, 3>
65 6669 : makeCCWTriangleHelper(const std::vector<Point> & nodes,
66 : const unsigned int a,
67 : const unsigned int b,
68 : const unsigned int c)
69 : {
70 6669 : if (orient2dHelper(nodes[a], nodes[b], nodes[c]) >= 0)
71 5661 : return {{a, b, c}};
72 1008 : return {{a, c, b}};
73 : }
74 :
75 : bool
76 6075 : pointInCircumcircleHelper(const Point & a, const Point & b, const Point & c, const Point & p)
77 : {
78 6075 : const auto ax = a(0) - p(0);
79 6075 : const auto ay = a(1) - p(1);
80 6075 : const auto bx = b(0) - p(0);
81 6075 : const auto by = b(1) - p(1);
82 6075 : const auto cx = c(0) - p(0);
83 6075 : const auto cy = c(1) - p(1);
84 6075 : const Real det = (ax * ax + ay * ay) * (bx * cy - by * cx) -
85 6075 : (bx * bx + by * by) * (ax * cy - ay * cx) +
86 6075 : (cx * cx + cy * cy) * (ax * by - ay * bx);
87 6075 : const Real orientation = orient2dHelper(a, b, c);
88 6075 : return orientation >= 0 ? det > TOLERANCE : det < -TOLERANCE;
89 : }
90 :
91 : void
92 4410 : performLocalDelaunayFlips(const std::vector<Point> & poly_nodes,
93 : const std::set<std::array<unsigned int, 2>> & constrained_edges,
94 : std::vector<std::array<unsigned int, 3>> & triangles)
95 : {
96 4410 : bool flipped = true;
97 9720 : while (flipped)
98 : {
99 5310 : flipped = false;
100 :
101 5310 : std::map<std::array<unsigned int, 2>, std::vector<unsigned int>> edge_to_triangles;
102 17667 : for (const auto tri_index : index_range(triangles))
103 : {
104 12357 : const auto & tri = triangles[tri_index];
105 12357 : edge_to_triangles[canonicalEdgeHelper(tri[0], tri[1])].push_back(tri_index);
106 12357 : edge_to_triangles[canonicalEdgeHelper(tri[1], tri[2])].push_back(tri_index);
107 12357 : edge_to_triangles[canonicalEdgeHelper(tri[2], tri[0])].push_back(tri_index);
108 : }
109 :
110 32022 : for (const auto & [edge, owning_triangles] : edge_to_triangles)
111 : {
112 27612 : if (owning_triangles.size() != 2 || constrained_edges.count(edge))
113 20844 : continue;
114 :
115 6768 : const auto first_tri_index = owning_triangles[0];
116 6768 : const auto second_tri_index = owning_triangles[1];
117 6768 : const auto & first_triangle = triangles[first_tri_index];
118 6768 : const auto & second_triangle = triangles[second_tri_index];
119 :
120 6768 : const auto a = edge[0];
121 6768 : const auto b = edge[1];
122 : const auto first_opposite =
123 6768 : *std::find_if(first_triangle.begin(),
124 : first_triangle.end(),
125 12069 : [a, b](const unsigned int vertex) { return vertex != a && vertex != b; });
126 : const auto second_opposite =
127 6768 : *std::find_if(second_triangle.begin(),
128 : second_triangle.end(),
129 14904 : [a, b](const unsigned int vertex) { return vertex != a && vertex != b; });
130 :
131 6768 : if (first_opposite == second_opposite)
132 0 : continue;
133 :
134 : const auto side_a =
135 6768 : orient2dHelper(poly_nodes[first_opposite], poly_nodes[second_opposite], poly_nodes[a]);
136 : const auto side_b =
137 6768 : orient2dHelper(poly_nodes[first_opposite], poly_nodes[second_opposite], poly_nodes[b]);
138 6768 : if (side_a * side_b >= -TOLERANCE)
139 693 : continue;
140 :
141 6075 : if (!pointInCircumcircleHelper(poly_nodes[first_triangle[0]],
142 6075 : poly_nodes[first_triangle[1]],
143 6075 : poly_nodes[first_triangle[2]],
144 6075 : poly_nodes[second_opposite]))
145 5175 : continue;
146 :
147 900 : triangles[first_tri_index] =
148 900 : makeCCWTriangleHelper(poly_nodes, first_opposite, second_opposite, b);
149 900 : triangles[second_tri_index] =
150 900 : makeCCWTriangleHelper(poly_nodes, second_opposite, first_opposite, a);
151 900 : flipped = true;
152 900 : break;
153 : }
154 5310 : }
155 4410 : }
156 :
157 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
158 : void
159 2205 : triangulateConstrainedDelaunayPolygon(std::vector<Point> & poly_nodes,
160 : const Real area_tol,
161 : const Real length_tol,
162 : std::vector<std::vector<unsigned int>> & tri_map)
163 : {
164 2205 : Parallel::Communicator comm_self;
165 2205 : ReplicatedMesh triangulation_mesh(comm_self, 2);
166 2205 : std::unordered_map<dof_id_type, unsigned int> node_id_to_local_index;
167 2205 : node_id_to_local_index.reserve(poly_nodes.size());
168 :
169 11484 : for (const auto i : index_range(poly_nodes))
170 9279 : triangulation_mesh.add_point(poly_nodes[i], i);
171 :
172 2205 : triangulation_mesh.set_mesh_dimension(2);
173 :
174 : #ifdef LIBMESH_HAVE_TRIANGLE
175 : TriangleInterface triangulator(triangulation_mesh);
176 : #else
177 2205 : Poly2TriTriangulator triangulator(triangulation_mesh);
178 2205 : triangulator.set_refine_boundary_allowed(false);
179 : #endif
180 :
181 2205 : triangulator.triangulation_type() = TriangulatorInterface::PSLG;
182 2205 : triangulator.elem_type() = TRI3;
183 2205 : triangulator.set_interpolate_boundary_points(0);
184 2205 : triangulator.set_verify_hole_boundaries(false);
185 2205 : triangulator.desired_area() = 0;
186 2205 : triangulator.minimum_angle() = 0;
187 2205 : triangulator.smooth_after_generating() = false;
188 2205 : triangulator.quiet() = true;
189 2205 : triangulator.segments.reserve(poly_nodes.size());
190 11484 : for (const auto i : index_range(poly_nodes))
191 9279 : triangulator.segments.emplace_back(i, (i + 1) % poly_nodes.size());
192 :
193 2205 : triangulator.triangulate();
194 :
195 : // node_ptr_range() and active_element_ptr_range() iterate in id order on this
196 : // serial ReplicatedMesh, so no explicit sort is needed.
197 11484 : for (const auto * const node : triangulation_mesh.node_ptr_range())
198 9279 : if (!node_id_to_local_index.count(node->id()))
199 : {
200 : // Node inherits from Point and the triangulator operates on a 2D plane, so
201 : // the libMesh node already lives at z = 0 and we can use it directly.
202 9279 : unsigned int matched_index = libMesh::invalid_uint;
203 9279 : Real best_distance = std::numeric_limits<Real>::max();
204 :
205 48690 : for (const auto i : index_range(poly_nodes))
206 : {
207 39411 : const Real distance = (*node - poly_nodes[i]).norm();
208 39411 : if (distance <= length_tol && distance < best_distance)
209 : {
210 9279 : matched_index = i;
211 9279 : best_distance = distance;
212 : }
213 : }
214 :
215 9279 : if (matched_index == libMesh::invalid_uint)
216 : {
217 0 : matched_index = cast_int<unsigned int>(poly_nodes.size());
218 0 : poly_nodes.push_back(*node);
219 : }
220 :
221 9279 : node_id_to_local_index.emplace(node->id(), matched_index);
222 2205 : }
223 :
224 2205 : std::vector<std::array<unsigned int, 3>> triangles;
225 2205 : triangles.reserve(triangulation_mesh.n_elem());
226 :
227 7074 : for (const auto * const elem : triangulation_mesh.active_element_ptr_range())
228 : {
229 : mooseAssert(elem->type() == TRI3,
230 : "The delaunay mortar triangulation backend produced a non-TRI3 element: "
231 : << static_cast<int>(elem->type()));
232 :
233 : std::array<unsigned int, 3> local_triangle;
234 19476 : for (const auto i : index_range(local_triangle))
235 14607 : local_triangle[i] = libmesh_map_find(node_id_to_local_index, elem->node_id(i));
236 :
237 4869 : const Real orientation = orient2dHelper(poly_nodes[local_triangle[0]],
238 4869 : poly_nodes[local_triangle[1]],
239 4869 : poly_nodes[local_triangle[2]]);
240 4869 : if (std::abs(orientation) <= 2. * area_tol)
241 0 : continue;
242 :
243 4869 : if (orientation < 0)
244 0 : std::swap(local_triangle[1], local_triangle[2]);
245 :
246 4869 : triangles.push_back(local_triangle);
247 2205 : }
248 :
249 2205 : std::set<std::array<unsigned int, 2>> constrained_edges;
250 11484 : for (const auto i : index_range(poly_nodes))
251 9279 : constrained_edges.insert(canonicalEdgeHelper(i, (i + 1) % poly_nodes.size()));
252 :
253 2205 : performLocalDelaunayFlips(poly_nodes, constrained_edges, triangles);
254 :
255 2205 : std::set<std::array<unsigned int, 3>> seen_triangles;
256 7074 : for (auto local_triangle : triangles)
257 : {
258 4869 : auto canonical_triangle = local_triangle;
259 4869 : std::sort(canonical_triangle.begin(), canonical_triangle.end());
260 4869 : if (!seen_triangles.insert(canonical_triangle).second)
261 0 : continue;
262 :
263 14607 : tri_map.push_back({local_triangle[0], local_triangle[1], local_triangle[2]});
264 : }
265 2205 : }
266 : #endif
267 :
268 : } // namespace
269 :
270 10134 : MortarSegmentHelper::MortarSegmentHelper(const std::vector<Point> secondary_nodes,
271 : const Point & center,
272 : const Point & normal,
273 : const MortarSegmentTriangulationMode triangulation_mode,
274 10134 : const bool triangulate_triangles)
275 10134 : : _center(center),
276 10134 : _normal(normal),
277 10134 : _debug(false),
278 10134 : _triangulation_mode(triangulation_mode),
279 10134 : _triangulate_triangles(triangulate_triangles)
280 : {
281 10134 : _secondary_poly.clear();
282 10134 : _secondary_poly.reserve(secondary_nodes.size());
283 :
284 : // Get orientation of secondary poly
285 10134 : const Point e1 = secondary_nodes[0] - secondary_nodes[1];
286 10134 : const Point e2 = secondary_nodes[2] - secondary_nodes[1];
287 10134 : const Real orient = e2.cross(e1) * _normal;
288 :
289 : // u and v define the tangent plane of the element (at center)
290 : // Note we embed orientation into our transformation to make 2D poly always
291 : // positively oriented
292 10134 : _u = _normal.cross(secondary_nodes[0] - center).unit();
293 10134 : _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
294 :
295 : // Transform problem to 2D plane spanned by u and v
296 45198 : for (const auto & node : secondary_nodes)
297 : {
298 35064 : Point pt = node - _center;
299 35064 : _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
300 : }
301 :
302 : // Initialize area of secondary polygon
303 10134 : _remaining_area_fraction = 1.0;
304 10134 : _secondary_area = area(_secondary_poly);
305 :
306 : // Tolerance for quantities with area dimensions
307 10134 : _area_tol = _tolerance * _secondary_area;
308 :
309 : // Tolerance for quantites with length dimensions
310 10134 : _length_tol = _tolerance * std::sqrt(_secondary_area);
311 10134 : }
312 :
313 : Point
314 185758 : MortarSegmentHelper::getIntersection(
315 : const Point & p1, const Point & p2, const Point & q1, const Point & q2, Real & s) const
316 : {
317 185758 : const Point dp = p2 - p1;
318 185758 : const Point dq = q2 - q1;
319 185758 : const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
320 185758 : const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
321 185758 : const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
322 185758 : const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
323 185758 : s = -alpha * (cp1q2 - cp1q1 - cq1q2);
324 :
325 : // Intersection should be between p1 and p2, if it's not (due to poor conditioning), simply
326 : // move it to one of the end points
327 185758 : s = s > 1 ? 1. : s;
328 185758 : s = s < 0 ? 0. : s;
329 185758 : return p1 + s * dp;
330 : }
331 :
332 : bool
333 0 : MortarSegmentHelper::isInsideSecondary(const Point & pt) const
334 : {
335 0 : for (auto i : index_range(_secondary_poly))
336 : {
337 0 : const Point & q1 = _secondary_poly[i];
338 0 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
339 :
340 0 : const Point e1 = q2 - q1;
341 0 : const Point e2 = pt - q1;
342 :
343 : // If point corresponds to one of the secondary vertices, skip
344 0 : if (e2.norm() < _tolerance)
345 0 : return true;
346 :
347 0 : const bool inside = (e1(0) * e2(1) - e1(1) * e2(0)) < _area_tol;
348 0 : if (!inside)
349 0 : return false;
350 : }
351 0 : return true;
352 : }
353 :
354 : bool
355 340981 : MortarSegmentHelper::isDisjoint(const std::vector<Point> & poly) const
356 : {
357 803443 : for (auto i : index_range(_secondary_poly))
358 : {
359 : // Get edge to check
360 752981 : const Point & q1 = _secondary_poly[i];
361 752981 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
362 752981 : const Point edg = q2 - q1;
363 752981 : const Real cp = q2(0) * q1(1) - q2(1) * q1(0);
364 :
365 : // If more optimization needed, could store these values for later
366 : // Check if point is to the left of (or on) clip_edge
367 2536080 : auto is_inside = [&edg, cp](Point & pt, Real tol)
368 2536080 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
369 :
370 752981 : bool all_outside = true;
371 3289061 : for (auto pt : poly)
372 2536080 : if (is_inside(pt, _area_tol))
373 1288835 : all_outside = false;
374 :
375 752981 : if (all_outside)
376 290519 : return true;
377 : }
378 50462 : return false;
379 : }
380 :
381 : std::vector<Point>
382 340981 : MortarSegmentHelper::projectPrimaryPoly(const std::vector<Point> & primary_nodes) const
383 : {
384 : // Check orientation of primary_poly
385 340981 : const Point e1 = primary_nodes[0] - primary_nodes[1];
386 340981 : const Point e2 = primary_nodes[2] - primary_nodes[1];
387 :
388 : // Note we use u x v here instead of normal because it may be flipped if secondary elem was
389 : // negatively oriented
390 340981 : const Real orient = e2.cross(e1) * _u.cross(_v);
391 :
392 : // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
393 340981 : std::vector<Point> primary_poly;
394 340981 : const int n_verts = primary_nodes.size();
395 340981 : primary_poly.reserve(primary_nodes.size());
396 1474041 : for (auto n : index_range(primary_nodes))
397 : {
398 1133060 : Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
399 1133060 : primary_poly.emplace_back(pt * _u, pt * _v, 0.);
400 : }
401 :
402 681962 : return primary_poly;
403 0 : }
404 :
405 : std::vector<Point>
406 340981 : MortarSegmentHelper::clipPoly(const std::vector<Point> & primary_nodes) const
407 : {
408 340981 : std::vector<Point> primary_poly = projectPrimaryPoly(primary_nodes);
409 :
410 340981 : if (isDisjoint(primary_poly))
411 : {
412 290519 : primary_poly.clear();
413 290519 : return primary_poly;
414 : }
415 :
416 : // Initialize clipped poly with secondary poly (secondary is target poly)
417 50462 : std::vector<Point> clipped_poly = _secondary_poly;
418 :
419 : // Loop through clipping edges
420 219947 : for (auto i : index_range(primary_poly))
421 : {
422 : // If clipped poly trivial, return
423 174191 : if (clipped_poly.size() < 3)
424 : {
425 4706 : clipped_poly.clear();
426 4706 : return clipped_poly;
427 : }
428 :
429 : // Set input poly to current clipped poly
430 169485 : std::vector<Point> input_poly(clipped_poly);
431 169485 : clipped_poly.clear();
432 :
433 : // Get clipping edge
434 169485 : const Point & clip_pt1 = primary_poly[i];
435 169485 : const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
436 169485 : const Point edg = clip_pt2 - clip_pt1;
437 169485 : const Real cp = clip_pt2(0) * clip_pt1(1) - clip_pt2(1) * clip_pt1(0);
438 :
439 : // Check if point is to the left of (or on) clip_edge
440 : /*
441 : * Note that use of tolerance here is to avoid degenerate case when lines are
442 : * essentially on top of each other (common when meshes match across interface)
443 : * since finding intersection is ill-conditioned in this case.
444 : */
445 1264296 : auto is_inside = [&edg, cp](const Point & pt, Real tol)
446 1264296 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < tol; };
447 :
448 : // Loop through edges of target polygon (with previous clippings already included)
449 801633 : for (auto j : index_range(input_poly))
450 : {
451 : // Get target edge
452 632148 : const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
453 632148 : const Point prev_pt = input_poly[j];
454 :
455 : // TODO: Don't need to calculate both each loop
456 632148 : const bool is_current_inside = is_inside(curr_pt, _area_tol);
457 632148 : const bool is_previous_inside = is_inside(prev_pt, _area_tol);
458 :
459 632148 : if (is_current_inside)
460 : {
461 437184 : if (!is_previous_inside)
462 : {
463 : Real s;
464 92879 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
465 :
466 : /*
467 : * s is the fraction of distance along clip poly edge that intersection lies
468 : * It is used here to avoid degenerate polygon cases. For example, consider a
469 : * case like:
470 : * o
471 : * | (inside)
472 : * ------|------
473 : * | (outside)
474 : * when the distance is small (< 1e-7) we don't want to to add both the point
475 : * and intersection. Also note that when distance on the scale of 1e-7,
476 : * area on scale of 1e-14 so is insignificant if this results in dropping
477 : * a tri (for example if next edge crosses again)
478 : */
479 92879 : if (s < (1 - _tolerance))
480 91800 : clipped_poly.push_back(intersect);
481 : }
482 437184 : clipped_poly.push_back(curr_pt);
483 : }
484 194964 : else if (is_previous_inside)
485 : {
486 : Real s;
487 92879 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
488 92879 : if (s > _tolerance)
489 91665 : clipped_poly.push_back(intersect);
490 : }
491 : }
492 169485 : }
493 :
494 : // Make sure final clipped poly is not trivial
495 45756 : if (clipped_poly.size() < 3)
496 : {
497 1674 : clipped_poly.clear();
498 1674 : return clipped_poly;
499 : }
500 :
501 : // Clean up result by removing any duplicate nodes
502 44082 : std::vector<Point> cleaned_poly;
503 44082 : cleaned_poly.push_back(clipped_poly.back());
504 167418 : for (auto i : make_range(clipped_poly.size() - 1))
505 : {
506 123336 : const Point prev_pt = cleaned_poly.back();
507 123336 : const Point curr_pt = clipped_poly[i];
508 :
509 : // If points are sufficiently distanced, add to output
510 123336 : if ((curr_pt - prev_pt).norm() > _length_tol)
511 123336 : cleaned_poly.push_back(curr_pt);
512 : }
513 :
514 : mooseAssert(
515 : cleaned_poly.size() <= 8,
516 : "Our distributed mesh numbering scheme assumes that we have at most 8 nodes resulting from "
517 : "clipping the projection of the primary sub-element onto the secondary sub-element");
518 44082 : return cleaned_poly;
519 340981 : }
520 :
521 : void
522 44082 : MortarSegmentHelper::triangulatePoly(std::vector<Point> & poly_nodes,
523 : std::vector<std::vector<unsigned int>> & tri_map) const
524 : {
525 : // tri_map is populated with triangle indices that are local to poly_nodes (starting at 0).
526 : // Callers are responsible for shifting these indices into a global node numbering.
527 2466 : const auto polygon_centroid = [](const std::vector<Point> & polygon_nodes)
528 : {
529 2466 : Point centroid(0);
530 2466 : Real double_area = 0;
531 12528 : for (const auto i : index_range(polygon_nodes))
532 : {
533 10062 : const auto & a = polygon_nodes[i];
534 10062 : const auto & b = polygon_nodes[(i + 1) % polygon_nodes.size()];
535 10062 : const Real cross = a(0) * b(1) - b(0) * a(1);
536 10062 : double_area += cross;
537 10062 : centroid(0) += (a(0) + b(0)) * cross;
538 10062 : centroid(1) += (a(1) + b(1)) * cross;
539 : }
540 :
541 2466 : if (std::abs(double_area) <= TOLERANCE)
542 : {
543 36 : for (const auto & node : polygon_nodes)
544 27 : centroid += node;
545 9 : centroid /= polygon_nodes.size();
546 9 : return centroid;
547 : }
548 :
549 2457 : centroid /= (3. * double_area);
550 2457 : centroid(2) = 0;
551 2457 : return centroid;
552 : };
553 :
554 10521 : const auto append_triangle = [this, &poly_nodes, &tri_map](
555 : const unsigned int a, const unsigned int b, const unsigned int c)
556 : {
557 10521 : if (triangleAreaHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) <= _area_tol)
558 27 : return false;
559 :
560 10494 : if (orient2dHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) >= 0)
561 31482 : tri_map.push_back({a, b, c});
562 : else
563 0 : tri_map.push_back({a, c, b});
564 :
565 10494 : return true;
566 44082 : };
567 :
568 : const auto point_in_triangle =
569 13410 : [this](const Point & p, const Point & a, const Point & b, const Point & c)
570 : {
571 13410 : const Real ab = orient2dHelper(a, b, p);
572 13410 : const Real bc = orient2dHelper(b, c, p);
573 13410 : const Real ca = orient2dHelper(c, a, p);
574 13410 : return ab >= -_area_tol && bc >= -_area_tol && ca >= -_area_tol;
575 44082 : };
576 :
577 11115 : const auto min_triangle_angle = [](const Point & a, const Point & b, const Point & c)
578 : {
579 33309 : const auto clamp_cos = [](Real value) { return std::max(-1., std::min(1., value)); };
580 : const auto angle_at =
581 33345 : [&clamp_cos](const Point & vertex, const Point & point_one, const Point & point_two)
582 : {
583 33345 : const Point edge_one = point_one - vertex;
584 33345 : const Point edge_two = point_two - vertex;
585 33345 : const Real denom = edge_one.norm() * edge_two.norm();
586 33345 : if (denom <= TOLERANCE)
587 36 : return 0.;
588 33309 : return std::acos(clamp_cos((edge_one * edge_two) / denom));
589 11115 : };
590 :
591 11115 : return std::min({angle_at(a, b, c), angle_at(b, c, a), angle_at(c, a, b)});
592 : };
593 :
594 9864 : const auto canonicalize_polygon = [this, &poly_nodes]()
595 : {
596 9864 : if (poly_nodes.size() < 3)
597 0 : return;
598 :
599 9864 : if (area(poly_nodes) < 0)
600 0 : std::reverse(poly_nodes.begin(), poly_nodes.end());
601 :
602 9864 : bool changed = true;
603 18684 : while (changed && poly_nodes.size() > 3)
604 : {
605 8820 : changed = false;
606 45936 : for (const auto i : index_range(poly_nodes))
607 : {
608 37116 : const auto prev = (i + poly_nodes.size() - 1) % poly_nodes.size();
609 37116 : const auto next = (i + 1) % poly_nodes.size();
610 37116 : if ((poly_nodes[i] - poly_nodes[prev]).norm() <= _length_tol ||
611 74232 : (poly_nodes[next] - poly_nodes[i]).norm() <= _length_tol ||
612 37116 : triangleAreaHelper(poly_nodes[prev], poly_nodes[i], poly_nodes[next]) <= _area_tol)
613 : {
614 0 : poly_nodes.erase(poly_nodes.begin() + i);
615 0 : changed = true;
616 0 : break;
617 : }
618 : }
619 : }
620 :
621 9864 : if (poly_nodes.size() >= 3 && area(poly_nodes) < 0)
622 0 : std::reverse(poly_nodes.begin(), poly_nodes.end());
623 44082 : };
624 :
625 : const auto triangulate_with_ear_clipping =
626 2205 : [this, &poly_nodes, &point_in_triangle, &min_triangle_angle](
627 : const bool perform_delaunay_flips)
628 : {
629 2205 : std::vector<std::array<unsigned int, 3>> triangles;
630 2205 : if (poly_nodes.size() < 3)
631 0 : return triangles;
632 :
633 2205 : if (poly_nodes.size() == 3)
634 : {
635 0 : triangles.push_back(makeCCWTriangleHelper(poly_nodes, 0, 1, 2));
636 0 : return triangles;
637 : }
638 :
639 2205 : std::vector<unsigned int> remaining_vertices(poly_nodes.size());
640 2205 : std::iota(remaining_vertices.begin(), remaining_vertices.end(), 0);
641 :
642 4869 : while (remaining_vertices.size() > 3)
643 : {
644 2664 : std::optional<std::size_t> best_position;
645 2664 : Real best_score = -std::numeric_limits<Real>::max();
646 2664 : Real best_area = -std::numeric_limits<Real>::max();
647 :
648 13779 : for (const auto position : index_range(remaining_vertices))
649 : {
650 : const auto prev_position =
651 11115 : (position + remaining_vertices.size() - 1) % remaining_vertices.size();
652 11115 : const auto next_position = (position + 1) % remaining_vertices.size();
653 11115 : const auto prev = remaining_vertices[prev_position];
654 11115 : const auto curr = remaining_vertices[position];
655 11115 : const auto next = remaining_vertices[next_position];
656 :
657 11115 : if (orient2dHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]) <= _area_tol)
658 0 : continue;
659 :
660 11115 : bool contains_other_vertex = false;
661 57870 : for (const auto other : remaining_vertices)
662 : {
663 46755 : if (other == prev || other == curr || other == next)
664 33345 : continue;
665 :
666 13410 : if (point_in_triangle(
667 13410 : poly_nodes[other], poly_nodes[prev], poly_nodes[curr], poly_nodes[next]))
668 : {
669 0 : contains_other_vertex = true;
670 0 : break;
671 : }
672 : }
673 :
674 11115 : if (contains_other_vertex)
675 0 : continue;
676 :
677 : const Real candidate_score =
678 11115 : min_triangle_angle(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
679 : const Real candidate_area =
680 11115 : triangleAreaHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
681 16209 : if (!best_position || candidate_score > best_score + TOLERANCE ||
682 5031 : (std::abs(candidate_score - best_score) <= TOLERANCE &&
683 63 : candidate_area > best_area + _area_tol))
684 : {
685 6120 : best_position = position;
686 6120 : best_score = candidate_score;
687 6120 : best_area = candidate_area;
688 : }
689 : }
690 :
691 2664 : if (!best_position)
692 : {
693 0 : std::vector<std::array<unsigned int, 3>> best_fan;
694 0 : Real best_fan_score = -std::numeric_limits<Real>::max();
695 0 : Real best_fan_area = -std::numeric_limits<Real>::max();
696 :
697 0 : for (const auto root_position : index_range(remaining_vertices))
698 : {
699 0 : std::vector<std::array<unsigned int, 3>> candidate_fan;
700 0 : Real candidate_score = std::numeric_limits<Real>::max();
701 0 : Real candidate_area = std::numeric_limits<Real>::max();
702 0 : bool valid_fan = true;
703 0 : const auto root = remaining_vertices[root_position];
704 :
705 0 : for (unsigned int step = 1; step + 1 < remaining_vertices.size(); ++step)
706 : {
707 0 : const auto next_position = (root_position + step) % remaining_vertices.size();
708 0 : const auto following_position = (root_position + step + 1) % remaining_vertices.size();
709 0 : const auto vertex_one = remaining_vertices[next_position];
710 0 : const auto vertex_two = remaining_vertices[following_position];
711 :
712 0 : if (orient2dHelper(poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]) <=
713 0 : _area_tol)
714 : {
715 0 : valid_fan = false;
716 0 : break;
717 : }
718 :
719 0 : candidate_fan.push_back(
720 0 : makeCCWTriangleHelper(poly_nodes, root, vertex_one, vertex_two));
721 0 : candidate_score =
722 0 : std::min(candidate_score,
723 0 : min_triangle_angle(
724 0 : poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
725 0 : candidate_area =
726 0 : std::min(candidate_area,
727 0 : triangleAreaHelper(
728 0 : poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
729 : }
730 :
731 0 : if (!valid_fan || candidate_fan.empty())
732 0 : continue;
733 :
734 0 : if (candidate_score > best_fan_score + TOLERANCE ||
735 0 : (std::abs(candidate_score - best_fan_score) <= TOLERANCE &&
736 0 : candidate_area > best_fan_area + _area_tol))
737 : {
738 0 : best_fan = std::move(candidate_fan);
739 0 : best_fan_score = candidate_score;
740 0 : best_fan_area = candidate_area;
741 : }
742 0 : }
743 :
744 0 : if (best_fan.empty())
745 0 : for (unsigned int i = 1; i + 1 < remaining_vertices.size(); ++i)
746 0 : best_fan.push_back(makeCCWTriangleHelper(poly_nodes,
747 0 : remaining_vertices[0],
748 0 : remaining_vertices[i],
749 0 : remaining_vertices[i + 1]));
750 :
751 0 : triangles.insert(triangles.end(), best_fan.begin(), best_fan.end());
752 0 : break;
753 0 : }
754 :
755 : const auto prev_position =
756 2664 : (*best_position + remaining_vertices.size() - 1) % remaining_vertices.size();
757 2664 : const auto next_position = (*best_position + 1) % remaining_vertices.size();
758 2664 : triangles.push_back(makeCCWTriangleHelper(poly_nodes,
759 2664 : remaining_vertices[prev_position],
760 2664 : remaining_vertices[*best_position],
761 2664 : remaining_vertices[next_position]));
762 2664 : remaining_vertices.erase(remaining_vertices.begin() + *best_position);
763 : }
764 :
765 2205 : if (remaining_vertices.size() == 3)
766 2205 : triangles.push_back(makeCCWTriangleHelper(
767 2205 : poly_nodes, remaining_vertices[0], remaining_vertices[1], remaining_vertices[2]));
768 :
769 2205 : if (!perform_delaunay_flips)
770 0 : return triangles;
771 :
772 2205 : std::set<std::array<unsigned int, 2>> boundary_edges;
773 11484 : for (const auto i : index_range(poly_nodes))
774 9279 : boundary_edges.insert(canonicalEdgeHelper(i, (i + 1) % poly_nodes.size()));
775 :
776 2205 : performLocalDelaunayFlips(poly_nodes, boundary_edges, triangles);
777 2205 : return triangles;
778 2205 : };
779 :
780 2205 : const auto is_convex_polygon = [this](const std::vector<Point> & polygon_nodes)
781 : {
782 2205 : if (polygon_nodes.size() <= 3)
783 0 : return true;
784 :
785 11484 : for (const auto i : index_range(polygon_nodes))
786 : {
787 9279 : const auto prev = (i + polygon_nodes.size() - 1) % polygon_nodes.size();
788 9279 : const auto next = (i + 1) % polygon_nodes.size();
789 9279 : if (orient2dHelper(polygon_nodes[prev], polygon_nodes[i], polygon_nodes[next]) <= _area_tol)
790 0 : return false;
791 : }
792 :
793 2205 : return true;
794 44082 : };
795 :
796 : // Fewer than 3 nodes can't be triangulated
797 44082 : if (poly_nodes.size() < 3)
798 0 : mooseError("Can't triangulate poly with fewer than 3 nodes");
799 :
800 : // Legacy centroid path: when the default triangulation (centroid) is selected
801 : // and triangle re-tessellation is not requested, reproduce the legacy
802 : // algorithm byte-for-byte so existing mortar baselines remain valid.
803 : // Uses the arithmetic mean of the vertices (not the area-weighted centroid),
804 : // emits one triangle per polygon edge without degeneracy filtering, and skips
805 : // the canonicalization pass which would drop near-degenerate vertices and
806 : // perturb integration weights in downstream test baselines.
807 44082 : if (_triangulation_mode == MortarSegmentTriangulationMode::Centroid && !_triangulate_triangles)
808 : {
809 34218 : if (poly_nodes.size() == 3)
810 : {
811 28350 : tri_map.push_back({0, 1, 2});
812 14175 : return;
813 : }
814 :
815 20043 : const unsigned int n_verts = poly_nodes.size();
816 20043 : Point poly_center;
817 104688 : for (const auto & node : poly_nodes)
818 84645 : poly_center += node;
819 20043 : poly_center /= n_verts;
820 :
821 104688 : for (const auto i : make_range(n_verts))
822 253935 : tri_map.push_back({i, (i + 1) % n_verts, n_verts});
823 :
824 20043 : poly_nodes.push_back(poly_center);
825 20043 : return;
826 : }
827 :
828 9864 : canonicalize_polygon();
829 9864 : if (poly_nodes.size() < 3)
830 0 : return;
831 :
832 9864 : if (poly_nodes.size() == 3 && !_triangulate_triangles)
833 : {
834 783 : append_triangle(0, 1, 2);
835 783 : return;
836 : }
837 :
838 9081 : const bool force_triangle_centroid_split = _triangulate_triangles && poly_nodes.size() == 3;
839 :
840 9081 : if (_triangulation_mode == MortarSegmentTriangulationMode::Vertex &&
841 2205 : !force_triangle_centroid_split)
842 : {
843 2205 : const unsigned int n_verts = poly_nodes.size();
844 7074 : for (unsigned int i = 1; i + 1 < n_verts; ++i)
845 4869 : append_triangle(0, i, i + 1);
846 2205 : return;
847 : }
848 :
849 6876 : if (_triangulation_mode == MortarSegmentTriangulationMode::Delaunay &&
850 2205 : !force_triangle_centroid_split)
851 : {
852 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
853 2205 : triangulateConstrainedDelaunayPolygon(poly_nodes, _area_tol, _length_tol, tri_map);
854 2205 : return;
855 : #else
856 : mooseError("The 'delaunay' mortar triangulation mode requires libMesh TriangleInterface or "
857 : "Poly2Tri support.");
858 : #endif
859 : }
860 :
861 4671 : if (_triangulation_mode == MortarSegmentTriangulationMode::EarClipping &&
862 2205 : !force_triangle_centroid_split)
863 : {
864 7074 : for (const auto & triangle : triangulate_with_ear_clipping(true))
865 7074 : append_triangle(triangle[0], triangle[1], triangle[2]);
866 2205 : return;
867 : }
868 :
869 2466 : if (!force_triangle_centroid_split && !is_convex_polygon(poly_nodes))
870 : {
871 0 : for (const auto & triangle : triangulate_with_ear_clipping(true))
872 0 : append_triangle(triangle[0], triangle[1], triangle[2]);
873 0 : return;
874 : }
875 :
876 2466 : const unsigned int n_verts = poly_nodes.size();
877 2466 : const Point poly_center = polygon_centroid(poly_nodes);
878 :
879 2466 : bool added_triangle = false;
880 12528 : for (const auto i : make_range(n_verts))
881 10062 : if (triangleAreaHelper(poly_nodes[i], poly_nodes[(i + 1) % n_verts], poly_center) > _area_tol)
882 : {
883 20070 : tri_map.push_back({i, (i + 1) % n_verts, n_verts});
884 10035 : added_triangle = true;
885 : }
886 :
887 2466 : if (added_triangle)
888 2457 : poly_nodes.push_back(poly_center);
889 : }
890 :
891 : void
892 340981 : MortarSegmentHelper::getMortarSegments(const std::vector<Point> & primary_nodes,
893 : std::vector<Point> & nodes,
894 : std::vector<std::vector<unsigned int>> & elem_to_nodes)
895 : {
896 : // Clip primary elem against secondary elem
897 340981 : std::vector<Point> clipped_poly = clipPoly(primary_nodes);
898 340981 : if (clipped_poly.size() < 3)
899 296899 : return;
900 :
901 44082 : if (_debug)
902 0 : for (auto pt : clipped_poly)
903 0 : if (!isInsideSecondary(pt))
904 0 : mooseError("Clipped polygon not inside linearized secondary element");
905 :
906 : // Compute area of clipped polygon, update remaining area fraction
907 44082 : _remaining_area_fraction -= area(clipped_poly) / _secondary_area;
908 :
909 : // Triangulate clip polygon. tri_map indices are local to clipped_poly (starting at 0); we
910 : // shift them into the global node numbering after appending the polygon nodes below.
911 44082 : std::vector<std::vector<unsigned int>> tri_map;
912 44082 : triangulatePoly(clipped_poly, tri_map);
913 44082 : if (tri_map.empty())
914 36 : return;
915 :
916 : // Transform clipped poly back to (linearized) 3d and append to list
917 44046 : const auto offset = cast_int<unsigned int>(nodes.size());
918 233856 : for (auto pt : clipped_poly)
919 189810 : nodes.emplace_back((pt(0) * _u) + (pt(1) * _v) + _center);
920 :
921 168264 : for (const auto & tri : tri_map)
922 : {
923 124218 : std::vector<unsigned int> shifted_tri;
924 124218 : shifted_tri.reserve(tri.size());
925 496872 : for (const auto local_index : tri)
926 372654 : shifted_tri.push_back(offset + local_index);
927 124218 : elem_to_nodes.push_back(std::move(shifted_tri));
928 124218 : }
929 341017 : }
930 :
931 : Real
932 73944 : MortarSegmentHelper::area(const std::vector<Point> & nodes) const
933 : {
934 73944 : Real poly_area = 0;
935 356922 : for (auto i : index_range(nodes))
936 282978 : poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
937 282978 : nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
938 73944 : poly_area *= 0.5;
939 73944 : return poly_area;
940 : }
|
