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