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_quad4.h"
21 : #include "libmesh/enum_io_package.h"
22 : #include "libmesh/enum_order.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 :
28 :
29 :
30 : // ------------------------------------------------------------
31 : // Quad4 class static member initialization
32 : const int Quad4::num_nodes;
33 : const int Quad4::nodes_per_side;
34 :
35 : const unsigned int Quad4::side_nodes_map[Quad4::num_sides][Quad4::nodes_per_side] =
36 : {
37 : {0, 1}, // Side 0
38 : {1, 2}, // Side 1
39 : {2, 3}, // Side 2
40 : {3, 0} // Side 3
41 : };
42 :
43 : #ifdef LIBMESH_ENABLE_AMR
44 :
45 : const Real Quad4::_embedding_matrix[Quad4::num_children][Quad4::num_nodes][Quad4::num_nodes] =
46 : {
47 : // embedding matrix for child 0
48 : {
49 : // 0 1 2 3
50 : {1.0, 0.0, 0.0, 0.0}, // 0
51 : {0.5, 0.5, 0.0, 0.0}, // 1
52 : {.25, .25, .25, .25}, // 2
53 : {0.5, 0.0, 0.0, 0.5} // 3
54 : },
55 :
56 : // embedding matrix for child 1
57 : {
58 : // 0 1 2 3
59 : {0.5, 0.5, 0.0, 0.0}, // 0
60 : {0.0, 1.0, 0.0, 0.0}, // 1
61 : {0.0, 0.5, 0.5, 0.0}, // 2
62 : {.25, .25, .25, .25} // 3
63 : },
64 :
65 : // embedding matrix for child 2
66 : {
67 : // 0 1 2 3
68 : {0.5, 0.0, 0.0, 0.5}, // 0
69 : {.25, .25, .25, .25}, // 1
70 : {0.0, 0.0, 0.5, 0.5}, // 2
71 : {0.0, 0.0, 0.0, 1.0} // 3
72 : },
73 :
74 : // embedding matrix for child 3
75 : {
76 : // 0 1 2 3
77 : {.25, .25, .25, .25}, // 0
78 : {0.0, 0.5, 0.5, 0.0}, // 1
79 : {0.0, 0.0, 1.0, 0.0}, // 2
80 : {0.0, 0.0, 0.5, 0.5} // 3
81 : }
82 : };
83 :
84 : #endif
85 :
86 :
87 :
88 :
89 :
90 : // ------------------------------------------------------------
91 : // Quad4 class member functions
92 :
93 23345172 : bool Quad4::is_vertex(const unsigned int libmesh_dbg_var(n)) const
94 : {
95 1861096 : libmesh_assert_not_equal_to (n, invalid_uint);
96 23345172 : return true;
97 : }
98 :
99 0 : bool Quad4::is_edge(const unsigned int) const
100 : {
101 0 : return false;
102 : }
103 :
104 0 : bool Quad4::is_face(const unsigned int) const
105 : {
106 0 : return false;
107 : }
108 :
109 347164 : bool Quad4::is_node_on_side(const unsigned int n,
110 : const unsigned int s) const
111 : {
112 54040 : libmesh_assert_less (s, n_sides());
113 54040 : return std::find(std::begin(side_nodes_map[s]),
114 347164 : std::end(side_nodes_map[s]),
115 347164 : n) != std::end(side_nodes_map[s]);
116 : }
117 :
118 : std::vector<unsigned>
119 112778 : Quad4::nodes_on_side(const unsigned int s) const
120 : {
121 3796 : libmesh_assert_less(s, n_sides());
122 116572 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
123 : }
124 :
125 : std::vector<unsigned>
126 1820 : Quad4::nodes_on_edge(const unsigned int e) const
127 : {
128 1820 : return nodes_on_side(e);
129 : }
130 :
131 54897578 : bool Quad4::has_affine_map() const
132 : {
133 10525741 : Point v = this->point(3) - this->point(0);
134 54897578 : return (v.relative_fuzzy_equals(this->point(2) - this->point(1), affine_tol));
135 : }
136 :
137 :
138 :
139 881 : bool Quad4::has_invertible_map(Real tol) const
140 : {
141 : // At the moment this only makes sense for Lagrange elements
142 46 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
143 :
144 : // Side vectors
145 92 : Point s0 = this->point(1) - this->point(0);
146 46 : Point s1 = this->point(2) - this->point(1);
147 46 : Point s2 = this->point(3) - this->point(2);
148 46 : Point s3 = this->point(0) - this->point(3);
149 :
150 : // Cross products of side vectors
151 835 : Point v1 = s3.cross(s0);
152 835 : Point v2 = s2.cross(s0);
153 835 : Point v3 = s3.cross(s1);
154 :
155 : // A unit vector in the direction of:
156 : // f(xi, eta) = (v1 + xi*v2 + eta*v3)
157 : // at the midpoint (xi, eta) = (1/2, 1/2) of the element. (Note that
158 : // we are using the [0,1]^2 reference element definition for the
159 : // Quad4 instead of the [-1,1]^2 reference element that is typically
160 : // used for FEM calculations.) We use this as a "reference" vector
161 : // and compare the sign of dot(n,f) at each vertex.
162 46 : Point n = v1 + Real(.5) * (v2 + v3);
163 835 : Real norm_n = n.norm();
164 :
165 : // If the Jacobian vector at the midpoint of the element is zero,
166 : // then it must be either zero for the entire element, or change
167 : // sign at the vertices. Either way the element is non-invertible.
168 881 : if (norm_n <= tol)
169 2 : return false;
170 :
171 44 : n /= norm_n;
172 :
173 : // Debugging
174 : // std::cout << "n=" << n << std::endl;
175 :
176 : // Compute scalar quantity n * (v1 + xi*v2 + eta*v3) at each
177 : // vertex. If it is non-zero and has the same sign at each
178 : // vertex, the the element is invertible, otherwise it is not.
179 : std::array<Real, 4> vertex_vals;
180 44 : unsigned int ctr = 0;
181 2430 : for (unsigned int i=0; i<2; ++i)
182 4860 : for (unsigned int j=0; j<2; ++j)
183 3592 : vertex_vals[ctr++] = n * (v1 + Real(i)*v2 + Real(j)*v3);
184 :
185 : // Debugging:
186 : // std::cout << "Vertex values: ";
187 : // for (const auto & val : vertex_vals)
188 : // std::cout << val << " ";
189 : // std::cout << std::endl;
190 :
191 44 : auto result = std::minmax_element(vertex_vals.begin(), vertex_vals.end());
192 810 : Real min_vertex = *(result.first);
193 810 : Real max_vertex = *(result.second);
194 :
195 : // Debugging
196 : // std::cout << "min_vertex=" << min_vertex << std::endl;
197 : // std::cout << "max_vertex=" << max_vertex << std::endl;
198 :
199 : // If max and min are both on the same side of 0, we are invertible, otherwise we are not.
200 810 : return ((max_vertex > 0 && min_vertex > 0) || (max_vertex < 0 && min_vertex < 0));
201 : }
202 :
203 :
204 :
205 67356787 : Order Quad4::default_order() const
206 : {
207 67356787 : return FIRST;
208 : }
209 :
210 :
211 :
212 18424170 : std::unique_ptr<Elem> Quad4::build_side_ptr (const unsigned int i)
213 : {
214 18424170 : return this->simple_build_side_ptr<Edge2, Quad4>(i);
215 : }
216 :
217 :
218 :
219 613917 : void Quad4::build_side_ptr (std::unique_ptr<Elem> & side,
220 : const unsigned int i)
221 : {
222 613917 : this->simple_build_side_ptr<Quad4>(side, i, EDGE2);
223 613917 : }
224 :
225 :
226 :
227 3568928 : void Quad4::connectivity(const unsigned int libmesh_dbg_var(sf),
228 : const IOPackage iop,
229 : std::vector<dof_id_type> & conn) const
230 : {
231 327606 : libmesh_assert_less (sf, this->n_sub_elem());
232 327606 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
233 :
234 : // Create storage.
235 3568928 : conn.resize(4);
236 :
237 3568928 : switch (iop)
238 : {
239 3568928 : case TECPLOT:
240 : {
241 3896534 : conn[0] = this->node_id(0)+1;
242 3568928 : conn[1] = this->node_id(1)+1;
243 3568928 : conn[2] = this->node_id(2)+1;
244 3568928 : conn[3] = this->node_id(3)+1;
245 3568928 : return;
246 : }
247 :
248 0 : case VTK:
249 : {
250 0 : conn[0] = this->node_id(0);
251 0 : conn[1] = this->node_id(1);
252 0 : conn[2] = this->node_id(2);
253 0 : conn[3] = this->node_id(3);
254 0 : return;
255 : }
256 :
257 0 : default:
258 0 : libmesh_error_msg("Unsupported IO package " << iop);
259 : }
260 : }
261 :
262 :
263 :
264 1190 : Point Quad4::true_centroid () const
265 : {
266 : // Convenient references to our points
267 : const Point
268 120 : &x0 = point(0), &x1 = point(1),
269 60 : &x2 = point(2), &x3 = point(3);
270 :
271 : // Construct "dx/d(xi)" and "dx/d(eta)" vectors which are columns of the Jacobian.
272 : // \vec{x}_{\xi} = \vec{a1}*eta + \vec{b1}
273 : // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}
274 : // This is copy-pasted directly from the output of a Python script. Note: we are off
275 : // by a factor of (1/4) here, but since the final result should have the same error
276 : // in the numerator and denominator, it should cancel out while saving us some math
277 : // operations.
278 : Point
279 60 : a1 = x0 - x1 + x2 - x3,
280 60 : b1 = -x0 + x1 + x2 - x3,
281 60 : a2 = a1,
282 60 : b2 = -x0 - x1 + x2 + x3;
283 :
284 : // Use 2x2 quadrature to compute the integral of each basis function
285 : // (as defined on the [-1,1]^2 reference domain). We use a 4-point
286 : // rule, which is exact for bi-cubics. The weights for this rule are
287 : // all equal to 1.
288 1190 : const Real q[2] = {-std::sqrt(Real(3))/3, std::sqrt(Real(3))/3.};
289 :
290 : // Nodal areas
291 60 : Real A0 = 0., A1 = 0., A2 = 0., A3 = 0.;
292 :
293 3570 : for (const auto & xi : q)
294 7140 : for (const auto & eta : q)
295 : {
296 5240 : Real jxw = cross_norm(eta*a1 + b1, xi*a2 + b2);
297 :
298 4760 : A0 += jxw * (1-xi) * (1-eta); // 4 * phi_0
299 4760 : A1 += jxw * (1+xi) * (1-eta); // 4 * phi_1
300 4760 : A2 += jxw * (1+xi) * (1+eta); // 4 * phi_2
301 4760 : A3 += jxw * (1-xi) * (1+eta); // 4 * phi_3
302 : }
303 :
304 : // Compute centroid
305 1250 : return Point(x0(0)*A0 + x1(0)*A1 + x2(0)*A2 + x3(0)*A3,
306 1190 : x0(1)*A0 + x1(1)*A1 + x2(1)*A2 + x3(1)*A3,
307 1310 : x0(2)*A0 + x1(2)*A1 + x2(2)*A2 + x3(2)*A3) / (A0 + A1 + A2 + A3);
308 : }
309 :
310 :
311 :
312 379213 : Real Quad4::volume () const
313 : {
314 : // Make copies of our points. It makes the subsequent calculations a bit
315 : // shorter and avoids dereferencing the same pointer multiple times.
316 : Point
317 443379 : x0 = point(0), x1 = point(1),
318 411296 : x2 = point(2), x3 = point(3);
319 :
320 : // Construct constant data vectors.
321 : // \vec{x}_{\xi} = \vec{a1}*eta + \vec{b1}
322 : // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}
323 : // This is copy-pasted directly from the output of a Python script.
324 : Point
325 32087 : a1 = x0/4 - x1/4 + x2/4 - x3/4,
326 32087 : b1 = -x0/4 + x1/4 + x2/4 - x3/4,
327 32087 : a2 = a1,
328 32087 : b2 = -x0/4 - x1/4 + x2/4 + x3/4;
329 :
330 : // Check for quick return for parallelogram QUAD4.
331 379213 : if (a1.relative_fuzzy_equals(Point(0,0,0)))
332 369454 : return 4. * b1.cross(b2).norm();
333 :
334 : // Otherwise, use 2x2 quadrature to approximate the surface area.
335 :
336 : // 4-point rule, exact for bi-cubics. The weights for this rule are
337 : // all equal to 1.
338 9759 : const Real q[2] = {-std::sqrt(Real(3))/3, std::sqrt(Real(3))/3.};
339 :
340 778 : Real vol=0.;
341 29277 : for (unsigned int i=0; i<2; ++i)
342 58554 : for (unsigned int j=0; j<2; ++j)
343 42148 : vol += cross_norm(q[j]*a1 + b1,
344 42148 : q[i]*a2 + b2);
345 :
346 778 : return vol;
347 : }
348 :
349 : BoundingBox
350 8218539 : Quad4::loose_bounding_box () const
351 : {
352 8218539 : return Elem::loose_bounding_box();
353 : }
354 :
355 :
356 1734 : void Quad4::permute(unsigned int perm_num)
357 : {
358 212 : libmesh_assert_less (perm_num, 4);
359 :
360 6360 : for (unsigned int i = 0; i != perm_num; ++i)
361 : {
362 4626 : swap4nodes(0,1,2,3);
363 4038 : swap4neighbors(0,1,2,3);
364 : }
365 1734 : }
366 :
367 :
368 288 : void Quad4::flip(BoundaryInfo * boundary_info)
369 : {
370 24 : libmesh_assert(boundary_info);
371 :
372 288 : swap2nodes(0,1);
373 288 : swap2nodes(2,3);
374 24 : swap2neighbors(1,3);
375 288 : swap2boundarysides(1,3,boundary_info);
376 288 : swap2boundaryedges(1,3,boundary_info);
377 288 : }
378 :
379 :
380 210191 : ElemType Quad4::side_type (const unsigned int libmesh_dbg_var(s)) const
381 : {
382 18040 : libmesh_assert_less (s, 4);
383 210191 : return EDGE2;
384 : }
385 :
386 : } // namespace libMesh
|
