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