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