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 : #include "MooseTypes.h"
12 :
13 : #include "libmesh/fe_interface.h"
14 : #include "libmesh/fe_map.h"
15 : #include "libmesh/face_quad4.h"
16 : #include "libmesh/face_tri3.h"
17 : #include "libmesh/int_range.h"
18 : #include "libmesh/node.h"
19 : #include "libmesh/utility.h"
20 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
21 : #include "libmesh/replicated_mesh.h"
22 : #include "libmesh/mesh_triangle_interface.h"
23 : #include "libmesh/poly2tri_triangulator.h"
24 : #endif
25 :
26 : #include <algorithm>
27 : #include <array>
28 : #include <cmath>
29 : #include <limits>
30 : #include <map>
31 : #include <numeric>
32 : #include <optional>
33 : #include <sstream>
34 : #include <set>
35 : #include <string>
36 : #include <unordered_map>
37 : #include <utility>
38 :
39 : using namespace libMesh;
40 :
41 : namespace
42 : {
43 :
44 : constexpr Real mortar_reference_mapping_tolerance = 1e-8;
45 :
46 : bool
47 368358 : isFinitePoint(const Point & point)
48 : {
49 1473432 : for (const auto component : make_range(Moose::dim))
50 1105074 : if (!std::isfinite(point(component)))
51 0 : return false;
52 :
53 368358 : return true;
54 : }
55 :
56 : // Signed-area test for the 2D triangle (a, b, c). Returns twice the signed area:
57 : // positive if a->b->c is counter-clockwise, negative if clockwise, zero if
58 : // collinear. Used as the building block for orientation, point-in-triangle, and
59 : // circumcircle predicates.
60 : Real
61 221600 : orient2dHelper(const Point & a, const Point & b, const Point & c)
62 : {
63 221600 : return (b(0) - a(0)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(0) - a(0));
64 : }
65 :
66 : Real
67 108918 : triangleAreaHelper(const Point & a, const Point & b, const Point & c)
68 : {
69 108918 : return 0.5 * std::abs(orient2dHelper(a, b, c));
70 : }
71 :
72 : // Canonical key for an undirected edge: the two endpoint indices sorted so that
73 : // (a, b) and (b, a) hash and compare equal. Used to dedupe / look up edges in
74 : // triangle-adjacency maps.
75 : std::array<unsigned int, 2>
76 55905 : canonicalEdgeHelper(const unsigned int a, const unsigned int b)
77 : {
78 55905 : return {{std::min(a, b), std::max(a, b)}};
79 : }
80 :
81 : // Reorder the three vertex indices (a, b, c) so the resulting triangle is wound
82 : // counter-clockwise (CCW) in the 2D plane spanned by \p nodes. Many of the
83 : // triangulation paths (orientation tests, area accumulation, ear-clipping
84 : // validity checks) assume CCW input, so we normalize before emitting triangles.
85 : std::array<unsigned int, 3>
86 6711 : makeCCWTriangleHelper(const std::vector<Point> & nodes,
87 : const unsigned int a,
88 : const unsigned int b,
89 : const unsigned int c)
90 : {
91 6711 : if (orient2dHelper(nodes[a], nodes[b], nodes[c]) >= 0)
92 5679 : return {{a, b, c}};
93 1032 : return {{a, c, b}};
94 : }
95 :
96 : bool
97 6069 : pointInCircumcircleHelper(const Point & a, const Point & b, const Point & c, const Point & p)
98 : {
99 6069 : const auto ax = a(0) - p(0);
100 6069 : const auto ay = a(1) - p(1);
101 6069 : const auto bx = b(0) - p(0);
102 6069 : const auto by = b(1) - p(1);
103 6069 : const auto cx = c(0) - p(0);
104 6069 : const auto cy = c(1) - p(1);
105 6069 : const Real det = (ax * ax + ay * ay) * (bx * cy - by * cx) -
106 6069 : (bx * bx + by * by) * (ax * cy - ay * cx) +
107 6069 : (cx * cx + cy * cy) * (ax * by - ay * bx);
108 6069 : const Real orientation = orient2dHelper(a, b, c);
109 6069 : return orientation >= 0 ? det > TOLERANCE : det < -TOLERANCE;
110 : }
111 :
112 : void
113 4422 : performLocalDelaunayFlips(const std::vector<Point> & poly_nodes,
114 : const std::set<std::array<unsigned int, 2>> & constrained_edges,
115 : std::vector<std::array<unsigned int, 3>> & triangles)
116 : {
117 4422 : bool flipped = true;
118 9756 : while (flipped)
119 : {
120 5334 : flipped = false;
121 :
122 5334 : std::map<std::array<unsigned int, 2>, std::vector<unsigned int>> edge_to_triangles;
123 17763 : for (const auto tri_index : index_range(triangles))
124 : {
125 12429 : const auto & tri = triangles[tri_index];
126 12429 : edge_to_triangles[canonicalEdgeHelper(tri[0], tri[1])].push_back(tri_index);
127 12429 : edge_to_triangles[canonicalEdgeHelper(tri[1], tri[2])].push_back(tri_index);
128 12429 : edge_to_triangles[canonicalEdgeHelper(tri[2], tri[0])].push_back(tri_index);
129 : }
130 :
131 32178 : for (const auto & [edge, owning_triangles] : edge_to_triangles)
132 : {
133 27756 : if (owning_triangles.size() != 2 || constrained_edges.count(edge))
134 20940 : continue;
135 :
136 6816 : const auto first_tri_index = owning_triangles[0];
137 6816 : const auto second_tri_index = owning_triangles[1];
138 6816 : const auto & first_triangle = triangles[first_tri_index];
139 6816 : const auto & second_triangle = triangles[second_tri_index];
140 :
141 6816 : const auto a = edge[0];
142 6816 : const auto b = edge[1];
143 : const auto first_opposite =
144 6816 : *std::find_if(first_triangle.begin(),
145 : first_triangle.end(),
146 12141 : [a, b](const unsigned int vertex) { return vertex != a && vertex != b; });
147 : const auto second_opposite =
148 6816 : *std::find_if(second_triangle.begin(),
149 : second_triangle.end(),
150 15024 : [a, b](const unsigned int vertex) { return vertex != a && vertex != b; });
151 :
152 6816 : if (first_opposite == second_opposite)
153 0 : continue;
154 :
155 : const auto side_a =
156 6816 : orient2dHelper(poly_nodes[first_opposite], poly_nodes[second_opposite], poly_nodes[a]);
157 : const auto side_b =
158 6816 : orient2dHelper(poly_nodes[first_opposite], poly_nodes[second_opposite], poly_nodes[b]);
159 6816 : if (side_a * side_b >= -TOLERANCE)
160 747 : continue;
161 :
162 6069 : if (!pointInCircumcircleHelper(poly_nodes[first_triangle[0]],
163 6069 : poly_nodes[first_triangle[1]],
164 6069 : poly_nodes[first_triangle[2]],
165 6069 : poly_nodes[second_opposite]))
166 5157 : continue;
167 :
168 912 : triangles[first_tri_index] =
169 912 : makeCCWTriangleHelper(poly_nodes, first_opposite, second_opposite, b);
170 912 : triangles[second_tri_index] =
171 912 : makeCCWTriangleHelper(poly_nodes, second_opposite, first_opposite, a);
172 912 : flipped = true;
173 912 : break;
174 : }
175 5334 : }
176 4422 : }
177 :
178 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
179 : void
180 2211 : triangulateConstrainedDelaunayPolygon(std::vector<Point> & poly_nodes,
181 : const Real area_tol,
182 : const Real length_tol,
183 : std::vector<std::vector<unsigned int>> & tri_map)
184 : {
185 2211 : Parallel::Communicator comm_self;
186 2211 : ReplicatedMesh triangulation_mesh(comm_self, 2);
187 2211 : std::unordered_map<dof_id_type, unsigned int> node_id_to_local_index;
188 2211 : node_id_to_local_index.reserve(poly_nodes.size());
189 :
190 11520 : for (const auto i : index_range(poly_nodes))
191 9309 : triangulation_mesh.add_point(poly_nodes[i], i);
192 :
193 2211 : triangulation_mesh.set_mesh_dimension(2);
194 :
195 : #ifdef LIBMESH_HAVE_TRIANGLE
196 : TriangleInterface triangulator(triangulation_mesh);
197 : #else
198 2211 : Poly2TriTriangulator triangulator(triangulation_mesh);
199 2211 : triangulator.set_refine_boundary_allowed(false);
200 : #endif
201 :
202 2211 : triangulator.triangulation_type() = TriangulatorInterface::PSLG;
203 2211 : triangulator.elem_type() = TRI3;
204 2211 : triangulator.set_interpolate_boundary_points(0);
205 2211 : triangulator.set_verify_hole_boundaries(false);
206 2211 : triangulator.desired_area() = 0;
207 2211 : triangulator.minimum_angle() = 0;
208 2211 : triangulator.smooth_after_generating() = false;
209 2211 : triangulator.quiet() = true;
210 2211 : triangulator.segments.reserve(poly_nodes.size());
211 11520 : for (const auto i : index_range(poly_nodes))
212 9309 : triangulator.segments.emplace_back(i, (i + 1) % poly_nodes.size());
213 :
214 2211 : triangulator.triangulate();
215 :
216 : // node_ptr_range() and active_element_ptr_range() iterate in id order on this
217 : // serial ReplicatedMesh, so no explicit sort is needed.
218 11520 : for (const auto * const node : triangulation_mesh.node_ptr_range())
219 9309 : if (!node_id_to_local_index.count(node->id()))
220 : {
221 : // Node inherits from Point and the triangulator operates on a 2D plane, so
222 : // the libMesh node already lives at z = 0 and we can use it directly.
223 9309 : unsigned int matched_index = libMesh::invalid_uint;
224 9309 : Real best_distance = std::numeric_limits<Real>::max();
225 :
226 48870 : for (const auto i : index_range(poly_nodes))
227 : {
228 39561 : const Real distance = (*node - poly_nodes[i]).norm();
229 39561 : if (distance <= length_tol && distance < best_distance)
230 : {
231 9309 : matched_index = i;
232 9309 : best_distance = distance;
233 : }
234 : }
235 :
236 9309 : if (matched_index == libMesh::invalid_uint)
237 : {
238 0 : matched_index = cast_int<unsigned int>(poly_nodes.size());
239 0 : poly_nodes.push_back(*node);
240 : }
241 :
242 9309 : node_id_to_local_index.emplace(node->id(), matched_index);
243 2211 : }
244 :
245 2211 : std::vector<std::array<unsigned int, 3>> triangles;
246 2211 : triangles.reserve(triangulation_mesh.n_elem());
247 :
248 7098 : for (const auto * const elem : triangulation_mesh.active_element_ptr_range())
249 : {
250 : mooseAssert(elem->type() == TRI3,
251 : "The delaunay mortar triangulation backend produced a non-TRI3 element: "
252 : << static_cast<int>(elem->type()));
253 :
254 : std::array<unsigned int, 3> local_triangle;
255 19548 : for (const auto i : index_range(local_triangle))
256 14661 : local_triangle[i] = libmesh_map_find(node_id_to_local_index, elem->node_id(i));
257 :
258 4887 : const Real orientation = orient2dHelper(poly_nodes[local_triangle[0]],
259 4887 : poly_nodes[local_triangle[1]],
260 4887 : poly_nodes[local_triangle[2]]);
261 4887 : if (std::abs(orientation) <= 2. * area_tol)
262 0 : continue;
263 :
264 4887 : if (orientation < 0)
265 0 : std::swap(local_triangle[1], local_triangle[2]);
266 :
267 4887 : triangles.push_back(local_triangle);
268 2211 : }
269 :
270 2211 : std::set<std::array<unsigned int, 2>> constrained_edges;
271 11520 : for (const auto i : index_range(poly_nodes))
272 9309 : constrained_edges.insert(canonicalEdgeHelper(i, (i + 1) % poly_nodes.size()));
273 :
274 2211 : performLocalDelaunayFlips(poly_nodes, constrained_edges, triangles);
275 :
276 2211 : std::set<std::array<unsigned int, 3>> seen_triangles;
277 7098 : for (auto local_triangle : triangles)
278 : {
279 4887 : auto canonical_triangle = local_triangle;
280 4887 : std::sort(canonical_triangle.begin(), canonical_triangle.end());
281 4887 : if (!seen_triangles.insert(canonical_triangle).second)
282 0 : continue;
283 :
284 14661 : tri_map.push_back({local_triangle[0], local_triangle[1], local_triangle[2]});
285 : }
286 2211 : }
287 : #endif
288 :
289 : } // namespace
290 :
291 11051 : MortarSegmentHelper::MortarSegmentHelper(std::vector<Point> secondary_nodes,
292 : const Point & center,
293 : const Point & normal,
294 : const MortarSegmentTriangulationMode triangulation_mode,
295 11051 : const bool triangulate_triangles)
296 : : MortarSegmentHelper(
297 11051 : std::move(secondary_nodes), {}, center, normal, triangulation_mode, triangulate_triangles)
298 : {
299 11051 : }
300 :
301 11781 : MortarSegmentHelper::MortarSegmentHelper(std::vector<Point> secondary_nodes,
302 : std::vector<Point> secondary_reference_points,
303 : const Point & center,
304 : const Point & normal,
305 : const MortarSegmentTriangulationMode triangulation_mode,
306 11781 : const bool triangulate_triangles)
307 11781 : : _center(center),
308 11781 : _normal(normal),
309 11781 : _debug(false),
310 11781 : _triangulation_mode(triangulation_mode),
311 11781 : _triangulate_triangles(triangulate_triangles),
312 11781 : _secondary_reference_points(std::move(secondary_reference_points))
313 : {
314 : mooseAssert(_secondary_reference_points.empty() ||
315 : secondary_nodes.size() == _secondary_reference_points.size(),
316 : "Each projected secondary node needs one parent reference point.");
317 :
318 11781 : _secondary_poly.clear();
319 11781 : _secondary_poly.reserve(secondary_nodes.size());
320 :
321 : // Get orientation of secondary poly
322 11781 : const Point e1 = secondary_nodes[0] - secondary_nodes[1];
323 11781 : const Point e2 = secondary_nodes[2] - secondary_nodes[1];
324 11781 : const Real orient = e2.cross(e1) * _normal;
325 :
326 : // u and v define the tangent plane of the element (at center)
327 : // Note we embed orientation into our transformation to make 2D poly always
328 : // positively oriented
329 11781 : _u = _normal.cross(secondary_nodes[0] - center).unit();
330 11781 : _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
331 :
332 : // Transform problem to 2D plane spanned by u and v
333 52695 : for (const auto & node : secondary_nodes)
334 : {
335 40914 : Point pt = node - _center;
336 40914 : _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
337 : }
338 :
339 : // Initialize area of secondary polygon
340 11781 : _remaining_area_fraction = 1.0;
341 11781 : _secondary_area = area(_secondary_poly);
342 :
343 : // Tolerance for quantities with area dimensions
344 11781 : _area_tol = _tolerance * _secondary_area;
345 :
346 : // Tolerance for quantites with length dimensions
347 11781 : _length_tol = _tolerance * std::sqrt(_secondary_area);
348 11781 : }
349 :
350 : Point
351 224006 : MortarSegmentHelper::getIntersection(
352 : const Point & p1, const Point & p2, const Point & q1, const Point & q2, Real & s) const
353 : {
354 224006 : const Point dp = p2 - p1;
355 224006 : const Point dq = q2 - q1;
356 224006 : const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
357 224006 : const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
358 224006 : const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
359 224006 : const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
360 224006 : s = -alpha * (cp1q2 - cp1q1 - cq1q2);
361 :
362 : // Intersection should be between p1 and p2, if it's not (due to poor conditioning), simply
363 : // move it to one of the end points
364 224006 : s = s > 1 ? 1. : s;
365 224006 : s = s < 0 ? 0. : s;
366 224006 : return p1 + s * dp;
367 : }
368 :
369 : bool
370 0 : MortarSegmentHelper::isInsideSecondary(const Point & pt) const
371 : {
372 0 : for (auto i : index_range(_secondary_poly))
373 : {
374 0 : const Point & q1 = _secondary_poly[i];
375 0 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
376 :
377 0 : const Point e1 = q2 - q1;
378 0 : const Point e2 = pt - q1;
379 :
380 : // If point corresponds to one of the secondary vertices, skip
381 0 : if (e2.norm() < _tolerance)
382 0 : return true;
383 :
384 0 : const bool inside = (e1(0) * e2(1) - e1(1) * e2(0)) < _area_tol;
385 0 : if (!inside)
386 0 : return false;
387 : }
388 0 : return true;
389 : }
390 :
391 : bool
392 393422 : MortarSegmentHelper::isDisjoint(const std::vector<Point> & poly) const
393 : {
394 935378 : for (auto i : index_range(_secondary_poly))
395 : {
396 : // Get edge to check
397 875031 : const Point & q1 = _secondary_poly[i];
398 875031 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
399 875031 : const Point edg = q2 - q1;
400 875031 : const Real cp = q2(0) * q1(1) - q2(1) * q1(0);
401 :
402 : // If more optimization needed, could store these values for later
403 : // Check if point is to the left of (or on) clip_edge
404 2956597 : auto is_inside = [&edg, cp](Point & pt, Real tol)
405 2956597 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
406 :
407 875031 : bool all_outside = true;
408 3831628 : for (auto pt : poly)
409 2956597 : if (is_inside(pt, _area_tol))
410 1518284 : all_outside = false;
411 :
412 875031 : if (all_outside)
413 333075 : return true;
414 : }
415 60347 : return false;
416 : }
417 :
418 : std::vector<Point>
419 393422 : MortarSegmentHelper::projectPrimaryPoly(const std::vector<Point> & primary_nodes) const
420 : {
421 : // Check orientation of primary_poly
422 393422 : const Point e1 = primary_nodes[0] - primary_nodes[1];
423 393422 : const Point e2 = primary_nodes[2] - primary_nodes[1];
424 :
425 : // Note we use u x v here instead of normal because it may be flipped if secondary elem was
426 : // negatively oriented
427 393422 : const Real orient = e2.cross(e1) * _u.cross(_v);
428 :
429 : // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
430 393422 : std::vector<Point> primary_poly;
431 393422 : const int n_verts = primary_nodes.size();
432 393422 : primary_poly.reserve(primary_nodes.size());
433 1704539 : for (auto n : index_range(primary_nodes))
434 : {
435 1311117 : Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
436 1311117 : primary_poly.emplace_back(pt * _u, pt * _v, 0.);
437 : }
438 :
439 786844 : return primary_poly;
440 0 : }
441 :
442 : std::vector<Point>
443 352108 : MortarSegmentHelper::clipPoly(const std::vector<Point> & primary_nodes) const
444 : {
445 352108 : return clipProjectedPoly(projectPrimaryPoly(primary_nodes));
446 : }
447 :
448 : std::vector<Point>
449 393422 : MortarSegmentHelper::clipProjectedPoly(const std::vector<Point> & primary_poly) const
450 : {
451 393422 : if (isDisjoint(primary_poly))
452 333075 : return {};
453 :
454 : // Initialize clipped poly with secondary poly (secondary is target poly)
455 60347 : std::vector<Point> clipped_poly = _secondary_poly;
456 :
457 : // Loop through clipping edges
458 263461 : for (auto i : index_range(primary_poly))
459 : {
460 : // If clipped poly trivial, return
461 208782 : if (clipped_poly.size() < 3)
462 : {
463 5668 : clipped_poly.clear();
464 5668 : return clipped_poly;
465 : }
466 :
467 : // Set input poly to current clipped poly
468 203114 : std::vector<Point> input_poly(clipped_poly);
469 203114 : clipped_poly.clear();
470 :
471 : // Get clipping edge
472 203114 : const Point & clip_pt1 = primary_poly[i];
473 203114 : const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
474 203114 : const Point edg = clip_pt2 - clip_pt1;
475 203114 : const Real cp = clip_pt2(0) * clip_pt1(1) - clip_pt2(1) * clip_pt1(0);
476 :
477 : // Check if point is to the left of (or on) clip_edge
478 : /*
479 : * Note that use of tolerance here is to avoid degenerate case when lines are
480 : * essentially on top of each other (common when meshes match across interface)
481 : * since finding intersection is ill-conditioned in this case.
482 : */
483 1518780 : auto is_inside = [&edg, cp](const Point & pt, Real tol)
484 1518780 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < tol; };
485 :
486 : // Loop through edges of target polygon (with previous clippings already included)
487 962504 : for (auto j : index_range(input_poly))
488 : {
489 : // Get target edge
490 759390 : const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
491 759390 : const Point prev_pt = input_poly[j];
492 :
493 : // TODO: Don't need to calculate both each loop
494 759390 : const bool is_current_inside = is_inside(curr_pt, _area_tol);
495 759390 : const bool is_previous_inside = is_inside(prev_pt, _area_tol);
496 :
497 759390 : if (is_current_inside)
498 : {
499 524367 : if (!is_previous_inside)
500 : {
501 : Real s;
502 112003 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
503 :
504 : /*
505 : * s is the fraction of distance along clip poly edge that intersection lies
506 : * It is used here to avoid degenerate polygon cases. For example, consider a
507 : * case like:
508 : * o
509 : * | (inside)
510 : * ------|------
511 : * | (outside)
512 : * when the distance is small (< 1e-7) we don't want to to add both the point
513 : * and intersection. Also note that when distance on the scale of 1e-7,
514 : * area on scale of 1e-14 so is insignificant if this results in dropping
515 : * a tri (for example if next edge crosses again)
516 : */
517 112003 : if (s < (1 - _tolerance))
518 110650 : clipped_poly.push_back(intersect);
519 : }
520 524367 : clipped_poly.push_back(curr_pt);
521 : }
522 235023 : else if (is_previous_inside)
523 : {
524 : Real s;
525 112003 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
526 112003 : if (s > _tolerance)
527 110379 : clipped_poly.push_back(intersect);
528 : }
529 : }
530 203114 : }
531 :
532 : // Make sure final clipped poly is not trivial
533 54679 : if (clipped_poly.size() < 3)
534 : {
535 1991 : clipped_poly.clear();
536 1991 : return clipped_poly;
537 : }
538 :
539 : // Clean up result by removing any duplicate nodes
540 52688 : std::vector<Point> cleaned_poly;
541 52688 : cleaned_poly.push_back(clipped_poly.back());
542 200526 : for (auto i : make_range(clipped_poly.size() - 1))
543 : {
544 147838 : const Point prev_pt = cleaned_poly.back();
545 147838 : const Point curr_pt = clipped_poly[i];
546 :
547 : // If points are sufficiently distanced, add to output
548 147838 : if ((curr_pt - prev_pt).norm() > _length_tol)
549 147838 : cleaned_poly.push_back(curr_pt);
550 : }
551 :
552 : mooseAssert(
553 : cleaned_poly.size() <= 8,
554 : "Our distributed mesh numbering scheme assumes that we have at most 8 nodes resulting from "
555 : "clipping the projection of the primary sub-element onto the secondary sub-element");
556 52688 : return cleaned_poly;
557 60347 : }
558 :
559 : void
560 52688 : MortarSegmentHelper::triangulatePoly(std::vector<Point> & poly_nodes,
561 : std::vector<std::vector<unsigned int>> & tri_map) const
562 : {
563 : // tri_map is populated with triangle indices that are local to poly_nodes (starting at 0).
564 : // Callers are responsible for shifting these indices into a global node numbering.
565 6522 : const auto polygon_centroid = [](const std::vector<Point> & polygon_nodes)
566 : {
567 6522 : Point centroid(0);
568 6522 : Real double_area = 0;
569 31622 : for (const auto i : index_range(polygon_nodes))
570 : {
571 25100 : const auto & a = polygon_nodes[i];
572 25100 : const auto & b = polygon_nodes[(i + 1) % polygon_nodes.size()];
573 25100 : const Real cross = a(0) * b(1) - b(0) * a(1);
574 25100 : double_area += cross;
575 25100 : centroid(0) += (a(0) + b(0)) * cross;
576 25100 : centroid(1) += (a(1) + b(1)) * cross;
577 : }
578 :
579 6522 : if (std::abs(double_area) <= TOLERANCE)
580 : {
581 410 : for (const auto & node : polygon_nodes)
582 309 : centroid += node;
583 101 : centroid /= polygon_nodes.size();
584 101 : return centroid;
585 : }
586 :
587 6421 : centroid /= (3. * double_area);
588 6421 : centroid(2) = 0;
589 6421 : return centroid;
590 : };
591 :
592 10557 : const auto append_triangle = [this, &poly_nodes, &tri_map](
593 : const unsigned int a, const unsigned int b, const unsigned int c)
594 : {
595 10557 : if (triangleAreaHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) <= _area_tol)
596 0 : return false;
597 :
598 10557 : if (orient2dHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) >= 0)
599 31671 : tri_map.push_back({a, b, c});
600 : else
601 0 : tri_map.push_back({a, c, b});
602 :
603 10557 : return true;
604 52688 : };
605 :
606 : const auto point_in_triangle =
607 13494 : [this](const Point & p, const Point & a, const Point & b, const Point & c)
608 : {
609 13494 : const Real ab = orient2dHelper(a, b, p);
610 13494 : const Real bc = orient2dHelper(b, c, p);
611 13494 : const Real ca = orient2dHelper(c, a, p);
612 13494 : return ab >= -_area_tol && bc >= -_area_tol && ca >= -_area_tol;
613 52688 : };
614 :
615 11169 : const auto min_triangle_angle = [](const Point & a, const Point & b, const Point & c)
616 : {
617 33471 : const auto clamp_cos = [](Real value) { return std::max(-1., std::min(1., value)); };
618 : const auto angle_at =
619 33507 : [&clamp_cos](const Point & vertex, const Point & point_one, const Point & point_two)
620 : {
621 33507 : const Point edge_one = point_one - vertex;
622 33507 : const Point edge_two = point_two - vertex;
623 33507 : const Real denom = edge_one.norm() * edge_two.norm();
624 33507 : if (denom <= TOLERANCE)
625 36 : return 0.;
626 33471 : return std::acos(clamp_cos((edge_one * edge_two) / denom));
627 11169 : };
628 :
629 11169 : return std::min({angle_at(a, b, c), angle_at(b, c, a), angle_at(c, a, b)});
630 : };
631 :
632 13938 : const auto canonicalize_polygon = [this, &poly_nodes]()
633 : {
634 13938 : if (poly_nodes.size() < 3)
635 0 : return;
636 :
637 13938 : if (area(poly_nodes) < 0)
638 0 : std::reverse(poly_nodes.begin(), poly_nodes.end());
639 :
640 13938 : bool changed = true;
641 25118 : while (changed && poly_nodes.size() > 3)
642 : {
643 11180 : changed = false;
644 58282 : for (const auto i : index_range(poly_nodes))
645 : {
646 47102 : const auto prev = (i + poly_nodes.size() - 1) % poly_nodes.size();
647 47102 : const auto next = (i + 1) % poly_nodes.size();
648 47102 : if ((poly_nodes[i] - poly_nodes[prev]).norm() <= _length_tol ||
649 94204 : (poly_nodes[next] - poly_nodes[i]).norm() <= _length_tol ||
650 47102 : triangleAreaHelper(poly_nodes[prev], poly_nodes[i], poly_nodes[next]) <= _area_tol)
651 : {
652 0 : poly_nodes.erase(poly_nodes.begin() + i);
653 0 : changed = true;
654 0 : break;
655 : }
656 : }
657 : }
658 :
659 13938 : if (poly_nodes.size() >= 3 && area(poly_nodes) < 0)
660 0 : std::reverse(poly_nodes.begin(), poly_nodes.end());
661 52688 : };
662 :
663 : const auto triangulate_with_ear_clipping =
664 2211 : [this, &poly_nodes, &point_in_triangle, &min_triangle_angle](
665 : const bool perform_delaunay_flips)
666 : {
667 2211 : std::vector<std::array<unsigned int, 3>> triangles;
668 2211 : if (poly_nodes.size() < 3)
669 0 : return triangles;
670 :
671 2211 : if (poly_nodes.size() == 3)
672 : {
673 0 : triangles.push_back(makeCCWTriangleHelper(poly_nodes, 0, 1, 2));
674 0 : return triangles;
675 : }
676 :
677 2211 : std::vector<unsigned int> remaining_vertices(poly_nodes.size());
678 2211 : std::iota(remaining_vertices.begin(), remaining_vertices.end(), 0);
679 :
680 4887 : while (remaining_vertices.size() > 3)
681 : {
682 2676 : std::optional<std::size_t> best_position;
683 2676 : Real best_score = -std::numeric_limits<Real>::max();
684 2676 : Real best_area = -std::numeric_limits<Real>::max();
685 :
686 13845 : for (const auto position : index_range(remaining_vertices))
687 : {
688 : const auto prev_position =
689 11169 : (position + remaining_vertices.size() - 1) % remaining_vertices.size();
690 11169 : const auto next_position = (position + 1) % remaining_vertices.size();
691 11169 : const auto prev = remaining_vertices[prev_position];
692 11169 : const auto curr = remaining_vertices[position];
693 11169 : const auto next = remaining_vertices[next_position];
694 :
695 11169 : if (orient2dHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]) <= _area_tol)
696 0 : continue;
697 :
698 11169 : bool contains_other_vertex = false;
699 58170 : for (const auto other : remaining_vertices)
700 : {
701 47001 : if (other == prev || other == curr || other == next)
702 33507 : continue;
703 :
704 13494 : if (point_in_triangle(
705 13494 : poly_nodes[other], poly_nodes[prev], poly_nodes[curr], poly_nodes[next]))
706 : {
707 0 : contains_other_vertex = true;
708 0 : break;
709 : }
710 : }
711 :
712 11169 : if (contains_other_vertex)
713 0 : continue;
714 :
715 : const Real candidate_score =
716 11169 : min_triangle_angle(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
717 : const Real candidate_area =
718 11169 : triangleAreaHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
719 16299 : if (!best_position || candidate_score > best_score + TOLERANCE ||
720 5058 : (std::abs(candidate_score - best_score) <= TOLERANCE &&
721 72 : candidate_area > best_area + _area_tol))
722 : {
723 6147 : best_position = position;
724 6147 : best_score = candidate_score;
725 6147 : best_area = candidate_area;
726 : }
727 : }
728 :
729 2676 : if (!best_position)
730 : {
731 0 : std::vector<std::array<unsigned int, 3>> best_fan;
732 0 : Real best_fan_score = -std::numeric_limits<Real>::max();
733 0 : Real best_fan_area = -std::numeric_limits<Real>::max();
734 :
735 0 : for (const auto root_position : index_range(remaining_vertices))
736 : {
737 0 : std::vector<std::array<unsigned int, 3>> candidate_fan;
738 0 : Real candidate_score = std::numeric_limits<Real>::max();
739 0 : Real candidate_area = std::numeric_limits<Real>::max();
740 0 : bool valid_fan = true;
741 0 : const auto root = remaining_vertices[root_position];
742 :
743 0 : for (unsigned int step = 1; step + 1 < remaining_vertices.size(); ++step)
744 : {
745 0 : const auto next_position = (root_position + step) % remaining_vertices.size();
746 0 : const auto following_position = (root_position + step + 1) % remaining_vertices.size();
747 0 : const auto vertex_one = remaining_vertices[next_position];
748 0 : const auto vertex_two = remaining_vertices[following_position];
749 :
750 0 : if (orient2dHelper(poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]) <=
751 0 : _area_tol)
752 : {
753 0 : valid_fan = false;
754 0 : break;
755 : }
756 :
757 0 : candidate_fan.push_back(
758 0 : makeCCWTriangleHelper(poly_nodes, root, vertex_one, vertex_two));
759 0 : candidate_score =
760 0 : std::min(candidate_score,
761 0 : min_triangle_angle(
762 0 : poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
763 0 : candidate_area =
764 0 : std::min(candidate_area,
765 0 : triangleAreaHelper(
766 0 : poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
767 : }
768 :
769 0 : if (!valid_fan || candidate_fan.empty())
770 0 : continue;
771 :
772 0 : if (candidate_score > best_fan_score + TOLERANCE ||
773 0 : (std::abs(candidate_score - best_fan_score) <= TOLERANCE &&
774 0 : candidate_area > best_fan_area + _area_tol))
775 : {
776 0 : best_fan = std::move(candidate_fan);
777 0 : best_fan_score = candidate_score;
778 0 : best_fan_area = candidate_area;
779 : }
780 0 : }
781 :
782 0 : if (best_fan.empty())
783 0 : for (unsigned int i = 1; i + 1 < remaining_vertices.size(); ++i)
784 0 : best_fan.push_back(makeCCWTriangleHelper(poly_nodes,
785 0 : remaining_vertices[0],
786 0 : remaining_vertices[i],
787 0 : remaining_vertices[i + 1]));
788 :
789 0 : triangles.insert(triangles.end(), best_fan.begin(), best_fan.end());
790 0 : break;
791 0 : }
792 :
793 : const auto prev_position =
794 2676 : (*best_position + remaining_vertices.size() - 1) % remaining_vertices.size();
795 2676 : const auto next_position = (*best_position + 1) % remaining_vertices.size();
796 2676 : triangles.push_back(makeCCWTriangleHelper(poly_nodes,
797 2676 : remaining_vertices[prev_position],
798 2676 : remaining_vertices[*best_position],
799 2676 : remaining_vertices[next_position]));
800 2676 : remaining_vertices.erase(remaining_vertices.begin() + *best_position);
801 : }
802 :
803 2211 : if (remaining_vertices.size() == 3)
804 2211 : triangles.push_back(makeCCWTriangleHelper(
805 2211 : poly_nodes, remaining_vertices[0], remaining_vertices[1], remaining_vertices[2]));
806 :
807 2211 : if (!perform_delaunay_flips)
808 0 : return triangles;
809 :
810 2211 : std::set<std::array<unsigned int, 2>> boundary_edges;
811 11520 : for (const auto i : index_range(poly_nodes))
812 9309 : boundary_edges.insert(canonicalEdgeHelper(i, (i + 1) % poly_nodes.size()));
813 :
814 2211 : performLocalDelaunayFlips(poly_nodes, boundary_edges, triangles);
815 2211 : return triangles;
816 2211 : };
817 :
818 4547 : const auto is_convex_polygon = [this](const std::vector<Point> & polygon_nodes)
819 : {
820 4547 : if (polygon_nodes.size() <= 3)
821 0 : return true;
822 :
823 23722 : for (const auto i : index_range(polygon_nodes))
824 : {
825 19175 : const auto prev = (i + polygon_nodes.size() - 1) % polygon_nodes.size();
826 19175 : const auto next = (i + 1) % polygon_nodes.size();
827 19175 : if (orient2dHelper(polygon_nodes[prev], polygon_nodes[i], polygon_nodes[next]) <= _area_tol)
828 0 : return false;
829 : }
830 :
831 4547 : return true;
832 52688 : };
833 :
834 : // Fewer than 3 nodes can't be triangulated
835 52688 : if (poly_nodes.size() < 3)
836 0 : mooseError("Can't triangulate poly with fewer than 3 nodes");
837 :
838 : // Legacy centroid path: when the default triangulation (centroid) is selected
839 : // and triangle re-tessellation is not requested, reproduce the legacy
840 : // algorithm byte-for-byte so existing mortar baselines remain valid.
841 : // Uses the arithmetic mean of the vertices (not the area-weighted centroid),
842 : // emits one triangle per polygon edge without degeneracy filtering, and skips
843 : // the canonicalization pass which would drop near-degenerate vertices and
844 : // perturb integration weights in downstream test baselines.
845 52688 : if (_triangulation_mode == MortarSegmentTriangulationMode::Centroid && !_triangulate_triangles)
846 : {
847 38750 : if (poly_nodes.size() == 3)
848 : {
849 30140 : tri_map.push_back({0, 1, 2});
850 15070 : return;
851 : }
852 :
853 23680 : const unsigned int n_verts = poly_nodes.size();
854 23680 : Point poly_center;
855 123620 : for (const auto & node : poly_nodes)
856 99940 : poly_center += node;
857 23680 : poly_center /= n_verts;
858 :
859 123620 : for (const auto i : make_range(n_verts))
860 299820 : tri_map.push_back({i, (i + 1) % n_verts, n_verts});
861 :
862 23680 : poly_nodes.push_back(poly_center);
863 23680 : return;
864 : }
865 :
866 13938 : canonicalize_polygon();
867 13938 : if (poly_nodes.size() < 3)
868 0 : return;
869 :
870 13938 : if (poly_nodes.size() == 3 && !_triangulate_triangles)
871 : {
872 783 : append_triangle(0, 1, 2);
873 783 : return;
874 : }
875 :
876 13155 : const bool force_triangle_centroid_split = _triangulate_triangles && poly_nodes.size() == 3;
877 :
878 13155 : if (_triangulation_mode == MortarSegmentTriangulationMode::Vertex &&
879 2215 : !force_triangle_centroid_split)
880 : {
881 2211 : const unsigned int n_verts = poly_nodes.size();
882 7098 : for (unsigned int i = 1; i + 1 < n_verts; ++i)
883 4887 : append_triangle(0, i, i + 1);
884 2211 : return;
885 : }
886 :
887 10944 : if (_triangulation_mode == MortarSegmentTriangulationMode::Delaunay &&
888 2215 : !force_triangle_centroid_split)
889 : {
890 : #if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
891 2211 : triangulateConstrainedDelaunayPolygon(poly_nodes, _area_tol, _length_tol, tri_map);
892 2211 : return;
893 : #else
894 : mooseError("The 'delaunay' mortar triangulation mode requires libMesh TriangleInterface or "
895 : "Poly2Tri support.");
896 : #endif
897 : }
898 :
899 8733 : if (_triangulation_mode == MortarSegmentTriangulationMode::EarClipping &&
900 2215 : !force_triangle_centroid_split)
901 : {
902 7098 : for (const auto & triangle : triangulate_with_ear_clipping(true))
903 7098 : append_triangle(triangle[0], triangle[1], triangle[2]);
904 2211 : return;
905 : }
906 :
907 6522 : if (!force_triangle_centroid_split && !is_convex_polygon(poly_nodes))
908 : {
909 0 : for (const auto & triangle : triangulate_with_ear_clipping(true))
910 0 : append_triangle(triangle[0], triangle[1], triangle[2]);
911 0 : return;
912 : }
913 :
914 6522 : const unsigned int n_verts = poly_nodes.size();
915 6522 : const Point poly_center = polygon_centroid(poly_nodes);
916 :
917 6522 : bool added_triangle = false;
918 31622 : for (const auto i : make_range(n_verts))
919 25100 : if (triangleAreaHelper(poly_nodes[i], poly_nodes[(i + 1) % n_verts], poly_center) > _area_tol)
920 : {
921 50200 : tri_map.push_back({i, (i + 1) % n_verts, n_verts});
922 25100 : added_triangle = true;
923 : }
924 :
925 6522 : if (added_triangle)
926 6522 : poly_nodes.push_back(poly_center);
927 : }
928 :
929 : void
930 352108 : MortarSegmentHelper::getMortarSegments(const std::vector<Point> & primary_nodes,
931 : std::vector<Point> & nodes,
932 : std::vector<std::vector<unsigned int>> & elem_to_nodes)
933 : {
934 352108 : getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, nullptr);
935 352108 : }
936 :
937 : void
938 41314 : MortarSegmentHelper::getMortarSegments(
939 : const std::vector<Point> & primary_nodes,
940 : const std::vector<Point> & primary_reference_points,
941 : std::vector<Point> & nodes,
942 : std::vector<std::vector<unsigned int>> & elem_to_nodes,
943 : std::vector<std::array<Point, 3>> & elem_to_secondary_reference_points,
944 : std::vector<std::array<Point, 3>> & elem_to_primary_reference_points,
945 : const Real minimum_segment_area)
946 : {
947 : ReferenceMappingData reference_mapping{primary_reference_points,
948 : elem_to_secondary_reference_points,
949 : elem_to_primary_reference_points,
950 41314 : minimum_segment_area};
951 41314 : getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, &reference_mapping);
952 41314 : }
953 :
954 : void
955 393422 : MortarSegmentHelper::getMortarSegmentsImpl(const std::vector<Point> & primary_nodes,
956 : std::vector<Point> & nodes,
957 : std::vector<std::vector<unsigned int>> & elem_to_nodes,
958 : ReferenceMappingData * const reference_mapping)
959 : {
960 393422 : std::vector<Point> primary_poly;
961 393422 : std::vector<Point> primary_poly_reference_points;
962 :
963 393422 : if (reference_mapping)
964 : {
965 41314 : if (primary_nodes.size() != reference_mapping->primary_reference_points.size())
966 0 : mooseError("Reference-interpolation mortar segment generation requires one primary "
967 : "reference point per primary sub-element node.");
968 41314 : if (_secondary_poly.size() != _secondary_reference_points.size())
969 0 : mooseError("Reference-interpolation mortar segment generation requires one secondary "
970 : "reference point per secondary sub-element node.");
971 82628 : if (reference_mapping->elem_to_secondary_reference_points.size() != elem_to_nodes.size() ||
972 41314 : reference_mapping->elem_to_primary_reference_points.size() != elem_to_nodes.size())
973 0 : mooseError("Reference-interpolation mortar segment outputs must be aligned before appending "
974 : "new segments.");
975 :
976 : // Keep reference points in the projected polygon's orientation.
977 41314 : const Point e1 = primary_nodes[0] - primary_nodes[1];
978 41314 : const Point e2 = primary_nodes[2] - primary_nodes[1];
979 41314 : const Real orient = e2.cross(e1) * _u.cross(_v);
980 41314 : const auto n_verts = primary_nodes.size();
981 :
982 41314 : primary_poly = projectPrimaryPoly(primary_nodes);
983 41314 : primary_poly_reference_points.reserve(reference_mapping->primary_reference_points.size());
984 176448 : for (const auto n : index_range(primary_nodes))
985 : {
986 135134 : const auto primary_node_index = (orient > 0) ? n : n_verts - 1 - n;
987 135134 : primary_poly_reference_points.push_back(
988 135134 : reference_mapping->primary_reference_points[primary_node_index]);
989 : }
990 : }
991 :
992 : // Clip primary elem against secondary elem. Reference mode preserves the projected primary
993 : // ordering so its reference points remain aligned.
994 : std::vector<Point> clipped_poly =
995 393422 : reference_mapping ? clipProjectedPoly(primary_poly) : clipPoly(primary_nodes);
996 393422 : if (clipped_poly.size() < 3)
997 340734 : return;
998 :
999 52688 : if (_debug)
1000 0 : for (const auto & point : clipped_poly)
1001 0 : if (!isInsideSecondary(point))
1002 0 : mooseError("Clipped polygon not inside linearized secondary element");
1003 :
1004 : // Compute area of clipped polygon, update remaining area fraction
1005 52688 : _remaining_area_fraction -= area(clipped_poly) / _secondary_area;
1006 :
1007 : // Triangulate clip polygon. tri_map indices are local to clipped_poly (starting at 0); we
1008 : // shift them into the global node numbering after appending the polygon nodes below.
1009 52688 : std::vector<std::vector<unsigned int>> tri_map;
1010 52688 : triangulatePoly(clipped_poly, tri_map);
1011 52688 : if (reference_mapping && reference_mapping->minimum_segment_area > 0.)
1012 8080 : tri_map.erase(
1013 4040 : std::remove_if(tri_map.begin(),
1014 : tri_map.end(),
1015 14990 : [&clipped_poly, reference_mapping](const std::vector<unsigned int> & tri)
1016 : {
1017 : mooseAssert(tri.size() == 3,
1018 : "Mortar segment triangulation should only produce TRI3 maps.");
1019 14990 : return triangleAreaHelper(clipped_poly[tri[0]],
1020 14990 : clipped_poly[tri[1]],
1021 14990 : clipped_poly[tri[2]]) <
1022 14990 : reference_mapping->minimum_segment_area;
1023 : }),
1024 8080 : tri_map.end());
1025 52688 : if (tri_map.empty())
1026 100 : return;
1027 :
1028 52588 : std::vector<Point> secondary_node_reference_points;
1029 52588 : std::vector<Point> primary_node_reference_points;
1030 52588 : if (reference_mapping)
1031 : {
1032 3964 : secondary_node_reference_points.reserve(clipped_poly.size());
1033 3964 : primary_node_reference_points.reserve(clipped_poly.size());
1034 :
1035 37476 : const auto recover_reference_point = [this](const Point & projected_point,
1036 : const std::vector<Point> & poly,
1037 : const std::vector<Point> & reference_points,
1038 : const char * const parent_name,
1039 : const std::size_t node_index)
1040 : {
1041 37476 : std::string failure_reason;
1042 : const auto reference_point =
1043 37476 : referencePoint(projected_point, poly, reference_points, &failure_reason);
1044 37476 : if (!reference_point)
1045 0 : mooseError("Unable to recover the ",
1046 : parent_name,
1047 : " parent reference point for retained 3D mortar overlap vertex ",
1048 : node_index,
1049 : " at projected point ",
1050 : projected_point,
1051 : ". Reason: ",
1052 : failure_reason,
1053 : ". Reference interpolation does not fall back to normal projection.");
1054 :
1055 74952 : return *reference_point;
1056 37476 : };
1057 :
1058 22702 : for (const auto node_index : index_range(clipped_poly))
1059 : {
1060 18738 : const auto & point = clipped_poly[node_index];
1061 18738 : secondary_node_reference_points.push_back(recover_reference_point(
1062 18738 : point, _secondary_poly, _secondary_reference_points, "secondary", node_index));
1063 18738 : primary_node_reference_points.push_back(recover_reference_point(
1064 : point, primary_poly, primary_poly_reference_points, "primary", node_index));
1065 : }
1066 : }
1067 :
1068 : // Transform clipped poly back to (linearized) 3d and append to list
1069 52588 : const auto offset = cast_int<unsigned int>(nodes.size());
1070 282908 : for (const auto & point : clipped_poly)
1071 230320 : nodes.emplace_back((point(0) * _u) + (point(1) * _v) + _center);
1072 :
1073 207830 : for (const auto & tri : tri_map)
1074 : {
1075 155242 : std::vector<unsigned int> shifted_tri;
1076 155242 : shifted_tri.reserve(tri.size());
1077 620968 : for (const auto local_index : tri)
1078 465726 : shifted_tri.push_back(offset + local_index);
1079 155242 : elem_to_nodes.push_back(std::move(shifted_tri));
1080 :
1081 155242 : if (reference_mapping)
1082 : {
1083 : mooseAssert(tri.size() == 3, "Mortar segment triangulation should only produce TRI3 maps.");
1084 14758 : std::array<Point, 3> elem_secondary_reference_points;
1085 14758 : std::array<Point, 3> elem_primary_reference_points;
1086 59032 : for (const auto n : index_range(tri))
1087 : {
1088 44274 : const auto local_node = tri[n];
1089 44274 : elem_secondary_reference_points[n] = secondary_node_reference_points[local_node];
1090 44274 : elem_primary_reference_points[n] = primary_node_reference_points[local_node];
1091 : }
1092 :
1093 14758 : reference_mapping->elem_to_secondary_reference_points.push_back(
1094 : elem_secondary_reference_points);
1095 14758 : reference_mapping->elem_to_primary_reference_points.push_back(elem_primary_reference_points);
1096 : }
1097 155242 : }
1098 1075190 : }
1099 :
1100 : std::optional<Point>
1101 37500 : MortarSegmentHelper::referencePoint(const Point & point,
1102 : const std::vector<Point> & poly,
1103 : const std::vector<Point> & reference_points,
1104 : std::string * const failure_reason) const
1105 : {
1106 : mooseAssert(poly.size() == reference_points.size(),
1107 : "Projected point and reference point containers should be the same size.");
1108 :
1109 37500 : if (failure_reason)
1110 37500 : failure_reason->clear();
1111 :
1112 12 : const auto fail = [failure_reason](const std::string & reason) -> std::optional<Point>
1113 : {
1114 12 : if (failure_reason)
1115 12 : *failure_reason = reason;
1116 12 : return std::nullopt;
1117 37500 : };
1118 :
1119 37500 : if (!isFinitePoint(point))
1120 0 : return fail("the projected target point contains a non-finite coordinate");
1121 :
1122 165440 : for (const auto i : index_range(poly))
1123 : {
1124 127940 : if (!isFinitePoint(poly[i]))
1125 0 : return fail("projected polygon vertex " + std::to_string(i) +
1126 0 : " contains a non-finite coordinate");
1127 127940 : if (!isFinitePoint(reference_points[i]))
1128 0 : return fail("parent reference vertex " + std::to_string(i) +
1129 0 : " contains a non-finite coordinate");
1130 : }
1131 :
1132 37500 : if (poly.size() != 3 && poly.size() != 4)
1133 0 : return fail("reference point recovery only supports triangular and quadrilateral mortar "
1134 0 : "sub-elements, but received " +
1135 0 : std::to_string(poly.size()) + " vertices");
1136 :
1137 37500 : Real minimum_edge_length = std::numeric_limits<Real>::max();
1138 37500 : Point local_origin;
1139 165440 : for (const auto & vertex : poly)
1140 127940 : local_origin += vertex;
1141 37500 : local_origin /= poly.size();
1142 :
1143 37500 : Real local_scale = 0.;
1144 165440 : for (const auto i : index_range(poly))
1145 : {
1146 127940 : minimum_edge_length =
1147 127940 : std::min(minimum_edge_length, (poly[(i + 1) % poly.size()] - poly[i]).norm());
1148 127940 : local_scale = std::max(local_scale, (poly[i] - local_origin).norm());
1149 : }
1150 :
1151 37500 : const Real singular_tolerance = 100. * std::numeric_limits<Real>::epsilon();
1152 37500 : if (!std::isfinite(local_scale) || local_scale <= singular_tolerance)
1153 0 : return fail("the projected polygon has a zero local length scale");
1154 75000 : if (!std::isfinite(minimum_edge_length) ||
1155 37500 : minimum_edge_length / local_scale <= singular_tolerance)
1156 0 : return fail("the projected polygon has a zero-length edge relative to its local scale");
1157 :
1158 37500 : std::vector<Point> normalized_poly;
1159 37500 : normalized_poly.reserve(poly.size());
1160 165440 : for (const auto & vertex : poly)
1161 127940 : normalized_poly.push_back((vertex - local_origin) / local_scale);
1162 37500 : const Point normalized_point = (point - local_origin) / local_scale;
1163 : const Real reference_tolerance =
1164 37500 : std::max(mortar_reference_mapping_tolerance, _area_tol / (minimum_edge_length * local_scale));
1165 37500 : std::array<Node, 4> element_nodes;
1166 37500 : const FEType fe_type(FIRST, LAGRANGE);
1167 37500 : const auto recover_with_libmesh = [&](auto & element) -> std::optional<Point>
1168 : {
1169 165440 : for (const auto i : index_range(normalized_poly))
1170 : {
1171 127940 : element_nodes[i] = normalized_poly[i];
1172 127940 : element_nodes[i].set_id(i);
1173 127940 : element.set_node(i, &element_nodes[i]);
1174 : }
1175 :
1176 37500 : if (!element.has_invertible_map(mortar_reference_mapping_tolerance))
1177 24 : return fail("the projected sub-element map is degenerate or non-invertible");
1178 :
1179 37492 : if (element.type() == QUAD4)
1180 77162 : for (const auto corner : make_range(element.n_vertices()))
1181 : {
1182 61730 : Point tangent_xi;
1183 61730 : Point tangent_eta;
1184 308650 : for (const auto node : make_range(element.n_nodes()))
1185 : {
1186 493840 : tangent_xi += FEInterface::shape_deriv(
1187 246920 : fe_type, 0, &element, node, 0, element.master_point(corner)) *
1188 246920 : element.point(node);
1189 493840 : tangent_eta += FEInterface::shape_deriv(
1190 493840 : fe_type, 0, &element, node, 1, element.master_point(corner)) *
1191 246920 : element.point(node);
1192 : }
1193 :
1194 61730 : const Real corner_jacobian = tangent_xi.cross(tangent_eta).norm();
1195 61730 : if (!std::isfinite(corner_jacobian) ||
1196 : corner_jacobian <= mortar_reference_mapping_tolerance)
1197 6 : return fail("the projected quadrilateral has a singular or ill-conditioned corner map");
1198 : }
1199 :
1200 37490 : Point local_reference = FEMap::inverse_map(
1201 : 2, &element, normalized_point, mortar_reference_mapping_tolerance, false, false);
1202 37490 : if (!isFinitePoint(local_reference))
1203 0 : return fail("libMesh inverse_map produced a non-finite reference point");
1204 :
1205 : const Real inverse_map_error =
1206 37490 : (FEMap::map(2, &element, local_reference) - normalized_point).norm();
1207 37490 : if (!std::isfinite(inverse_map_error) || inverse_map_error > mortar_reference_mapping_tolerance)
1208 : {
1209 0 : std::ostringstream reason;
1210 0 : reason << "the normalized inverse-map error " << inverse_map_error << " exceeds "
1211 0 : << mortar_reference_mapping_tolerance;
1212 0 : return fail(reason.str());
1213 0 : }
1214 :
1215 37490 : if (element.type() == TRI3)
1216 : {
1217 : std::array<Real, 3> weights;
1218 88232 : for (const auto i : index_range(weights))
1219 66174 : weights[i] = FEInterface::shape(fe_type, &element, i, local_reference, false);
1220 :
1221 88232 : for (auto & weight : weights)
1222 66174 : weight = std::clamp(weight, 0., 1.);
1223 22058 : const Real weight_sum = std::accumulate(weights.begin(), weights.end(), 0.);
1224 22058 : if (!std::isfinite(weight_sum) || weight_sum <= singular_tolerance)
1225 0 : return fail("clamped triangle barycentric coordinates have a zero or non-finite sum");
1226 88232 : for (auto & weight : weights)
1227 66174 : weight /= weight_sum;
1228 :
1229 : // Preserve partition of unity after clamping tolerance-sized violations.
1230 : const auto corrected_weight =
1231 22058 : std::distance(weights.begin(), std::max_element(weights.begin(), weights.end()));
1232 22058 : weights[corrected_weight] = 1.;
1233 88232 : for (const auto i : index_range(weights))
1234 66174 : if (i != static_cast<unsigned int>(corrected_weight))
1235 44116 : weights[corrected_weight] -= weights[i];
1236 :
1237 22058 : local_reference = Point();
1238 88232 : for (const auto i : index_range(weights))
1239 66174 : local_reference += weights[i] * element.master_point(i);
1240 : }
1241 : else
1242 : {
1243 15432 : local_reference(0) = std::clamp(local_reference(0), -1., 1.);
1244 15432 : local_reference(1) = std::clamp(local_reference(1), -1., 1.);
1245 15432 : local_reference(2) = 0.;
1246 : }
1247 :
1248 37490 : if (!element.on_reference_element(local_reference, mortar_reference_mapping_tolerance))
1249 0 : return fail("the clamped inverse-map result is outside the reference element");
1250 :
1251 : const Real round_trip_error =
1252 37490 : (FEMap::map(2, &element, local_reference) - normalized_point).norm();
1253 37490 : if (!std::isfinite(round_trip_error) || round_trip_error > reference_tolerance)
1254 : {
1255 2 : std::ostringstream reason;
1256 2 : reason << "the normalized inverse-map round-trip error " << round_trip_error
1257 2 : << " exceeds the clipping-consistent tolerance " << reference_tolerance;
1258 2 : return fail(reason.str());
1259 2 : }
1260 :
1261 37488 : Point parent_reference;
1262 165384 : for (const auto i : index_range(normalized_poly))
1263 127896 : parent_reference +=
1264 255792 : FEInterface::shape(fe_type, &element, i, local_reference, false) * reference_points[i];
1265 :
1266 37488 : if (!isFinitePoint(parent_reference))
1267 0 : return fail("reference interpolation produced a non-finite parent reference point");
1268 :
1269 37488 : return parent_reference;
1270 37500 : };
1271 :
1272 37500 : if (poly.size() == 3)
1273 : {
1274 22060 : Tri3 element;
1275 22060 : return recover_with_libmesh(element);
1276 22060 : }
1277 :
1278 15440 : Quad4 element;
1279 15440 : return recover_with_libmesh(element);
1280 37500 : }
1281 :
1282 : Real
1283 92345 : MortarSegmentHelper::area(const std::vector<Point> & nodes) const
1284 : {
1285 92345 : Real poly_area = 0;
1286 444537 : for (auto i : index_range(nodes))
1287 352192 : poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
1288 352192 : nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
1289 92345 : poly_area *= 0.5;
1290 92345 : return poly_area;
1291 : }
|
