Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 31210 : C0Polygon::C0Polygon (const unsigned int num_sides, Elem * p) :
33 31210 : Polygon(num_sides, num_sides, p)
34 : {
35 : // A default triangulation is better than nothing
36 100177 : for (int i : make_range(num_sides-2))
37 68967 : this->_triangulation.push_back({0, i+1, i+2});
38 31210 : }
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 8813 : bool C0Polygon::is_vertex(const unsigned int libmesh_dbg_var(i)) const
66 : {
67 836 : libmesh_assert (i < this->n_nodes());
68 :
69 8813 : return true;
70 : }
71 :
72 3053 : bool C0Polygon::is_edge(const unsigned int libmesh_dbg_var(i)) const
73 : {
74 86 : libmesh_assert (i < this->n_nodes());
75 :
76 3053 : 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 12778 : C0Polygon::nodes_on_side(const unsigned int s) const
99 : {
100 916 : const auto ns = this->n_sides();
101 916 : libmesh_assert(!(this->n_nodes() % ns));
102 916 : libmesh_assert_less(s, ns);
103 :
104 12778 : std::vector<unsigned int> returnval(2);
105 12778 : returnval[0] = s;
106 12778 : returnval[1] = (s+1)%ns;
107 :
108 13694 : 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 356750 : Order C0Polygon::default_order() const
163 : {
164 356750 : return FIRST;
165 : }
166 :
167 :
168 :
169 4195 : std::unique_ptr<Elem> C0Polygon::build_side_ptr (const unsigned int i)
170 : {
171 131 : const auto ns = this->n_sides();
172 131 : libmesh_assert_less (i, ns);
173 :
174 4195 : std::unique_ptr<Elem> sidep = std::make_unique<Edge2>();
175 4326 : sidep->set_node(0, this->node_ptr(i));
176 4326 : sidep->set_node(1, this->node_ptr((i+1)%ns));
177 :
178 4195 : sidep->set_interior_parent(this);
179 4064 : sidep->inherit_data_from(*this);
180 :
181 4326 : return sidep;
182 0 : }
183 :
184 :
185 :
186 95964 : void C0Polygon::build_side_ptr (std::unique_ptr<Elem> & side,
187 : const unsigned int i)
188 : {
189 2751 : const auto ns = this->n_sides();
190 2751 : libmesh_assert_less (i, ns);
191 :
192 95964 : if (!side.get() || side->type() != EDGE2)
193 : {
194 7914 : side = this->build_side_ptr(i);
195 : }
196 : else
197 : {
198 89314 : side->inherit_data_from(*this);
199 :
200 94584 : side->set_node(0, this->node_ptr(i));
201 94584 : side->set_node(1, this->node_ptr((i+1)%ns));
202 : }
203 95964 : }
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 284 : Real C0Polygon::volume () const
217 : {
218 : // This specialization is good for Lagrange mappings only
219 284 : if (this->mapping_type() != LAGRANGE_MAP)
220 0 : return this->Elem::volume();
221 :
222 : // We use a triangulation to calculate here; assert that it's as
223 : // consistent as possible.
224 8 : libmesh_assert_equal_to (this->_triangulation.size(),
225 : this->n_nodes() - 2);
226 :
227 8 : Real double_area = 0;
228 1065 : for (const auto & triangle : this->_triangulation)
229 : {
230 781 : Point v01 = this->point(triangle[1]) -
231 803 : this->point(triangle[0]);
232 781 : Point v02 = this->point(triangle[2]) -
233 44 : this->point(triangle[0]);
234 781 : double_area += std::sqrt(cross_norm_sq(v01, v02));
235 : }
236 :
237 284 : return double_area/2;
238 : }
239 :
240 :
241 :
242 72 : Point C0Polygon::true_centroid () const
243 : {
244 : // This specialization is good for Lagrange mappings only
245 72 : if (this->mapping_type() != LAGRANGE_MAP)
246 0 : return this->Elem::true_centroid();
247 :
248 : // We use a triangulation to calculate here; assert that it's as
249 : // consistent as possible.
250 6 : libmesh_assert_equal_to (this->_triangulation.size(),
251 : this->n_nodes() - 2);
252 :
253 6 : Real double_area = 0;
254 6 : Point double_area_weighted_centroid;
255 288 : for (const auto & triangle : this->_triangulation)
256 : {
257 216 : Point v01 = this->point(triangle[1]) -
258 234 : this->point(triangle[0]);
259 216 : Point v02 = this->point(triangle[2]) -
260 36 : this->point(triangle[0]);
261 :
262 216 : const Real double_tri_area = std::sqrt(cross_norm_sq(v01, v02));
263 216 : const Point tri_centroid = (this->point(triangle[0]) +
264 234 : this->point(triangle[1]) +
265 234 : this->point(triangle[2]))/3;
266 :
267 216 : double_area += double_tri_area;
268 :
269 18 : double_area_weighted_centroid += double_tri_area * tri_centroid;
270 : }
271 :
272 6 : return double_area_weighted_centroid / double_area;
273 : }
274 :
275 :
276 :
277 : std::pair<unsigned short int, unsigned short int>
278 0 : C0Polygon::second_order_child_vertex (const unsigned int /*n*/) const
279 : {
280 0 : libmesh_not_implemented();
281 : return std::pair<unsigned short int, unsigned short int> (0, 0);
282 : }
283 :
284 :
285 600 : void C0Polygon::permute(unsigned int perm_num)
286 : {
287 77 : const auto ns = this->n_sides();
288 77 : libmesh_assert_less (perm_num, ns);
289 :
290 : // This is mostly for unit testing so I'll just make it O(N).
291 2880 : for (unsigned int p : make_range(perm_num))
292 : {
293 298 : libmesh_ignore(p);
294 :
295 596 : Node * tempnode = this->node_ptr(0);
296 596 : Elem * tempneigh = this->neighbor_ptr(0);
297 11400 : for (unsigned int s : make_range(ns-1))
298 : {
299 10312 : this->set_node(s, this->node_ptr((s+1)%ns));
300 2384 : this->set_neighbor(s, this->neighbor_ptr(s+1));
301 : }
302 2280 : this->set_node(ns-1, tempnode);
303 596 : this->set_neighbor(ns-1, tempneigh);
304 :
305 : // Keep the same triangulation, just with permuted node order,
306 : // so we can really expect this to act like the *same* polygon
307 9120 : for (auto & triangle : this->_triangulation)
308 27360 : for (auto i : make_range(3))
309 20520 : triangle[i] = (triangle[i]+1)%ns;
310 : }
311 600 : }
312 :
313 :
314 580 : void C0Polygon::flip(BoundaryInfo * boundary_info)
315 : {
316 17 : libmesh_assert(boundary_info);
317 :
318 17 : const auto ns = this->n_sides();
319 :
320 : // Reorder nodes in such a way as to keep side ns-1 the same
321 1882 : for (auto s : make_range(0u, ns/2))
322 : {
323 1302 : swap2nodes(s, ns-1-s);
324 1264 : swap2neighbors(s, ns-2-s);
325 1302 : swap2boundarysides(s, ns-2-s, boundary_info);
326 1302 : swap2boundaryedges(s, ns-2-s, boundary_info);
327 : }
328 580 : }
329 :
330 :
331 :
332 180 : ElemType C0Polygon::side_type (const unsigned int libmesh_dbg_var(s)) const
333 : {
334 15 : libmesh_assert_less (s, this->n_sides());
335 180 : return EDGE2;
336 : }
337 :
338 :
339 : Point
340 710 : C0Polygon::side_vertex_average_normal(const unsigned int s) const
341 : {
342 20 : const auto n_sides = this->n_sides();
343 20 : libmesh_assert_less (s, n_sides);
344 20 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
345 750 : const Point side_t = this->point((s+1) % n_sides) -
346 40 : this->point(s);
347 20 : Point plane_normal(0, 0, 1);
348 : // Find out what plane we're in, if we're in 3D
349 : // Note we don't support non-planar C0-polygon at this time
350 : #if LIBMESH_DIM > 2
351 710 : plane_normal(2) = 0;
352 710 : const auto vavg = this->vertex_average();
353 710 : for (auto i : make_range(n_sides))
354 : {
355 40 : const Point vi = this->point(i) - vavg;
356 710 : const Point viplus = this->point((i+1)%n_sides) - vavg;
357 690 : plane_normal += vi.cross(viplus);
358 : // Since we know the polygon is planar
359 710 : if (plane_normal.norm_sq() > TOLERANCE)
360 20 : break;
361 : }
362 710 : plane_normal = plane_normal.unit();
363 : #endif
364 690 : const Point v = side_t.cross(plane_normal);
365 710 : return v.unit();
366 : }
367 :
368 :
369 0 : void C0Polygon::retriangulate()
370 : {
371 0 : this->_triangulation.clear();
372 :
373 : // Start with the full polygon
374 0 : std::vector<int> remaining_nodes(this->n_nodes());
375 0 : std::iota(remaining_nodes.begin(), remaining_nodes.end(), 0);
376 :
377 0 : const auto ns = this->n_sides();
378 :
379 0 : const Point vavg = this->vertex_average();
380 :
381 : // Find out what plane we're in, if we're in 3D
382 : #if LIBMESH_DIM > 2
383 0 : Point plane_normal;
384 0 : for (auto i : make_range(ns))
385 : {
386 0 : const Point vi = this->point(i) - vavg;
387 0 : const Point viplus = this->point((i+1)%ns) - vavg;
388 0 : plane_normal += vi.cross(viplus);
389 : }
390 0 : plane_normal = plane_normal.unit();
391 : #endif
392 :
393 : // Greedy heuristic: find the node with the most acutely convex
394 : // angle and strip it out as a triangle.
395 0 : while (remaining_nodes.size() > 2)
396 : {
397 0 : Real min_cos_angle = 1;
398 0 : int best_vertex = -1;
399 0 : for (auto n : index_range(remaining_nodes))
400 : {
401 0 : const Point & pn = this->point(remaining_nodes[n]);
402 0 : const Point & pnext = this->point(remaining_nodes[(n+1)%ns]);
403 0 : const Point & pprev = this->point(remaining_nodes[(n+ns-1)%ns]);
404 0 : const Point vprev = (pn - pprev).unit();
405 0 : const Point vnext = (pnext - pn).unit();
406 :
407 0 : const Real sign_check = (vprev.cross(vnext)) * plane_normal;
408 0 : if (sign_check <= 0)
409 0 : continue;
410 :
411 0 : const Real cos_angle = vprev * vnext;
412 0 : if (cos_angle < min_cos_angle)
413 : {
414 0 : min_cos_angle = cos_angle;
415 0 : best_vertex = n;
416 : }
417 : }
418 :
419 0 : libmesh_assert(best_vertex >= 0);
420 :
421 0 : this->_triangulation.push_back({remaining_nodes[(best_vertex+ns-1)%ns],
422 0 : remaining_nodes[best_vertex],
423 0 : remaining_nodes[(best_vertex+1)%ns]});
424 0 : remaining_nodes.erase(remaining_nodes.begin()+best_vertex);
425 : }
426 0 : }
427 :
428 :
429 :
430 : } // namespace libMesh
|
