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_tri.h"
20 : #include "libmesh/edge_edge2.h"
21 : #include "libmesh/face_tri3.h"
22 : #include "libmesh/enum_elem_quality.h"
23 :
24 : // C++ includes
25 : #include <array>
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 : // ------------------------------------------------------------
32 : // Tri class static member initializations
33 : const int Tri::num_sides;
34 : const int Tri::num_children;
35 :
36 : // Note: we can omit initialization of the third entry of each row because
37 : // static variables are automatically zero-initialized.
38 : const Real Tri::_master_points[7][3] =
39 : {
40 : {0, 0},
41 : {1, 0},
42 : {0, 1},
43 : {0.5, 0},
44 : {0.5, 0.5},
45 : {0, 0.5},
46 : {Real(1)/3, Real(1)/3}
47 : };
48 :
49 : const unsigned int Tri::adjacent_sides_map[/*num_vertices*/3][/*n_adjacent_sides*/2] =
50 : {
51 : {0, 2}, // Sides adjacent to node 0
52 : {0, 1}, // Sides adjacent to node 1
53 : {1, 2} // Sides adjacent to node 2
54 : };
55 :
56 :
57 :
58 : // ------------------------------------------------------------
59 : // Tri class member functions
60 0 : dof_id_type Tri::key (const unsigned int s) const
61 : {
62 0 : libmesh_assert_less (s, this->n_sides());
63 :
64 0 : return this->compute_key(this->node_id(Tri3::side_nodes_map[s][0]),
65 0 : this->node_id(Tri3::side_nodes_map[s][1]));
66 : }
67 :
68 :
69 :
70 33552771 : dof_id_type Tri::low_order_key (const unsigned int s) const
71 : {
72 846084 : libmesh_assert_less (s, this->n_sides());
73 :
74 34398477 : return this->compute_key(this->node_id(Tri3::side_nodes_map[s][0]),
75 34398855 : this->node_id(Tri3::side_nodes_map[s][1]));
76 : }
77 :
78 :
79 :
80 4689204 : unsigned int Tri::local_side_node(unsigned int side,
81 : unsigned int side_node) const
82 : {
83 390844 : libmesh_assert_less (side, this->n_sides());
84 390844 : libmesh_assert_less (side_node, Tri3::nodes_per_side);
85 :
86 4689204 : return Tri3::side_nodes_map[side][side_node];
87 : }
88 :
89 :
90 :
91 111038504 : unsigned int Tri::local_edge_node(unsigned int edge,
92 : unsigned int edge_node) const
93 : {
94 111038504 : return local_side_node(edge, edge_node);
95 : }
96 :
97 :
98 :
99 48 : dof_id_type Tri::key () const
100 : {
101 52 : return this->compute_key(this->node_id(0),
102 : this->node_id(1),
103 48 : this->node_id(2));
104 : }
105 :
106 :
107 :
108 213395 : std::unique_ptr<Elem> Tri::side_ptr (const unsigned int i)
109 : {
110 6300 : libmesh_assert_less (i, this->n_sides());
111 :
112 213395 : std::unique_ptr<Elem> edge = std::make_unique<Edge2>();
113 :
114 640185 : for (auto n : edge->node_index_range())
115 439326 : edge->set_node(n, this->node_ptr(Tri3::side_nodes_map[i][n]));
116 :
117 213395 : return edge;
118 0 : }
119 :
120 :
121 :
122 31617703 : void Tri::side_ptr (std::unique_ptr<Elem> & side,
123 : const unsigned int i)
124 : {
125 31617703 : this->simple_side_ptr<Tri,Tri3>(side, i, EDGE2);
126 31617703 : }
127 :
128 :
129 :
130 11175298 : bool Tri::is_child_on_side(const unsigned int c,
131 : const unsigned int s) const
132 : {
133 3230222 : libmesh_assert_less (c, this->n_children());
134 3230222 : libmesh_assert_less (s, this->n_sides());
135 :
136 11175298 : return (c == s || c == (s+1)%3);
137 : }
138 :
139 :
140 46054 : bool Tri::is_flipped() const
141 : {
142 : return (
143 : #if LIBMESH_DIM > 2
144 : // Don't bother outside the XY plane
145 87854 : !this->point(0)(2) && !this->point(1)(2) &&
146 133908 : !this->point(2)(2) &&
147 : #endif
148 46310 : ((this->point(1)(0)-this->point(0)(0))*
149 46054 : (this->point(2)(1)-this->point(0)(1)) <
150 46310 : (this->point(2)(0)-this->point(0)(0))*
151 87854 : (this->point(1)(1)-this->point(0)(1))));
152 : }
153 :
154 :
155 : std::vector<unsigned int>
156 9120 : Tri::edges_adjacent_to_node(const unsigned int n) const
157 : {
158 760 : libmesh_assert_less(n, this->n_nodes());
159 :
160 : // For vertices, we use the Tri::adjacent_sides_map, otherwise each
161 : // of the mid-edge nodes is adjacent only to the edge it is on, and the
162 : // center node is not adjacent to any edge.
163 9120 : if (this->is_vertex(n))
164 6240 : return {std::begin(adjacent_sides_map[n]), std::end(adjacent_sides_map[n])};
165 3360 : else if (this->is_edge(n))
166 2880 : return {n - this->n_vertices()};
167 :
168 40 : libmesh_assert(this->is_face(n));
169 440 : return {};
170 : }
171 :
172 :
173 2133 : Real Tri::quality (const ElemQuality q) const
174 : {
175 2133 : switch (q)
176 : {
177 213 : case ASPECT_RATIO:
178 : {
179 : // Aspect Ratio definition from Ansys Theory Manual.
180 : // Reference: Ansys, Inc. Theory Reference, Ansys Release 9.0, 2004 (Chapter: 13.7.3)
181 :
182 : // Compute midpoint positions along each edge
183 : Point m[3] = {
184 12 : Real(0.5) * (this->point(0) + this->point(1)), // side opposite vertex 2
185 6 : Real(0.5) * (this->point(1) + this->point(2)), // side opposite vertex 0
186 18 : Real(0.5) * (this->point(2) + this->point(0))}; // side opposite vertex 1
187 :
188 : // opposite[i] is the side index which is "opposite" vertex i
189 : static const unsigned int opposite[3] = {1, 2, 0};
190 :
191 : // other[i] is the side index which is _not_ i and _not_ opposite[i]
192 : static const unsigned int other[3] = {2, 0, 1};
193 :
194 : // Input is vertex index, i = 0, 1, 2
195 639 : auto vertex_aspect_ratio = [&](unsigned int i) -> Real
196 : {
197 : // Compute vectors:
198 : // v0: (vertex v, opposite midpoint)
199 : // v1: (midpoint[i], last midpoint)
200 657 : Point v0 = m[opposite[i]] - this->point(i);
201 639 : Point v1 = m[other[i]] - m[i];
202 :
203 : // Compute the length of the midlines
204 639 : Real v0_norm = v0.norm();
205 639 : Real v1_norm = v1.norm();
206 :
207 : // Instead of dividing by zero in the next step, just return
208 : // 0. The optimal aspect ratio is 1.0, and "high" aspect
209 : // ratios are bad, but an aspect ratio of 0 should also be
210 : // considered bad.
211 639 : if (v0_norm == 0. || v1_norm == 0.)
212 0 : return 0.;
213 :
214 : // Compute sine of the angle between v0, v1.
215 639 : Real sin_theta = cross_norm(v0, v1) / v0_norm / v1_norm;
216 639 : Real v0s = v0_norm*sin_theta;
217 639 : Real v1s = v1_norm*sin_theta;
218 :
219 : // Determine the min, max of each midline length and its
220 : // projection.
221 18 : auto [min0, max0] = std::minmax(v0_norm, v1s);
222 18 : auto [min1, max1] = std::minmax(v0s, v1_norm);
223 :
224 : // Return the max of the two quotients
225 777 : return std::max(max0/min0, max1/min1);
226 207 : };
227 :
228 420 : return std::max(std::max(vertex_aspect_ratio(0), vertex_aspect_ratio(1)), vertex_aspect_ratio(2)) / std::sqrt(3);
229 : }
230 :
231 : /**
232 : * Source: Netgen, meshtool.cpp, TriangleQualityInst
233 : */
234 0 : case DISTORTION:
235 : case STRETCH:
236 : {
237 0 : const Point & p1 = this->point(0);
238 0 : const Point & p2 = this->point(1);
239 0 : const Point & p3 = this->point(2);
240 :
241 0 : Point v1 = p2 - p1;
242 0 : Point v2 = p3 - p1;
243 0 : Point v3 = p3 - p2;
244 0 : const Real l1 = v1.norm();
245 0 : const Real l2 = v2.norm();
246 0 : const Real l3 = v3.norm();
247 :
248 : // if one length is 0, quality is quite bad!
249 0 : if ((l1 <=0.) || (l2 <= 0.) || (l3 <= 0.))
250 0 : return 0.;
251 :
252 0 : const Real s1 = std::sin(std::acos(v1*v2/l1/l2)/2.);
253 0 : v1 *= -1;
254 0 : const Real s2 = std::sin(std::acos(v1*v3/l1/l3)/2.);
255 0 : const Real s3 = std::sin(std::acos(v2*v3/l2/l3)/2.);
256 :
257 0 : return 8. * s1 * s2 * s3;
258 :
259 : }
260 :
261 : // From: P. Knupp, "Algebraic mesh quality metrics for
262 : // unstructured initial meshes," Finite Elements in Analysis
263 : // and Design 39, 2003, p. 217-241, Section 3.2.
264 0 : case SHAPE:
265 : {
266 : // Unlike Quads, the Tri SHAPE metric is independent of the
267 : // node at which it is computed, we choose to compute it for
268 : // node 0.
269 :
270 : // The nodal Jacobian matrix A is a 3x2 matrix, hence we
271 : // represent it by a std:array with 6 entries.
272 : Point
273 0 : d01 = point(1) - point(0),
274 0 : d02 = point(2) - point(0);
275 :
276 : std::array<Real, 6> A =
277 0 : {{d01(0), d02(0),
278 0 : d01(1), d02(1),
279 0 : d01(2), d02(2)}};
280 :
281 : // Compute metric tensor entries, T = A^T * A.
282 : // This is a symmetric 2x2 matrix so we only
283 : // compute one of the off-diagonal entries.
284 : // As in the paper, we define lambda_ij := T_ij.
285 : Real
286 0 : lambda11 = A[0]*A[0] + A[2]*A[2] + A[4]*A[4],
287 0 : lambda12 = A[0]*A[1] + A[2]*A[3] + A[4]*A[5],
288 0 : lambda22 = A[1]*A[1] + A[3]*A[3] + A[5]*A[5];
289 :
290 : // Compute the denominator of the metric. If it is exactly
291 : // zero then return 0 (lowest quality) for this metric.
292 0 : Real den = lambda11 + lambda22 - lambda12;
293 0 : if (den == 0.0)
294 0 : return 0.;
295 :
296 : // Compute the nodal area
297 0 : Real alpha = std::sqrt(lambda11 * lambda22 - lambda12 * lambda12);
298 :
299 : // Finally, compute and return the metric.
300 0 : return std::sqrt(3) * alpha / den;
301 : }
302 :
303 1920 : default:
304 1920 : return Elem::quality(q);
305 : }
306 :
307 : // We won't get here.
308 : return Elem::quality(q);
309 : }
310 :
311 :
312 :
313 :
314 :
315 :
316 0 : std::pair<Real, Real> Tri::qual_bounds (const ElemQuality q) const
317 : {
318 0 : std::pair<Real, Real> bounds;
319 :
320 0 : switch (q)
321 : {
322 : // A recent copy of the cubit manual [0] does not list bounds
323 : // for EDGE_LENGTH_RATIO or ASPECT_RATIO quality metrics, so we
324 : // have arbitrarily adopted the same values used for Quads here.
325 : // I'm open to suggestions of other appropriate values.
326 : //
327 : // [0]: https://cubit.sandia.gov/files/cubit/16.08/help_manual/WebHelp/mesh_generation/mesh_quality_assessment/triangular_metrics.htm
328 0 : case EDGE_LENGTH_RATIO:
329 : case ASPECT_RATIO:
330 0 : bounds.first = 1.;
331 0 : bounds.second = 4.;
332 0 : break;
333 :
334 0 : case MAX_ANGLE:
335 0 : bounds.first = 60.;
336 0 : bounds.second = 90.;
337 0 : break;
338 :
339 0 : case MIN_ANGLE:
340 0 : bounds.first = 30.;
341 0 : bounds.second = 60.;
342 0 : break;
343 :
344 0 : case CONDITION:
345 0 : bounds.first = 1.;
346 0 : bounds.second = 1.3;
347 0 : break;
348 :
349 0 : case JACOBIAN:
350 : case SCALED_JACOBIAN:
351 0 : bounds.first = 0.5;
352 0 : bounds.second = 1.155;
353 0 : break;
354 :
355 0 : case SIZE:
356 : case SHAPE:
357 0 : bounds.first = 0.25;
358 0 : bounds.second = 1.;
359 0 : break;
360 :
361 0 : case DISTORTION:
362 0 : bounds.first = 0.6;
363 0 : bounds.second = 1.;
364 0 : break;
365 :
366 0 : default:
367 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
368 0 : bounds.first = -1;
369 0 : bounds.second = -1;
370 : }
371 :
372 0 : return bounds;
373 : }
374 :
375 :
376 14153942 : bool Tri::on_reference_element(const Point & p,
377 : const Real eps) const
378 : {
379 2941960 : const Real & xi = p(0);
380 2941960 : const Real & eta = p(1);
381 27566558 : return ((xi >= 0.-eps) &&
382 17033758 : (eta >= 0.-eps) &&
383 15776148 : ((xi + eta) <= 1.+eps));
384 : }
385 :
386 :
387 : } // namespace libMesh
|
