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 5446 : MortarSegmentHelper::MortarSegmentHelper(const std::vector<Point> secondary_nodes,
17 : const Point & center,
18 5446 : const Point & normal)
19 5446 : : _center(center), _normal(normal), _debug(false)
20 : {
21 5446 : _secondary_poly.clear();
22 5446 : _secondary_poly.reserve(secondary_nodes.size());
23 :
24 : // Get orientation of secondary poly
25 5446 : const Point e1 = secondary_nodes[0] - secondary_nodes[1];
26 5446 : const Point e2 = secondary_nodes[2] - secondary_nodes[1];
27 5446 : 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 5446 : _u = _normal.cross(secondary_nodes[0] - center).unit();
33 5446 : _v = (orient > 0) ? _normal.cross(_u).unit() : _u.cross(_normal).unit();
34 :
35 : // Transform problem to 2D plane spanned by u and v
36 24142 : for (const auto & node : secondary_nodes)
37 : {
38 18696 : Point pt = node - _center;
39 18696 : _secondary_poly.emplace_back(pt * _u, pt * _v, 0);
40 : }
41 :
42 : // Initialize area of secondary polygon
43 5446 : _remaining_area_fraction = 1.0;
44 5446 : _secondary_area = area(_secondary_poly);
45 :
46 : // Tolerance for quantities with area dimensions
47 5446 : _area_tol = _tolerance * _secondary_area;
48 :
49 : // Tolerance for quantites with length dimensions
50 5446 : _length_tol = _tolerance * std::sqrt(_secondary_area);
51 5446 : }
52 :
53 : Point
54 76118 : MortarSegmentHelper::getIntersection(
55 : const Point & p1, const Point & p2, const Point & q1, const Point & q2, Real & s) const
56 : {
57 76118 : const Point dp = p2 - p1;
58 76118 : const Point dq = q2 - q1;
59 76118 : const Real cp1q1 = p1(0) * q1(1) - p1(1) * q1(0);
60 76118 : const Real cp1q2 = p1(0) * q2(1) - p1(1) * q2(0);
61 76118 : const Real cq1q2 = q1(0) * q2(1) - q1(1) * q2(0);
62 76118 : const Real alpha = 1. / (dp(0) * dq(1) - dp(1) * dq(0));
63 76118 : 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 76118 : s = s > 1 ? 1. : s;
68 76118 : s = s < 0 ? 0. : s;
69 76118 : 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 164150 : MortarSegmentHelper::isDisjoint(const std::vector<Point> & poly) const
96 : {
97 379307 : for (auto i : index_range(_secondary_poly))
98 : {
99 : // Get edge to check
100 356948 : const Point & q1 = _secondary_poly[i];
101 356948 : const Point & q2 = _secondary_poly[(i + 1) % _secondary_poly.size()];
102 356948 : const Point edg = q2 - q1;
103 356948 : 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 4782608 : auto is_inside = [&edg, cp](Point & pt, Real tol)
108 1195652 : { return pt(0) * edg(1) - pt(1) * edg(0) + cp < -tol; };
109 :
110 356948 : bool all_outside = true;
111 1552600 : for (auto pt : poly)
112 1195652 : if (is_inside(pt, _area_tol))
113 587260 : all_outside = false;
114 :
115 356948 : if (all_outside)
116 141791 : return true;
117 : }
118 22359 : return false;
119 : }
120 :
121 : std::vector<Point>
122 164150 : MortarSegmentHelper::clipPoly(const std::vector<Point> & primary_nodes) const
123 : {
124 : // Check orientation of primary_poly
125 164150 : const Point e1 = primary_nodes[0] - primary_nodes[1];
126 164150 : 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 164150 : const Real orient = e2.cross(e1) * _u.cross(_v);
131 :
132 : // Get primary_poly (primary is clipping poly). If negatively oriented, reverse
133 164150 : std::vector<Point> primary_poly;
134 164150 : const int n_verts = primary_nodes.size();
135 707367 : for (auto n : index_range(primary_nodes))
136 : {
137 543217 : Point pt = (orient > 0) ? primary_nodes[n] - _center : primary_nodes[n_verts - 1 - n] - _center;
138 543217 : primary_poly.emplace_back(pt * _u, pt * _v, 0.);
139 : }
140 :
141 164150 : if (isDisjoint(primary_poly))
142 : {
143 141791 : primary_poly.clear();
144 141791 : return primary_poly;
145 : }
146 :
147 : // Initialize clipped poly with secondary poly (secondary is target poly)
148 22359 : std::vector<Point> clipped_poly = _secondary_poly;
149 :
150 : // Loop through clipping edges
151 96431 : for (auto i : index_range(primary_poly))
152 : {
153 : // If clipped poly trivial, return
154 76213 : 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 74072 : std::vector<Point> input_poly(clipped_poly);
162 74072 : clipped_poly.clear();
163 :
164 : // Get clipping edge
165 74072 : const Point & clip_pt1 = primary_poly[i];
166 74072 : const Point & clip_pt2 = primary_poly[(i + 1) % primary_poly.size()];
167 74072 : const Point edg = clip_pt2 - clip_pt1;
168 74072 : 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 2186704 : auto is_inside = [&edg, cp](const Point & pt, Real tol)
177 546676 : { 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 347410 : for (auto j : index_range(input_poly))
181 : {
182 : // Get target edge
183 273338 : const Point curr_pt = input_poly[(j + 1) % input_poly.size()];
184 273338 : const Point prev_pt = input_poly[j];
185 :
186 : // TODO: Don't need to calculate both each loop
187 273338 : const bool is_current_inside = is_inside(curr_pt, _area_tol);
188 273338 : const bool is_previous_inside = is_inside(prev_pt, _area_tol);
189 :
190 273338 : if (is_current_inside)
191 : {
192 193214 : if (!is_previous_inside)
193 : {
194 : Real s;
195 38059 : 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 38059 : if (s < (1 - _tolerance))
211 37498 : clipped_poly.push_back(intersect);
212 : }
213 193214 : clipped_poly.push_back(curr_pt);
214 : }
215 80124 : else if (is_previous_inside)
216 : {
217 : Real s;
218 38059 : Point intersect = getIntersection(prev_pt, curr_pt, clip_pt1, clip_pt2, s);
219 38059 : if (s > _tolerance)
220 37507 : clipped_poly.push_back(intersect);
221 : }
222 : }
223 74072 : }
224 :
225 : // return clipped_poly;
226 : // Make sure final clipped poly is not trivial
227 20218 : 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 19437 : std::vector<Point> cleaned_poly;
235 :
236 19437 : cleaned_poly.push_back(clipped_poly.back());
237 73192 : for (auto i : make_range(clipped_poly.size() - 1))
238 : {
239 53755 : const Point prev_pt = cleaned_poly.back();
240 53755 : const Point curr_pt = clipped_poly[i];
241 :
242 : // If points are sufficiently distanced, add to output
243 53755 : if ((curr_pt - prev_pt).norm() > _length_tol)
244 53755 : cleaned_poly.push_back(curr_pt);
245 : }
246 :
247 19437 : return cleaned_poly;
248 164150 : }
249 :
250 : void
251 19437 : MortarSegmentHelper::triangulatePoly(std::vector<Point> & poly_nodes,
252 : const unsigned int offset,
253 : std::vector<std::vector<unsigned int>> & tri_map) const
254 : {
255 : // Note offset is important because it converts from list of nodes on local poly to list of
256 : // nodes on entire element. This is a sloppy and error-prone way to do this.
257 : // Could be redesigned by just building and adding elements inside helper, but leaving for now
258 : // to keep buildMortarSegmentMesh similar for 2D and 3D
259 :
260 : // Fewer than 3 nodes can't be triangulated
261 19437 : if (poly_nodes.size() < 3)
262 0 : mooseError("Can't triangulate poly with fewer than 3 nodes");
263 : // If 3 nodes, already a triangle, simply pass back map
264 19437 : else if (poly_nodes.size() == 3)
265 : {
266 7256 : tri_map.push_back({0 + offset, 1 + offset, 2 + offset});
267 7256 : return;
268 : }
269 : // Otherwise use simple center point triangulation
270 : // Note: Could use better algorithm to reduce number of triangles
271 : // - Delaunay produces good quality triangulation but might be a bit overkill
272 : // - Monotone polygon triangulation algorithm is O(N) but quality is not guaranteed
273 : else
274 : {
275 12181 : Point poly_center;
276 12181 : const unsigned int n_verts = poly_nodes.size();
277 :
278 : // Get geometric center of polygon
279 63605 : for (auto node : poly_nodes)
280 51424 : poly_center += node;
281 12181 : poly_center /= n_verts;
282 :
283 : // Add triangles formed by outer edge and center point
284 63605 : for (auto i : make_range(n_verts))
285 51424 : tri_map.push_back({i + offset, (i + 1) % n_verts + offset, n_verts + offset});
286 :
287 : // Add center point to list of polygon nodes
288 12181 : poly_nodes.push_back(poly_center);
289 12181 : return;
290 : }
291 : }
292 :
293 : void
294 164150 : MortarSegmentHelper::getMortarSegments(const std::vector<Point> & primary_nodes,
295 : std::vector<Point> & nodes,
296 : std::vector<std::vector<unsigned int>> & elem_to_nodes)
297 : {
298 : // Clip primary elem against secondary elem
299 164150 : std::vector<Point> clipped_poly = clipPoly(primary_nodes);
300 164150 : if (clipped_poly.size() < 3)
301 144713 : return;
302 :
303 19437 : if (_debug)
304 0 : for (auto pt : clipped_poly)
305 0 : if (!isInsideSecondary(pt))
306 0 : mooseError("Clipped polygon not inside linearized secondary element");
307 :
308 : // Compute area of clipped polygon, update remaining area fraction
309 19437 : _remaining_area_fraction -= area(clipped_poly) / _secondary_area;
310 :
311 : // Triangulate clip polygon
312 19437 : triangulatePoly(clipped_poly, nodes.size(), elem_to_nodes);
313 :
314 : // Transform clipped poly back to (linearized) 3d and append to list
315 104810 : for (auto pt : clipped_poly)
316 85373 : nodes.emplace_back((pt(0) * _u) + (pt(1) * _v) + _center);
317 164150 : }
318 :
319 : Real
320 24883 : MortarSegmentHelper::area(const std::vector<Point> & nodes) const
321 : {
322 24883 : Real poly_area = 0;
323 116771 : for (auto i : index_range(nodes))
324 91888 : poly_area += nodes[i](0) * nodes[(i + 1) % nodes.size()](1) -
325 91888 : nodes[i](1) * nodes[(i + 1) % nodes.size()](0);
326 24883 : poly_area *= 0.5;
327 24883 : return poly_area;
328 : }
|
