Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 : #include "MortarSegmentHelper.h"
10 : #include "MooseError.h"
11 :
12 : #include "libmesh/int_range.h"
13 :
14 : using namespace libMesh;
15 :
16 5590 : MortarSegmentHelper::MortarSegmentHelper(const std::vector<Point> secondary_nodes,
17 : const Point & center,
18 5590 : const Point & normal)
19 5590 : : _center(center), _normal(normal), _debug(false)
20 : {
21 5590 : _secondary_poly.clear();
22 5590 : _secondary_poly.reserve(secondary_nodes.size());
23 :
24 : // Get orientation of secondary poly
25 5590 : const Point e1 = secondary_nodes[0] - secondary_nodes[1];
26 5590 : const Point e2 = secondary_nodes[2] - secondary_nodes[1];
27 5590 : const Real orient = e2.cross(e1) * _normal;
28 :
29 : // u and v define the tangent plane of the element (at center)
30 : // Note we embed orientation into our transformation to make 2D poly always
31 : // positively oriented
32 5590 : _u = _normal.cross(secondary_nodes[0] - center).unit();
33 5590 : _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
34 :
35 : // Transform problem to 2D plane spanned by u and v
36 24718 : for (const auto & node : secondary_nodes)
37 : {
38 19128 : Point pt = node - _center;
39 19128 : _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
40 : }
41 :
42 : // Initialize area of secondary polygon
43 5590 : _remaining_area_fraction = 1.0;
44 5590 : _secondary_area = area(_secondary_poly);
45 :
46 : // Tolerance for quantities with area dimensions
47 5590 : _area_tol = _tolerance * _secondary_area;
48 :
49 : // Tolerance for quantites with length dimensions
50 5590 : _length_tol = _tolerance * std::sqrt(_secondary_area);
51 5590 : }
52 :
53 : Point
54 78278 : MortarSegmentHelper::getIntersection(
55 : const Point & p1, const Point & p2, const Point & q1, const Point & q2, Real & s) const
56 : {
57 78278 : const Point dp = p2 - p1;
58 78278 : const Point dq = q2 - q1;
59 78278 : const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
60 78278 : const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
61 78278 : const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
62 78278 : const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
63 78278 : s = -alpha * (cp1q2 - cp1q1 - cq1q2);
64 :
65 : // Intersection should be between p1 and p2, if it's not (due to poor conditioning), simply
66 : // move it to one of the end points
67 78278 : s = s > 1 ? 1. : s;
68 78278 : s = s < 0 ? 0. : s;
69 78278 : return p1 + s * dp;
70 : }
71 :
72 : bool
73 0 : MortarSegmentHelper::isInsideSecondary(const Point & pt) const
74 : {
75 0 : for (auto i : index_range(_secondary_poly))
76 : {
77 0 : const Point & q1 = _secondary_poly[i];
78 0 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
79 :
80 0 : const Point e1 = q2 - q1;
81 0 : const Point e2 = pt - q1;
82 :
83 : // If point corresponds to one of the secondary vertices, skip
84 0 : if (e2.norm() < _tolerance)
85 0 : return true;
86 :
87 0 : bool inside = (e1(0) * e2(1) - e1(1) * e2(0)) < _area_tol;
88 0 : if (!inside)
89 0 : return false;
90 : }
91 0 : return true;
92 : }
93 :
94 : bool
95 165487 : MortarSegmentHelper::isDisjoint(const std::vector<Point> & poly) const
96 : {
97 383398 : for (auto i : index_range(_secondary_poly))
98 : {
99 : // Get edge to check
100 360391 : const Point & q1 = _secondary_poly[i];
101 360391 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
102 360391 : const Point edg = q2 - q1;
103 360391 : const Real cp = q2(0) * q1(1) - q2(1) * q1(0);
104 :
105 : // If more optimization needed, could store these values for later
106 : // Check if point is to the left of (or on) clip_edge
107 1205981 : auto is_inside = [&edg, cp](Point & pt, Real tol)
108 1205981 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
109 :
110 360391 : bool all_outside = true;
111 1566372 : for (auto pt : poly)
112 1205981 : if (is_inside(pt, _area_tol))
113 592706 : all_outside = false;
114 :
115 360391 : if (all_outside)
116 142480 : return true;
117 : }
118 23007 : return false;
119 : }
120 :
121 : std::vector<Point>
122 165487 : MortarSegmentHelper::clipPoly(const std::vector<Point> & primary_nodes) const
123 : {
124 : // Check orientation of primary_poly
125 165487 : const Point e1 = primary_nodes[0] - primary_nodes[1];
126 165487 : const Point e2 = primary_nodes[2] - primary_nodes[1];
127 :
128 : // Note we use u x v here instead of normal because it may be flipped if secondary elem was
129 : // negatively oriented
130 165487 : const Real orient = e2.cross(e1) * _u.cross(_v);
131 :
132 : // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
133 165487 : std::vector<Point> primary_poly;
134 165487 : const int n_verts = primary_nodes.size();
135 712715 : for (auto n : index_range(primary_nodes))
136 : {
137 547228 : Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
138 547228 : primary_poly.emplace_back(pt * _u, pt * _v, 0.);
139 : }
140 :
141 165487 : if (isDisjoint(primary_poly))
142 : {
143 142480 : primary_poly.clear();
144 142480 : return primary_poly;
145 : }
146 :
147 : // Initialize clipped poly with secondary poly (secondary is target poly)
148 23007 : std::vector<Point> clipped_poly = _secondary_poly;
149 :
150 : // Loop through clipping edges
151 99023 : for (auto i : index_range(primary_poly))
152 : {
153 : // If clipped poly trivial, return
154 78157 : if (clipped_poly.size() < 3)
155 : {
156 2141 : clipped_poly.clear();
157 2141 : return clipped_poly;
158 : }
159 :
160 : // Set input poly to current clipped poly
161 76016 : std::vector<Point> input_poly(clipped_poly);
162 76016 : clipped_poly.clear();
163 :
164 : // Get clipping edge
165 76016 : const Point & clip_pt1 = primary_poly[i];
166 76016 : const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
167 76016 : const Point edg = clip_pt2 - clip_pt1;
168 76016 : const Real cp = clip_pt2(0) * clip_pt1(1) - clip_pt2(1) * clip_pt1(0);
169 :
170 : // Check if point is to the left of (or on) clip_edge
171 : /*
172 : * Note that use of tolerance here is to avoid degenerate case when lines are
173 : * essentially on top of each other (common when meshes match across interface)
174 : * since finding intersection is ill-conditioned in this case.
175 : */
176 559348 : auto is_inside = [&edg, cp](const Point & pt, Real tol)
177 559348 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < tol; };
178 :
179 : // Loop through edges of target polygon (with previous clippings already included)
180 355690 : for (auto j : index_range(input_poly))
181 : {
182 : // Get target edge
183 279674 : const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
184 279674 : const Point prev_pt = input_poly[j];
185 :
186 : // TODO: Don't need to calculate both each loop
187 279674 : const bool is_current_inside = is_inside(curr_pt, _area_tol);
188 279674 : const bool is_previous_inside = is_inside(prev_pt, _area_tol);
189 :
190 279674 : if (is_current_inside)
191 : {
192 197966 : if (!is_previous_inside)
193 : {
194 : Real s;
195 39139 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
196 :
197 : /*
198 : * s is the fraction of distance along clip poly edge that intersection lies
199 : * It is used here to avoid degenerate polygon cases. For example, consider a
200 : * case like:
201 : * o
202 : * | (inside)
203 : * ------|------
204 : * | (outside)
205 : * when the distance is small (< 1e-7) we don't want to to add both the point
206 : * and intersection. Also note that when distance on the scale of 1e-7,
207 : * area on scale of 1e-14 so is insignificant if this results in dropping
208 : * a tri (for example if next edge crosses again)
209 : */
210 39139 : if (s < (1 - _tolerance))
211 38434 : clipped_poly.push_back(intersect);
212 : }
213 197966 : clipped_poly.push_back(curr_pt);
214 : }
215 81708 : else if (is_previous_inside)
216 : {
217 : Real s;
218 39139 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
219 39139 : if (s > _tolerance)
220 38299 : clipped_poly.push_back(intersect);
221 : }
222 : }
223 76016 : }
224 :
225 : // return clipped_poly;
226 : // Make sure final clipped poly is not trivial
227 20866 : if (clipped_poly.size() < 3)
228 : {
229 781 : clipped_poly.clear();
230 781 : return clipped_poly;
231 : }
232 :
233 : // Clean up result by removing any duplicate nodes
234 20085 : std::vector<Point> cleaned_poly;
235 :
236 20085 : cleaned_poly.push_back(clipped_poly.back());
237 75280 : for (auto i : make_range(clipped_poly.size() - 1))
238 : {
239 55195 : const Point prev_pt = cleaned_poly.back();
240 55195 : const Point curr_pt = clipped_poly[i];
241 :
242 : // If points are sufficiently distanced, add to output
243 55195 : if ((curr_pt - prev_pt).norm() > _length_tol)
244 55195 : cleaned_poly.push_back(curr_pt);
245 : }
246 :
247 : mooseAssert(
248 : cleaned_poly.size() <= 8,
249 : "Our distributed mesh numbering scheme assumes that we have at most 8 nodes resulting from "
250 : "clipping the projection of the primary sub-element onto the secondary sub-element");
251 20085 : return cleaned_poly;
252 165487 : }
253 :
254 : void
255 20085 : MortarSegmentHelper::triangulatePoly(std::vector<Point> & poly_nodes,
256 : const unsigned int offset,
257 : std::vector<std::vector<unsigned int>> & tri_map) const
258 : {
259 : // Note offset is important because it converts from list of nodes on local poly to list of
260 : // nodes on entire element. This is a sloppy and error-prone way to do this.
261 : // Could be redesigned by just building and adding elements inside helper, but leaving for now
262 : // to keep buildMortarSegmentMesh similar for 2D and 3D
263 :
264 : // Fewer than 3 nodes can't be triangulated
265 20085 : if (poly_nodes.size() < 3)
266 0 : mooseError("Can't triangulate poly with fewer than 3 nodes");
267 : // If 3 nodes, already a triangle, simply pass back map
268 20085 : else if (poly_nodes.size() == 3)
269 : {
270 15520 : tri_map.push_back({0 + offset, 1 + offset, 2 + offset});
271 7760 : return;
272 : }
273 : // Otherwise use simple center point triangulation
274 : // Note: Could use better algorithm to reduce number of triangles
275 : // - Delaunay produces good quality triangulation but might be a bit overkill
276 : // - Monotone polygon triangulation algorithm is O(N) but quality is not guaranteed
277 : else
278 : {
279 12325 : Point poly_center;
280 12325 : const unsigned int n_verts = poly_nodes.size();
281 :
282 : // Get geometric center of polygon
283 64325 : for (auto node : poly_nodes)
284 52000 : poly_center += node;
285 12325 : poly_center /= n_verts;
286 :
287 : // Add triangles formed by outer edge and center point
288 64325 : for (auto i : make_range(n_verts))
289 156000 : tri_map.push_back({i + offset, (i + 1) % n_verts + offset, n_verts + offset});
290 :
291 : // Add center point to list of polygon nodes
292 12325 : poly_nodes.push_back(poly_center);
293 12325 : return;
294 : }
295 : }
296 :
297 : void
298 165487 : MortarSegmentHelper::getMortarSegments(const std::vector<Point> & primary_nodes,
299 : std::vector<Point> & nodes,
300 : std::vector<std::vector<unsigned int>> & elem_to_nodes)
301 : {
302 : // Clip primary elem against secondary elem
303 165487 : std::vector<Point> clipped_poly = clipPoly(primary_nodes);
304 165487 : if (clipped_poly.size() < 3)
305 145402 : return;
306 :
307 20085 : if (_debug)
308 0 : for (auto pt : clipped_poly)
309 0 : if (!isInsideSecondary(pt))
310 0 : mooseError("Clipped polygon not inside linearized secondary element");
311 :
312 : // Compute area of clipped polygon, update remaining area fraction
313 20085 : _remaining_area_fraction -= area(clipped_poly) / _secondary_area;
314 :
315 : // Triangulate clip polygon
316 20085 : triangulatePoly(clipped_poly, nodes.size(), elem_to_nodes);
317 :
318 : // Transform clipped poly back to (linearized) 3d and append to list
319 107690 : for (auto pt : clipped_poly)
320 87605 : nodes.emplace_back((pt(0) * _u) + (pt(1) * _v) + _center);
321 165487 : }
322 :
323 : Real
324 25675 : MortarSegmentHelper::area(const std::vector<Point> & nodes) const
325 : {
326 25675 : Real poly_area = 0;
327 120083 : for (auto i : index_range(nodes))
328 94408 : poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
329 94408 : nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
330 25675 : poly_area *= 0.5;
331 25675 : return poly_area;
332 : }
|
