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 :
19 :
20 : // Local includes
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/enum_io_package.h"
23 : #include "libmesh/enum_order.h"
24 :
25 : namespace libMesh
26 : {
27 :
28 : // ------------------------------------------------------------
29 : // Edge3 class static member initializations
30 : const int Edge3::num_nodes;
31 :
32 : #ifdef LIBMESH_ENABLE_AMR
33 :
34 : const Real Edge3::_embedding_matrix[2][3][3] =
35 : {
36 : // embedding matrix for child 0
37 : {
38 : // 0 1 2
39 : {1.0, 0.0, 0.0}, // left
40 : {0.0, 0.0, 1.0}, // right
41 : {0.375,-0.125,0.75} // middle
42 : },
43 :
44 : // embedding matrix for child 1
45 : {
46 : // 0 1 2
47 : {0.0, 0.0, 1.0}, // left
48 : {0.0, 1.0, 0.0}, // right
49 : {-0.125,0.375,0.75} // middle
50 : }
51 : };
52 :
53 : #endif
54 :
55 1296197 : bool Edge3::is_vertex(const unsigned int i) const
56 : {
57 1296197 : if (i < 2)
58 866039 : return true;
59 38737 : return false;
60 : }
61 :
62 0 : bool Edge3::is_edge(const unsigned int i) const
63 : {
64 0 : if (i < 2)
65 0 : return false;
66 0 : return true;
67 : }
68 :
69 0 : bool Edge3::is_face(const unsigned int ) const
70 : {
71 0 : return false;
72 : }
73 :
74 16211 : bool Edge3::is_node_on_side(const unsigned int n,
75 : const unsigned int s) const
76 : {
77 254 : libmesh_assert_less (s, 2);
78 254 : libmesh_assert_less (n, Edge3::num_nodes);
79 16211 : return (s == n);
80 : }
81 :
82 0 : bool Edge3::is_node_on_edge(const unsigned int,
83 : const unsigned int libmesh_dbg_var(e)) const
84 : {
85 0 : libmesh_assert_equal_to (e, 0);
86 0 : return true;
87 : }
88 :
89 :
90 :
91 232049 : bool Edge3::has_affine_map() const
92 : {
93 38326 : Point v = this->point(1) - this->point(0);
94 : return (v.relative_fuzzy_equals
95 232049 : ((this->point(2) - this->point(0))*2, affine_tol));
96 : }
97 :
98 :
99 :
100 213 : bool Edge3::has_invertible_map(Real tol) const
101 : {
102 : // At the moment this only makes sense for Lagrange elements
103 6 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
104 :
105 : // The "Jacobian vector" (dx/dxi, dy/dxi, dz/dxi) is:
106 : // j(xi) := a*xi + b, where
107 12 : Point a = this->point(0) + this->point(1) - 2 * this->point(2);
108 6 : Point b = Real(.5) * (this->point(1) - this->point(0));
109 :
110 : // Now we solve for the point xi_m where j(xi_m) \cdot j(0) = 0.
111 : // If this occurs somewhere on the reference element, then the
112 : // element is not invertible.
113 : // j(xi_m) . j(0) = 0
114 : // <=>
115 : // (a*xi_m + b) . b = 0
116 : // <=>
117 : // (a.b)*xi_m + b.b = 0
118 : // <=>
119 : // xi_m = -(b.b) / (a.b)
120 :
121 : // 1.) If b.b==0, then the endpoints of the Edge3 are at the same
122 : // location, and the map is therefore not invertible.
123 6 : Real b_norm2 = b.norm_sq();
124 213 : if (b_norm2 <= tol*tol)
125 0 : return false;
126 :
127 : // 2.) If a.b==0, but b != 0 (see above), then the element is
128 : // invertible but we don't want to divide by zero in the
129 : // formula, so simply return true.
130 6 : Real ab = a * b;
131 213 : if (std::abs(ab) <= tol*tol)
132 2 : return true;
133 :
134 142 : Real xi_m = -b_norm2 / ab;
135 142 : return (xi_m < -1.) || (xi_m > 1.);
136 : }
137 :
138 :
139 :
140 33947428 : Order Edge3::default_order() const
141 : {
142 33947428 : return SECOND;
143 : }
144 :
145 :
146 :
147 1728 : void Edge3::connectivity(const unsigned int sc,
148 : const IOPackage iop,
149 : std::vector<dof_id_type> & conn) const
150 : {
151 144 : libmesh_assert_less_equal (sc, 1);
152 144 : libmesh_assert_less (sc, this->n_sub_elem());
153 144 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
154 :
155 : // Create storage
156 1728 : conn.resize(2);
157 :
158 1728 : switch (iop)
159 : {
160 1728 : case TECPLOT:
161 : {
162 1728 : switch (sc)
163 : {
164 864 : case 0:
165 936 : conn[0] = this->node_id(0)+1;
166 864 : conn[1] = this->node_id(2)+1;
167 864 : return;
168 :
169 864 : case 1:
170 936 : conn[0] = this->node_id(2)+1;
171 864 : conn[1] = this->node_id(1)+1;
172 864 : return;
173 :
174 0 : default:
175 0 : libmesh_error_msg("Invalid sc = " << sc);
176 : }
177 : }
178 :
179 :
180 0 : case VTK:
181 : {
182 0 : conn.resize(3);
183 0 : conn[0] = this->node_id(0);
184 0 : conn[1] = this->node_id(1);
185 0 : conn[2] = this->node_id(2);
186 0 : return;
187 :
188 : /*
189 : switch (sc)
190 : {
191 : case 0:
192 : conn[0] = this->node_id(0);
193 : conn[1] = this->node_id(2);
194 :
195 : return;
196 :
197 : case 1:
198 : conn[0] = this->node_id(2);
199 : conn[1] = this->node_id(1);
200 :
201 : return;
202 :
203 : default:
204 : libmesh_error_msg("Invalid sc = " << sc);
205 : }
206 : */
207 : }
208 :
209 0 : default:
210 0 : libmesh_error_msg("Unsupported IO package " << iop);
211 : }
212 : }
213 :
214 :
215 :
216 : std::pair<unsigned short int, unsigned short int>
217 0 : Edge3::second_order_child_vertex (const unsigned int) const
218 : {
219 0 : return std::pair<unsigned short int, unsigned short int>(0,0);
220 : }
221 :
222 :
223 :
224 281 : Real Edge3::volume () const
225 : {
226 : // Finding the (exact) length of a general quadratic element
227 : // is a surprisingly complicated formula.
228 52 : Point A = this->point(0) + this->point(1) - 2*this->point(2);
229 28 : Point B = (this->point(1) - this->point(0))/2;
230 :
231 28 : const Real a = A.norm_sq();
232 28 : const Real c = B.norm_sq();
233 :
234 : // Degenerate straight line case
235 281 : if (a == 0.)
236 221 : return 2. * std::sqrt(c);
237 :
238 60 : const Real b = -2.*std::abs(A*B);
239 :
240 60 : const Real sqrt_term1 = a - b + c;
241 60 : const Real sqrt_term2 = a + b + c;
242 :
243 : // Fall back on straight line case instead of computing nan
244 : // Note: b can be positive or negative so we have to check both cases.
245 60 : if (sqrt_term1 < 0. || sqrt_term2 < 0.)
246 0 : return 2. * std::sqrt(c);
247 :
248 60 : const Real r1 = std::sqrt(sqrt_term1);
249 60 : const Real r2 = std::sqrt(sqrt_term2);
250 60 : const Real rsum = r1 + r2;
251 60 : const Real root_a = std::sqrt(a);
252 60 : const Real b_over_root_a = b / root_a;
253 :
254 : // Pre-compute the denominator of the log term. If it's zero, fall
255 : // back on the straight line case.
256 60 : const Real log_term_denom = -root_a - 0.5*b_over_root_a + r2;
257 60 : if (log_term_denom == 0. || rsum == 0.)
258 0 : return 2. * std::sqrt(c);
259 :
260 60 : Real log1p_arg = 2*(root_a - b/rsum) / log_term_denom;
261 :
262 : return
263 76 : 0.5*(rsum + b_over_root_a*b_over_root_a/rsum +
264 60 : (c-0.25*b_over_root_a*b_over_root_a)*std::log1p(log1p_arg)/root_a);
265 : }
266 :
267 :
268 :
269 132892 : BoundingBox Edge3::loose_bounding_box () const
270 : {
271 : // This might be a curved line through 2-space or 3-space, in which
272 : // case the full bounding box can be larger than the bounding box of
273 : // just the nodes.
274 5077 : Point pmin, pmax;
275 :
276 531568 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
277 : {
278 413913 : Real center = this->point(2)(d);
279 413913 : Real hd = std::max(std::abs(center - this->point(0)(d)),
280 429144 : std::abs(center - this->point(1)(d)));
281 :
282 398676 : pmin(d) = center - hd;
283 398676 : pmax(d) = center + hd;
284 : }
285 :
286 137969 : return BoundingBox(pmin, pmax);
287 : }
288 :
289 :
290 :
291 1776 : dof_id_type Edge3::key () const
292 : {
293 1924 : return this->compute_key(this->node_id(2));
294 : }
295 :
296 :
297 622 : void Edge3::flip(BoundaryInfo * boundary_info)
298 : {
299 56 : libmesh_assert(boundary_info);
300 :
301 622 : swap2nodes(0,1);
302 56 : swap2neighbors(0,1);
303 622 : swap2boundarysides(0,1,boundary_info);
304 622 : }
305 :
306 :
307 : } // namespace libMesh
|
