https://mooseframework.inl.gov
Loading...
Searching...
No Matches
Classes | Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
MortarSegmentHelper Class Reference

This class supports defining mortar segment mesh elements in 3D by projecting secondary and primary elements onto a linearized plane, computing the overlapping polygon formed by their projections, and triangulating the resulting nodes. More...

#include <MortarSegmentHelper.h>

Classes

struct  ReferenceMappingData
 Output containers and filtering data used while generating reference-coordinate mappings. More...
 

Public Member Functions

 MortarSegmentHelper (std::vector< Point > secondary_nodes, const Point &center, const Point &normal, const MortarSegmentTriangulationMode triangulation_mode, const bool triangulate_triangles)
 Construct a helper that generates mortar segment geometry only.
 
 MortarSegmentHelper (std::vector< Point > secondary_nodes, std::vector< Point > secondary_reference_points, const Point &center, const Point &normal, const MortarSegmentTriangulationMode triangulation_mode, const bool triangulate_triangles)
 Construct a helper that also tracks secondary parent-reference coordinates.
 
Point getIntersection (const Point &p1, const Point &p2, const Point &q1, const Point &q2, Real &s) const
 Computes the intersection between line segments defined by point pairs (p1,p2) and (q1,q2) Also computes s, the ratio of distance between (p1,p2) that the intersection falls, quantity s is useful in avoiding adding nearly degenerate nodes.
 
bool isInsideSecondary (const Point &pt) const
 Check that a point is inside the secondary polygon (for verification only)
 
bool isDisjoint (const std::vector< Point > &poly) const
 Checks whether polygons are disjoint for an easy out.
 
std::vector< PointprojectPrimaryPoly (const std::vector< Point > &primary_nodes) const
 Project a primary polygon into the helper plane while preserving the clipping orientation.
 
std::vector< PointclipPoly (const std::vector< Point > &primary_nodes) const
 Clip secondary element (defined in instantiation) against given primary polygon result is a set of 2D nodes defining clipped polygon.
 
void triangulatePoly (std::vector< Point > &poly_nodes, std::vector< std::vector< unsigned int > > &tri_map) const
 Triangulate a polygon according to the configured mortar-segment triangulation mode.
 
void getMortarSegments (const std::vector< Point > &primary_nodes, std::vector< Point > &nodes, std::vector< std::vector< unsigned int > > &elem_to_nodes)
 Get mortar segments generated by a secondary and primary element pair.
 
void getMortarSegments (const std::vector< Point > &primary_nodes, const std::vector< Point > &primary_reference_points, std::vector< Point > &nodes, std::vector< std::vector< unsigned int > > &elem_to_nodes, std::vector< std::array< Point, 3 > > &elem_to_secondary_reference_points, std::vector< std::array< Point, 3 > > &elem_to_primary_reference_points, Real minimum_segment_area=0.)
 Get mortar segments and aligned parent-reference coordinates by appending to the output containers.
 
Real area (const std::vector< Point > &nodes) const
 Compute area of polygon.
 
const Pointcenter () const
 Get center point of secondary element.
 
const Pointnormal () const
 Get the unit normal of the projection plane.
 
Real remainder () const
 Get area fraction remaining after clipping against primary elements.
 
Point point (unsigned int i) const
 Get 3D position of node of linearized secondary element.
 

Private Member Functions

void getMortarSegmentsImpl (const std::vector< Point > &primary_nodes, std::vector< Point > &nodes, std::vector< std::vector< unsigned int > > &elem_to_nodes, ReferenceMappingData *reference_mapping)
 
std::vector< PointclipProjectedPoly (const std::vector< Point > &primary_poly) const
 Clip an already projected primary polygon against the secondary polygon.
 
std::optional< PointreferencePoint (const Point &point, const std::vector< Point > &poly, const std::vector< Point > &reference_points, std::string *failure_reason=nullptr) const
 Recover a parent-reference point from a projected sub-element map.
 

Private Attributes

Point _center
 Geometric center of secondary element.
 
Point _normal
 Unit normal of the plane used to project and clip the linearized secondary subpatch.
 
Point _u
 Vectors orthogonal to normal that span the plane projection will be performed on.
 
Point _v
 
Real _secondary_area
 Area of projected secondary element.
 
Real _remaining_area_fraction
 Fraction of area remaining after overlapping primary polygons clipped.
 
bool _debug
 
Real _tolerance = 1e-8
 Tolerance for intersection and clipping.
 
const MortarSegmentTriangulationMode _triangulation_mode
 Triangulation mode used for clipped polygons.
 
const bool _triangulate_triangles
 Whether already-triangular polygons should still be centroid-subdivided.
 
Real _area_tol
 Tolerance times secondary area for dimensional consistency.
 
Real _length_tol
 Tolerance times secondary area for dimensional consistency.
 
std::vector< Point_secondary_poly
 List of projected points on the linearized secondary element.
 
std::vector< Point_secondary_reference_points
 Parent reference points corresponding to _secondary_poly.
 

Friends

class MortarSegmentHelperTest
 

Detailed Description

This class supports defining mortar segment mesh elements in 3D by projecting secondary and primary elements onto a linearized plane, computing the overlapping polygon formed by their projections, and triangulating the resulting nodes.

Definition at line 35 of file MortarSegmentHelper.h.

Constructor & Destructor Documentation

◆ MortarSegmentHelper() [1/2]

MortarSegmentHelper::MortarSegmentHelper ( std::vector< Point secondary_nodes,
const Point center,
const Point normal,
const MortarSegmentTriangulationMode  triangulation_mode,
const bool  triangulate_triangles 
)

Construct a helper that generates mortar segment geometry only.

Definition at line 324 of file MortarSegmentHelper.C.

330 std::move(secondary_nodes), {}, center, normal, triangulation_mode, triangulate_triangles)
331{
332}
This class supports defining mortar segment mesh elements in 3D by projecting secondary and primary e...
const Point & normal() const
Get the unit normal of the projection plane.
const Point & center() const
Get center point of secondary element.

◆ MortarSegmentHelper() [2/2]

MortarSegmentHelper::MortarSegmentHelper ( std::vector< Point secondary_nodes,
std::vector< Point secondary_reference_points,
const Point center,
const Point normal,
const MortarSegmentTriangulationMode  triangulation_mode,
const bool  triangulate_triangles 
)

Construct a helper that also tracks secondary parent-reference coordinates.

Definition at line 334 of file MortarSegmentHelper.C.

