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/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 20134308 : bool Tri3::is_vertex(const unsigned int libmesh_dbg_var(n)) const
86 : {
87 703056 : libmesh_assert_not_equal_to (n, invalid_uint);
88 20134308 : 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 572024 : bool Tri3::is_node_on_side(const unsigned int n,
102 : const unsigned int s) const
103 : {
104 16304 : libmesh_assert_less (s, n_sides());
105 16304 : return std::find(std::begin(side_nodes_map[s]),
106 572024 : std::end(side_nodes_map[s]),
107 572024 : n) != std::end(side_nodes_map[s]);
108 : }
109 :
110 : std::vector<unsigned>
111 35586 : Tri3::nodes_on_side(const unsigned int s) const
112 : {
113 1692 : libmesh_assert_less(s, n_sides());
114 37278 : 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 37743698 : Order Tri3::default_order() const
124 : {
125 37743698 : return FIRST;
126 : }
127 :
128 576 : bool Tri3::has_invertible_map(Real tol) const
129 : {
130 576 : return this->volume() > tol;
131 : }
132 :
133 8087149 : std::unique_ptr<Elem> Tri3::build_side_ptr (const unsigned int i)
134 : {
135 8087149 : return this->simple_build_side_ptr<Edge2, Tri3>(i);
136 : }
137 :
138 :
139 :
140 595551 : void Tri3::build_side_ptr (std::unique_ptr<Elem> & side,
141 : const unsigned int i)
142 : {
143 595551 : this->simple_build_side_ptr<Tri3>(side, i, EDGE2);
144 595551 : }
145 :
146 :
147 :
148 287418 : void Tri3::connectivity(const unsigned int libmesh_dbg_var(sf),
149 : const IOPackage iop,
150 : std::vector<dof_id_type> & conn) const
151 : {
152 23952 : libmesh_assert_less (sf, this->n_sub_elem());
153 23952 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
154 :
155 287418 : switch (iop)
156 : {
157 287418 : case TECPLOT:
158 : {
159 287418 : conn.resize(4);
160 311370 : conn[0] = this->node_id(0)+1;
161 287418 : conn[1] = this->node_id(1)+1;
162 287418 : conn[2] = this->node_id(2)+1;
163 287418 : conn[3] = this->node_id(2)+1;
164 287418 : 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 849408 : Real Tri3::volume () const
192 : {
193 : // 3-node triangles have the following formula for computing the area
194 849408 : return 0.5 * cross_norm(point(1) - point(0),
195 880733 : point(2) - point(0));
196 : }
197 :
198 :
199 :
200 0 : std::pair<Real, Real> Tri3::min_and_max_angle() const
201 : {
202 0 : Point v10 ( this->point(1) - this->point(0) );
203 0 : Point v20 ( this->point(2) - this->point(0) );
204 0 : Point v21 ( this->point(2) - this->point(1) );
205 :
206 : const Real
207 0 : len_10=v10.norm(),
208 0 : len_20=v20.norm(),
209 0 : len_21=v21.norm();
210 :
211 : const Real
212 0 : theta0=std::acos(( v10*v20)/len_10/len_20),
213 0 : theta1=std::acos((-v10*v21)/len_10/len_21),
214 0 : theta2=libMesh::pi - theta0 - theta1
215 : ;
216 :
217 0 : libmesh_assert_greater (theta0, 0.);
218 0 : libmesh_assert_greater (theta1, 0.);
219 0 : libmesh_assert_greater (theta2, 0.);
220 :
221 : return std::make_pair(std::min(theta0, std::min(theta1,theta2)),
222 0 : std::max(theta0, std::max(theta1,theta2)));
223 : }
224 :
225 2114516 : bool Tri3::contains_point (const Point & p, Real tol) const
226 : {
227 : // See "Barycentric Technique" section at
228 : // http://www.blackpawn.com/texts/pointinpoly for details.
229 :
230 : // Compute vectors
231 376348 : Point v0 = this->point(1) - this->point(0);
232 188174 : Point v1 = this->point(2) - this->point(0);
233 188174 : Point v2 = p - this->point(0);
234 :
235 : // Compute dot products
236 2114516 : Real dot00 = v0 * v0;
237 188174 : Real dot01 = v0 * v1;
238 188174 : Real dot02 = v0 * v2;
239 2114516 : Real dot11 = v1 * v1;
240 188174 : Real dot12 = v1 * v2;
241 :
242 : // Out of plane check
243 2768808 : if (std::abs(triple_product(v2, v0, v1)) / std::max(dot00, dot11) > tol)
244 6 : return false;
245 :
246 : // Compute barycentric coordinates
247 2114303 : Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
248 2114303 : Real u = (dot11 * dot02 - dot01 * dot12) * invDenom;
249 2114303 : Real v = (dot00 * dot12 - dot01 * dot02) * invDenom;
250 :
251 : // Check if point is in triangle
252 2114303 : return (u > -tol) && (v > -tol) && (u + v < 1 + tol);
253 : }
254 :
255 : BoundingBox
256 14438118 : Tri3::loose_bounding_box () const
257 : {
258 14438118 : return Elem::loose_bounding_box();
259 : }
260 :
261 :
262 4536 : void Tri3::permute(unsigned int perm_num)
263 : {
264 408 : libmesh_assert_less (perm_num, 3);
265 :
266 9072 : for (unsigned int i = 0; i != perm_num; ++i)
267 : {
268 4536 : swap3nodes(0,1,2);
269 4128 : swap3neighbors(0,1,2);
270 : }
271 4536 : }
272 :
273 :
274 576 : void Tri3::flip(BoundaryInfo * boundary_info)
275 : {
276 48 : libmesh_assert(boundary_info);
277 :
278 576 : swap2nodes(0,1);
279 48 : swap2neighbors(1,2);
280 576 : swap2boundarysides(1,2,boundary_info);
281 576 : swap2boundaryedges(1,2,boundary_info);
282 576 : }
283 :
284 :
285 : ElemType
286 401608 : Tri3::side_type (const unsigned int libmesh_dbg_var(s)) const
287 : {
288 34760 : libmesh_assert_less (s, 3);
289 401608 : return EDGE2;
290 : }
291 :
292 : } // namespace libMesh
|
