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 : #include "libmesh/face_polygon.h"
19 :
20 : // Local includes
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/enum_elem_quality.h"
23 : #include "libmesh/hashword.h"
24 :
25 : // C++ includes
26 : #include <array>
27 :
28 :
29 : namespace libMesh
30 : {
31 :
32 : // ------------------------------------------------------------
33 : // Polygon class static member initializations
34 :
35 : const int Polygon::num_children;
36 :
37 : // ------------------------------------------------------------
38 : // Polygon class member functions
39 :
40 :
41 : /**
42 : * We heap-allocate our element and node links, but we can't do that
43 : * until *after* we've intitialized our parent class, which means
44 : * we need to do a little initialization of those links manually too.
45 : */
46 29293 : Polygon::Polygon (const unsigned int nn,
47 : const unsigned int ns,
48 29293 : Elem * p) :
49 : Face(nn, ns, p, nullptr, nullptr),
50 1728 : _elemlinks_data(ns+2), // neighbors + parent + interior_parent
51 29293 : _nodelinks_data(nn)
52 : {
53 : // Do the manual initialization that Elem::Elem couldn't. No need
54 : // to manually set nullptr, though, since std::vector does that.
55 29293 : this->_elemlinks = _elemlinks_data.data();
56 29293 : this->_nodes = _nodelinks_data.data();
57 29293 : this->_elemlinks[0] = p;
58 :
59 : // Is this likely to ever be used? We may do refinement with
60 : // polygons but it's probably not going to have a hierarchy...
61 29293 : if (p)
62 : {
63 0 : this->subdomain_id() = p->subdomain_id();
64 0 : this->processor_id() = p->processor_id();
65 0 : _map_type = p->mapping_type();
66 0 : _map_data = p->mapping_data();
67 :
68 : #ifdef LIBMESH_ENABLE_AMR
69 0 : this->set_p_level(p->p_level());
70 : #endif
71 : }
72 :
73 : // Make sure the interior parent isn't undefined
74 : if (LIBMESH_DIM > 2)
75 29293 : this->set_interior_parent(nullptr);
76 29293 : }
77 :
78 :
79 :
80 11966161 : Point Polygon::master_point (const unsigned int i) const
81 : {
82 1032614 : const unsigned int ns = this->n_sides();
83 :
84 : // Non-vertex points need to be handled by a subclass overriding
85 : // this.
86 11966161 : if (i > ns)
87 0 : libmesh_not_implemented();
88 :
89 : // Center-to-midside to center-to-vertex angle
90 11966161 : const Real pi_over_ns = libMesh::pi / ns;
91 :
92 : // Center-to-midside distance
93 11966161 : const Real min_r = 0.5/tan(pi_over_ns);
94 :
95 : // Center-to-vertex distance
96 11966161 : const Real max_r = std::sqrt(min_r*min_r + 0.25);
97 :
98 1032614 : const Point center(0.5, min_r);
99 :
100 1032614 : libmesh_assert_less(i, this->n_nodes());
101 :
102 11966161 : return center + Point(max_r*sin((int(i)*2-1)*pi_over_ns),
103 13977493 : -max_r*cos((int(i)*2-1)*pi_over_ns));
104 :
105 : }
106 :
107 :
108 :
109 25880 : bool Polygon::on_reference_element(const Point & p,
110 : const Real eps) const
111 : {
112 3825 : const unsigned int ns = this->n_sides();
113 :
114 : // Center-to-midside to center-to-vertex angle
115 25880 : const Real pi_over_ns = libMesh::pi / ns;
116 :
117 : // Center-to-midside distance
118 25880 : const Real min_r = 0.5/tan(pi_over_ns);
119 :
120 : // Center-to-vertex distance
121 25880 : const Real max_r = std::sqrt(min_r*min_r + 0.25);
122 :
123 3825 : const Point center(0.5, min_r);
124 :
125 : // Check that the point is on the same side of all the faces by
126 : // testing whether:
127 : //
128 : // n_i.(p - x_i) <= 0
129 : //
130 : // for each i, where:
131 : // n_i is the outward normal of face i,
132 : // x_i is a point on face i.
133 :
134 146640 : for (auto i : make_range(ns))
135 : {
136 : const Point x_i =
137 141505 : center + Point(max_r*sin((int(i)*2-1)*pi_over_ns),
138 122920 : -max_r*cos((int(i)*2-1)*pi_over_ns));
139 : const Point n_i =
140 18585 : Point(sin(int(i)*2*pi_over_ns),
141 122920 : -cos(int(i)*2*pi_over_ns));
142 :
143 122920 : if (n_i * (p - x_i) > eps)
144 180 : return false;
145 : }
146 :
147 3645 : return true;
148 : }
149 :
150 :
151 :
152 0 : dof_id_type Polygon::key (const unsigned int s) const
153 : {
154 0 : const int ns = this->n_sides();
155 0 : libmesh_assert_less (s, ns);
156 :
157 0 : return this->compute_key(this->node_id(s),
158 0 : this->node_id((s+1)%ns));
159 : }
160 :
161 :
162 :
163 18815 : dof_id_type Polygon::low_order_key (const unsigned int s) const
164 : {
165 530 : const int ns = this->n_sides();
166 530 : libmesh_assert_less (s, ns);
167 :
168 18815 : return this->compute_key(this->node_id(s),
169 19875 : this->node_id((s+1)%ns));
170 : }
171 :
172 :
173 :
174 10220 : unsigned int Polygon::local_side_node(unsigned int side,
175 : unsigned int side_node) const
176 : {
177 554 : const auto ns = this->n_sides();
178 554 : libmesh_assert_less (side, ns);
179 554 : libmesh_assert_less (side_node, this->n_nodes_per_side());
180 :
181 10220 : if (!side_node)
182 277 : return side;
183 :
184 5110 : if (side_node == 1)
185 5110 : return (side + 1) % ns;
186 :
187 0 : return (side_node - 1) * ns + side;
188 : }
189 :
190 :
191 :
192 7522 : unsigned int Polygon::local_edge_node(unsigned int edge,
193 : unsigned int edge_node) const
194 : {
195 7522 : return local_side_node(edge, edge_node);
196 : }
197 :
198 :
199 :
200 0 : dof_id_type Polygon::key () const
201 : {
202 0 : std::vector<dof_id_type> node_ids;
203 0 : for (const auto & n : this->node_ref_range())
204 0 : node_ids.push_back(n.id());
205 :
206 0 : return Utility::hashword(node_ids);
207 : }
208 :
209 :
210 :
211 1349 : std::unique_ptr<Elem> Polygon::side_ptr (const unsigned int i)
212 : {
213 38 : const auto ns = this->n_sides();
214 38 : libmesh_assert_less (i, ns);
215 :
216 1349 : std::unique_ptr<Elem> edge = std::make_unique<Edge2>();
217 :
218 1387 : edge->set_node(0, this->node_ptr(i));
219 1387 : edge->set_node(1, this->node_ptr((i+1)%ns));
220 :
221 1387 : return edge;
222 0 : }
223 :
224 :
225 :
226 0 : void Polygon::side_ptr (std::unique_ptr<Elem> & side,
227 : const unsigned int i)
228 : {
229 0 : const auto ns = this->n_sides();
230 0 : libmesh_assert_less (i, ns);
231 :
232 0 : if (!side.get() || side->type() != EDGE2)
233 0 : side = std::make_unique<Edge2>();
234 :
235 0 : side->subdomain_id() = this->subdomain_id();
236 :
237 0 : side->set_node(0, this->node_ptr(i));
238 0 : side->set_node(1, this->node_ptr((i+1)%ns));
239 0 : }
240 :
241 :
242 :
243 0 : bool Polygon::is_child_on_side(const unsigned int /*c*/,
244 : const unsigned int /*s*/) const
245 : {
246 0 : libmesh_not_implemented();
247 : return false;
248 : }
249 :
250 :
251 :
252 994 : unsigned int Polygon::opposite_side(const unsigned int side_in) const
253 : {
254 28 : const auto ns = this->n_sides();
255 994 : if (ns % 2)
256 0 : libmesh_error();
257 :
258 994 : return (ns / 2 + side_in) % ns;
259 : }
260 :
261 :
262 :
263 631 : bool Polygon::is_flipped() const
264 : {
265 : // Just check based on vertices, the first ns nodes
266 22 : const auto ns = this->n_sides();
267 :
268 : #if LIBMESH_DIM > 2
269 : // Don't bother outside the XY plane
270 3644 : for (auto n : make_range(ns))
271 3119 : if (this->point(n)(2))
272 0 : return false;
273 : #endif
274 :
275 22 : Real twice_signed_area = 0;
276 3644 : for (auto n : make_range(ns))
277 : {
278 212 : const Point & p1 = this->point(n);
279 3013 : const Point & p2 = this->point((n+1)%ns);
280 3013 : twice_signed_area += p1(0)*p2(1)-p1(1)*p2(0);
281 : }
282 :
283 631 : return (twice_signed_area < 0);
284 : }
285 :
286 :
287 : std::vector<unsigned int>
288 300 : Polygon::edges_adjacent_to_node(const unsigned int n) const
289 : {
290 25 : libmesh_assert_less(n, this->n_nodes());
291 :
292 : // For mid-face nodes, the subclass had better have overridden this.
293 25 : libmesh_assert_less(n, this->n_sides() * (this->n_nodes_per_side() + 1));
294 :
295 25 : const unsigned int ns = this->n_sides();
296 :
297 : // For vertices, we have two adjacent edges; otherwise each of the
298 : // mid-edge nodes is adjacent only to the edge it is on.
299 300 : if (this->is_vertex(n))
300 300 : return {n, (n+ns-1)%ns};
301 :
302 0 : return {n % ns};
303 : }
304 :
305 :
306 0 : std::pair<Real, Real> Polygon::qual_bounds (const ElemQuality q) const
307 : {
308 0 : std::pair<Real, Real> bounds;
309 :
310 0 : switch (q)
311 : {
312 0 : case EDGE_LENGTH_RATIO:
313 0 : bounds.first = 1.;
314 0 : bounds.second = 4.;
315 0 : break;
316 :
317 0 : case MIN_ANGLE:
318 0 : bounds.first = 90. - (180./this->n_sides());
319 0 : bounds.second = 180. - (360./this->n_sides());
320 0 : break;
321 :
322 0 : case MAX_ANGLE:
323 0 : bounds.first = 180. - (360./this->n_sides());
324 0 : bounds.second = 180. - (180./this->n_sides());
325 0 : break;
326 :
327 0 : case JACOBIAN:
328 : case SCALED_JACOBIAN:
329 0 : bounds.first = 0.5;
330 0 : bounds.second = 1.;
331 0 : break;
332 :
333 0 : default:
334 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
335 0 : bounds.first = -1;
336 0 : bounds.second = -1;
337 : }
338 :
339 0 : return bounds;
340 : }
341 :
342 :
343 3972201 : std::array<Point, 3> Polygon::master_subtriangle (unsigned int i) const
344 : {
345 341676 : libmesh_assert_less(i, this->_triangulation.size());
346 :
347 3972201 : const auto & tri = this->_triangulation[i];
348 :
349 3972201 : return { this->master_point(tri[0]),
350 3972201 : this->master_point(tri[1]),
351 3972201 : this->master_point(tri[2]) };
352 : }
353 :
354 :
355 : std::tuple<unsigned int, Real, Real>
356 2003000 : Polygon::subtriangle_coordinates (const Point & p, Real tol) const
357 : {
358 173535 : std::tuple<unsigned int, Real, Real> returnval = {libMesh::invalid_uint, -1, -1};
359 :
360 173535 : Real best_bad_coord = -1;
361 :
362 3241892 : for (auto s : make_range(this->n_subtriangles()))
363 : {
364 : const std::array<Point, 3> subtri =
365 3190592 : this->master_subtriangle (s);
366 :
367 : // Find barycentric coordinates in subtri.
368 : // Hand coding these since we don't need a full 3D cross product
369 : // in 2D
370 3190592 : const Real twice_area = (- subtri[1](1) * subtri[2](0)
371 3190592 : - subtri[0](1) * subtri[1](0)
372 3190592 : + subtri[0](1) * subtri[2](0)
373 3190592 : + subtri[0](0) * subtri[1](1)
374 3190592 : - subtri[0](0) * subtri[2](1)
375 3190592 : + subtri[1](0) * subtri[2](1));
376 :
377 277619 : const Real xi = ( subtri[0](1)*subtri[2](0)
378 3190592 : - subtri[0](0)*subtri[2](1)
379 3190592 : + (subtri[2](1)-subtri[0](1))*p(0)
380 3190592 : + (subtri[0](0)-subtri[2](0))*p(1)) / twice_area;
381 :
382 277619 : const Real eta = ( subtri[0](0)*subtri[1](1)
383 3190592 : - subtri[0](1)*subtri[1](0)
384 3190592 : + (subtri[0](1)-subtri[1](1))*p(0)
385 3190592 : + (subtri[1](0)-subtri[0](0))*p(1)) / twice_area;
386 :
387 3190592 : if (xi>=0 && eta>=0 && xi+eta<=1)
388 1951700 : return { s, xi, eta };
389 :
390 1238892 : Real my_best_bad_coord = std::min(std::min(xi, eta), 1-xi-eta);
391 :
392 1238892 : if (my_best_bad_coord > best_bad_coord)
393 : {
394 83001 : best_bad_coord = my_best_bad_coord;
395 83001 : returnval = { s, xi, eta };
396 : }
397 : }
398 :
399 51300 : if (best_bad_coord > -tol)
400 4275 : return returnval;
401 :
402 0 : return {libMesh::invalid_uint, -1, -1};
403 : }
404 :
405 :
406 : } // namespace libMesh
|
