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 : // libMesh includes
20 : #include "libmesh/quadrature_conical.h"
21 : #include "libmesh/quadrature_gauss.h"
22 : #include "libmesh/quadrature_jacobi.h"
23 : #include "libmesh/enum_quadrature_type.h"
24 :
25 : namespace libMesh
26 : {
27 :
28 : // See also the files:
29 : // quadrature_conical_2D.C
30 : // quadrature_conical_3D.C
31 : // for additional implementation.
32 :
33 0 : QuadratureType QConical::type() const
34 : {
35 0 : return QCONICAL;
36 : }
37 :
38 0 : std::unique_ptr<QBase> QConical::clone() const
39 : {
40 0 : return std::make_unique<QConical>(*this);
41 : }
42 :
43 0 : void QConical::init_1D()
44 : {
45 0 : QGauss gauss1D(1, get_order());
46 0 : gauss1D.init(*this);
47 :
48 : // Swap points and weights with the about-to-be destroyed rule.
49 0 : _points.swap(gauss1D.get_points());
50 0 : _weights.swap(gauss1D.get_weights());
51 0 : }
52 :
53 :
54 :
55 : // Builds and scales a Gauss rule and a Jacobi rule.
56 : // Then combines them to compute points and weights
57 : // of a 2D conical product rule.
58 710 : void QConical::conical_product_tri()
59 : {
60 : // Be sure the underlying rule object was built with the same dimension as the
61 : // rule we are about to construct.
62 20 : libmesh_assert_equal_to (this->get_dim(), 2);
63 :
64 750 : QGauss gauss1D(1, get_order());
65 750 : QJacobi jac1D(1, get_order(), 1, 0);
66 :
67 : // The Gauss rule needs to be scaled to [0,1]
68 710 : std::pair<Real, Real> old_range(-1, 1);
69 710 : std::pair<Real, Real> new_range( 0, 1);
70 710 : gauss1D.scale(old_range,
71 : new_range);
72 :
73 : // Now construct the points and weights for the conical product rule.
74 :
75 : // Both rules should have the same number of points.
76 20 : libmesh_assert_equal_to (gauss1D.n_points(), jac1D.n_points());
77 :
78 : // Save the number of points as a convenient variable
79 20 : const unsigned int np = gauss1D.n_points();
80 :
81 : // Both rules should be between x=0 and x=1
82 20 : libmesh_assert_greater_equal (gauss1D.qp(0)(0), 0.0);
83 20 : libmesh_assert_less_equal (gauss1D.qp(np-1)(0), 1.0);
84 20 : libmesh_assert_greater_equal (jac1D.qp(0)(0), 0.0);
85 20 : libmesh_assert_less_equal (jac1D.qp(np-1)(0), 1.0);
86 :
87 : // Resize the points and weights vectors
88 710 : _points.resize(np * np);
89 710 : _weights.resize(np * np);
90 :
91 : // Compute the conical product
92 20 : unsigned int gp = 0;
93 2840 : for (unsigned int i=0; i<np; i++)
94 9940 : for (unsigned int j=0; j<np; j++)
95 : {
96 8030 : _points[gp](0) = jac1D.qp(j)(0); //s[j];
97 7810 : _points[gp](1) = gauss1D.qp(i)(0) * (1.-jac1D.qp(j)(0)); //r[i]*(1.-s[j]);
98 7810 : _weights[gp] = gauss1D.w(i) * jac1D.w(j); //A[i]*B[j];
99 7810 : gp++;
100 : }
101 710 : }
102 :
103 :
104 :
105 :
106 : // Builds and scales a Gauss rule and a Jacobi rule.
107 : // Then combines them to compute points and weights
108 : // of a 3D conical product rule for the Tet.
109 639 : void QConical::conical_product_tet()
110 : {
111 : // Be sure the underlying rule object was built with the same dimension as the
112 : // rule we are about to construct.
113 18 : libmesh_assert_equal_to (this->get_dim(), 3);
114 :
115 675 : QGauss gauss1D(1, get_order());
116 675 : QJacobi jacA1D(1, get_order(), /*alpha=*/1, /*beta=*/0);
117 675 : QJacobi jacB1D(1, get_order(), /*alpha=*/2, /*beta=*/0);
118 :
119 : // The Gauss rule needs to be scaled to [0,1]
120 639 : std::pair<Real, Real> old_range(-1, 1);
121 639 : std::pair<Real, Real> new_range( 0, 1);
122 639 : gauss1D.scale(old_range,
123 : new_range);
124 :
125 : // Now construct the points and weights for the conical product rule.
126 :
127 : // All rules should have the same number of points
128 18 : libmesh_assert_equal_to (gauss1D.n_points(), jacA1D.n_points());
129 18 : libmesh_assert_equal_to (jacA1D.n_points(), jacB1D.n_points());
130 :
131 : // Save the number of points as a convenient variable
132 18 : const unsigned int np = gauss1D.n_points();
133 :
134 : // All rules should be between x=0 and x=1
135 18 : libmesh_assert_greater_equal (gauss1D.qp(0)(0), 0.0);
136 18 : libmesh_assert_less_equal (gauss1D.qp(np-1)(0), 1.0);
137 18 : libmesh_assert_greater_equal (jacA1D.qp(0)(0), 0.0);
138 18 : libmesh_assert_less_equal (jacA1D.qp(np-1)(0), 1.0);
139 18 : libmesh_assert_greater_equal (jacB1D.qp(0)(0), 0.0);
140 18 : libmesh_assert_less_equal (jacB1D.qp(np-1)(0), 1.0);
141 :
142 : // Resize the points and weights vectors
143 639 : _points.resize(np * np * np);
144 639 : _weights.resize(np * np * np);
145 :
146 : // Compute the conical product
147 18 : unsigned int gp = 0;
148 2059 : for (unsigned int i=0; i<np; i++)
149 5112 : for (unsigned int j=0; j<np; j++)
150 14484 : for (unsigned int k=0; k<np; k++)
151 : {
152 11096 : _points[gp](0) = jacB1D.qp(k)(0); //t[k];
153 10792 : _points[gp](1) = jacA1D.qp(j)(0) * (1.-jacB1D.qp(k)(0)); //s[j]*(1.-t[k]);
154 11096 : _points[gp](2) = gauss1D.qp(i)(0) * (1.-jacA1D.qp(j)(0)) * (1.-jacB1D.qp(k)(0)); //r[i]*(1.-s[j])*(1.-t[k]);
155 11096 : _weights[gp] = gauss1D.w(i) * jacA1D.w(j) * jacB1D.w(k); //A[i]*B[j]*C[k];
156 10792 : gp++;
157 : }
158 639 : }
159 :
160 :
161 :
162 :
163 :
164 : // Builds and scales a Gauss rule and a Jacobi rule.
165 : // Then combines them to compute points and weights
166 : // of a 3D conical product rule for the Pyramid. The integral
167 : // over the reference Pyramid can be written (in LaTeX notation) as:
168 : //
169 : // If := \int_0^1 dz \int_{-(1-z)}^{(1-z)} dy \int_{-(1-z)}^{(1-z)} f(x,y,z) dx (1)
170 : //
171 : // (Imagine a stack of infinitely thin squares which decrease in size as
172 : // you approach the apex.) Under the transformation of variables:
173 : //
174 : // z=w
175 : // y=(1-z)v
176 : // x=(1-z)u,
177 : //
178 : // The Jacobian determinant of this transformation is |J|=(1-w)^2, and
179 : // the integral itself is transformed to:
180 : //
181 : // If = \int_0^1 (1-w)^2 dw \int_{-1}^{1} dv \int_{-1}^{1} f((1-w)u, (1-w)v, w) du (2)
182 : //
183 : // The integral can now be approximated by the product of three 1D quadrature rules:
184 : // A Jacobi rule with alpha==2, beta==0 in w, and Gauss rules in v and u. In this way
185 : // we can obtain 3D rules to any order for which the 1D rules exist.
186 18819 : void QConical::conical_product_pyramid()
187 : {
188 : // Be sure the underlying rule object was built with the same dimension as the
189 : // rule we are about to construct.
190 1018 : libmesh_assert_equal_to (this->get_dim(), 3);
191 :
192 20855 : QGauss gauss1D(1, get_order());
193 20855 : QJacobi jac1D(1, get_order(), 2, 0);
194 :
195 : // These rules should have the same number of points
196 1018 : libmesh_assert_equal_to (gauss1D.n_points(), jac1D.n_points());
197 :
198 : // Save the number of points as a convenient variable
199 1018 : const unsigned int np = gauss1D.n_points();
200 :
201 : // Resize the points and weights vectors
202 18819 : _points.resize(np * np * np);
203 18819 : _weights.resize(np * np * np);
204 :
205 : // Compute the conical product
206 1018 : unsigned int q = 0;
207 77730 : for (unsigned int i=0; i<np; ++i)
208 288996 : for (unsigned int j=0; j<np; ++j)
209 1407408 : for (unsigned int k=0; k<np; ++k, ++q)
210 : {
211 52266 : const Real xi=gauss1D.qp(i)(0);
212 52266 : const Real yj=gauss1D.qp(j)(0);
213 52266 : const Real zk=jac1D.qp(k)(0);
214 :
215 1177323 : _points[q](0) = (1.-zk) * xi;
216 1177323 : _points[q](1) = (1.-zk) * yj;
217 1177323 : _points[q](2) = zk;
218 1281855 : _weights[q] = gauss1D.w(i) * gauss1D.w(j) * jac1D.w(k);
219 : }
220 18819 : }
221 :
222 : } // namespace libMesh
|
