libMesh
quadrature_conical.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 #include "libmesh/quadrature_conical.h"
20 #include "libmesh/quadrature_gauss.h"
21 #include "libmesh/quadrature_jacobi.h"
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
43 QConical::~QConical()
44 {
45 }
46 
47 
48 
49 void QConical::init_1D(const ElemType /*type_in*/,
50  unsigned int p)
51 {
52  QGauss gauss1D(1, static_cast<Order>(_order+2*p));
53 
54  // Swap points and weights with the about-to-be destroyed rule.
55  _points.swap(gauss1D.get_points());
56  _weights.swap(gauss1D.get_weights());
57 }
58 
59 
60 
61 // Builds and scales a Gauss rule and a Jacobi rule.
62 // Then combines them to compute points and weights
63 // of a 2D conical product rule.
64 void QConical::conical_product_tri(unsigned int p)
65 {
66  // Be sure the underlying rule object was built with the same dimension as the
67  // rule we are about to construct.
68  libmesh_assert_equal_to (this->get_dim(), 2);
69 
70  QGauss gauss1D(1,static_cast<Order>(_order+2*p));
71  QJacobi jac1D(1,static_cast<Order>(_order+2*p),1,0);
72 
73  // The Gauss rule needs to be scaled to [0,1]
74  std::pair<Real, Real> old_range(-1.0L, 1.0L);
75  std::pair<Real, Real> new_range( 0.0L, 1.0L);
76  gauss1D.scale(old_range,
77  new_range);
78 
79  // Now construct the points and weights for the conical product rule.
80 
81  // Both rules should have the same number of points.
82  libmesh_assert_equal_to (gauss1D.n_points(), jac1D.n_points());
83 
84  // Save the number of points as a convenient variable
85  const unsigned int np = gauss1D.n_points();
86 
87  // Both rules should be between x=0 and x=1
88  libmesh_assert_greater_equal (gauss1D.qp(0)(0), 0.0);
89  libmesh_assert_less_equal (gauss1D.qp(np-1)(0), 1.0);
90  libmesh_assert_greater_equal (jac1D.qp(0)(0), 0.0);
91  libmesh_assert_less_equal (jac1D.qp(np-1)(0), 1.0);
92 
93  // Resize the points and weights vectors
94  _points.resize(np * np);
95  _weights.resize(np * np);
96 
97  // Compute the conical product
98  unsigned int gp = 0;
99  for (unsigned int i=0; i<np; i++)
100  for (unsigned int j=0; j<np; j++)
101  {
102  _points[gp](0) = jac1D.qp(j)(0); //s[j];
103  _points[gp](1) = gauss1D.qp(i)(0) * (1.-jac1D.qp(j)(0)); //r[i]*(1.-s[j]);
104  _weights[gp] = gauss1D.w(i) * jac1D.w(j); //A[i]*B[j];
105  gp++;
106  }
107 }
108 
109 
110 
111 
112 // Builds and scales a Gauss rule and a Jacobi rule.
113 // Then combines them to compute points and weights
114 // of a 3D conical product rule for the Tet.
115 void QConical::conical_product_tet(unsigned int p)
116 {
117  // Be sure the underlying rule object was built with the same dimension as the
118  // rule we are about to construct.
119  libmesh_assert_equal_to (this->get_dim(), 3);
120 
121  QGauss gauss1D(1,static_cast<Order>(_order+2*p));
122  QJacobi jacA1D(1,static_cast<Order>(_order+2*p),1,0);
123  QJacobi jacB1D(1,static_cast<Order>(_order+2*p),2,0);
124 
125  // The Gauss rule needs to be scaled to [0,1]
126  std::pair<Real, Real> old_range(-1.0L, 1.0L);
127  std::pair<Real, Real> new_range( 0.0L, 1.0L);
128  gauss1D.scale(old_range,
129  new_range);
130 
131  // Now construct the points and weights for the conical product rule.
132 
133  // All rules should have the same number of points
134  libmesh_assert_equal_to (gauss1D.n_points(), jacA1D.n_points());
135  libmesh_assert_equal_to (jacA1D.n_points(), jacB1D.n_points());
136 
137  // Save the number of points as a convenient variable
138  const unsigned int np = gauss1D.n_points();
139 
140  // All rules should be between x=0 and x=1
141  libmesh_assert_greater_equal (gauss1D.qp(0)(0), 0.0);
142  libmesh_assert_less_equal (gauss1D.qp(np-1)(0), 1.0);
143  libmesh_assert_greater_equal (jacA1D.qp(0)(0), 0.0);
144  libmesh_assert_less_equal (jacA1D.qp(np-1)(0), 1.0);
145  libmesh_assert_greater_equal (jacB1D.qp(0)(0), 0.0);
146  libmesh_assert_less_equal (jacB1D.qp(np-1)(0), 1.0);
147 
148  // Resize the points and weights vectors
149  _points.resize(np * np * np);
150  _weights.resize(np * np * np);
151 
152  // Compute the conical product
153  unsigned int gp = 0;
154  for (unsigned int i=0; i<np; i++)
155  for (unsigned int j=0; j<np; j++)
156  for (unsigned int k=0; k<np; k++)
157  {
158  _points[gp](0) = jacB1D.qp(k)(0); //t[k];
159  _points[gp](1) = jacA1D.qp(j)(0) * (1.-jacB1D.qp(k)(0)); //s[j]*(1.-t[k]);
160  _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]);
161  _weights[gp] = gauss1D.w(i) * jacA1D.w(j) * jacB1D.w(k); //A[i]*B[j]*C[k];
162  gp++;
163  }
164 }
165 
166 
167 
168 
169 
170 // Builds and scales a Gauss rule and a Jacobi rule.
171 // Then combines them to compute points and weights
172 // of a 3D conical product rule for the Pyramid. The integral
173 // over the reference Pyramid can be written (in LaTeX notation) as:
174 //
175 // If := \int_0^1 dz \int_{-(1-z)}^{(1-z)} dy \int_{-(1-z)}^{(1-z)} f(x,y,z) dx (1)
176 //
177 // (Imagine a stack of infinitely thin squares which decrease in size as
178 // you approach the apex.) Under the transformation of variables:
179 //
180 // z=w
181 // y=(1-z)v
182 // x=(1-z)u,
183 //
184 // The Jacobian determinant of this transformation is |J|=(1-w)^2, and
185 // the integral itself is transformed to:
186 //
187 // 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)
188 //
189 // The integral can now be approximated by the product of three 1D quadrature rules:
190 // A Jacobi rule with alpha==2, beta==0 in w, and Gauss rules in v and u. In this way
191 // we can obtain 3D rules to any order for which the 1D rules exist.
192 void QConical::conical_product_pyramid(unsigned int p)
193 {
194  // Be sure the underlying rule object was built with the same dimension as the
195  // rule we are about to construct.
196  libmesh_assert_equal_to (this->get_dim(), 3);
197 
198  QGauss gauss1D(1,static_cast<Order>(_order+2*p));
199  QJacobi jac1D(1,static_cast<Order>(_order+2*p),2,0);
200 
201  // These rules should have the same number of points
202  libmesh_assert_equal_to (gauss1D.n_points(), jac1D.n_points());
203 
204  // Save the number of points as a convenient variable
205  const unsigned int np = gauss1D.n_points();
206 
207  // Resize the points and weights vectors
208  _points.resize(np * np * np);
209  _weights.resize(np * np * np);
210 
211  // Compute the conical product
212  unsigned int q = 0;
213  for (unsigned int i=0; i<np; ++i)
214  for (unsigned int j=0; j<np; ++j)
215  for (unsigned int k=0; k<np; ++k, ++q)
216  {
217  const Real xi=gauss1D.qp(i)(0);
218  const Real yj=gauss1D.qp(j)(0);
219  const Real zk=jac1D.qp(k)(0);
220 
221  _points[q](0) = (1.-zk) * xi;
222  _points[q](1) = (1.-zk) * yj;
223  _points[q](2) = zk;
224  _weights[q] = gauss1D.w(i) * gauss1D.w(j) * jac1D.w(k);
225  }
226 }
227 
228 } // namespace libMesh
ElemType
Defines an enum for geometric element types.
The libMesh namespace provides an interface to certain functionality in the library.
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int get_dim() const
Definition: quadrature.h:122
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
const Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:317