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< 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 > 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.
 
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 Point & center () const
 Get center point of secondary element.
 
const Point & normal () 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< Point > clipProjectedPoly (const std::vector< Point > &primary_poly) const
 Clip an already projected primary polygon against the secondary polygon.
 
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.
 

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 32 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 322 of file MortarSegmentHelper.C.

328 std::move(secondary_nodes), {}, center, normal, triangulation_mode, triangulate_triangles)
329{
330}
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 332 of file MortarSegmentHelper.C.

338 : _center(center),
340 _debug(false),
341 _triangulation_mode(triangulation_mode),
342 _triangulate_triangles(triangulate_triangles),
343 _secondary_reference_points(std::move(secondary_reference_points))
344{
345 mooseAssert(_secondary_reference_points.empty() ||
346 secondary_nodes.size() == _secondary_reference_points.size(),
347 "Each projected secondary node needs one parent reference point.");
348
349 _secondary_poly.clear();
350 _secondary_poly.reserve(secondary_nodes.size());
351
352 // Get orientation of secondary poly
353 const Point e1 = secondary_nodes[0] - secondary_nodes[1];
354 const Point e2 = secondary_nodes[2] - secondary_nodes[1];
355 const Real orient = e2.cross(e1) * _normal;
356
357 // u and v define the tangent plane of the element (at center)
358 // Note we embed orientation into our transformation to make 2D poly always
359 // positively oriented
360 _u = _normal.cross(secondary_nodes[0] - center).unit();
361 _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
362
363 // Transform problem to 2D plane spanned by u and v
364 for (const auto & node : secondary_nodes)
365 {
366 Point pt = node - _center;
367 _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
368 }
369
370 // Half-plane clipping assumes that the projected secondary polygon has a unique orientation.
371 validateProjectedQuadrilateral(_secondary_poly, "secondary");
372
373 // Initialize area of secondary polygon
376
377 // Tolerance for quantities with area dimensions
379
380 // Tolerance for quantites with length dimensions
382}
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.
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 1320 of file MortarSegmentHelper.C.

1321{
1322 Real poly_area = 0;
1323 for (auto i : index_range(nodes))
1324 poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
1325 nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
1326 poly_area *= 0.5;
1327 return poly_area;
1328}
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 127 of file MortarSegmentHelper.h.

127{ 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 480 of file MortarSegmentHelper.C.

481{
482 return clipProjectedPoly(projectPrimaryPoly(primary_nodes));
483}
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 486 of file MortarSegmentHelper.C.

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

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

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 975 of file MortarSegmentHelper.C.

983{
984 ReferenceMappingData reference_mapping{primary_reference_points,
985 elem_to_secondary_reference_points,
986 elem_to_primary_reference_points,
987 minimum_segment_area};
988 getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, &reference_mapping);
989}
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 967 of file MortarSegmentHelper.C.

970{
971 getMortarSegmentsImpl(primary_nodes, nodes, elem_to_nodes, nullptr);
972}

◆ 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 992 of file MortarSegmentHelper.C.

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

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 426 of file MortarSegmentHelper.C.

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

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 404 of file MortarSegmentHelper.C.

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

Referenced by getMortarSegmentsImpl().

◆ normal()

const Point & MortarSegmentHelper::normal ( ) const
inline

Get the unit normal of the projection plane.

Definition at line 132 of file MortarSegmentHelper.h.

132{ return _normal; }

◆ point()

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

Get 3D position of node of linearized secondary element.

Definition at line 142 of file MortarSegmentHelper.h.

143 {
144 return (_secondary_poly[i](0) * _u) + (_secondary_poly[i](1) * _v) + _center;
145 }

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 453 of file MortarSegmentHelper.C.

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

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 1138 of file MortarSegmentHelper.C.

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

137{ 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 597 of file MortarSegmentHelper.C.

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

Member Data Documentation

◆ _area_tol

Real MortarSegmentHelper::_area_tol
private

Tolerance times secondary area for dimensional consistency.

Definition at line 230 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 186 of file MortarSegmentHelper.h.

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

◆ _debug

bool MortarSegmentHelper::_debug
private

Definition at line 210 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl().

◆ _length_tol

Real MortarSegmentHelper::_length_tol
private

Tolerance times secondary area for dimensional consistency.

Definition at line 235 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 191 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 208 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 203 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 240 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 243 of file MortarSegmentHelper.h.

Referenced by getMortarSegmentsImpl(), and MortarSegmentHelper().

◆ _tolerance

Real MortarSegmentHelper::_tolerance = 1e-8
private

Tolerance for intersection and clipping.

Definition at line 215 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 225 of file MortarSegmentHelper.h.

◆ _triangulation_mode

const MortarSegmentTriangulationMode MortarSegmentHelper::_triangulation_mode
private

Triangulation mode used for clipped polygons.

Definition at line 220 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 198 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: