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>

Public Member Functions
	MortarSegmentHelper (const std::vector< Point > secondary_nodes, const Point &center, const Point &normal)

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. More...

bool	isInsideSecondary (const Point &pt) const
	Check that a point is inside the secondary polygon (for verification only) More...

bool	isDisjoint (const std::vector< Point > &poly) const
	Checks whether polygons are disjoint for an easy out. More...

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. More...

void	triangulatePoly (std::vector< Point > &poly_nodes, const unsigned int offset, std::vector< std::vector< unsigned int >> &tri_map) const
	Triangulate a polygon (currently uses center of polygon to define triangulation) More...

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. More...

Real	area (const std::vector< Point > &nodes) const
	Compute area of polygon. More...

const Point &	center () const
	Get center point of secondary element. More...

Real	remainder () const
	Get area fraction remaining after clipping against primary elements. More...

Point	point (unsigned int i) const
	Get 3D position of node of linearized secondary element. More...

Private Attributes
Point	_center
	Geometric center of secondary element. More...

Point	_normal
	Normal at geometric center of secondary element. More...

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

Point	_v

Real	_secondary_area
	Area of projected secondary element. More...

Real	_remaining_area_fraction
	Fraction of area remaining after overlapping primary polygons clipped. More...

bool	_debug

Real	_tolerance = 1e-8
	Tolerance for intersection and clipping. More...

Real	_area_tol
	Tolerance times secondary area for dimensional consistency. More...

Real	_length_tol
	Tolerance times secondary area for dimensional consistency. More...

std::vector< Point >	_secondary_poly
	List of projected points on the linearized secondary element. More...

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 26 of file MortarSegmentHelper.h.

Constructor & Destructor Documentation

◆ MortarSegmentHelper()

MortarSegmentHelper::MortarSegmentHelper	(	const std::vector< Point >	secondary_nodes,
		const Point &	center,
		const Point &	normal
	)

Definition at line 16 of file MortarSegmentHelper.C.

   : _center(center), _normal(normal), _debug(false)
 {
   _secondary_poly.clear();
   _secondary_poly.reserve(secondary_nodes.size());
 
   // Get orientation of secondary poly
   const Point e1 = secondary_nodes[0] - secondary_nodes[1];
   const Point e2 = secondary_nodes[2] - secondary_nodes[1];
   const Real orient = e2.cross(e1) * _normal;
 
   // u and v define the tangent plane of the element (at center)
   // Note we embed orientation into our transformation to make 2D poly always
   // positively oriented
   _u = _normal.cross(secondary_nodes[0] - center).unit();
   _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
 
   // Transform problem to 2D plane spanned by u and v
   for (const auto & node : secondary_nodes)
   {
     Point pt = node - _center;
     _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
   }
 
   // Initialize area of secondary polygon
   _remaining_area_fraction = 1.0;
   _secondary_area = area(_secondary_poly);
 
   // Tolerance for quantities with area dimensions
   _area_tol = _tolerance * _secondary_area;
 
   // Tolerance for quantites with length dimensions
   _length_tol = _tolerance * std::sqrt(_secondary_area);
 }

Member Function Documentation

◆ area()

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

Compute area of polygon.

Definition at line 320 of file MortarSegmentHelper.C.

Referenced by getMortarSegments(), and MortarSegmentHelper().

 {
   Real poly_area = 0;
   for (auto i : index_range(nodes))
     poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
                  nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
   poly_area *= 0.5;
   return poly_area;
 }

◆ center()

const Point& MortarSegmentHelper::center ( ) const

inline

Get center point of secondary element.

Definition at line 85 of file MortarSegmentHelper.h.

Referenced by MortarSegmentHelper().

85 { return _center; }

MortarSegmentHelper::_center

Point _center

Geometric center of secondary element.

Definition: MortarSegmentHelper.h:104

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

Referenced by getMortarSegments().

 {
   // Check orientation of primary_poly
   const Point e1 = primary_nodes[0] - primary_nodes[1];
   const Point e2 = primary_nodes[2] - primary_nodes[1];
 
   // Note we use u x v here instead of normal because it may be flipped if secondary elem was
   // negatively oriented
   const Real orient = e2.cross(e1) * _u.cross(_v);
 
   // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
   std::vector<Point> primary_poly;
   const int n_verts = primary_nodes.size();
   for (auto n : index_range(primary_nodes))
   {
     Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
     primary_poly.emplace_back(pt * _u, pt * _v, 0.);
   }
 
   if (isDisjoint(primary_poly))
   {
     primary_poly.clear();
     return primary_poly;
   }
 
   // Initialize clipped poly with secondary poly (secondary is target poly)
   std::vector<Point> clipped_poly = _secondary_poly;
 
   // Loop through clipping edges
   for (auto i : index_range(primary_poly))
   {
     // If clipped poly trivial, return
     if (clipped_poly.size() < 3)
     {
       clipped_poly.clear();
       return clipped_poly;
     }
 
     // Set input poly to current clipped poly
     std::vector<Point> input_poly(clipped_poly);
     clipped_poly.clear();
 
     // Get clipping edge
     const Point & clip_pt1 = primary_poly[i];
     const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
     const Point edg = clip_pt2 - clip_pt1;
     const Real cp = clip_pt2(0) * clip_pt1(1) - clip_pt2(1) * clip_pt1(0);
 
     // Check if point is to the left of (or on) clip_edge
     /*
      * Note that use of tolerance here is to avoid degenerate case when lines are
      * essentially on top of each other (common when meshes match across interface)
      * since finding intersection is ill-conditioned in this case.
      */
     auto is_inside = [&edg, cp](const Point & pt, Real tol)
     { return pt(0) * edg(1) - pt(1) * edg(0) + cp < tol; };
 
     // Loop through edges of target polygon (with previous clippings already included)
     for (auto j : index_range(input_poly))
     {
       // Get target edge
       const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
       const Point prev_pt = input_poly[j];
 
       // TODO: Don't need to calculate both each loop
       const bool is_current_inside = is_inside(curr_pt, _area_tol);
       const bool is_previous_inside = is_inside(prev_pt, _area_tol);
 
       if (is_current_inside)
       {
         if (!is_previous_inside)
         {
           Real s;
           Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
 
           /*
            * s is the fraction of distance along clip poly edge that intersection lies
            * It is used here to avoid degenerate polygon cases. For example, consider a
            * case like:
            *          o
            *          |    (inside)
            *    ------|------
            *          |    (outside)
            * when the distance is small (< 1e-7) we don't want to to add both the point
            * and intersection. Also note that when distance on the scale of 1e-7,
            * area on scale of 1e-14 so is insignificant if this results in dropping
            * a tri (for example if next edge crosses again)
            */
           if (s < (1 - _tolerance))
             clipped_poly.push_back(intersect);
         }
         clipped_poly.push_back(curr_pt);
       }
       else if (is_previous_inside)
       {
         Real s;
         Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
         if (s > _tolerance)
           clipped_poly.push_back(intersect);
       }
     }
   }
 
   // return clipped_poly;
   // Make sure final clipped poly is not trivial
   if (clipped_poly.size() < 3)
   {
     clipped_poly.clear();
     return clipped_poly;
   }
 
   // Clean up result by removing any duplicate nodes
   std::vector<Point> cleaned_poly;
 
   cleaned_poly.push_back(clipped_poly.back());
   for (auto i : make_range(clipped_poly.size() - 1))
   {
     const Point prev_pt = cleaned_poly.back();
     const Point curr_pt = clipped_poly[i];
 
     // If points are sufficiently distanced, add to output
     if ((curr_pt - prev_pt).norm() > _length_tol)
       cleaned_poly.push_back(curr_pt);
   }
 
   return cleaned_poly;
 }

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

Referenced by clipPoly().

 {
   const Point dp = p2 - p1;
   const Point dq = q2 - q1;
   const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
   const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
   const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
   const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
   s = -alpha * (cp1q2 - cp1q1 - cq1q2);
 
   // Intersection should be between p1 and p2, if it's not (due to poor conditioning), simply
   // move it to one of the end points
   s = s > 1 ? 1. : s;
   s = s < 0 ? 0. : s;
   return p1 + s * dp;
 }