340 : _center(center),
342 _debug(false),
343 _triangulation_mode(triangulation_mode),
344 _triangulate_triangles(triangulate_triangles),
345 _secondary_reference_points(std::move(secondary_reference_points))
346{
347 mooseAssert(_secondary_reference_points.empty() ||
348 secondary_nodes.size() == _secondary_reference_points.size(),
349 "Each projected secondary node needs one parent reference point.");
350
351 _secondary_poly.clear();
352 _secondary_poly.reserve(secondary_nodes.size());
353
354 // Get orientation of secondary poly
355 const Point e1 = secondary_nodes[0] - secondary_nodes[1];
356 const Point e2 = secondary_nodes[2] - secondary_nodes[1];
357 const Real orient = e2.cross(e1) * _normal;
358
359 // u and v define the tangent plane of the element (at center)
360 // Note we embed orientation into our transformation to make 2D poly always
361 // positively oriented
362 _u = _normal.cross(secondary_nodes[0] - center).unit();
363 _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
364
365 // Transform problem to 2D plane spanned by u and v
366 for (const auto & node : secondary_nodes)
367 {
368 Point pt = node - _center;
369 _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
370 }
371
372 // Half-plane clipping assumes that the projected secondary polygon has a unique orientation.
373 validateProjectedQuadrilateral(_secondary_poly, "secondary");
374
375 // Initialize area of secondary polygon
378
379 // Tolerance for quantities with area dimensions
381
382 // Tolerance for quantites with length dimensions
384}
Real _area_tol
Tolerance times secondary area for dimensional consistency.
Point _u
Vectors orthogonal to normal that span the plane projection will be performed on.
Point _center
Geometric center of secondary element.
Real _secondary_area
Area of projected secondary element.
std::vector< Point > _secondary_poly
List of projected points on the linearized secondary element.
const MortarSegmentTriangulationMode _triangulation_mode
Triangulation mode used for clipped polygons.
const bool _triangulate_triangles
Whether already-triangular polygons should still be centroid-subdivided.
Point _normal
Unit normal of the plane used to project and clip the linearized secondary subpatch.
Real _remaining_area_fraction
Fraction of area remaining after overlapping primary polygons clipped.
Real _length_tol
Tolerance times secondary area for dimensional consistency.
Real area(const std::vector< Point > &nodes) const
Compute area of polygon.
Real _tolerance
Tolerance for intersection and clipping.
std::vector< Point > _secondary_reference_points
Parent reference points corresponding to _secondary_poly.
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

Member Function Documentation

◆ area()

Real MortarSegmentHelper::area ( const std::vector< Point > &  nodes) const

Compute area of polygon.

Definition at line 1322 of file MortarSegmentHelper.C.

1323{
1324 Real poly_area = 0;
1325 for (auto i : index_range(nodes))
1326 poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
1327 nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
1328 poly_area *= 0.5;
1329 return poly_area;
1330}
auto index_range(const T &sizable)

Referenced by getMortarSegmentsImpl(), and MortarSegmentHelper().

◆ center()

const Point & MortarSegmentHelper::center ( ) const
inline

Get center point of secondary element.

Definition at line 130 of file MortarSegmentHelper.h.

130{ return _center; }

Referenced by MortarSegmentHelper().

◆ clipPoly()

std::vector< Point > MortarSegmentHelper::clipPoly ( const std::vector< Point > &  primary_nodes) const

Clip secondary element (defined in instantiation) against given primary polygon result is a set of 2D nodes defining clipped polygon.

Definition at line 482 of file MortarSegmentHelper.C.

483{
484 return clipProjectedPoly(projectPrimaryPoly(primary_nodes));
485}
std::vector< Point > projectPrimaryPoly(const std::vector< Point > &primary_nodes) const
Project a primary polygon into the helper plane while preserving the clipping orientation.
std::vector< Point > clipProjectedPoly(const std::vector< Point > &primary_poly) const
Clip an already projected primary polygon against the secondary polygon.

Referenced by getMortarSegmentsImpl().

◆ clipProjectedPoly()

std::vector< Point > MortarSegmentHelper::clipProjectedPoly ( const std::vector< Point > &  primary_poly) const
private

Clip an already projected primary polygon against the secondary polygon.

Keeping projection separate preserves the ordering shared with primary reference points.

Definition at line 488 of file MortarSegmentHelper.C.

