libMesh
face_tri.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 
34 
35 // We need to require C++11...
36 const Real Tri::_master_points[6][3] =
37  {
38  {0, 0},
39  {1, 0},
40  {0, 1},
41  {0.5, 0},
42  {0.5, 0.5},
43  {0, 0.5}
44  };
45 
46 
47 
48 
49 
50 
51 
52 // ------------------------------------------------------------
53 // Tri class member functions
54 dof_id_type Tri::key (const unsigned int s) const
55 {
56  libmesh_assert_less (s, this->n_sides());
57 
58  return this->compute_key(this->node_id(Tri3::side_nodes_map[s][0]),
59  this->node_id(Tri3::side_nodes_map[s][1]));
60 }
61 
62 
63 
64 dof_id_type Tri::low_order_key (const unsigned int s) const
65 {
66  libmesh_assert_less (s, this->n_sides());
67 
68  return this->compute_key(this->node_id(Tri3::side_nodes_map[s][0]),
69  this->node_id(Tri3::side_nodes_map[s][1]));
70 }
71 
72 
73 
74 unsigned int Tri::local_side_node(unsigned int side,
75  unsigned int side_node) const
76 {
77  libmesh_assert_less (side, this->n_sides());
78  libmesh_assert_less (side_node, Tri3::nodes_per_side);
79 
80  return Tri3::side_nodes_map[side][side_node];
81 }
82 
83 
84 
85 unsigned int Tri::local_edge_node(unsigned int edge,
86  unsigned int edge_node) const
87 {
88  return local_side_node(edge, edge_node);
89 }
90 
91 
92 
94 {
95  return this->compute_key(this->node_id(0),
96  this->node_id(1),
97  this->node_id(2));
98 }
99 
100 
101 
102 std::unique_ptr<Elem> Tri::side_ptr (const unsigned int i)
103 {
104  libmesh_assert_less (i, this->n_sides());
105 
106  std::unique_ptr<Elem> edge = std::make_unique<Edge2>();
107 
108  for (auto n : edge->node_index_range())
109  edge->set_node(n) = this->node_ptr(Tri3::side_nodes_map[i][n]);
110 
111  return edge;
112 }
113 
114 
115 
116 void Tri::side_ptr (std::unique_ptr<Elem> & side,
117  const unsigned int i)
118 {
119  this->simple_side_ptr<Tri,Tri3>(side, i, EDGE2);
120 }
121 
122 
123 
124 bool Tri::is_child_on_side(const unsigned int c,
125  const unsigned int s) const
126 {
127  libmesh_assert_less (c, this->n_children());
128  libmesh_assert_less (s, this->n_sides());
129 
130  return (c == s || c == (s+1)%3);
131 }
132 
133 
134 bool Tri::is_flipped() const
135 {
136  return (
137 #if LIBMESH_DIM > 2
138  // Don't bother outside the XY plane
139  !this->point(0)(2) && !this->point(1)(2) &&
140  !this->point(2)(2) &&
141 #endif
142  ((this->point(1)(0)-this->point(0)(0))*
143  (this->point(2)(1)-this->point(0)(1)) <
144  (this->point(2)(0)-this->point(0)(0))*
145  (this->point(1)(1)-this->point(0)(1))));
146 }
147 
148 
149 
151 {
152  switch (q)
153  {
154 
158  case DISTORTION:
159  case STRETCH:
160  {
161  const Point & p1 = this->point(0);
162  const Point & p2 = this->point(1);
163  const Point & p3 = this->point(2);
164 
165  Point v1 = p2 - p1;
166  Point v2 = p3 - p1;
167  Point v3 = p3 - p2;
168  const Real l1 = v1.norm();
169  const Real l2 = v2.norm();
170  const Real l3 = v3.norm();
171 
172  // if one length is 0, quality is quite bad!
173  if ((l1 <=0.) || (l2 <= 0.) || (l3 <= 0.))
174  return 0.;
175 
176  const Real s1 = std::sin(std::acos(v1*v2/l1/l2)/2.);
177  v1 *= -1;
178  const Real s2 = std::sin(std::acos(v1*v3/l1/l3)/2.);
179  const Real s3 = std::sin(std::acos(v2*v3/l2/l3)/2.);
180 
181  return 8. * s1 * s2 * s3;
182 
183  }
184 
185  // From: P. Knupp, "Algebraic mesh quality metrics for
186  // unstructured initial meshes," Finite Elements in Analysis
187  // and Design 39, 2003, p. 217-241, Section 3.2.
188  case SHAPE:
189  {
190  // Unlike Quads, the Tri SHAPE metric is independent of the
191  // node at which it is computed, we choose to compute it for
192  // node 0.
193 
194  // The nodal Jacobian matrix A is a 3x2 matrix, hence we
195  // represent it by a std:array with 6 entries.
196  Point
197  d01 = point(1) - point(0),
198  d02 = point(2) - point(0);
199 
200  std::array<Real, 6> A =
201  {{d01(0), d02(0),
202  d01(1), d02(1),
203  d01(2), d02(2)}};
204 
205  // Compute metric tensor entries, T = A^T * A.
206  // This is a symmetric 2x2 matrix so we only
207  // compute one of the off-diagonal entries.
208  // As in the paper, we define lambda_ij := T_ij.
209  Real
210  lambda11 = A[0]*A[0] + A[2]*A[2] + A[4]*A[4],
211  lambda12 = A[0]*A[1] + A[2]*A[3] + A[4]*A[5],
212  lambda22 = A[1]*A[1] + A[3]*A[3] + A[5]*A[5];
213 
214  // Compute the denominator of the metric. If it is exactly
215  // zero then return 0 (lowest quality) for this metric.
216  Real den = lambda11 + lambda22 - lambda12;
217  if (den == 0.0)
218  return 0.;
219 
220  // Compute the nodal area
221  Real alpha = std::sqrt(lambda11 * lambda22 - lambda12 * lambda12);
222 
223  // Finally, compute and return the metric.
224  return std::sqrt(3) * alpha / den;
225  }
226 
227  default:
228  return Elem::quality(q);
229  }
230 
231  // We won't get here.
232  return Elem::quality(q);
233 }
234 
235 
236 
237 
238 
239 
240 std::pair<Real, Real> Tri::qual_bounds (const ElemQuality q) const
241 {
242  std::pair<Real, Real> bounds;
243 
244  switch (q)
245  {
246 
247  case MAX_ANGLE:
248  bounds.first = 60.;
249  bounds.second = 90.;
250  break;
251 
252  case MIN_ANGLE:
253  bounds.first = 30.;
254  bounds.second = 60.;
255  break;
256 
257  case CONDITION:
258  bounds.first = 1.;
259  bounds.second = 1.3;
260  break;
261 
262  case JACOBIAN:
263  bounds.first = 0.5;
264  bounds.second = 1.155;
265  break;
266 
267  case SIZE:
268  case SHAPE:
269  bounds.first = 0.25;
270  bounds.second = 1.;
271  break;
272 
273  case DISTORTION:
274  bounds.first = 0.6;
275  bounds.second = 1.;
276  break;
277 
278  default:
279  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
280  bounds.first = -1;
281  bounds.second = -1;
282  }
283 
284  return bounds;
285 }
286 
287 } // namespace libMesh
virtual std::unique_ptr< Elem > side_ptr(const unsigned int i) override final
Definition: face_tri.C:102
static const Real _master_points[6][3]
Master element node locations.
Definition: face_tri.h:197
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:929
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Calls local_side_node(edge, edge_node).
Definition: face_tri.C:85
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition: face_tri3.h:169
virtual bool is_flipped() const override final
Definition: face_tri.C:134
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int n_sides() const override final
Definition: face_tri.h:92
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const override final
Definition: face_tri.C:124
virtual dof_id_type key() const override final
Definition: face_tri.C:93
virtual dof_id_type low_order_key(const unsigned int s) const override
Definition: face_tri.C:64
ElemQuality
Defines an enum for element quality metrics.
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition: face_tri.C:74
static const int nodes_per_side
Definition: face_tri3.h:163
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2331
OStreamProxy out
virtual Real quality(const ElemQuality q) const
Definition: elem.C:1582
virtual std::pair< Real, Real > qual_bounds(const ElemQuality q) const override
Definition: face_tri.C:240
virtual unsigned int n_children() const override final
Definition: face_tri.h:107
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:3131
virtual Real quality(const ElemQuality q) const override
Definition: face_tri.C:150
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2299
const Point & point(const unsigned int i) const
Definition: elem.h:2277
uint8_t dof_id_type
Definition: id_types.h:67