◆ getMortarSegments()

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_nodes List 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 294 of file MortarSegmentHelper.C.

 {
   // Clip primary elem against secondary elem
   std::vector<Point> clipped_poly = clipPoly(primary_nodes);
   if (clipped_poly.size() < 3)
     return;
 
   if (_debug)
     for (auto pt : clipped_poly)
       if (!isInsideSecondary(pt))
         mooseError("Clipped polygon not inside linearized secondary element");
 
   // Compute area of clipped polygon, update remaining area fraction
   _remaining_area_fraction -= area(clipped_poly) / _secondary_area;
 
   // Triangulate clip polygon
   triangulatePoly(clipped_poly, nodes.size(), elem_to_nodes);
 
   // Transform clipped poly back to (linearized) 3d and append to list
   for (auto pt : clipped_poly)
     nodes.emplace_back((pt(0) * _u) + (pt(1) * _v) + _center);
 }

◆ isDisjoint()

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

Checks whether polygons are disjoint for an easy out.

Definition at line 95 of file MortarSegmentHelper.C.

Referenced by clipPoly().

 {
   for (auto i : index_range(_secondary_poly))
   {
     // Get edge to check
     const Point & q1 = _secondary_poly[i];
     const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
     const Point edg = q2 - q1;
     const Real cp = q2(0) * q1(1) - q2(1) * q1(0);
 
     // If more optimization needed, could store these values for later
     // Check if point is to the left of (or on) clip_edge
     auto is_inside = [&edg, cp](Point & pt, Real tol)
     { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
 
     bool all_outside = true;
     for (auto pt : poly)
       if (is_inside(pt, _area_tol))
         all_outside = false;
 
     if (all_outside)
       return true;
   }
   return false;
 }

◆ isInsideSecondary()

bool MortarSegmentHelper::isInsideSecondary ( const Point & pt ) const

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

Definition at line 73 of file MortarSegmentHelper.C.

Referenced by getMortarSegments().

 {
   for (auto i : index_range(_secondary_poly))
   {
     const Point & q1 = _secondary_poly[i];
     const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
 
     const Point e1 = q2 - q1;
     const Point e2 = pt - q1;
 
     // If point corresponds to one of the secondary vertices, skip
     if (e2.norm() < _tolerance)
       return true;
 
     bool inside = (e1(0) * e2(1) - e1(1) * e2(0)) < _area_tol;
     if (!inside)
       return false;
   }
   return true;
 }

◆ point()

Point MortarSegmentHelper::point ( unsigned int i ) const

inline

Get 3D position of node of linearized secondary element.

Definition at line 95 of file MortarSegmentHelper.h.

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

◆ remainder()

Real MortarSegmentHelper::remainder ( ) const

inline

Get area fraction remaining after clipping against primary elements.

Definition at line 90 of file MortarSegmentHelper.h.

90 { return _remaining_area_fraction; }

MortarSegmentHelper::_remaining_area_fraction

Real _remaining_area_fraction

Fraction of area remaining after overlapping primary polygons clipped.

Definition: MortarSegmentHelper.h:126

◆ triangulatePoly()

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

Triangulate a polygon (currently uses center of polygon to define triangulation)

Parameters

poly_nodes	List of 2D nodes defining polygon
offset	Current size of 3D nodes array (not poly_nodes)

Returns: tri_map List of integer arrays defining which nodes belong to each triangle

Definition at line 251 of file MortarSegmentHelper.C.

Referenced by getMortarSegments().

 {
   // Note offset is important because it converts from list of nodes on local poly to list of
   // nodes on entire element. This is a sloppy and error-prone way to do this.
   // Could be redesigned by just building and adding elements inside helper, but leaving for now
   // to keep buildMortarSegmentMesh similar for 2D and 3D
 
   // Fewer than 3 nodes can't be triangulated
   if (poly_nodes.size() < 3)
     mooseError("Can't triangulate poly with fewer than 3 nodes");
   // If 3 nodes, already a triangle, simply pass back map
   else if (poly_nodes.size() == 3)
   {
     tri_map.push_back({0 + offset, 1 + offset, 2 + offset});
     return;
   }
   // Otherwise use simple center point triangulation
   // Note: Could use better algorithm to reduce number of triangles
   //    - Delaunay produces good quality triangulation but might be a bit overkill
   //    - Monotone polygon triangulation algorithm is O(N) but quality is not guaranteed
   else
   {
     Point poly_center;
     const unsigned int n_verts = poly_nodes.size();
 
     // Get geometric center of polygon
     for (auto node : poly_nodes)
       poly_center += node;
     poly_center /= n_verts;
 
     // Add triangles formed by outer edge and center point
     for (auto i : make_range(n_verts))
       tri_map.push_back({i + offset, (i + 1) % n_verts + offset, n_verts + offset});
 
     // Add center point to list of polygon nodes
     poly_nodes.push_back(poly_center);
     return;
   }
 }

Member Data Documentation

◆ _area_tol

Real MortarSegmentHelper::_area_tol

private

Tolerance times secondary area for dimensional consistency.

Definition at line 138 of file MortarSegmentHelper.h.

Referenced by clipPoly(), isDisjoint(), isInsideSecondary(), and MortarSegmentHelper().

◆ _center

Point MortarSegmentHelper::_center

private

Geometric center of secondary element.

Definition at line 104 of file MortarSegmentHelper.h.

Referenced by center(), clipPoly(), getMortarSegments(), MortarSegmentHelper(), and point().

◆ _debug

bool MortarSegmentHelper::_debug

private

Definition at line 128 of file MortarSegmentHelper.h.

Referenced by getMortarSegments().

◆ _length_tol

Real MortarSegmentHelper::_length_tol

private

Tolerance times secondary area for dimensional consistency.

Definition at line 143 of file MortarSegmentHelper.h.

Referenced by clipPoly(), and MortarSegmentHelper().

◆ _normal

Point MortarSegmentHelper::_normal

private

Normal at geometric center of secondary element.

Definition at line 109 of file MortarSegmentHelper.h.

Referenced by MortarSegmentHelper().

◆ _remaining_area_fraction

Real MortarSegmentHelper::_remaining_area_fraction

private

Fraction of area remaining after overlapping primary polygons clipped.

Definition at line 126 of file MortarSegmentHelper.h.

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

◆ _secondary_area

Real MortarSegmentHelper::_secondary_area

private

Area of projected secondary element.

Definition at line 121 of file MortarSegmentHelper.h.

Referenced by getMortarSegments(), and MortarSegmentHelper().

◆ _secondary_poly

std::vector<Point> MortarSegmentHelper::_secondary_poly

private

List of projected points on the linearized secondary element.

Definition at line 148 of file MortarSegmentHelper.h.

Referenced by clipPoly(), isDisjoint(), isInsideSecondary(), MortarSegmentHelper(), and point().

◆ _tolerance

Real MortarSegmentHelper::_tolerance = 1e-8

private

Tolerance for intersection and clipping.

Definition at line 133 of file MortarSegmentHelper.h.

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

◆ _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 116 of file MortarSegmentHelper.h.

Referenced by clipPoly(), getMortarSegments(), MortarSegmentHelper(), and point().

◆ _v

Point MortarSegmentHelper::_v

private

Definition at line 116 of file MortarSegmentHelper.h.

Referenced by clipPoly(), getMortarSegments(), MortarSegmentHelper(), and point().

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

include/constraints/MortarSegmentHelper.h
src/constraints/MortarSegmentHelper.C

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MortarSegmentHelper()

Member Function Documentation

◆ area()

◆ center()

◆ clipPoly()

◆ getIntersection()

◆ getMortarSegments()

◆ isDisjoint()

◆ isInsideSecondary()

◆ point()

◆ remainder()

◆ triangulatePoly()

Member Data Documentation

◆ _area_tol

◆ _center

◆ _debug

◆ _length_tol

◆ _normal

◆ _remaining_area_fraction

◆ _secondary_area

◆ _secondary_poly

◆ _tolerance

◆ _u

◆ _v