489{
490 if (isDisjoint(primary_poly))
491 return {};
492
493 // Initialize clipped poly with secondary poly (secondary is target poly)
494 std::vector<Point> clipped_poly = _secondary_poly;
495
496 // Loop through clipping edges
497 for (auto i : index_range(primary_poly))
498 {
499 // If clipped poly trivial, return
500 if (clipped_poly.size() < 3)
501 {
502 clipped_poly.clear();
503 return clipped_poly;
504 }
505
506 // Set input poly to current clipped poly
507 std::vector<Point> input_poly(clipped_poly);
508 clipped_poly.clear();
509
510 // Get clipping edge
511 const Point & clip_pt1 = primary_poly[i];
512 const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
513 const Point edg = clip_pt2 - clip_pt1;
514 const Real cp = clip_pt2(0) * clip_pt1(1) - clip_pt2(1) * clip_pt1(0);
515
516 // Check if point is to the left of (or on) clip_edge
517 /*
518 * Note that use of tolerance here is to avoid degenerate case when lines are
519 * essentially on top of each other (common when meshes match across interface)
520 * since finding intersection is ill-conditioned in this case.
521 */
522 auto is_inside = [&edg, cp](const Point & pt, Real tol)
523 { return pt(0) * edg(1) - pt(1) * edg(0) + cp < tol; };
524
525 // Loop through edges of target polygon (with previous clippings already included)
526 for (auto j : index_range(input_poly))
527 {
528 // Get target edge
529 const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
530 const Point prev_pt = input_poly[j];
531
532 // TODO: Don't need to calculate both each loop
533 const bool is_current_inside = is_inside(curr_pt, _area_tol);
534 const bool is_previous_inside = is_inside(prev_pt, _area_tol);
535
536 if (is_current_inside)
537 {
538 if (!is_previous_inside)
539 {
540 Real s;
541 Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
542
543 /*
544 * s is the fraction of distance along clip poly edge that intersection lies
545 * It is used here to avoid degenerate polygon cases. For example, consider a
546 * case like:
547 * o
548 * | (inside)
549 * ------|------
550 * | (outside)
551 * when the distance is small (< 1e-7) we don't want to to add both the point
552 * and intersection. Also note that when distance on the scale of 1e-7,
553 * area on scale of 1e-14 so is insignificant if this results in dropping
554 * a tri (for example if next edge crosses again)
555 */
556 if (s < (1 - _tolerance))
557 clipped_poly.push_back(intersect);
558 }
559 clipped_poly.push_back(curr_pt);
560 }
561 else if (is_previous_inside)
562 {
563 Real s;
564 Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
565 if (s > _tolerance)
566 clipped_poly.push_back(intersect);
567 }
568 }
569 }
570
571 // Make sure final clipped poly is not trivial
572 if (clipped_poly.size() < 3)
573 {
574 clipped_poly.clear();
575 return clipped_poly;
576 }
577
578 // Clean up result by removing any duplicate nodes
579 std::vector<Point> cleaned_poly;
580 cleaned_poly.push_back(clipped_poly.back());
581 for (auto i : make_range(clipped_poly.size() - 1))
582 {
583 const Point prev_pt = cleaned_poly.back();
584 const Point curr_pt = clipped_poly[i];
585
586 // If points are sufficiently distanced, add to output
587 if ((curr_pt - prev_pt).norm() > _length_tol)
588 cleaned_poly.push_back(curr_pt);
589 }
590
591 mooseAssert(
592 cleaned_poly.size() <= 8,
593 "Our distributed mesh numbering scheme assumes that we have at most 8 nodes resulting from "
594 "clipping the projection of the primary sub-element onto the secondary sub-element");
595 return cleaned_poly;
596}
Point getIntersection(const Point &p1, const Point &p2, const Point &q1, const Point &q2, Real &s) const
Computes the intersection between line segments defined by point pairs (p1,p2) and (q1,...
bool isDisjoint(const std::vector< Point > &poly) const
Checks whether polygons are disjoint for an easy out.
auto norm(const T &a)
IntRange< T > make_range(T beg, T end)

Referenced by clipPoly(), and getMortarSegmentsImpl().

◆ getIntersection()

Point MortarSegmentHelper::getIntersection ( const Point p1,
const Point p2,
const Point q1,
const Point q2,
Real &  s 
) const

Computes the intersection between line segments defined by point pairs (p1,p2) and (q1,q2) Also computes s, the ratio of distance between (p1,p2) that the intersection falls, quantity s is useful in avoiding adding nearly degenerate nodes.

Definition at line 387 of file MortarSegmentHelper.C.

389{
390 const Point dp = p2 - p1;
391 const Point dq = q2 - q1;
392 const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
393 const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
394 const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
395 const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
396 s = -alpha * (cp1q2 - cp1q1 - cq1q2);
397
398 // Intersection should be between p1 and p2, if it's not (due to poor conditioning), simply
399 // move it to one of the end points
400 s = s > 1 ? 1. : s;
401 s = s < 0 ? 0. : s;
402 return p1 + s * dp;
403}

Referenced by clipProjectedPoly().

◆ getMortarSegments() [1/2]

void MortarSegmentHelper::getMortarSegments ( const std::vector< Point > &  primary_nodes,
const std::vector< Point > &  primary_reference_points,
std::vector< Point > &  nodes,
std::vector< std::vector< unsigned int > > &  elem_to_nodes,
std::vector< std::array< Point, 3 > > &  elem_to_secondary_reference_points,
std::vector< std::array< Point, 3 > > &  elem_to_primary_reference_points,
Real  minimum_segment_area = 0. 
)

Get mortar segments and aligned parent-reference coordinates by appending to the output containers.

Segments below minimum_segment_area are removed before mapping; retained vertices must map uniquely.

Definition at line 977 of file MortarSegmentHelper.C.

985{
986 ReferenceMappingData reference_mapping{primary_reference_points,
987 elem_to_secondary_reference_points,
988 elem_to_primary_reference_points,
989 minimum_segment_area};
990 getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, &reference_mapping);
991}
void getMortarSegmentsImpl(const std::vector< Point > &primary_nodes, std::vector< Point > &nodes, std::vector< std::vector< unsigned int > > &elem_to_nodes, ReferenceMappingData *reference_mapping)

◆ getMortarSegments() [2/2]

void MortarSegmentHelper::getMortarSegments ( const std::vector< Point > &  primary_nodes,
std::vector< Point > &  nodes,
std::vector< std::vector< unsigned int > > &  elem_to_nodes 
)

Get mortar segments generated by a secondary and primary element pair.

Parameters
primary_nodesList of primary element 3D nodes
Returns
nodes List of 3D mortar segment nodes
tri_map List of integer arrays defining which nodes belong to each mortar segment

Definition at line 969 of file MortarSegmentHelper.C.

972{
973 getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, nullptr);
974}

◆ getMortarSegmentsImpl()

void MortarSegmentHelper::getMortarSegmentsImpl ( const std::vector< Point > &  primary_nodes,
std::vector< Point > &  nodes,
std::vector< std::vector< unsigned int > > &  elem_to_nodes,
ReferenceMappingData reference_mapping 
)
private

Definition at line 994 of file MortarSegmentHelper.C.

998{
999 std::vector<Point> primary_poly;
1000 std::vector<Point> primary_poly_reference_points;
1001
1002 if (reference_mapping)
1003 {
1004 if (primary_nodes.size() != reference_mapping->primary_reference_points.size())
1005 mooseError("Reference-interpolation mortar segment generation requires one primary "
1006 "reference point per primary sub-element node.");
1007 if (_secondary_poly.size() != _secondary_reference_points.size())
1008 mooseError("Reference-interpolation mortar segment generation requires one secondary "
1009 "reference point per secondary sub-element node.");
1010 if (reference_mapping->elem_to_secondary_reference_points.size() != elem_to_nodes.size() ||
1011 reference_mapping->elem_to_primary_reference_points.size() != elem_to_nodes.size())
1012 mooseError("Reference-interpolation mortar segment outputs must be aligned before appending "
1013 "new segments.");
1014
1015 // Keep reference points in the projected polygon's orientation.
1016 const Point e1 = primary_nodes[0] - primary_nodes[1];
1017 const Point e2 = primary_nodes[2] - primary_nodes[1];
1018 const Real orient = e2.cross(e1) * _u.cross(_v);
1019 const auto n_verts = primary_nodes.size();
1020
1021 primary_poly = projectPrimaryPoly(primary_nodes);
1022 primary_poly_reference_points.reserve(reference_mapping->primary_reference_points.size());
1023 for (const auto n : index_range(primary_nodes))
1024 {
1025 const auto primary_node_index = (orient > 0) ? n : n_verts - 1 - n;
1026 primary_poly_reference_points.push_back(
1027 reference_mapping->primary_reference_points[primary_node_index]);
1028 }
1029 }
1030
1031 // Clip primary elem against secondary elem. Reference mode preserves the projected primary
1032 // ordering so its reference points remain aligned.
1033 std::vector<Point> clipped_poly =
1034 reference_mapping ? clipProjectedPoly(primary_poly) : clipPoly(primary_nodes);
1035 if (clipped_poly.size() < 3)
1036 return;
1037
1038 if (_debug)
1039 for (const auto & point : clipped_poly)
1041 mooseError("Clipped polygon not inside linearized secondary element");
1042
1043 // Compute area of clipped polygon, update remaining area fraction
1045
1046 // Triangulate clip polygon. tri_map indices are local to clipped_poly (starting at 0); we
1047 // shift them into the global node numbering after appending the polygon nodes below.
1048 std::vector<std::vector<unsigned int>> tri_map;
1049 triangulatePoly(clipped_poly, tri_map);
1050 if (reference_mapping && reference_mapping->minimum_segment_area > 0.)
1051 tri_map.erase(
1052 std::remove_if(tri_map.begin(),
1053 tri_map.end(),
1054 [&clipped_poly, reference_mapping](const std::vector<unsigned int> & tri)
1055 {
1056 mooseAssert(tri.size() == 3,
1057 "Mortar segment triangulation should only produce TRI3 maps.");
1058 return triangleAreaHelper(clipped_poly[tri[0]],
1059 clipped_poly[tri[1]],
1060 clipped_poly[tri[2]]) <
1061 reference_mapping->minimum_segment_area;
1062 }),
1063 tri_map.end());
1064 if (tri_map.empty())
1065 return;
1066
1067 std::vector<Point> secondary_node_reference_points;
1068 std::vector<Point> primary_node_reference_points;
1069 if (reference_mapping)
1070 {
1071 secondary_node_reference_points.reserve(clipped_poly.size());
1072 primary_node_reference_points.reserve(clipped_poly.size());
1073
1074 const auto recover_reference_point = [this](const Point & projected_point,
1075 const std::vector<Point> & poly,
1076 const std::vector<Point> & reference_points,
1077 const char * const parent_name,
1078 const std::size_t node_index)
1079 {
1080 std::string failure_reason;
1081 const auto reference_point =
1082 referencePoint(projected_point, poly, reference_points, &failure_reason);
1083 if (!reference_point)
1084 mooseError("Unable to recover the ",
1085 parent_name,
1086 " parent reference point for retained 3D mortar overlap vertex ",
1087 node_index,
1088 " at projected point ",
1089 projected_point,
1090 ". Reason: ",
1091 failure_reason,
1092 ". Reference interpolation does not fall back to normal projection.");
1093
1094 return *reference_point;
1095 };
1096
1097 for (const auto node_index : index_range(clipped_poly))
1098 {
1099 const auto & point = clipped_poly[node_index];
1100 secondary_node_reference_points.push_back(recover_reference_point(
1101 point, _secondary_poly, _secondary_reference_points, "secondary", node_index));
1102 primary_node_reference_points.push_back(recover_reference_point(
1103 point, primary_poly, primary_poly_reference_points, "primary", node_index));
1104 }
1105 }
1106
1107 // Transform clipped poly back to (linearized) 3d and append to list
1108 const auto offset = cast_int<unsigned int>(nodes.size());
1109 for (const auto & point : clipped_poly)
1110 nodes.emplace_back((point(0) * _u) + (point(1) * _v) + _center);
1111
1112 for (const auto & tri : tri_map)
1113 {
1114 std::vector<unsigned int> shifted_tri;
1115 shifted_tri.reserve(tri.size());
1116 for (const auto local_index : tri)
1117 shifted_tri.push_back(offset + local_index);
1118 elem_to_nodes.push_back(std::move(shifted_tri));
1119
1120 if (reference_mapping)
1121 {
1122 mooseAssert(tri.size() == 3, "Mortar segment triangulation should only produce TRI3 maps.");
1123 std::array<Point, 3> elem_secondary_reference_points;
1124 std::array<Point, 3> elem_primary_reference_points;
1125 for (const auto n : index_range(tri))
1126 {
1127 const auto local_node = tri[n];
1128 elem_secondary_reference_points[n] = secondary_node_reference_points[local_node];
1129 elem_primary_reference_points[n] = primary_node_reference_points[local_node];
1130 }
1131
1132 reference_mapping->elem_to_secondary_reference_points.push_back(
1133 elem_secondary_reference_points);
1134 reference_mapping->elem_to_primary_reference_points.push_back(elem_primary_reference_points);
1135 }
1136 }
1137}
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application.
Definition MooseError.h:311
std::vector< Point > clipPoly(const std::vector< Point > &primary_nodes) const
Clip secondary element (defined in instantiation) against given primary polygon result is a set of 2D...
Point point(unsigned int i) const
Get 3D position of node of linearized secondary element.
bool isInsideSecondary(const Point &pt) const
Check that a point is inside the secondary polygon (for verification only)
void triangulatePoly(std::vector< Point > &poly_nodes, std::vector< std::vector< unsigned int > > &tri_map) const
Triangulate a polygon according to the configured mortar-segment triangulation mode.
std::optional< Point > referencePoint(const Point &point, const std::vector< Point > &poly, const std::vector< Point > &reference_points, std::string *failure_reason=nullptr) const
Recover a parent-reference point from a projected sub-element map.
R poly(const C &c, const T x, const bool derivative=false)
Evaluate a polynomial with the coefficients c at x.
Definition MathUtils.h:238
if(subdm)

Referenced by getMortarSegments(), and getMortarSegments().

◆ isDisjoint()

bool MortarSegmentHelper::isDisjoint ( const std::vector< Point > &  poly) const

Checks whether polygons are disjoint for an easy out.

Definition at line 428 of file MortarSegmentHelper.C.

