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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo