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/edge_edge2.h"
20 : #include "libmesh/face_tri3.h"
21 : #include "libmesh/enum_io_package.h"
22 : #include "libmesh/enum_order.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 :
28 :
29 : // ------------------------------------------------------------
30 : // Tri3 class static member initializations
31 : const int Tri3::num_nodes;
32 : const int Tri3::nodes_per_side;
33 :
34 : const unsigned int Tri3::side_nodes_map[Tri3::num_sides][Tri3::nodes_per_side] =
35 : {
36 : {0, 1}, // Side 0
37 : {1, 2}, // Side 1
38 : {2, 0} // Side 2
39 : };
40 :
41 : #ifdef LIBMESH_ENABLE_AMR
42 :
43 : const Real Tri3::_embedding_matrix[Tri3::num_children][Tri3::num_nodes][Tri3::num_nodes] =
44 : {
45 : // embedding matrix for child 0
46 : {
47 : // 0 1 2
48 : {1.0, 0.0, 0.0}, // 0
49 : {0.5, 0.5, 0.0}, // 1
50 : {0.5, 0.0, 0.5} // 2
51 : },
52 :
53 : // embedding matrix for child 1
54 : {
55 : // 0 1 2
56 : {0.5, 0.5, 0.0}, // 0
57 : {0.0, 1.0, 0.0}, // 1
58 : {0.0, 0.5, 0.5} // 2
59 : },
60 :
61 : // embedding matrix for child 2
62 : {
63 : // 0 1 2
64 : {0.5, 0.0, 0.5}, // 0
65 : {0.0, 0.5, 0.5}, // 1
66 : {0.0, 0.0, 1.0} // 2
67 : },
68 :
69 : // embedding matrix for child 3
70 : {
71 : // 0 1 2
72 : {0.5, 0.5, 0.0}, // 0
73 : {0.0, 0.5, 0.5}, // 1
74 : {0.5, 0.0, 0.5} // 2
75 : }
76 : };
77 :
78 : #endif
79 :
80 :
81 :
82 : // ------------------------------------------------------------
83 : // Tri3 class member functions
84 :
85 21065196 : bool Tri3::is_vertex(const unsigned int libmesh_dbg_var(n)) const
86 : {
87 730920 : libmesh_assert_not_equal_to (n, invalid_uint);
88 21065196 : return true;
89 : }
90 :
91 0 : bool Tri3::is_edge(const unsigned int) const
92 : {
93 0 : return false;
94 : }
95 :
96 0 : bool Tri3::is_face(const unsigned int) const
97 : {
98 0 : return false;
99 : }
100 :
101 708738 : bool Tri3::is_node_on_side(const unsigned int n,
102 : const unsigned int s) const
103 : {
104 21344 : libmesh_assert_less (s, n_sides());
105 21344 : return std::find(std::begin(side_nodes_map[s]),
106 708738 : std::end(side_nodes_map[s]),
107 708738 : n) != std::end(side_nodes_map[s]);
108 : }
109 :
110 : std::vector<unsigned>
111 22278 : Tri3::nodes_on_side(const unsigned int s) const
112 : {
113 1526 : libmesh_assert_less(s, n_sides());
114 23804 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
115 : }
116 :
117 : std::vector<unsigned>
118 2517 : Tri3::nodes_on_edge(const unsigned int e) const
119 : {
120 2517 : return nodes_on_side(e);
121 : }
122 :
123 37847602 : Order Tri3::default_order() const
124 : {
125 37847602 : return FIRST;
126 : }
127 :
128 576 : bool Tri3::has_invertible_map(Real tol) const
129 : {
130 576 : return this->volume() > tol;
131 : }
132 :
133 8087432 : std::unique_ptr<Elem> Tri3::build_side_ptr (const unsigned int i)
134 : {
135 8087432 : return this->simple_build_side_ptr<Edge2, Tri3>(i);
136 : }
137 :
138 :
139 :
140 618009 : void Tri3::build_side_ptr (std::unique_ptr<Elem> & side,
141 : const unsigned int i)
142 : {
143 618009 : this->simple_build_side_ptr<Tri3>(side, i, EDGE2);
144 618009 : }
145 :
146 :
147 :
148 135744 : void Tri3::connectivity(const unsigned int libmesh_dbg_var(sf),
149 : const IOPackage iop,
150 : std::vector<dof_id_type> & conn) const
151 : {
152 11312 : libmesh_assert_less (sf, this->n_sub_elem());
153 11312 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
154 :
155 135744 : switch (iop)
156 : {
157 135744 : case TECPLOT:
158 : {
159 135744 : conn.resize(4);
160 147056 : conn[0] = this->node_id(0)+1;
161 135744 : conn[1] = this->node_id(1)+1;
162 135744 : conn[2] = this->node_id(2)+1;
163 135744 : conn[3] = this->node_id(2)+1;
164 135744 : return;
165 : }
166 :
167 0 : case VTK:
168 : {
169 0 : conn.resize(3);
170 0 : conn[0] = this->node_id(0);
171 0 : conn[1] = this->node_id(1);
172 0 : conn[2] = this->node_id(2);
173 0 : return;
174 : }
175 :
176 0 : default:
177 0 : libmesh_error_msg("Unsupported IO package " << iop);
178 : }
179 : }
180 :
181 :
182 :
183 :
184 :
185 :
186 839 : Point Tri3::true_centroid () const
187 : {
188 839 : return Elem::vertex_average();
189 : }
190 :
191 849139 : Real Tri3::volume () const
192 : {
193 : // This specialization is good for non-Lagrange mappings too! Even
194 : // with differing vertex weights we've still got a p=1 planar face.
195 :
196 : // 3-node triangles have the following formula for computing the area
197 849139 : return 0.5 * cross_norm(point(1) - point(0),
198 880696 : point(2) - point(0));
199 : }
200 :
201 :
202 :
203 0 : std::pair<Real, Real> Tri3::min_and_max_angle() const
204 : {
205 0 : Point v10 ( this->point(1) - this->point(0) );
206 0 : Point v20 ( this->point(2) - this->point(0) );
207 0 : Point v21 ( this->point(2) - this->point(1) );
208 :
209 : const Real
210 0 : len_10=v10.norm(),
211 0 : len_20=v20.norm(),
212 0 : len_21=v21.norm();
213 :
214 : const Real
215 0 : theta0=std::acos(( v10*v20)/len_10/len_20),
216 0 : theta1=std::acos((-v10*v21)/len_10/len_21),
217 0 : theta2=libMesh::pi - theta0 - theta1
218 : ;
219 :
220 0 : libmesh_assert_greater (theta0, 0.);
221 0 : libmesh_assert_greater (theta1, 0.);
222 0 : libmesh_assert_greater (theta2, 0.);
223 :
224 : return std::make_pair(std::min(theta0, std::min(theta1,theta2)),
225 0 : std::max(theta0, std::max(theta1,theta2)));
226 : }
227 :
228 2154954 : bool Tri3::contains_point (const Point & p, Real tol) const
229 : {
230 : // See "Barycentric Technique" section at
231 : // http://www.blackpawn.com/texts/pointinpoly for details.
232 :
233 : // Compute vectors
234 378828 : Point v0 = this->point(1) - this->point(0);
235 189414 : Point v1 = this->point(2) - this->point(0);
236 189414 : Point v2 = p - this->point(0);
237 :
238 : // Compute dot products
239 2154954 : Real dot00 = v0 * v0;
240 189414 : Real dot01 = v0 * v1;
241 189414 : Real dot02 = v0 * v2;
242 2154954 : Real dot11 = v1 * v1;
243 189414 : Real dot12 = v1 * v2;
244 :
245 : // Out of plane check
246 2828622 : if (std::abs(triple_product(v2, v0, v1)) / std::max(dot00, dot11) > tol)
247 6 : return false;
248 :
249 : // Compute barycentric coordinates
250 2154741 : Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
251 2154741 : Real u = (dot11 * dot02 - dot01 * dot12) * invDenom;
252 2154741 : Real v = (dot00 * dot12 - dot01 * dot02) * invDenom;
253 :
254 : // Check if point is in triangle
255 2154741 : return (u > -tol) && (v > -tol) && (u + v < 1 + tol);
256 : }
257 :
258 : BoundingBox
259 16178262 : Tri3::loose_bounding_box () const
260 : {
261 16178262 : return Elem::loose_bounding_box();
262 : }
263 :
264 :
265 4536 : void Tri3::permute(unsigned int perm_num)
266 : {
267 408 : libmesh_assert_less (perm_num, 3);
268 :
269 9072 : for (unsigned int i = 0; i != perm_num; ++i)
270 : {
271 4536 : swap3nodes(0,1,2);
272 4128 : swap3neighbors(0,1,2);
273 : }
274 4536 : }
275 :
276 :
277 3132 : void Tri3::flip(BoundaryInfo * boundary_info)
278 : {
279 120 : libmesh_assert(boundary_info);
280 :
281 3132 : swap2nodes(0,1);
282 120 : swap2neighbors(1,2);
283 3132 : swap2boundarysides(1,2,boundary_info);
284 3132 : swap2boundaryedges(1,2,boundary_info);
285 3132 : }
286 :
287 :
288 : ElemType
289 401644 : Tri3::side_type (const unsigned int libmesh_dbg_var(s)) const
290 : {
291 34760 : libmesh_assert_less (s, 3);
292 401644 : return EDGE2;
293 : }
294 :
295 :
296 : Point
297 426 : Tri3::side_vertex_average_normal(const unsigned int s) const
298 : {
299 12 : libmesh_assert_less (s, 3);
300 450 : const Point side_t = this->point(side_nodes_map[s][1]) -
301 426 : this->point(side_nodes_map[s][0]);
302 426 : const unsigned int another_side = s > 1 ? 0 : s + 1;
303 450 : const Point other_side_t = this->point(side_nodes_map[another_side][1]) -
304 426 : this->point(side_nodes_map[another_side][0]);
305 414 : const Point pointing_up = side_t.cross(other_side_t);
306 414 : const Point v = side_t.cross(pointing_up);
307 426 : return v.unit();
308 : }
309 :
310 :
311 : } // namespace libMesh
|