429{
430 for (auto i : index_range(_secondary_poly))
431 {
432 // Get edge to check
433 const Point & q1 = _secondary_poly[i];
434 const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
435 const Point edg = q2 - q1;
436 const Real cp = q2(0) * q1(1) - q2(1) * q1(0);
437
438 // If more optimization needed, could store these values for later
439 // Check if point is to the left of (or on) clip_edge
440 auto is_inside = [&edg, cp](Point & pt, Real tol)
441 { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
442
443 bool all_outside = true;
444 for (auto pt : poly)
445 if (is_inside(pt, _area_tol))
446 all_outside = false;
447
448 if (all_outside)
449 return true;
450 }
451 return false;
452}

Referenced by clipProjectedPoly().

◆ isInsideSecondary()

bool MortarSegmentHelper::isInsideSecondary ( const Point pt) const

Check that a point is inside the secondary polygon (for verification only)

Definition at line 406 of file MortarSegmentHelper.C.

407{
408 for (auto i : index_range(_secondary_poly))
409 {
410 const Point & q1 = _secondary_poly[i];
411 const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
412
413 const Point e1 = q2 - q1;
414 const Point e2 = pt - q1;
415
416 // If point corresponds to one of the secondary vertices, skip
417 if (e2.norm() < _tolerance)
418 return true;
419
420 const bool inside = (e1(0) * e2(1) - e1(1) * e2(0)) < _area_tol;
421 if (!inside)
422 return false;
423 }
424 return true;
425}
auto norm() const

Referenced by getMortarSegmentsImpl().

◆ normal()

const Point & MortarSegmentHelper::normal ( ) const
inline

Get the unit normal of the projection plane.

Definition at line 135 of file MortarSegmentHelper.h.

135{ return _normal; }

◆ point()

Point MortarSegmentHelper::point ( unsigned int  i) const
inline

Get 3D position of node of linearized secondary element.

Definition at line 145 of file MortarSegmentHelper.h.

146 {
147 return (_secondary_poly[i](0) * _u) + (_secondary_poly[i](1) * _v) + _center;
148 }

Referenced by getMortarSegmentsImpl(), and referencePoint().

◆ projectPrimaryPoly()

std::vector< Point > MortarSegmentHelper::projectPrimaryPoly ( const std::vector< Point > &  primary_nodes) const

Project a primary polygon into the helper plane while preserving the clipping orientation.

Definition at line 455 of file MortarSegmentHelper.C.

456{
457 // Check orientation of primary_poly
458 const Point e1 = primary_nodes[0] - primary_nodes[1];
459 const Point e2 = primary_nodes[2] - primary_nodes[1];
460
461 // Note we use u x v here instead of normal because it may be flipped if secondary elem was
462 // negatively oriented
463 const Real orient = e2.cross(e1) * _u.cross(_v);
464
465 // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
466 std::vector<Point> primary_poly;
467 const int n_verts = primary_nodes.size();
468 primary_poly.reserve(primary_nodes.size());
469 for (auto n : index_range(primary_nodes))
470 {
471 Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
472 primary_poly.emplace_back(pt * _u, pt * _v, 0.);
473 }
474
475 // Validate the primary clipping polygon once before its edges are used as half planes.
476 validateProjectedQuadrilateral(primary_poly, "primary");
477
478 return primary_poly;
479}

Referenced by clipPoly(), and getMortarSegmentsImpl().

◆ referencePoint()

std::optional< Point > MortarSegmentHelper::referencePoint ( const Point point,
const std::vector< Point > &  poly,
const std::vector< Point > &  reference_points,
std::string *  failure_reason = nullptr 
) const
private

Recover a parent-reference point from a projected sub-element map.

Returns
No value for invalid or out-of-tolerance maps.

Definition at line 1140 of file MortarSegmentHelper.C.

1144{
1145 mooseAssert(poly.size() == reference_points.size(),
1146 "Projected point and reference point containers should be the same size.");
1147
1148 if (failure_reason)
1149 failure_reason->clear();
1150
1151 const auto fail = [failure_reason](const std::string & reason) -> std::optional<Point>
1152 {
1153 if (failure_reason)
1154 *failure_reason = reason;
1155 return std::nullopt;
1156 };
1157
1159 return fail("the projected target point contains a non-finite coordinate");
1160
1161 for (const auto i : index_range(poly))
1162 {
1163 if (!MooseUtils::isFinitePoint(poly[i]))
1164 return fail("projected polygon vertex " + std::to_string(i) +
1165 " contains a non-finite coordinate");
1166 if (!MooseUtils::isFinitePoint(reference_points[i]))
1167 return fail("parent reference vertex " + std::to_string(i) +
1168 " contains a non-finite coordinate");
1169 }
1170
1171 if (poly.size() != 3 && poly.size() != 4)
1172 return fail("reference point recovery only supports triangular and quadrilateral mortar "
1173 "sub-elements, but received " +
1174 std::to_string(poly.size()) + " vertices");
1175
1176 Real minimum_edge_length = std::numeric_limits<Real>::max();
1177 Point local_origin;
1178 for (const auto & vertex : poly)
1179 local_origin += vertex;
1180 local_origin /= poly.size();
1181
1182 Real local_scale = 0.;
1183 for (const auto i : index_range(poly))
1184 {
1185 minimum_edge_length =
1186 std::min(minimum_edge_length, (poly[(i + 1) % poly.size()] - poly[i]).norm());
1187 local_scale = std::max(local_scale, (poly[i] - local_origin).norm());
1188 }
1189
1190 const Real singular_tolerance = 100. * std::numeric_limits<Real>::epsilon();
1191 if (!std::isfinite(local_scale) || local_scale <= singular_tolerance)
1192 return fail("the projected polygon has a zero local length scale");
1193 if (!std::isfinite(minimum_edge_length) ||
1194 minimum_edge_length / local_scale <= singular_tolerance)
1195 return fail("the projected polygon has a zero-length edge relative to its local scale");
1196
1197 std::vector<Point> normalized_poly;
1198 normalized_poly.reserve(poly.size());
1199 for (const auto & vertex : poly)
1200 normalized_poly.push_back((vertex - local_origin) / local_scale);
1201 const Point normalized_point = (point - local_origin) / local_scale;
1202 const Real reference_tolerance =
1203 std::max(mortar_reference_mapping_tolerance, _area_tol / (minimum_edge_length * local_scale));
1204 std::array<Node, 4> element_nodes;
1205 const FEType fe_type(FIRST, LAGRANGE);
1206 const auto recover_with_libmesh = [&](auto & element) -> std::optional<Point>
1207 {
1208 for (const auto i : index_range(normalized_poly))
1209 {
1210 element_nodes[i] = normalized_poly[i];
1211 element_nodes[i].set_id(i);
1212 element.set_node(i, &element_nodes[i]);
1213 }
1214
1215 if (!element.has_invertible_map(mortar_reference_mapping_tolerance))
1216 return fail("the projected sub-element map is degenerate or non-invertible");
1217
1218 if (element.type() == QUAD4)
1219 for (const auto corner : make_range(element.n_vertices()))
1220 {
1221 Point tangent_xi;
1222 Point tangent_eta;
1223 for (const auto node : make_range(element.n_nodes()))
1224 {
1225 tangent_xi += FEInterface::shape_deriv(
1226 fe_type, 0, &element, node, 0, element.master_point(corner)) *
1227 element.point(node);
1228 tangent_eta += FEInterface::shape_deriv(
1229 fe_type, 0, &element, node, 1, element.master_point(corner)) *
1230 element.point(node);
1231 }
1232
1233 const Real corner_jacobian = tangent_xi.cross(tangent_eta).norm();
1234 if (!std::isfinite(corner_jacobian) ||
1235 corner_jacobian <= mortar_reference_mapping_tolerance)
1236 return fail("the projected quadrilateral has a singular or ill-conditioned corner map");
1237 }
1238
1239 Point local_reference = FEMap::inverse_map(
1240 2, &element, normalized_point, mortar_reference_mapping_tolerance, false, false);
1241 if (!MooseUtils::isFinitePoint(local_reference))
1242 return fail("libMesh inverse_map produced a non-finite reference point");
1243
1244 const Real inverse_map_error =
1245 (FEMap::map(2, &element, local_reference) - normalized_point).norm();
1246 if (!std::isfinite(inverse_map_error) || inverse_map_error > mortar_reference_mapping_tolerance)
1247 {
1248 std::ostringstream reason;
1249 reason << "the normalized inverse-map error " << inverse_map_error << " exceeds "
1250 << mortar_reference_mapping_tolerance;
1251 return fail(reason.str());
1252 }
1253
1254 if (element.type() == TRI3)
1255 {
1256 std::array<Real, 3> weights;
1257 for (const auto i : index_range(weights))
1258 weights[i] = FEInterface::shape(fe_type, &element, i, local_reference, false);
1259
1260 for (auto & weight : weights)
1261 weight = std::clamp(weight, 0., 1.);
1262 const Real weight_sum = std::accumulate(weights.begin(), weights.end(), 0.);
1263 if (!std::isfinite(weight_sum) || weight_sum <= singular_tolerance)
1264 return fail("clamped triangle barycentric coordinates have a zero or non-finite sum");
1265 for (auto & weight : weights)
1266 weight /= weight_sum;
1267
1268 // Preserve partition of unity after clamping tolerance-sized violations.
1269 const auto corrected_weight =
1270 std::distance(weights.begin(), std::max_element(weights.begin(), weights.end()));
1271 weights[corrected_weight] = 1.;
1272 for (const auto i : index_range(weights))
1273 if (i != static_cast<unsigned int>(corrected_weight))
1274 weights[corrected_weight] -= weights[i];
1275
1276 local_reference = Point();
1277 for (const auto i : index_range(weights))
1278 local_reference += weights[i] * element.master_point(i);
1279 }
1280 else
1281 {
1282 local_reference(0) = std::clamp(local_reference(0), -1., 1.);
1283 local_reference(1) = std::clamp(local_reference(1), -1., 1.);
1284 local_reference(2) = 0.;
1285 }
1286
1287 if (!element.on_reference_element(local_reference, mortar_reference_mapping_tolerance))
1288 return fail("the clamped inverse-map result is outside the reference element");
1289
1290 const Real round_trip_error =
1291 (FEMap::map(2, &element, local_reference) - normalized_point).norm();
1292 if (!std::isfinite(round_trip_error) || round_trip_error > reference_tolerance)
1293 {
1294 std::ostringstream reason;
1295 reason << "the normalized inverse-map round-trip error " << round_trip_error
1296 << " exceeds the clipping-consistent tolerance " << reference_tolerance;
1297 return fail(reason.str());
1298 }
1299
1300 Point parent_reference;
1301 for (const auto i : index_range(normalized_poly))
1302 parent_reference +=
1303 FEInterface::shape(fe_type, &element, i, local_reference, false) * reference_points[i];
1304
1305 if (!MooseUtils::isFinitePoint(parent_reference))
1306 return fail("reference interpolation produced a non-finite parent reference point");
1307
1308 return parent_reference;
1309 };
1310
1311 if (poly.size() == 3)
1312 {
1313 Tri3 element;
1314 return recover_with_libmesh(element);
1315 }
1316
1317 Quad4 element;
1318 return recover_with_libmesh(element);
1319}
void ErrorVector unsigned int
static Real shape_deriv(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const unsigned int j, const Point &p)
static Point map(const unsigned int dim, const Elem *elem, const Point &reference_point)
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
T clamp(const T &x, T2 lowerlimit, T2 upperlimit)
Definition MathUtils.h:310
bool isFinitePoint(const Point &point)
Definition MooseUtils.C:60
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
const dof_id_type n_nodes

Referenced by getMortarSegmentsImpl().

◆ remainder()

Real MortarSegmentHelper::remainder ( ) const
inline

Get area fraction remaining after clipping against primary elements.

Definition at line 140 of file MortarSegmentHelper.h.

140{ return _remaining_area_fraction; }

◆ triangulatePoly()

void MortarSegmentHelper::triangulatePoly ( std::vector< Point > &  poly_nodes,
std::vector< std::vector< unsigned int > > &  tri_map 
) const

Triangulate a polygon according to the configured mortar-segment triangulation mode.

Parameters
poly_nodesList of 2D nodes defining polygon. May be augmented with extra interior nodes (e.g. centroid) by triangulation modes that require them; callers should append the result to their nodes list before applying the offset to tri_map.
tri_mapOutput triangle list expressed in indices local to poly_nodes (i.e. starting at 0). Callers are responsible for shifting these indices into the global node numbering.

Definition at line 599 of file MortarSegmentHelper.C.

601{
602 // tri_map is populated with triangle indices that are local to poly_nodes (starting at 0).
603 // Callers are responsible for shifting these indices into a global node numbering.
604 const auto polygon_centroid = [](const std::vector<Point> & polygon_nodes)
605 {
606 Point centroid(0);
607 Real double_area = 0;
608 for (const auto i : index_range(polygon_nodes))
609 {
610 const auto & a = polygon_nodes[i];
611 const auto & b = polygon_nodes[(i + 1) % polygon_nodes.size()];
612 const Real cross = a(0) * b(1) - b(0) * a(1);
613 double_area += cross;
614 centroid(0) += (a(0) + b(0)) * cross;
615 centroid(1) += (a(1) + b(1)) * cross;
616 }
617
618 if (std::abs(double_area) <= TOLERANCE)
619 {
620 for (const auto & node : polygon_nodes)
621 centroid += node;
622 centroid /= polygon_nodes.size();
623 return centroid;
624 }
625
626 centroid /= (3. * double_area);
627 centroid(2) = 0;
628 return centroid;
629 };
630
631 const auto append_triangle = [this, &poly_nodes, &tri_map](
632 const unsigned int a, const unsigned int b, const unsigned int c)
633 {
634 if (triangleAreaHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) <= _area_tol)
635 return false;
636
637 if (orient2dHelper(poly_nodes[a], poly_nodes[b], poly_nodes[c]) >= 0)
638 tri_map.push_back({a, b, c});
639 else
640 tri_map.push_back({a, c, b});
641
642 return true;
643 };
644
645 const auto point_in_triangle =
646 [this](const Point & p, const Point & a, const Point & b, const Point & c)
647 {
648 const Real ab = orient2dHelper(a, b, p);
649 const Real bc = orient2dHelper(b, c, p);
650 const Real ca = orient2dHelper(c, a, p);
651 return ab >= -_area_tol && bc >= -_area_tol && ca >= -_area_tol;
652 };
653
654 const auto min_triangle_angle = [](const Point & a, const Point & b, const Point & c)
655 {
656 const auto clamp_cos = [](Real value) { return std::max(-1., std::min(1., value)); };
657 const auto angle_at =
658 [&clamp_cos](const Point & vertex, const Point & point_one, const Point & point_two)
659 {
660 const Point edge_one = point_one - vertex;
661 const Point edge_two = point_two - vertex;
662 const Real denom = edge_one.norm() * edge_two.norm();
663 if (denom <= TOLERANCE)
664 return 0.;
665 return std::acos(clamp_cos((edge_one * edge_two) / denom));
666 };
667
668 return std::min({angle_at(a, b, c), angle_at(b, c, a), angle_at(c, a, b)});
669 };
670
671 const auto canonicalize_polygon = [this, &poly_nodes]()
672 {
673 if (poly_nodes.size() < 3)
674 return;
675
676 if (area(poly_nodes) < 0)
677 std::reverse(poly_nodes.begin(), poly_nodes.end());
678
679 bool changed = true;
680 while (changed && poly_nodes.size() > 3)
681 {
682 changed = false;
683 for (const auto i : index_range(poly_nodes))
684 {
685 const auto prev = (i + poly_nodes.size() - 1) % poly_nodes.size();
686 const auto next = (i + 1) % poly_nodes.size();
687 if ((poly_nodes[i] - poly_nodes[prev]).norm() <= _length_tol ||
688 (poly_nodes[next] - poly_nodes[i]).norm() <= _length_tol ||
689 triangleAreaHelper(poly_nodes[prev], poly_nodes[i], poly_nodes[next]) <= _area_tol)
690 {
691 poly_nodes.erase(poly_nodes.begin() + i);
692 changed = true;
693 break;
694 }
695 }
696 }
697
698 if (poly_nodes.size() >= 3 && area(poly_nodes) < 0)
699 std::reverse(poly_nodes.begin(), poly_nodes.end());
700 };
701
702 const auto triangulate_with_ear_clipping =
703 [this, &poly_nodes, &point_in_triangle, &min_triangle_angle](
704 const bool perform_delaunay_flips)
705 {
706 std::vector<std::array<unsigned int, 3>> triangles;
707 if (poly_nodes.size() < 3)
708 return triangles;
709
710 if (poly_nodes.size() == 3)
711 {
712 triangles.push_back(makeCCWTriangleHelper(poly_nodes, 0, 1, 2));
713 return triangles;
714 }
715
716 std::vector<unsigned int> remaining_vertices(poly_nodes.size());
717 std::iota(remaining_vertices.begin(), remaining_vertices.end(), 0);
718
719 while (remaining_vertices.size() > 3)
720 {
721 std::optional<std::size_t> best_position;
722 Real best_score = -std::numeric_limits<Real>::max();
723 Real best_area = -std::numeric_limits<Real>::max();
724
725 for (const auto position : index_range(remaining_vertices))
726 {
727 const auto prev_position =
728 (position + remaining_vertices.size() - 1) % remaining_vertices.size();
729 const auto next_position = (position + 1) % remaining_vertices.size();
730 const auto prev = remaining_vertices[prev_position];
731 const auto curr = remaining_vertices[position];
732 const auto next = remaining_vertices[next_position];
733
734 if (orient2dHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]) <= _area_tol)
735 continue;
736
737 bool contains_other_vertex = false;
738 for (const auto other : remaining_vertices)
739 {
740 if (other == prev || other == curr || other == next)
741 continue;
742
743 if (point_in_triangle(
744 poly_nodes[other], poly_nodes[prev], poly_nodes[curr], poly_nodes[next]))
745 {
746 contains_other_vertex = true;
747 break;
748 }
749 }
750
751 if (contains_other_vertex)
752 continue;
753
754 const Real candidate_score =
755 min_triangle_angle(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
756 const Real candidate_area =
757 triangleAreaHelper(poly_nodes[prev], poly_nodes[curr], poly_nodes[next]);
758 if (!best_position || candidate_score > best_score + TOLERANCE ||
759 (std::abs(candidate_score - best_score) <= TOLERANCE &&
760 candidate_area > best_area + _area_tol))
761 {
762 best_position = position;
763 best_score = candidate_score;
764 best_area = candidate_area;
765 }
766 }
767
768 if (!best_position)
769 {
770 std::vector<std::array<unsigned int, 3>> best_fan;
771 Real best_fan_score = -std::numeric_limits<Real>::max();
772 Real best_fan_area = -std::numeric_limits<Real>::max();
773
774 for (const auto root_position : index_range(remaining_vertices))
775 {
776 std::vector<std::array<unsigned int, 3>> candidate_fan;
777 Real candidate_score = std::numeric_limits<Real>::max();
778 Real candidate_area = std::numeric_limits<Real>::max();
779 bool valid_fan = true;
780 const auto root = remaining_vertices[root_position];
781
782 for (unsigned int step = 1; step + 1 < remaining_vertices.size(); ++step)
783 {
784 const auto next_position = (root_position + step) % remaining_vertices.size();
785 const auto following_position = (root_position + step + 1) % remaining_vertices.size();
786 const auto vertex_one = remaining_vertices[next_position];
787 const auto vertex_two = remaining_vertices[following_position];
788
789 if (orient2dHelper(poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]) <=
790 _area_tol)
791 {
792 valid_fan = false;
793 break;
794 }
795
796 candidate_fan.push_back(
797 makeCCWTriangleHelper(poly_nodes, root, vertex_one, vertex_two));
798 candidate_score =
799 std::min(candidate_score,
800 min_triangle_angle(
801 poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
802 candidate_area =
803 std::min(candidate_area,
804 triangleAreaHelper(
805 poly_nodes[root], poly_nodes[vertex_one], poly_nodes[vertex_two]));
806 }
807
808 if (!valid_fan || candidate_fan.empty())
809 continue;
810
811 if (candidate_score > best_fan_score + TOLERANCE ||
812 (std::abs(candidate_score - best_fan_score) <= TOLERANCE &&
813 candidate_area > best_fan_area + _area_tol))
814 {
815 best_fan = std::move(candidate_fan);
816 best_fan_score = candidate_score;
817 best_fan_area = candidate_area;
818 }
819 }
820
821 if (best_fan.empty())
822 for (unsigned int i = 1; i + 1 < remaining_vertices.size(); ++i)
823 best_fan.push_back(makeCCWTriangleHelper(poly_nodes,
824 remaining_vertices[0],
825 remaining_vertices[i],
826 remaining_vertices[i + 1]));
827
828 triangles.insert(triangles.end(), best_fan.begin(), best_fan.end());
829 break;
830 }
831
832 const auto prev_position =
833 (*best_position + remaining_vertices.size() - 1) % remaining_vertices.size();
834 const auto next_position = (*best_position + 1) % remaining_vertices.size();
835 triangles.push_back(makeCCWTriangleHelper(poly_nodes,
836 remaining_vertices[prev_position],
837 remaining_vertices[*best_position],
838 remaining_vertices[next_position]));
839 remaining_vertices.erase(remaining_vertices.begin() + *best_position);
840 }
841
842 if (remaining_vertices.size() == 3)
843 triangles.push_back(makeCCWTriangleHelper(
844 poly_nodes, remaining_vertices[0], remaining_vertices[1], remaining_vertices[2]));
845
846 if (!perform_delaunay_flips)
847 return triangles;
848
849 std::set<std::array<unsigned int, 2>> boundary_edges;
850 for (const auto i : index_range(poly_nodes))
851 boundary_edges.insert(canonicalEdgeHelper(i, (i + 1) % poly_nodes.size()));
852
853 performLocalDelaunayFlips(poly_nodes, boundary_edges, triangles);
854 return triangles;
855 };
856
857 const auto is_convex_polygon = [this](const std::vector<Point> & polygon_nodes)
858 {
859 if (polygon_nodes.size() <= 3)
860 return true;
861
862 for (const auto i : index_range(polygon_nodes))
863 {
864 const auto prev = (i + polygon_nodes.size() - 1) % polygon_nodes.size();
865 const auto next = (i + 1) % polygon_nodes.size();
866 if (orient2dHelper(polygon_nodes[prev], polygon_nodes[i], polygon_nodes[next]) <= _area_tol)
867 return false;
868 }
869
870 return true;
871 };
872
873 // Fewer than 3 nodes can't be triangulated
874 if (poly_nodes.size() < 3)
875 mooseError("Can't triangulate poly with fewer than 3 nodes");
876
877 // Legacy centroid path: when the default triangulation (centroid) is selected
878 // and triangle re-tessellation is not requested, reproduce the legacy
879 // algorithm byte-for-byte so existing mortar baselines remain valid.
880 // Uses the arithmetic mean of the vertices (not the area-weighted centroid),
881 // emits one triangle per polygon edge without degeneracy filtering, and skips
882 // the canonicalization pass which would drop near-degenerate vertices and
883 // perturb integration weights in downstream test baselines.
885 {
886 if (poly_nodes.size() == 3)
887 {
888 tri_map.push_back({0, 1, 2});
889 return;
890 }
891
892 const unsigned int n_verts = poly_nodes.size();
893 Point poly_center;
894 for (const auto & node : poly_nodes)
895 poly_center += node;
896 poly_center /= n_verts;
897
898 for (const auto i : make_range(n_verts))
899 tri_map.push_back({i, (i + 1) % n_verts, n_verts});
900
901 poly_nodes.push_back(poly_center);
902 return;
903 }
904
905 canonicalize_polygon();
906 if (poly_nodes.size() < 3)
907 return;
908
909 if (poly_nodes.size() == 3 && !_triangulate_triangles)
910 {
911 append_triangle(0, 1, 2);
912 return;
913 }
914
915 const bool force_triangle_centroid_split = _triangulate_triangles && poly_nodes.size() == 3;
916
918 !force_triangle_centroid_split)
919 {
920 const unsigned int n_verts = poly_nodes.size();
921 for (unsigned int i = 1; i + 1 < n_verts; ++i)
922 append_triangle(0, i, i + 1);
923 return;
924 }
925
927 !force_triangle_centroid_split)
928 {
929#if defined(LIBMESH_HAVE_TRIANGLE) || defined(LIBMESH_HAVE_POLY2TRI)
930 triangulateConstrainedDelaunayPolygon(poly_nodes, _area_tol, _length_tol, tri_map);
931 return;
932#else
933 mooseError("The 'delaunay' mortar triangulation mode requires libMesh TriangleInterface or "
934 "Poly2Tri support.");
935#endif
936 }
937
939 !force_triangle_centroid_split)
940 {
941 for (const auto & triangle : triangulate_with_ear_clipping(true))
942 append_triangle(triangle[0], triangle[1], triangle[2]);
943 return;
944 }
945
946 if (!force_triangle_centroid_split && !is_convex_polygon(poly_nodes))
947 {
948 for (const auto & triangle : triangulate_with_ear_clipping(true))
949 append_triangle(triangle[0], triangle[1], triangle[2]);
950 return;
951 }
952
953 const unsigned int n_verts = poly_nodes.size();
954 const Point poly_center = polygon_centroid(poly_nodes);
955
956 bool added_triangle = false;
957 for (const auto i : make_range(n_verts))
958 if (triangleAreaHelper(poly_nodes[i], poly_nodes[(i + 1) % n_verts], poly_center) > _area_tol)
959 {
960 tri_map.push_back({i, (i + 1) % n_verts, n_verts});
961 added_triangle = true;
962 }
963
964 if (added_triangle)
965 poly_nodes.push_back(poly_center);
966}
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)

