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

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

Hosted By:
SourceForge.net Logo