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_quad.h"
20 : #include "libmesh/edge_edge2.h"
21 : #include "libmesh/face_quad4.h"
22 : #include "libmesh/enum_elem_quality.h"
23 :
24 : // C++ includes
25 : #include <array>
26 :
27 :
28 : namespace libMesh
29 : {
30 :
31 : // ------------------------------------------------------------
32 : // Quad class static member initializations
33 : const int Quad::num_sides;
34 : const int Quad::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 Quad::_master_points[9][3] =
39 : {
40 : {-1, -1},
41 : {1, -1},
42 : {1, 1},
43 : {-1, 1},
44 : {0, -1},
45 : {1, 0},
46 : {0, 1},
47 : {-1, 0},
48 : {0, 0}
49 : };
50 :
51 : const unsigned int Quad::adjacent_sides_map[/*num_vertices*/4][/*n_adjacent_sides*/2] =
52 : {
53 : {0, 3}, // Sides adjacent to node 0
54 : {0, 1}, // Sides adjacent to node 1
55 : {1, 2}, // Sides adjacent to node 2
56 : {2, 3} // Sides adjacent to node 3
57 : };
58 :
59 :
60 :
61 : // ------------------------------------------------------------
62 : // Quad class member functions
63 0 : dof_id_type Quad::key (const unsigned int s) const
64 : {
65 0 : libmesh_assert_less (s, this->n_sides());
66 :
67 0 : return this->compute_key(this->node_id(Quad4::side_nodes_map[s][0]),
68 0 : this->node_id(Quad4::side_nodes_map[s][1]));
69 : }
70 :
71 :
72 :
73 146924694 : dof_id_type Quad::low_order_key (const unsigned int s) const
74 : {
75 1890688 : libmesh_assert_less (s, this->n_sides());
76 :
77 148815382 : return this->compute_key(this->node_id(Quad4::side_nodes_map[s][0]),
78 148815382 : this->node_id(Quad4::side_nodes_map[s][1]));
79 : }
80 :
81 :
82 :
83 65513642 : unsigned int Quad::local_side_node(unsigned int side,
84 : unsigned int side_node) const
85 : {
86 5782800 : libmesh_assert_less (side, this->n_sides());
87 5782800 : libmesh_assert_less (side_node, Quad4::nodes_per_side);
88 :
89 65513642 : return Quad4::side_nodes_map[side][side_node];
90 : }
91 :
92 :
93 :
94 153171132 : unsigned int Quad::local_edge_node(unsigned int edge,
95 : unsigned int edge_node) const
96 : {
97 153171132 : return local_side_node(edge, edge_node);
98 : }
99 :
100 :
101 :
102 0 : dof_id_type Quad::key () const
103 : {
104 0 : return this->compute_key(this->node_id(0),
105 : this->node_id(1),
106 : this->node_id(2),
107 0 : this->node_id(3));
108 : }
109 :
110 :
111 :
112 239886 : std::unique_ptr<Elem> Quad::side_ptr (const unsigned int i)
113 : {
114 6674 : libmesh_assert_less (i, this->n_sides());
115 :
116 239886 : std::unique_ptr<Elem> edge = std::make_unique<Edge2>();
117 :
118 719658 : for (auto n : edge->node_index_range())
119 493120 : edge->set_node(n, this->node_ptr(Quad4::side_nodes_map[i][n]));
120 :
121 239886 : return edge;
122 0 : }
123 :
124 :
125 :
126 140057957 : void Quad::side_ptr (std::unique_ptr<Elem> & side,
127 : const unsigned int i)
128 : {
129 140057957 : this->simple_side_ptr<Quad,Quad4>(side, i, EDGE2);
130 140057957 : }
131 :
132 :
133 :
134 12451391 : bool Quad::is_child_on_side(const unsigned int c,
135 : const unsigned int s) const
136 : {
137 8579529 : libmesh_assert_less (c, this->n_children());
138 8579529 : libmesh_assert_less (s, this->n_sides());
139 :
140 : // A quad's children and nodes don't share the same ordering:
141 : // child 2 and 3 are swapped;
142 12451391 : unsigned int n = (c < 2) ? c : 5-c;
143 12451391 : return (n == s || n == (s+1)%4);
144 : }
145 :
146 :
147 :
148 0 : unsigned int Quad::opposite_side(const unsigned int side_in) const
149 : {
150 0 : libmesh_assert_less (side_in, 4);
151 :
152 0 : return (side_in + 2) % 4;
153 : }
154 :
155 :
156 :
157 0 : unsigned int Quad::opposite_node(const unsigned int node_in,
158 : const unsigned int side_in) const
159 : {
160 0 : libmesh_assert_less (node_in, 8);
161 0 : libmesh_assert_less (node_in, this->n_nodes());
162 0 : libmesh_assert_less (side_in, this->n_sides());
163 0 : libmesh_assert(this->is_node_on_side(node_in, side_in));
164 :
165 : static const unsigned char side02_nodes_map[] =
166 : {3, 2, 1, 0, 6, 255, 4, 255};
167 : static const unsigned char side13_nodes_map[] =
168 : {1, 0, 3, 2, 255, 7, 255, 5};
169 :
170 0 : switch (side_in)
171 : {
172 0 : case 0:
173 : case 2:
174 0 : return side02_nodes_map[node_in];
175 0 : case 1:
176 : case 3:
177 0 : return side13_nodes_map[node_in];
178 0 : default:
179 0 : libmesh_error_msg("Unsupported side_in = " << side_in);
180 : }
181 : }
182 :
183 :
184 2280 : bool Quad::is_flipped() const
185 : {
186 : return (
187 : #if LIBMESH_DIM > 2
188 : // Don't bother outside the XY plane
189 2448 : !this->point(0)(2) && !this->point(1)(2) &&
190 4728 : !this->point(2)(2) && !this->point(3)(2) &&
191 : #endif
192 2448 : ((this->point(1)(0)-this->point(0)(0))*
193 2280 : (this->point(3)(1)-this->point(0)(1)) <
194 2448 : (this->point(3)(0)-this->point(0)(0))*
195 2448 : (this->point(1)(1)-this->point(0)(1))));
196 : }
197 :
198 :
199 : std::vector<unsigned int>
200 14624 : Quad::edges_adjacent_to_node(const unsigned int n) const
201 : {
202 968 : libmesh_assert_less(n, this->n_nodes());
203 :
204 : // For vertices, we use the Quad::adjacent_sides_map, otherwise each
205 : // of the mid-edge nodes is adjacent only to the edge it is on, and the
206 : // center node is not adjacent to any edge.
207 14624 : if (this->is_vertex(n))
208 10912 : return {std::begin(adjacent_sides_map[n]), std::end(adjacent_sides_map[n])};
209 4320 : else if (this->is_edge(n))
210 3840 : return {n - this->n_vertices()};
211 :
212 40 : libmesh_assert(this->is_face(n));
213 440 : return {};
214 : }
215 :
216 3783 : Real Quad::quality (const ElemQuality q) const
217 : {
218 3783 : switch (q)
219 : {
220 : // The EDGE_LENGTH_RATIO metric is handled by the base class.
221 : //
222 : // The CUBIT 15.1 User Documentation refers to the
223 : // EDGE_LENGTH_RATIO as the "Aspect Ratio" metric for Quads,
224 : // however, we prefer the slightly more robust aspect ratio
225 : // formula employed by Ansys, and have therefore designated that
226 : // as our ASPECT_RATIO metric here.
227 : //
228 : // As a counter-example for employing the EDGE_LENGTH_RATIO in
229 : // Quads, consider a "rhombus". All four side lengths are equal,
230 : // so the EDGE_LENGTH_RATIO of this element is always (the
231 : // optimal value) 1.0, regardless of the internal angle "theta"
232 : // of the rhombus. A more sensitive shape quality metric should
233 : // take this internal angle into account, so that very thin
234 : // rhombus-shaped Quads are not considered optimal.
235 :
236 781 : case ASPECT_RATIO:
237 : {
238 : // Aspect Ratio definition from Ansys Theory Manual.
239 : // Reference: Ansys, Inc. Theory Reference, Ansys Release 9.0, 2004 (Chapter: 13.7.3)
240 : //
241 : // For the rhombus case described above, the aspect ratio of
242 : // the element using the Ansys Aspect Ratio metric is
243 : // 1/sin(theta), where theta is the acute interior angle of
244 : // the rhombus. The aspect ratio therefore has a minimum value
245 : // of 1.0 (when theta = pi/2), and blows up to infinity as
246 : // sin(theta) -> 0, rather than remaining equal to 1.0 as in
247 : // the CUBIT formula.
248 :
249 : // Compute midpoint positions along each edge
250 44 : Point m0 = Real(0.5) * (this->point(0) + this->point(1));
251 22 : Point m1 = Real(0.5) * (this->point(1) + this->point(2));
252 22 : Point m2 = Real(0.5) * (this->point(2) + this->point(3));
253 22 : Point m3 = Real(0.5) * (this->point(3) + this->point(0));
254 :
255 : // Compute vectors adjoining opposite side midpoints
256 22 : Point v0 = m2 - m0;
257 22 : Point v1 = m1 - m3;
258 :
259 : // Compute the length of the midlines
260 781 : Real v0_norm = v0.norm();
261 781 : Real v1_norm = v1.norm();
262 :
263 : // Instead of dividing by zero in the next step, just return
264 : // 0. The optimal aspect ratio is 1.0, and "high" aspect
265 : // ratios are bad, but an aspect ratio of 0 should also be
266 : // considered bad.
267 781 : if (v0_norm == 0. || v1_norm == 0.)
268 0 : return 0.;
269 :
270 : // Compute sine of the angle between v0, v1. For rectangular
271 : // elements, sin_theta == 1.
272 781 : Real sin_theta = cross_norm(v0, v1) / v0_norm / v1_norm;
273 781 : Real v0s = v0_norm*sin_theta;
274 781 : Real v1s = v1_norm*sin_theta;
275 :
276 : // Determine the min, max of each midline length and its
277 : // projection.
278 : //
279 : // Note: The return values min{0,1} and max{0,1} here are
280 : // *references* and we cannot pass a temporary to
281 : // std::minmax() since the reference returned becomes a
282 : // dangling reference at "the end of the full expression that
283 : // contains the call to std::minmax()"
284 : // https://en.cppreference.com/w/cpp/algorithm/minmax
285 22 : auto [min0, max0] = std::minmax(v0_norm, v1s);
286 22 : auto [min1, max1] = std::minmax(v0s, v1_norm);
287 :
288 : // Return the max of the two quotients
289 919 : return std::max(max0/min0, max1/min1);
290 : }
291 :
292 : // Compute the min/max diagonal ratio.
293 : // This is modeled after the Hex element
294 0 : case DISTORTION:
295 : case DIAGONAL:
296 : {
297 : // Diagonal between node 0 and node 2
298 0 : const Real d02 = this->length(0,2);
299 :
300 : // Diagonal between node 1 and node 3
301 0 : const Real d13 = this->length(1,3);
302 :
303 : // Find the biggest and smallest diagonals
304 0 : if ((d02 > 0.) && (d13 >0.))
305 0 : if (d02 < d13) return d02 / d13;
306 0 : else return d13 / d02;
307 : else
308 0 : return 0.;
309 : break;
310 : }
311 :
312 : // CUBIT 15.1 User Documentation:
313 : // Stretch: Sqrt(2) * minimum edge length / maximum diagonal length
314 0 : case STRETCH:
315 : {
316 0 : Real lengths[4] = {this->length(0,1), this->length(1,2), this->length(2,3), this->length(3,0)};
317 0 : Real min_edge = *std::min_element(lengths, lengths+4);
318 0 : Real d_max = std::max(this->length(0,2), this->length(1,3));
319 :
320 : // Return 0. instead of dividing by zero.
321 0 : if (d_max == 0.)
322 0 : return 0.;
323 : else
324 0 : return std::sqrt(2) * min_edge / d_max;
325 : }
326 :
327 0 : case SHAPE:
328 : case SKEW:
329 : {
330 : // From: P. Knupp, "Algebraic mesh quality metrics for
331 : // unstructured initial meshes," Finite Elements in Analysis
332 : // and Design 39, 2003, p. 217-241, Sections 5.2 and 5.3.
333 : typedef std::array<Real, 4> Array4;
334 : typedef std::array<Real, 6> Array6;
335 :
336 : // x, y, z node coordinates.
337 : std::vector<Real>
338 0 : x = {this->point(0)(0), this->point(1)(0), this->point(2)(0), this->point(3)(0)},
339 0 : y = {this->point(0)(1), this->point(1)(1), this->point(2)(1), this->point(3)(1)},
340 0 : z = {this->point(0)(2), this->point(1)(2), this->point(2)(2), this->point(3)(2)};
341 :
342 : // Nodal Jacobians. These are 3x2 matrices, hence we represent
343 : // them by Array6.
344 : std::array<Array6, 4> A;
345 0 : for (unsigned int k=0; k<4; ++k)
346 : {
347 : unsigned int
348 0 : kp1 = k+1 > 3 ? k+1-4 : k+1,
349 0 : kp3 = k+3 > 3 ? k+3-4 : k+3;
350 :
351 : // To initialize std::array we need double curly braces in
352 : // C++11 but not C++14 apparently.
353 0 : A[k] = {{x[kp1] - x[k], x[kp3] - x[k],
354 0 : y[kp1] - y[k], y[kp3] - y[k],
355 0 : z[kp1] - z[k], z[kp3] - z[k]}};
356 : }
357 :
358 : // Compute metric tensors, T_k = A_k^T * A_k. These are 2x2
359 : // square matrices, hence we represent them by Array4.
360 : std::array<Array4, 4> T;
361 0 : for (unsigned int k=0; k<4; ++k)
362 : {
363 : Real
364 0 : top_left = A[k][0]*A[k][0] + A[k][2]*A[k][2] + A[k][4]*A[k][4],
365 0 : off_diag = A[k][0]*A[k][1] + A[k][2]*A[k][3] + A[k][4]*A[k][5],
366 0 : bot_rigt = A[k][1]*A[k][1] + A[k][3]*A[k][3] + A[k][5]*A[k][5];
367 :
368 0 : T[k] = {{top_left, off_diag,
369 : off_diag, bot_rigt}};
370 : }
371 :
372 :
373 : // Nodal areas. These are approximated as sqrt(det(A^T * A))
374 : // to handle the general case of a 2D element living in 3D.
375 : Array4 alpha;
376 0 : for (unsigned int k=0; k<4; ++k)
377 0 : alpha[k] = std::sqrt(T[k][0]*T[k][3] - T[k][1]*T[k][2]);
378 :
379 : // All nodal areas must be strictly positive. Return 0 (the
380 : // lowest quality) otherwise. We know they can't be negative
381 : // because they are the result of a sqrt, but they might be
382 : // identically 0.
383 0 : if (*std::min_element(alpha.begin(), alpha.end()) == 0.)
384 0 : return 0.;
385 :
386 : // Compute and return the shape metric. These only use the
387 : // diagonal entries of the T_k.
388 0 : Real den = 0.;
389 0 : if (q == SHAPE)
390 : {
391 0 : for (unsigned int k=0; k<4; ++k)
392 0 : den += (T[k][0] + T[k][3]) / alpha[k];
393 0 : return (den == 0.) ? 0 : (8. / den);
394 : }
395 : else
396 : {
397 0 : for (unsigned int k=0; k<4; ++k)
398 0 : den += std::sqrt(T[k][0] * T[k][3]) / alpha[k];
399 0 : return (den == 0.) ? 0 : (4. / den);
400 : }
401 : }
402 :
403 : // This test returns 0 if a Quad:
404 : // 1) Is "twisted" (i.e. has an invalid numbering)
405 : // 2) Has nearly parallel adjacent edges
406 : // 3) Has a nearly zero length edge
407 : // Otherwise, it returns 1.
408 0 : case TWIST:
409 : {
410 : // Compute the cross product induced by each "corner" of the QUAD
411 : // and check that:
412 : // 1) None of the cross products are zero (within a tolernace), and
413 : // 2) all of the cross products point in the same direction.
414 :
415 : // Corner 0
416 0 : Point vec_01 = point(1) - point(0);
417 0 : Point vec_03 = point(3) - point(0);
418 0 : Point corner_0_vec = vec_01.cross(vec_03);
419 :
420 : // Corner 1
421 0 : Point vec_12 = point(2) - point(1);
422 0 : Point vec_10 = point(0) - point(1);
423 0 : Point corner_1_vec = vec_12.cross(vec_10);
424 :
425 : // Corner 2
426 0 : Point vec_23 = point(3) - point(2);
427 0 : Point vec_21 = point(1) - point(2);
428 0 : Point corner_2_vec = vec_23.cross(vec_21);
429 :
430 : // Corner 3
431 0 : Point vec_30 = point(0) - point(3);
432 0 : Point vec_32 = point(2) - point(3);
433 0 : Point corner_3_vec = vec_30.cross(vec_32);
434 :
435 : // If any of these cross products is nearly zero, then either
436 : // we have nearly parallel adjacent edges or a nearly zero
437 : // length edge. We return 0 in this case.
438 0 : if ((corner_0_vec.norm() < TOLERANCE*TOLERANCE) ||
439 0 : (corner_1_vec.norm() < TOLERANCE*TOLERANCE) ||
440 0 : (corner_2_vec.norm() < TOLERANCE*TOLERANCE) ||
441 0 : (corner_3_vec.norm() < TOLERANCE*TOLERANCE))
442 : {
443 0 : return 0.;
444 : }
445 :
446 : // Now check whether the element is twisted.
447 0 : Real dot_01 = corner_0_vec * corner_1_vec;
448 0 : Real dot_02 = corner_0_vec * corner_2_vec;
449 0 : Real dot_03 = corner_0_vec * corner_3_vec;
450 :
451 0 : if ((dot_01 <= 0.) || (dot_02 <= 0.) || (dot_03 <= 0.))
452 0 : return 0.;
453 :
454 : // If we made it here, then Elem is not twisted.
455 0 : return 1.;
456 : }
457 :
458 426 : case WARP:
459 : {
460 : // We compute the angle between the two planes formed by
461 : // drawing an imaginary diagonal separating opposite vertices
462 : // A and B, where (A,B) = (0,2) and (1,3). The element warpage
463 : // is then taken to be the smaller of these two angles. For a
464 : // flat quadrilateral, the planes are parallel, thus the angle
465 : // between them is zero, and the element warpage is therefore
466 : // 1.
467 : //
468 : // Reference: https://cubit.sandia.gov/files/cubit/16.08/help_manual/WebHelp/mesh_generation/mesh_quality_assessment/quadrilateral_metrics.htm
469 : Point
470 438 : n0 = (this->point(1) - this->point(0)).cross(this->point(3) - this->point(0)).unit(),
471 438 : n1 = (this->point(2) - this->point(1)).cross(this->point(0) - this->point(1)).unit(),
472 438 : n2 = (this->point(3) - this->point(2)).cross(this->point(1) - this->point(2)).unit(),
473 438 : n3 = (this->point(0) - this->point(3)).cross(this->point(2) - this->point(3)).unit();
474 :
475 576 : return std::min(n0*n2, n1*n3);
476 : }
477 :
478 2576 : default:
479 2576 : return Elem::quality(q);
480 : }
481 : }
482 :
483 :
484 :
485 :
486 :
487 :
488 0 : std::pair<Real, Real> Quad::qual_bounds (const ElemQuality q) const
489 : {
490 0 : std::pair<Real, Real> bounds;
491 :
492 0 : switch (q)
493 : {
494 0 : case EDGE_LENGTH_RATIO:
495 : case ASPECT_RATIO:
496 0 : bounds.first = 1.;
497 0 : bounds.second = 4.;
498 0 : break;
499 :
500 0 : case TAPER:
501 0 : bounds.first = 0.;
502 0 : bounds.second = 0.7;
503 0 : break;
504 :
505 0 : case WARP:
506 0 : bounds.first = 0.9;
507 0 : bounds.second = 1.;
508 0 : break;
509 :
510 0 : case STRETCH:
511 0 : bounds.first = 0.25;
512 0 : bounds.second = 1.;
513 0 : break;
514 :
515 0 : case MIN_ANGLE:
516 0 : bounds.first = 45.;
517 0 : bounds.second = 90.;
518 0 : break;
519 :
520 0 : case MAX_ANGLE:
521 0 : bounds.first = 90.;
522 0 : bounds.second = 135.;
523 0 : break;
524 :
525 0 : case CONDITION:
526 0 : bounds.first = 1.;
527 0 : bounds.second = 4.;
528 0 : break;
529 :
530 0 : case JACOBIAN:
531 : case SCALED_JACOBIAN:
532 0 : bounds.first = 0.5;
533 0 : bounds.second = 1.;
534 0 : break;
535 :
536 0 : case SHEAR:
537 : case SHAPE:
538 : case SKEW:
539 : case SIZE:
540 0 : bounds.first = 0.3;
541 0 : bounds.second = 1.;
542 0 : break;
543 :
544 0 : case DISTORTION:
545 0 : bounds.first = 0.6;
546 0 : bounds.second = 1.;
547 0 : break;
548 :
549 0 : case TWIST:
550 0 : bounds.first = 0.;
551 0 : bounds.second = 1.;
552 0 : break;
553 :
554 0 : default:
555 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
556 0 : bounds.first = -1;
557 0 : bounds.second = -1;
558 : }
559 :
560 0 : return bounds;
561 : }
562 :
563 :
564 :
565 12432592 : bool Quad::on_reference_element(const Point & p,
566 : const Real eps) const
567 : {
568 1874831 : const Real & xi = p(0);
569 1874831 : const Real & eta = p(1);
570 :
571 : // The reference quadrilateral element is [-1,1]^2.
572 22801869 : return ((xi >= -1.-eps) &&
573 10369277 : (xi <= 1.+eps) &&
574 22362737 : (eta >= -1.-eps) &&
575 14063208 : (eta <= 1.+eps));
576 : }
577 :
578 :
579 :
580 : const unsigned short int Quad::_second_order_adjacent_vertices[4][2] =
581 : {
582 : {0, 1}, // vertices adjacent to node 4
583 : {1, 2}, // vertices adjacent to node 5
584 : {2, 3}, // vertices adjacent to node 6
585 : {0, 3} // vertices adjacent to node 7
586 : };
587 :
588 :
589 :
590 : const unsigned short int Quad::_second_order_vertex_child_number[9] =
591 : {
592 : 99,99,99,99, // Vertices
593 : 0,1,2,0, // Edges
594 : 0 // Interior
595 : };
596 :
597 :
598 :
599 : const unsigned short int Quad::_second_order_vertex_child_index[9] =
600 : {
601 : 99,99,99,99, // Vertices
602 : 1,2,3,3, // Edges
603 : 2 // Interior
604 : };
605 :
606 :
607 :
608 : #ifdef LIBMESH_ENABLE_AMR
609 :
610 : // We number 25 "possible node locations" for a 2x2 refinement of
611 : // quads with up to 3x3 nodes each
612 : const int Quad::_child_node_lookup[4][9] =
613 : {
614 : // node lookup for child 0 (near node 0)
615 : { 0, 2, 12, 10, 1, 7, 11, 5, 6},
616 :
617 : // node lookup for child 1 (near node 1)
618 : { 2, 4, 14, 12, 3, 9, 13, 7, 8},
619 :
620 : // node lookup for child 2 (near node 3)
621 : { 10, 12, 22, 20, 11, 17, 21, 15, 16},
622 :
623 : // node lookup for child 3 (near node 2)
624 : { 12, 14, 24, 22, 13, 19, 23, 17, 18}
625 : };
626 :
627 :
628 : #endif
629 :
630 : } // namespace libMesh
|