Referenced by getMortarSegmentsImpl().

Friends And Related Symbol Documentation

◆ MortarSegmentHelperTest

friend class MortarSegmentHelperTest
friend

Definition at line 183 of file MortarSegmentHelper.h.

Member Data Documentation

◆ _area_tol

Real MortarSegmentHelper::_area_tol
private

Tolerance times secondary area for dimensional consistency.

Definition at line 233 of file MortarSegmentHelper.h.

Referenced by clipProjectedPoly(), isDisjoint(), isInsideSecondary(), MortarSegmentHelper(), referencePoint(), and triangulatePoly().

◆ _center

Point MortarSegmentHelper::_center
private

Geometric center of secondary element.

Definition at line 189 of file MortarSegmentHelper.h.

Referenced by center(), getMortarSegmentsImpl(), MortarSegmentHelper(), point(), and projectPrimaryPoly().

◆ _debug

bool MortarSegmentHelper::_debug
private

Definition at line 213 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl().

◆ _length_tol

Real MortarSegmentHelper::_length_tol
private

Tolerance times secondary area for dimensional consistency.

Definition at line 238 of file MortarSegmentHelper.h.

Referenced by clipProjectedPoly(), and MortarSegmentHelper().

◆ _normal

Point MortarSegmentHelper::_normal
private

