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_edge4.h"
22 : #include "libmesh/enum_io_package.h"
23 : #include "libmesh/enum_order.h"
24 :
25 : namespace libMesh
26 : {
27 :
28 : // ------------------------------------------------------------
29 : // Edge4 class static member initializations
30 : const int Edge4::num_nodes;
31 :
32 : #ifdef LIBMESH_ENABLE_AMR
33 :
34 : const Real Edge4::_embedding_matrix[Edge4::num_children][Edge4::num_nodes][Edge4::num_nodes] =
35 : {
36 : // embedding matrix for child 0
37 :
38 : {
39 : // 0 1 2 3 // Shape function index
40 : {1.0, 0.0, 0.0, 0.0}, // left, xi = -1
41 : {-0.0625, -0.0625, 0.5625, 0.5625}, // right, xi = 0
42 : {0.3125, 0.0625, 0.9375, -0.3125}, // middle left, xi = -2/3
43 : {0.0, 0.0, 1.0, 0.0} // middle right, xi = -1/3
44 : },
45 :
46 : // embedding matrix for child 1
47 : {
48 : // 0 1 2 3 // Shape function index
49 : {-0.0625, -0.0625, 0.5625, 0.5625}, // left, xi = 0
50 : {0.0, 1.0, 0.0, 0.0}, // right, xi = 1
51 : {0.0, 0.0, 0.0, 1.0}, // middle left, xi = 1/3
52 : {0.0625, 0.3125, -0.3125, 0.9375} // middle right, xi = 2/3
53 : }
54 : };
55 :
56 : #endif
57 :
58 96 : bool Edge4::is_vertex(const unsigned int i) const
59 : {
60 96 : return (i==0) || (i==1);
61 : }
62 :
63 0 : bool Edge4::is_edge(const unsigned int i) const
64 : {
65 0 : return (i==2) || (i==3);
66 : }
67 :
68 0 : bool Edge4::is_face(const unsigned int ) const
69 : {
70 0 : return false;
71 : }
72 :
73 616 : bool Edge4::is_node_on_side(const unsigned int n,
74 : const unsigned int s) const
75 : {
76 20 : libmesh_assert_less (s, 2);
77 20 : libmesh_assert_less (n, Edge4::num_nodes);
78 616 : return (s == n);
79 : }
80 :
81 0 : bool Edge4::is_node_on_edge(const unsigned int,
82 : const unsigned int libmesh_dbg_var(e)) const
83 : {
84 0 : libmesh_assert_equal_to (e, 0);
85 0 : return true;
86 : }
87 :
88 :
89 :
90 1018 : bool Edge4::has_affine_map() const
91 : {
92 80 : Point v = this->point(1) - this->point(0);
93 1058 : if (!v.relative_fuzzy_equals
94 1018 : ((this->point(2) - this->point(0))*3, affine_tol))
95 0 : return false;
96 1058 : if (!v.relative_fuzzy_equals
97 1018 : ((this->point(3) - this->point(0))*1.5, affine_tol))
98 0 : return false;
99 40 : return true;
100 : }
101 :
102 :
103 :
104 355 : bool Edge4::has_invertible_map(Real tol) const
105 : {
106 : // At the moment this only makes sense for Lagrange elements
107 10 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
108 :
109 : // dx/dxi = a*xi^2 + b*xi + c,
110 : // where a, b, and c are vector quantities that depend on the
111 : // nodal positions as follows:
112 20 : const Point & x0 = this->point(0);
113 10 : const Point & x1 = this->point(1);
114 10 : const Point & x2 = this->point(2);
115 10 : const Point & x3 = this->point(3);
116 :
117 10 : Point a = Real(27)/16 * (x1 - x0 + 3*(x2 - x3));
118 10 : Point b = Real(18)/16 * (x0 + x1 - x2 - x3);
119 10 : Point c = Real(1)/16 * (x0 - x1 + 27*(x3 - x2));
120 :
121 : // Normalized midpoint value of Jacobian. If c==0, then dx/dxi(0) =
122 : // 0 and the Jacobian is either zero on the whole element, or
123 : // changes sign within the element. Either way, the element is not
124 : // invertible. Note: use <= tol, so that if someone passes 0 for tol
125 : // it still works.
126 345 : Real c_norm = c.norm();
127 355 : if (c_norm <= tol)
128 0 : return false;
129 :
130 : // Coefficients of the quadratic scalar function
131 : // j(xi) := dot(c.unit(), dx/dxi)
132 : // = alpha*xi^2 + beta*xi + gamma
133 355 : Real alpha = (a * c) / c_norm;
134 355 : Real beta = (b * c) / c_norm;
135 10 : Real gamma = c_norm;
136 :
137 : // If alpha and beta are both (approximately) zero, then the
138 : // Jacobian is actually constant but it's not zero (as we already
139 : // checked) so the element is invertible!
140 355 : if ((std::abs(alpha) <= tol) && (std::abs(beta) <= tol))
141 2 : return true;
142 :
143 : // If alpha is approximately zero, but beta is not, then j(xi) is
144 : // actually linear and the Jacobian changes sign at the point xi =
145 : // -gamma/beta, so we need to check whether that's in the
146 : // element.
147 284 : if (std::abs(alpha) <= tol)
148 : {
149 : // Debugging:
150 : // libMesh::out << "alpha = " << alpha << ", std::abs(alpha) <= tol" << std::endl;
151 :
152 0 : Real xi_0 = -gamma / beta;
153 0 : return ((xi_0 < -1.) || (xi_0 > 1.));
154 : }
155 :
156 : // If alpha is not (approximately) zero, then j(xi) is quadratic
157 : // and we need to solve for the roots.
158 284 : Real sqrt_term = beta*beta - 4*alpha*gamma;
159 :
160 : // Debugging:
161 : // libMesh::out << "sqrt_term = " << sqrt_term << std::endl;
162 :
163 : // If the term under the square root is negative, then the roots are
164 : // imaginary and the element is invertible.
165 284 : if (sqrt_term < 0)
166 2 : return true;
167 :
168 213 : sqrt_term = std::sqrt(sqrt_term);
169 : Real
170 213 : xi_1 = 0.5 * (-beta + sqrt_term) / alpha,
171 213 : xi_2 = 0.5 * (-beta - sqrt_term) / alpha;
172 :
173 : // libMesh::out << "xi_1 = " << xi_1 << std::endl;
174 : // libMesh::out << "xi_2 = " << xi_2 << std::endl;
175 :
176 : // If a root is outside the reference element, it is "OK".
177 : bool
178 213 : xi_1_ok = ((xi_1 < -1.) || (xi_1 > 1.)),
179 213 : xi_2_ok = ((xi_2 < -1.) || (xi_2 > 1.));
180 :
181 : // If both roots are OK, the element is invertible.
182 213 : return xi_1_ok && xi_2_ok;
183 : }
184 :
185 :
186 :
187 25893 : Order Edge4::default_order() const
188 : {
189 25893 : return THIRD;
190 : }
191 :
192 :
193 :
194 0 : void Edge4::connectivity(const unsigned int sc,
195 : const IOPackage iop,
196 : std::vector<dof_id_type> & conn) const
197 : {
198 0 : libmesh_assert_less_equal (sc, 2);
199 0 : libmesh_assert_less (sc, this->n_sub_elem());
200 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
201 :
202 : // Create storage
203 0 : conn.resize(2);
204 :
205 0 : switch(iop)
206 : {
207 0 : case TECPLOT:
208 : {
209 0 : switch (sc)
210 : {
211 0 : case 0:
212 0 : conn[0] = this->node_id(0)+1;
213 0 : conn[1] = this->node_id(2)+1;
214 0 : return;
215 :
216 0 : case 1:
217 0 : conn[0] = this->node_id(2)+1;
218 0 : conn[1] = this->node_id(3)+1;
219 0 : return;
220 :
221 0 : case 2:
222 0 : conn[0] = this->node_id(3)+1;
223 0 : conn[1] = this->node_id(1)+1;
224 0 : return;
225 :
226 0 : default:
227 0 : libmesh_error_msg("Invalid sc = " << sc);
228 : }
229 :
230 : }
231 :
232 0 : case VTK:
233 : {
234 :
235 0 : switch (sc)
236 : {
237 0 : case 0:
238 0 : conn[0] = this->node_id(0);
239 0 : conn[1] = this->node_id(2);
240 0 : return;
241 :
242 0 : case 1:
243 0 : conn[0] = this->node_id(2);
244 0 : conn[1] = this->node_id(3);
245 0 : return;
246 :
247 0 : case 2:
248 0 : conn[0] = this->node_id(3);
249 0 : conn[1] = this->node_id(1);
250 0 : return;
251 :
252 0 : default:
253 0 : libmesh_error_msg("Invalid sc = " << sc);
254 : }
255 : }
256 :
257 0 : default:
258 0 : libmesh_error_msg("Unsupported IO package " << iop);
259 : }
260 : }
261 :
262 :
263 :
264 8147 : BoundingBox Edge4::loose_bounding_box () const
265 : {
266 : // This might be a curved line through 2-space or 3-space, in which
267 : // case the full bounding box can be larger than the bounding box of
268 : // just the nodes.
269 : //
270 : // FIXME - I haven't yet proven the formula below to be correct for
271 : // cubics - RHS
272 299 : Point pmin, pmax;
273 :
274 32588 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
275 : {
276 26235 : Real center = (this->point(2)(d) + this->point(3)(d))/2;
277 25338 : Real hd = std::max(std::abs(center - this->point(0)(d)),
278 26235 : std::abs(center - this->point(1)(d)));
279 :
280 24441 : pmin(d) = center - hd;
281 24441 : pmax(d) = center + hd;
282 : }
283 :
284 8446 : return BoundingBox(pmin, pmax);
285 : }
286 :
287 :
288 :
289 0 : dof_id_type Edge4::key () const
290 : {
291 0 : return this->compute_key(this->node_id(0),
292 : this->node_id(1),
293 : this->node_id(2),
294 0 : this->node_id(3));
295 : }
296 :
297 :
298 :
299 0 : Real Edge4::volume () const
300 : {
301 : // Make copies of our points. It makes the subsequent calculations a bit
302 : // shorter and avoids dereferencing the same pointer multiple times.
303 : Point
304 0 : x0 = point(0),
305 0 : x1 = point(1),
306 0 : x2 = point(2),
307 0 : x3 = point(3);
308 :
309 : // Construct constant data vectors.
310 : // \vec{x}_{\xi} = \vec{a1}*xi**2 + \vec{b1}*xi + \vec{c1}
311 : // This is copy-pasted directly from the output of a Python script.
312 : Point
313 0 : a1 = -27*x0/16 + 27*x1/16 + 81*x2/16 - 81*x3/16,
314 0 : b1 = 9*x0/8 + 9*x1/8 - 9*x2/8 - 9*x3/8,
315 0 : c1 = x0/16 - x1/16 - 27*x2/16 + 27*x3/16;
316 :
317 : // 4 point quadrature, 7th-order accurate
318 0 : const unsigned int N = 4;
319 0 : const Real q[N] = {-std::sqrt(525 + 70*std::sqrt(30.)) / 35,
320 : -std::sqrt(525 - 70*std::sqrt(30.)) / 35,
321 : std::sqrt(525 - 70*std::sqrt(30.)) / 35,
322 : std::sqrt(525 + 70*std::sqrt(30.)) / 35};
323 0 : const Real w[N] = {(18 - std::sqrt(30.)) / 36,
324 : (18 + std::sqrt(30.)) / 36,
325 : (18 + std::sqrt(30.)) / 36,
326 : (18 - std::sqrt(30.)) / 36};
327 :
328 0 : Real vol=0.;
329 0 : for (unsigned int i=0; i<N; ++i)
330 0 : vol += w[i] * (q[i]*q[i]*a1 + q[i]*b1 + c1).norm();
331 :
332 0 : return vol;
333 : }
334 :
335 :
336 72 : void Edge4::flip(BoundaryInfo * boundary_info)
337 : {
338 6 : libmesh_assert(boundary_info);
339 :
340 72 : swap2nodes(0,1);
341 72 : swap2nodes(2,3);
342 6 : swap2neighbors(0,1);
343 72 : swap2boundarysides(0,1,boundary_info);
344 72 : }
345 :
346 :
347 : } // namespace libMesh
|
