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/quadrature_nodal.h"
22 :
23 : #include "libmesh/enum_to_string.h"
24 : #include "libmesh/face_c0polygon.h"
25 : #include "libmesh/quadrature_trap.h"
26 : #include "libmesh/quadrature_simpson.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 1972 : void QNodal::init_2D()
32 : {
33 : #if LIBMESH_DIM > 1
34 :
35 1972 : switch (_type)
36 : {
37 :
38 142 : case QUAD4:
39 : case QUADSHELL4:
40 : case TRI3:
41 : case TRISHELL3:
42 : {
43 146 : QTrap rule(/*dim=*/2, /*ignored*/_order);
44 142 : rule.init(*this);
45 4 : _points.swap (rule.get_points());
46 4 : _weights.swap(rule.get_weights());
47 4 : return;
48 : }
49 :
50 36 : case QUAD8:
51 : case QUADSHELL8:
52 : {
53 : // A rule with 8 points which is exact for linears, and
54 : // naturally produces a lumped approximation to the mass
55 : // matrix. The quadrature points are numbered the same way as
56 : // the reference element nodes.
57 980 : _points =
58 : {
59 : Point(-1,-1), Point(+1,-1), Point(+1,+1), Point(-1,+1),
60 : Point(0.,-1), Point(+1,0.), Point(0.,+1), Point(-1,0.)
61 526 : };
62 :
63 : // The weights for the Quad8 nodal quadrature rule are
64 : // obtained from the following specific steps. Other
65 : // "serendipity" type rules are obtained similarly.
66 : //
67 : // 1.) Due to the symmetry of the bi-unit square domain, we
68 : // first note that there are only two "classes" of node in the
69 : // Quad8: vertices (with associated weight wv) and edges (with
70 : // associated weight we).
71 : //
72 : // 2.) In order for such a nodal quadrature rule to be exact
73 : // for constants, the weights must sum to the area of the
74 : // reference element, and therefore we must have:
75 : // 4*wv + 4*we = 4 --> wv + we = 1
76 : //
77 : // 3.) Due to symmetry (once again), such a rule is then
78 : // automatically exact for the linear polynomials "x" and "y",
79 : // regardless of the values of wv and we.
80 : //
81 : // 4.) We therefore still have two unknowns and one equation.
82 : // Attempting to additionally make the nodal quadrature rule
83 : // exact for the quadratic polynomials "x^2" and "y^2" leads
84 : // to a rule with negative weights, namely:
85 : // wv = -1/3
86 : // we = 4/3
87 : // Note: these are the same values one obtains by integrating
88 : // the corresponding Quad8 Lagrange basis functions over the
89 : // reference element.
90 : //
91 : // Since the weights appear on the diagonal of the nodal
92 : // quadrature rule's approximation to the mass matrix, rules
93 : // with negative weights yield an indefinite mass matrix,
94 : // i.e. one with both positive and negative eigenvalues. Rules
95 : // with negative weights can also produce a negative integral
96 : // approximation to a strictly positive integrand, which may
97 : // be unacceptable in some situations.
98 : //
99 : // 5.) Requiring all positive weights therefore restricts the
100 : // nodal quadrature rule to only be exact for linears. But
101 : // there is still one free parameter remaining, and thus we
102 : // need some other criterion in order to complete the
103 : // description of the rule.
104 : //
105 : // Here, we have decided to choose the quadrature weights in
106 : // such a way that the resulting nodal quadrature mass matrix
107 : // approximation is "proportional" to the true mass matrix's
108 : // diagonal entries for the reference element. We therefore
109 : // pose the following constrained optimization problem:
110 : //
111 : // { min_{wv, we} |diag(M) - C*diag(W)|^2
112 : // (O) {
113 : // { subject to wv + we = 1
114 : //
115 : // where:
116 : // * M is the true mass matrix
117 : // * W is the nodal quadrature approximation to the mass
118 : // matrix. In this particular case:
119 : // diag(W) = [wv,wv,wv,wv,we,we,we,we]
120 : // * C = tr(M) / vol(E) is the ratio between the trace of the
121 : // true mass matrix and the volume of the reference
122 : // element. For all Lagrange finite elements, we have C<1.
123 : //
124 : // 6.) The optimization problem (O) is solved directly by
125 : // substituting the algebraic constraint into the functional
126 : // which is to be minimized, then setting the derivative with
127 : // respect to the remaining parameter equal to zero and
128 : // solving. In the Quad8 and Hex20 cases, there is only one
129 : // free parameter, while in the Prism15 case there are two
130 : // free parameters, so a 2x2 system of linear equations must
131 : // be solved.
132 36 : Real wv = Real(12) / 79;
133 36 : Real we = Real(67) / 79;
134 :
135 526 : _weights = {wv, wv, wv, wv, we, we, we, we};
136 :
137 526 : return;
138 : }
139 :
140 885 : case QUAD9:
141 : case QUADSHELL9:
142 : case TRI6:
143 : {
144 947 : QSimpson rule(/*dim=*/2, /*ignored*/_order);
145 885 : rule.init(*this);
146 62 : _points.swap (rule.get_points());
147 62 : _weights.swap(rule.get_weights());
148 62 : return;
149 : }
150 :
151 26 : case TRI7:
152 : {
153 : // We can't exactly represent cubics with only seven nodes,
154 : // but with w_i = integral(phi_i) for Lagrange shape functions
155 : // phi_i, we not only get exact integrals of every Lagrange
156 : // shape function, including the cubic bubble, we also get
157 : // exact integrals of the rest of P^3 too.
158 666 : _points =
159 : {
160 : Point(0.,0.), Point(+1,0.), Point(0.,+1), Point(.5,0.),
161 : Point(.5,.5), Point(0.,.5), Point(1/Real(3),1/Real(3))
162 359 : };
163 :
164 26 : Real wv = Real(1)/40;
165 26 : Real we = Real(1)/15;
166 359 : _weights = {wv, wv, wv, we, we, we, Real(9)/40};
167 359 : return;
168 : }
169 :
170 : //---------------------------------------------
171 : // Arbitrary polygon quadrature rules
172 60 : case C0POLYGON:
173 : {
174 : // C0Polygon requires the newer Quadrature API
175 60 : if (!_elem)
176 0 : libmesh_error();
177 :
178 5 : libmesh_assert(_elem->type() == C0POLYGON);
179 :
180 5 : const C0Polygon & poly = *cast_ptr<const C0Polygon *>(_elem);
181 :
182 5 : const unsigned int ns = poly.n_sides();
183 60 : const Real pi_over_ns = libMesh::pi / ns;
184 60 : const Real master_poly_area = ns * 0.25 / tan(pi_over_ns);
185 :
186 60 : const unsigned int nn = poly.n_nodes();
187 60 : const Real weight = master_poly_area / nn;
188 60 : _weights.resize(nn, weight);
189 :
190 60 : _points.resize(poly.n_nodes());
191 360 : for (auto n : make_range(nn))
192 300 : _points[n] = poly.master_point(n);
193 :
194 5 : return;
195 : }
196 :
197 0 : default:
198 0 : libmesh_error_msg("Element type not supported!:" << Utility::enum_to_string(_type));
199 : }
200 : #endif
201 : }
202 :
203 : } // namespace libMesh
|
