libMesh
quadrature_gm.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_gm.h"
20 #include "libmesh/quadrature_gauss.h"
21 
22 namespace libMesh
23 {
24 
25 // See also the file:
26 // quadrature_gm_3D.C
27 // for additional implementation.
28 
29 
30 
31 
32 // Constructor
33 QGrundmann_Moller::QGrundmann_Moller(const unsigned int d,
34  const Order o) : QBase(d,o)
35 {
36 }
37 
38 
39 // Destructor
40 QGrundmann_Moller::~QGrundmann_Moller()
41 {
42 }
43 
44 
45 
46 void QGrundmann_Moller::init_1D(const ElemType /*type_in*/,
47  unsigned int p)
48 {
49  QGauss gauss1D(1, static_cast<Order>(_order+2*p));
50 
51  // Swap points and weights with the about-to-be destroyed rule.
52  _points.swap(gauss1D.get_points());
53  _weights.swap(gauss1D.get_weights());
54 }
55 
56 
57 
58 
59 
60 void QGrundmann_Moller::gm_rule(unsigned int s, unsigned int dim)
61 {
62  // Only dim==2 or dim==3 is allowed
63  libmesh_assert(dim==2 || dim==3);
64 
65  // A GM rule of index s can integrate polynomials of degree 2*s+1 exactly
66  const unsigned int degree = 2*s+1;
67 
68  // The number of points for rule of index s is
69  // (dim+1+s)! / (dim+1)! / s!
70  // In 3D, this is = 1/24 * Pi_{i=1}^4 (s+i)
71  // In 2D, this is = 1/6 * Pi_{i=1}^3 (s+i)
72  const unsigned int n_pts = dim==2 ? (s+3)*(s+2)*(s+1) / 6 : (s+4)*(s+3)*(s+2)*(s+1) / 24;
73  //libMesh::out << "n_pts=" << n_pts << std::endl;
74 
75  // Allocate space for points and weights
76  _points.resize(n_pts);
77  _weights.resize(n_pts);
78 
79  // (-1)^i -> This one flips sign at each iteration of the i-loop below.
80  int one_pm=1;
81 
82  // Where we store all the integer point compositions/permutations
83  std::vector<std::vector<unsigned int>> permutations;
84 
85  // Index into the vector where we should start adding the next round of points/weights
86  std::size_t offset=0;
87 
88  // Implement the GM formula 4.1 on page 286 of the paper
89  for (unsigned int i=0; i<=s; ++i)
90  {
91  // Get all the ordered compositions (and their permutations)
92  // of |beta| = s-i into dim+1 parts
93  compose_all(s-i, dim+1, permutations);
94  //libMesh::out << "n. permutations=" << permutations.size() << std::endl;
95 
96  for (std::size_t p=0; p<permutations.size(); ++p)
97  {
98  // We use the first dim entries of each permutation to
99  // construct an integration point.
100  for (unsigned int j=0; j<dim; ++j)
101  _points[offset+p](j) =
102  static_cast<Real>(2.*permutations[p][j] + 1.) /
103  static_cast<Real>( degree + dim - 2.*i );
104  }
105 
106  // Compute the weight for this i, being careful to avoid overflow.
107  // This technique is borrowed from Burkardt's code as well.
108  // Use once for each of the points obtained from the permutations array.
109  Real weight = one_pm;
110 
111  // This for loop needs to run for dim, degree, or dim+degree-i iterations,
112  // whichever is largest.
113  const unsigned int weight_loop_index =
114  std::max(dim, std::max(degree, degree+dim-i))+1;
115 
116  for (unsigned int j=1; j<weight_loop_index; ++j)
117  {
118  if (j <= degree) // Accumulate (d+n-2i)^d term
119  weight *= static_cast<Real>(degree+dim-2*i);
120 
121  if (j <= 2*s) // Accumulate 2^{-2s}
122  weight *= 0.5;
123 
124  if (j <= i) // Accumulate (i!)^{-1}
125  weight /= static_cast<Real>(j);
126 
127  if (j <= degree+dim-i) // Accumulate ( (d+n-i)! )^{-1}
128  weight /= static_cast<Real>(j);
129  }
130 
131  // This is the weight for each of the points computed previously
132  for (std::size_t j=0; j<permutations.size(); ++j)
133  _weights[offset+j] = weight;
134 
135  // Change sign for next iteration
136  one_pm = -one_pm;
137 
138  // Update offset for the next set of points
139  offset += permutations.size();
140  }
141 }
142 
143 
144 
145 
146 // This routine for computing compositions and their permutations is
147 // originally due to:
148 //
149 // Albert Nijenhuis, Herbert Wilf,
150 // Combinatorial Algorithms for Computers and Calculators,
151 // Second Edition,
152 // Academic Press, 1978,
153 // ISBN: 0-12-519260-6,
154 // LC: QA164.N54.
155 //
156 // The routine is deceptively simple: I still don't understand exactly
157 // why it works, but it does.
158 void QGrundmann_Moller::compose_all(unsigned int s, // number to be compositioned
159  unsigned int p, // # of partitions
160  std::vector<std::vector<unsigned int>> & result)
161 {
162  // Clear out results remaining from previous calls
163  result.clear();
164 
165  // Allocate storage for a workspace. The workspace will periodically
166  // be copied into the result container.
167  std::vector<unsigned int> workspace(p);
168 
169  // The first result is always (s,0,...,0)
170  workspace[0] = s;
171  result.push_back(workspace);
172 
173  // the value of the first non-zero entry
174  unsigned int head_value=s;
175 
176  // When head_index=-1, it refers to "off the front" of the array. Therefore,
177  // this needs to be a regular int rather than unsigned. I initially tried to
178  // do this with head_index unsigned and an else statement below, but then there
179  // is the special case: (1,0,...,0) which does not work correctly.
180  int head_index = -1;
181 
182  // At the end, all the entries will be in the final slot of workspace
183  while (workspace.back() != s)
184  {
185  // If the previous head value is still larger than 1, reset the index
186  // to "off the front" of the array
187  if (head_value > 1)
188  head_index = -1;
189 
190  // Either move the index onto the front of the array or on to
191  // the next value.
192  head_index++;
193 
194  // Get current value of the head entry
195  head_value = workspace[head_index];
196 
197  // Put a zero into the head_index of the array. If head_index==0,
198  // this will be overwritten in the next line with head_value-1.
199  workspace[head_index] = 0;
200 
201  // The initial entry gets the current head value, minus 1.
202  // If head_value > 1, the next loop iteration will start back
203  // at workspace[0] again.
204  libmesh_assert_greater (head_value, 0);
205  workspace[0] = head_value - 1;
206 
207  // Increment the head+1 value
208  workspace[head_index+1] += 1;
209 
210  // Save this composition in the results
211  result.push_back(workspace);
212  }
213 }
214 
215 } // namespace libMesh
unsigned int dim
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
long double max(long double a, double b)
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
libmesh_assert(j)
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
QGauss(const unsigned int _dim, const Order _order=INVALID_ORDER)
Constructor.
const Order _order
The polynomial order which the quadrature rule is capable of integrating exactly. ...
Definition: quadrature.h:317