Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // Local includes
19 : #include "libmesh/face_c0polygon.h"
20 :
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/enum_order.h"
23 : #include "libmesh/tensor_value.h"
24 :
25 : // C++ headers
26 : #include <numeric> // std::iota
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 72461 : C0Polygon::C0Polygon (const unsigned int num_sides, Elem * p) :
33 72461 : Polygon(num_sides, num_sides, p)
34 : {
35 : // A default triangulation is better than nothing
36 250129 : for (int i : make_range(num_sides-2))
37 177668 : this->_triangulation.push_back({0, i+1, i+2});
38 72461 : }
39 :
40 :
41 0 : unsigned int C0Polygon::opposite_node(const unsigned int node_in,
42 : const unsigned int side_in) const
43 : {
44 0 : const auto ns = this->n_sides();
45 0 : if (ns % 2)
46 0 : libmesh_error();
47 :
48 0 : libmesh_assert_less (node_in, ns);
49 0 : libmesh_assert_less (node_in, this->n_nodes());
50 0 : libmesh_assert_less (side_in, this->n_sides());
51 0 : libmesh_assert(this->is_node_on_side(node_in, side_in));
52 :
53 0 : if (node_in == side_in)
54 0 : return (node_in + ns/2 + 1) % ns;
55 :
56 0 : libmesh_assert(node_in == side_in + 1);
57 0 : return (node_in + ns/2 - 1) % ns;
58 : }
59 :
60 :
61 :
62 : // ------------------------------------------------------------
63 : // C0Polygon class member functions
64 :
65 11369 : bool C0Polygon::is_vertex(const unsigned int libmesh_dbg_var(i)) const
66 : {
67 908 : libmesh_assert (i < this->n_nodes());
68 :
69 11369 : return true;
70 : }
71 :
72 5609 : bool C0Polygon::is_edge(const unsigned int libmesh_dbg_var(i)) const
73 : {
74 158 : libmesh_assert (i < this->n_nodes());
75 :
76 5609 : return false;
77 : }
78 :
79 1349 : bool C0Polygon::is_face(const unsigned int libmesh_dbg_var(i)) const
80 : {
81 38 : libmesh_assert (i < this->n_nodes());
82 :
83 1349 : return false;
84 : }
85 :
86 5756 : bool C0Polygon::is_node_on_side(const unsigned int n,
87 : const unsigned int s) const
88 : {
89 182 : const auto ns = this->n_sides();
90 182 : libmesh_assert_less (s, ns);
91 182 : libmesh_assert_less (n, this->n_nodes());
92 :
93 5847 : return ((n % ns) == s) ||
94 3060 : ((n < ns) && ((s+1)%ns) == n);
95 : }
96 :
97 : std::vector<unsigned int>
98 12742 : C0Polygon::nodes_on_side(const unsigned int s) const
99 : {
100 913 : const auto ns = this->n_sides();
101 913 : libmesh_assert(!(this->n_nodes() % ns));
102 913 : libmesh_assert_less(s, ns);
103 :
104 12742 : std::vector<unsigned int> returnval(2);
105 12742 : returnval[0] = s;
106 12742 : returnval[1] = (s+1)%ns;
107 :
108 13655 : return returnval;
109 : }
110 :
111 : std::vector<unsigned int>
112 1589 : C0Polygon::nodes_on_edge(const unsigned int e) const
113 : {
114 1589 : return this->nodes_on_side(e);
115 : }
116 :
117 22964 : bool C0Polygon::has_affine_map() const
118 : {
119 1898 : const unsigned int ns = this->n_sides();
120 :
121 : // Get a good tolerance to use
122 1898 : Real perimeter_l1 = 0;
123 137713 : for (auto s : make_range(ns))
124 229498 : perimeter_l1 += (this->point((s+1)%ns) - this->point(s)).l1_norm();
125 22964 : const Real tol = perimeter_l1 * affine_tol;
126 :
127 : // Find the map we expect to have
128 3796 : const Point veci = this->point(1) - this->point(0);
129 1898 : const Point vec12 = this->point(2) - this->point(1);
130 :
131 : // Exterior angle
132 22964 : const Real theta = 2 * libMesh::pi / ns;
133 22964 : const Real costheta = cos(theta);
134 22964 : const Real sintheta = sin(theta);
135 :
136 3796 : RealTensor map;
137 22964 : map(0, 0) = veci(0);
138 22964 : map(1, 0) = veci(1);
139 22964 : map(2, 0) = veci(2);
140 22964 : map(0, 1) = (vec12(0) - veci(0)*costheta)/sintheta;
141 22964 : map(1, 1) = (vec12(1) - veci(1)*costheta)/sintheta;
142 22964 : map(2, 1) = (vec12(2) - veci(2)*costheta)/sintheta;
143 :
144 1898 : libmesh_assert_less((map * this->master_point(1) -
145 : (this->point(1) - this->point(0))).l1_norm(),
146 : tol);
147 1898 : libmesh_assert_less((map * this->master_point(2) -
148 : (this->point(2) - this->point(0))).l1_norm(),
149 : tol);
150 :
151 67960 : for (auto i : make_range(3u, ns))
152 87159 : if ((map * this->master_point(i) -
153 11295 : (this->point(i) - this->point(0))).l1_norm() >
154 : tol)
155 31 : return false;
156 :
157 1867 : return true;
158 : }
159 :
160 :
161 :
162 409970 : Order C0Polygon::default_order() const
163 : {
164 409970 : return FIRST;
165 : }
166 :
167 :
168 :
169 6467 : std::unique_ptr<Elem> C0Polygon::build_side_ptr (const unsigned int i)
170 : {
171 195 : const auto ns = this->n_sides();
172 195 : libmesh_assert_less (i, ns);
173 :
174 6467 : std::unique_ptr<Elem> sidep = std::make_unique<Edge2>();
175 6662 : sidep->set_node(0, this->node_ptr(i));
176 6662 : sidep->set_node(1, this->node_ptr((i+1)%ns));
177 :
178 6467 : sidep->set_interior_parent(this);
179 6272 : sidep->inherit_data_from(*this);
180 :
181 6662 : return sidep;
182 0 : }
183 :
184 :
185 :
186 150137 : void C0Polygon::build_side_ptr (std::unique_ptr<Elem> & side,
187 : const unsigned int i)
188 : {
189 4277 : const auto ns = this->n_sides();
190 4277 : libmesh_assert_less (i, ns);
191 :
192 150137 : if (!side.get() || side->type() != EDGE2)
193 : {
194 12394 : side = this->build_side_ptr(i);
195 : }
196 : else
197 : {
198 139753 : side->inherit_data_from(*this);
199 :
200 147947 : side->set_node(0, this->node_ptr(i));
201 147947 : side->set_node(1, this->node_ptr((i+1)%ns));
202 : }
203 150137 : }
204 :
205 :
206 :
207 0 : void C0Polygon::connectivity(const unsigned int /*sf*/,
208 : const IOPackage /*iop*/,
209 : std::vector<dof_id_type> & /*conn*/) const
210 : {
211 0 : libmesh_not_implemented();
212 : }
213 :
214 :
215 :
216 355 : Real C0Polygon::volume () const
217 : {
218 : // This specialization is good for non-Lagrange mappings, but
219 : // remember to be more careful when we extend it to non-planar or
220 : // non p=1 polygon types.
221 : // if (this->mapping_type() != LAGRANGE_MAP)
222 : // return this->Elem::volume();
223 :
224 : // We use a triangulation to calculate here; assert that it's as
225 : // consistent as possible.
226 10 : libmesh_assert_equal_to (this->_triangulation.size(),
227 : this->n_nodes() - 2);
228 :
229 10 : Real double_area = 0;
230 1562 : for (const auto & triangle : this->_triangulation)
231 : {
232 1207 : Point v01 = this->point(triangle[1]) -
233 1241 : this->point(triangle[0]);
234 1207 : Point v02 = this->point(triangle[2]) -
235 68 : this->point(triangle[0]);
236 1207 : double_area += std::sqrt(cross_norm_sq(v01, v02));
237 : }
238 :
239 355 : return double_area/2;
240 : }
241 :
242 :
243 :
244 72 : Point C0Polygon::true_centroid () const
245 : {
246 : // This specialization is good for Lagrange mappings only
247 72 : if (this->mapping_type() != LAGRANGE_MAP)
248 0 : return this->Elem::true_centroid();
249 :
250 : // We use a triangulation to calculate here; assert that it's as
251 : // consistent as possible.
252 6 : libmesh_assert_equal_to (this->_triangulation.size(),
253 : this->n_nodes() - 2);
254 :
255 6 : Real double_area = 0;
256 6 : Point double_area_weighted_centroid;
257 288 : for (const auto & triangle : this->_triangulation)
258 : {
259 216 : Point v01 = this->point(triangle[1]) -
260 234 : this->point(triangle[0]);
261 216 : Point v02 = this->point(triangle[2]) -
262 36 : this->point(triangle[0]);
263 :
264 216 : const Real double_tri_area = std::sqrt(cross_norm_sq(v01, v02));
265 216 : const Point tri_centroid = (this->point(triangle[0]) +
266 234 : this->point(triangle[1]) +
267 234 : this->point(triangle[2]))/3;
268 :
269 216 : double_area += double_tri_area;
270 :
271 18 : double_area_weighted_centroid += double_tri_area * tri_centroid;
272 : }
273 :
274 6 : return double_area_weighted_centroid / double_area;
275 : }
276 :
277 :
278 :
279 : std::pair<unsigned short int, unsigned short int>
280 0 : C0Polygon::second_order_child_vertex (const unsigned int /*n*/) const
281 : {
282 0 : libmesh_not_implemented();
283 : return std::pair<unsigned short int, unsigned short int> (0, 0);
284 : }
285 :
286 :
287 600 : void C0Polygon::permute(unsigned int perm_num)
288 : {
289 77 : const auto ns = this->n_sides();
290 77 : libmesh_assert_less (perm_num, ns);
291 :
292 : // This is mostly for unit testing so I'll just make it O(N).
293 2880 : for (unsigned int p : make_range(perm_num))
294 : {
295 298 : libmesh_ignore(p);
296 :
297 596 : Node * tempnode = this->node_ptr(0);
298 596 : Elem * tempneigh = this->neighbor_ptr(0);
299 11400 : for (unsigned int s : make_range(ns-1))
300 : {
301 10312 : this->set_node(s, this->node_ptr((s+1)%ns));
302 2384 : this->set_neighbor(s, this->neighbor_ptr(s+1));
303 : }
304 2280 : this->set_node(ns-1, tempnode);
305 596 : this->set_neighbor(ns-1, tempneigh);
306 :
307 : // Keep the same triangulation, just with permuted node order,
308 : // so we can really expect this to act like the *same* polygon
309 9120 : for (auto & triangle : this->_triangulation)
310 27360 : for (auto i : make_range(3))
311 20520 : triangle[i] = (triangle[i]+1)%ns;
312 : }
313 600 : }
314 :
315 :
316 580 : void C0Polygon::flip(BoundaryInfo * boundary_info)
317 : {
318 17 : libmesh_assert(boundary_info);
319 :
320 17 : const auto ns = this->n_sides();
321 :
322 : // Reorder nodes in such a way as to keep side ns-1 the same
323 1882 : for (auto s : make_range(0u, ns/2))
324 : {
325 1302 : swap2nodes(s, ns-1-s);
326 1264 : swap2neighbors(s, ns-2-s);
327 1302 : swap2boundarysides(s, ns-2-s, boundary_info);
328 1302 : swap2boundaryedges(s, ns-2-s, boundary_info);
329 : }
330 580 : }
331 :
332 :
333 :
334 180 : ElemType C0Polygon::side_type (const unsigned int libmesh_dbg_var(s)) const
335 : {
336 15 : libmesh_assert_less (s, this->n_sides());
337 180 : return EDGE2;
338 : }
339 :
340 :
341 : Point
342 710 : C0Polygon::side_vertex_average_normal(const unsigned int s) const
343 : {
344 20 : const auto n_sides = this->n_sides();
345 20 : libmesh_assert_less (s, n_sides);
346 20 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
347 750 : const Point side_t = this->point((s+1) % n_sides) -
348 40 : this->point(s);
349 20 : Point plane_normal(0, 0, 1);
350 : // Find out what plane we're in, if we're in 3D
351 : // Note we don't support non-planar C0-polygon at this time
352 : #if LIBMESH_DIM > 2
353 710 : plane_normal(2) = 0;
354 710 : const auto vavg = this->vertex_average();
355 710 : for (auto i : make_range(n_sides))
356 : {
357 40 : const Point vi = this->point(i) - vavg;
358 710 : const Point viplus = this->point((i+1)%n_sides) - vavg;
359 690 : plane_normal += vi.cross(viplus);
360 : // Since we know the polygon is planar
361 710 : if (plane_normal.norm_sq() > TOLERANCE)
362 20 : break;
363 : }
364 710 : plane_normal = plane_normal.unit();
365 : #endif
366 690 : const Point v = side_t.cross(plane_normal);
367 710 : return v.unit();
368 : }
369 :
370 :
371 0 : void C0Polygon::retriangulate()
372 : {
373 0 : this->_triangulation.clear();
374 :
375 : // Start with the full polygon
376 0 : std::vector<int> remaining_nodes(this->n_nodes());
377 0 : std::iota(remaining_nodes.begin(), remaining_nodes.end(), 0);
378 :
379 0 : const auto ns = this->n_sides();
380 :
381 0 : const Point vavg = this->vertex_average();
382 :
383 : // Find out what plane we're in, if we're in 3D
384 : #if LIBMESH_DIM > 2
385 0 : Point plane_normal;
386 0 : for (auto i : make_range(ns))
387 : {
388 0 : const Point vi = this->point(i) - vavg;
389 0 : const Point viplus = this->point((i+1)%ns) - vavg;
390 0 : plane_normal += vi.cross(viplus);
391 : }
392 0 : plane_normal = plane_normal.unit();
393 : #endif
394 :
395 : // Greedy heuristic: find the node with the most acutely convex
396 : // angle and strip it out as a triangle.
397 0 : while (remaining_nodes.size() > 2)
398 : {
399 0 : Real min_cos_angle = 1;
400 0 : int best_vertex = -1;
401 0 : for (auto n : index_range(remaining_nodes))
402 : {
403 0 : const Point & pn = this->point(remaining_nodes[n]);
404 0 : const Point & pnext = this->point(remaining_nodes[(n+1)%ns]);
405 0 : const Point & pprev = this->point(remaining_nodes[(n+ns-1)%ns]);
406 0 : const Point vprev = (pn - pprev).unit();
407 0 : const Point vnext = (pnext - pn).unit();
408 :
409 0 : const Real sign_check = (vprev.cross(vnext)) * plane_normal;
410 0 : if (sign_check <= 0)
411 0 : continue;
412 :
413 0 : const Real cos_angle = vprev * vnext;
414 0 : if (cos_angle < min_cos_angle)
415 : {
416 0 : min_cos_angle = cos_angle;
417 0 : best_vertex = n;
418 : }
419 : }
420 :
421 0 : libmesh_assert(best_vertex >= 0);
422 :
423 0 : this->_triangulation.push_back({remaining_nodes[(best_vertex+ns-1)%ns],
424 0 : remaining_nodes[best_vertex],
425 0 : remaining_nodes[(best_vertex+1)%ns]});
426 0 : remaining_nodes.erase(remaining_nodes.begin()+best_vertex);
427 : }
428 0 : }
429 :
430 :
431 :
432 : } // namespace libMesh
|
