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