Unit normal of the plane used to project and clip the linearized secondary subpatch.

Definition at line 194 of file MortarSegmentHelper.h.

Referenced by MortarSegmentHelper(), and normal().

◆ _remaining_area_fraction

Real MortarSegmentHelper::_remaining_area_fraction
private

Fraction of area remaining after overlapping primary polygons clipped.

Definition at line 211 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl(), MortarSegmentHelper(), and remainder().

◆ _secondary_area

Real MortarSegmentHelper::_secondary_area
private

Area of projected secondary element.

Definition at line 206 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl(), and MortarSegmentHelper().

◆ _secondary_poly

std::vector<Point> MortarSegmentHelper::_secondary_poly
private

List of projected points on the linearized secondary element.

Definition at line 243 of file MortarSegmentHelper.h.

Referenced by clipProjectedPoly(), getMortarSegmentsImpl(), isDisjoint(), isInsideSecondary(), MortarSegmentHelper(), and point().

◆ _secondary_reference_points

std::vector<Point> MortarSegmentHelper::_secondary_reference_points
private

Parent reference points corresponding to _secondary_poly.

Definition at line 246 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl(), and MortarSegmentHelper().

◆ _tolerance

Real MortarSegmentHelper::_tolerance = 1e-8
private

Tolerance for intersection and clipping.

Definition at line 218 of file MortarSegmentHelper.h.

Referenced by clipProjectedPoly(), isInsideSecondary(), and MortarSegmentHelper().

◆ _triangulate_triangles

const bool MortarSegmentHelper::_triangulate_triangles
private

Whether already-triangular polygons should still be centroid-subdivided.

Definition at line 228 of file MortarSegmentHelper.h.

◆ _triangulation_mode

const MortarSegmentTriangulationMode MortarSegmentHelper::_triangulation_mode
private

Triangulation mode used for clipped polygons.

Definition at line 223 of file MortarSegmentHelper.h.

◆ _u

Point MortarSegmentHelper::_u
private

Vectors orthogonal to normal that span the plane projection will be performed on.

These vectors are used to project the polygon clipping problem on a 2D plane, they are defined so the nodes of the projected polygon are listed with positive orientation

Definition at line 201 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl(), MortarSegmentHelper(), point(), and projectPrimaryPoly().

◆ _v

Point MortarSegmentHelper::_v
private

The documentation for this class was generated from the following files: