libMesh
quadrature_monomial.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // libMesh includes
20 #include "libmesh/quadrature_monomial.h"
21 #include "libmesh/enum_quadrature_type.h"
22 
23 namespace libMesh
24 {
25 
26 
27 // See the files:
28 // quadrature_monomial_2D.C
29 // quadrature_monomial_3D.C
30 // for implementation of specific element types.
31 
32 QuadratureType QMonomial::type() const
33 {
34  return QMONOMIAL;
35 }
36 
37 std::unique_ptr<QBase> QMonomial::clone() const
38 {
39  return std::make_unique<QMonomial>(*this);
40 }
41 
42 void QMonomial::wissmann_rule(const Real rule_data[][3],
43  const unsigned int n_pts)
44 {
45  for (unsigned int i=0, c=0; i<n_pts; ++i)
46  {
47  _points[c] = Point( rule_data[i][0], rule_data[i][1] );
48  _weights[c++] = rule_data[i][2];
49 
50  // This may be an (x1,x2) -> (-x1,x2) point, in which case
51  // we will also generate the mirror point using the same weight.
52  if (rule_data[i][0] != static_cast<Real>(0.0))
53  {
54  _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
55  _weights[c++] = rule_data[i][2];
56  }
57  }
58 }
59 
60 
61 
62 void QMonomial::stroud_rule(const Real rule_data[][3],
63  const unsigned int * rule_symmetry,
64  const unsigned int n_pts)
65 {
66  for (unsigned int i=0, c=0; i<n_pts; ++i)
67  {
68  const Real
69  x=rule_data[i][0],
70  y=rule_data[i][1],
71  wt=rule_data[i][2];
72 
73  switch(rule_symmetry[i])
74  {
75  case 0: // Single point (no symmetry)
76  {
77  _points[c] = Point( x, y);
78  _weights[c++] = wt;
79 
80  break;
81  }
82  case 1: // Fully-symmetric (x,y)
83  {
84  _points[c] = Point( x, y);
85  _weights[c++] = wt;
86 
87  _points[c] = Point(-x, y);
88  _weights[c++] = wt;
89 
90  _points[c] = Point( x,-y);
91  _weights[c++] = wt;
92 
93  _points[c] = Point(-x,-y);
94  _weights[c++] = wt;
95 
96  _points[c] = Point( y, x);
97  _weights[c++] = wt;
98 
99  _points[c] = Point(-y, x);
100  _weights[c++] = wt;
101 
102  _points[c] = Point( y,-x);
103  _weights[c++] = wt;
104 
105  _points[c] = Point(-y,-x);
106  _weights[c++] = wt;
107 
108  break;
109  }
110  case 2: // Fully-symmetric (x,x)
111  {
112  _points[c] = Point( x, x);
113  _weights[c++] = wt;
114 
115  _points[c] = Point(-x, x);
116  _weights[c++] = wt;
117 
118  _points[c] = Point( x,-x);
119  _weights[c++] = wt;
120 
121  _points[c] = Point(-x,-x);
122  _weights[c++] = wt;
123 
124  break;
125  }
126  case 3: // Fully-symmetric (x,0)
127  {
128  libmesh_assert_equal_to (y, 0.0);
129 
130  _points[c] = Point( x,0.);
131  _weights[c++] = wt;
132 
133  _points[c] = Point(-x,0.);
134  _weights[c++] = wt;
135 
136  _points[c] = Point(0., x);
137  _weights[c++] = wt;
138 
139  _points[c] = Point(0.,-x);
140  _weights[c++] = wt;
141 
142  break;
143  }
144  case 4: // Rotational invariant
145  {
146  _points[c] = Point( x, y);
147  _weights[c++] = wt;
148 
149  _points[c] = Point(-x,-y);
150  _weights[c++] = wt;
151 
152  _points[c] = Point(-y, x);
153  _weights[c++] = wt;
154 
155  _points[c] = Point( y,-x);
156  _weights[c++] = wt;
157 
158  break;
159  }
160  case 5: // Partial symmetry (Wissman's rules)
161  {
162  libmesh_assert_not_equal_to (x, 0.0);
163 
164  _points[c] = Point( x, y);
165  _weights[c++] = wt;
166 
167  _points[c] = Point(-x, y);
168  _weights[c++] = wt;
169 
170  break;
171  }
172  case 6: // Rectangular symmetry
173  {
174  _points[c] = Point( x, y);
175  _weights[c++] = wt;
176 
177  _points[c] = Point(-x, y);
178  _weights[c++] = wt;
179 
180  _points[c] = Point(-x,-y);
181  _weights[c++] = wt;
182 
183  _points[c] = Point( x,-y);
184  _weights[c++] = wt;
185 
186  break;
187  }
188  case 7: // Central symmetry
189  {
190  libmesh_assert_equal_to (x, 0.0);
191  libmesh_assert_not_equal_to (y, 0.0);
192 
193  _points[c] = Point(0., y);
194  _weights[c++] = wt;
195 
196  _points[c] = Point(0.,-y);
197  _weights[c++] = wt;
198 
199  break;
200  }
201  default:
202  libmesh_error_msg("Unknown symmetry!");
203  } // end switch(rule_symmetry[i])
204  }
205 }
206 
207 
208 
209 
210 void QMonomial::kim_rule(const Real rule_data[][4],
211  const unsigned int * rule_id,
212  const unsigned int n_pts)
213 {
214  for (unsigned int i=0, c=0; i<n_pts; ++i)
215  {
216  const Real
217  x=rule_data[i][0],
218  y=rule_data[i][1],
219  z=rule_data[i][2],
220  wt=rule_data[i][3];
221 
222  switch(rule_id[i])
223  {
224  case 0: // (0,0,0) 1 permutation
225  {
226  _points[c] = Point( x, y, z); _weights[c++] = wt;
227 
228  break;
229  }
230  case 1: // (x,0,0) 6 permutations
231  {
232  _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
233  _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
234  _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
235  _points[c] = Point(0., x, 0.); _weights[c++] = wt;
236  _points[c] = Point(0., 0., -x); _weights[c++] = wt;
237  _points[c] = Point(0., 0., x); _weights[c++] = wt;
238 
239  break;
240  }
241  case 2: // (x,x,0) 12 permutations
242  {
243  _points[c] = Point( x, x, 0.); _weights[c++] = wt;
244  _points[c] = Point( x, -x, 0.); _weights[c++] = wt;
245  _points[c] = Point(-x, -x, 0.); _weights[c++] = wt;
246  _points[c] = Point(-x, x, 0.); _weights[c++] = wt;
247  _points[c] = Point( x, 0., -x); _weights[c++] = wt;
248  _points[c] = Point( x, 0., x); _weights[c++] = wt;
249  _points[c] = Point(0., x, -x); _weights[c++] = wt;
250  _points[c] = Point(0., x, x); _weights[c++] = wt;
251  _points[c] = Point(0., -x, -x); _weights[c++] = wt;
252  _points[c] = Point(-x, 0., -x); _weights[c++] = wt;
253  _points[c] = Point(0., -x, x); _weights[c++] = wt;
254  _points[c] = Point(-x, 0., x); _weights[c++] = wt;
255 
256  break;
257  }
258  case 3: // (x,y,0) 24 permutations
259  {
260  _points[c] = Point( x, y, 0.); _weights[c++] = wt;
261  _points[c] = Point( y, -x, 0.); _weights[c++] = wt;
262  _points[c] = Point(-x, -y, 0.); _weights[c++] = wt;
263  _points[c] = Point(-y, x, 0.); _weights[c++] = wt;
264  _points[c] = Point( x, 0., -y); _weights[c++] = wt;
265  _points[c] = Point( x, -y, 0.); _weights[c++] = wt;
266  _points[c] = Point( x, 0., y); _weights[c++] = wt;
267  _points[c] = Point(0., y, -x); _weights[c++] = wt;
268  _points[c] = Point(-x, y, 0.); _weights[c++] = wt;
269  _points[c] = Point(0., y, x); _weights[c++] = wt;
270  _points[c] = Point( y, 0., -x); _weights[c++] = wt;
271  _points[c] = Point(0., -y, -x); _weights[c++] = wt;
272  _points[c] = Point(-y, 0., -x); _weights[c++] = wt;
273  _points[c] = Point( y, x, 0.); _weights[c++] = wt;
274  _points[c] = Point(-y, -x, 0.); _weights[c++] = wt;
275  _points[c] = Point( y, 0., x); _weights[c++] = wt;
276  _points[c] = Point(0., -y, x); _weights[c++] = wt;
277  _points[c] = Point(-y, 0., x); _weights[c++] = wt;
278  _points[c] = Point(-x, 0., y); _weights[c++] = wt;
279  _points[c] = Point(0., -x, -y); _weights[c++] = wt;
280  _points[c] = Point(0., -x, y); _weights[c++] = wt;
281  _points[c] = Point(-x, 0., -y); _weights[c++] = wt;
282  _points[c] = Point(0., x, y); _weights[c++] = wt;
283  _points[c] = Point(0., x, -y); _weights[c++] = wt;
284 
285  break;
286  }
287  case 4: // (x,x,x) 8 permutations
288  {
289  _points[c] = Point( x, x, x); _weights[c++] = wt;
290  _points[c] = Point( x, -x, x); _weights[c++] = wt;
291  _points[c] = Point(-x, -x, x); _weights[c++] = wt;
292  _points[c] = Point(-x, x, x); _weights[c++] = wt;
293  _points[c] = Point( x, x, -x); _weights[c++] = wt;
294  _points[c] = Point( x, -x, -x); _weights[c++] = wt;
295  _points[c] = Point(-x, x, -x); _weights[c++] = wt;
296  _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
297 
298  break;
299  }
300  case 5: // (x,x,z) 24 permutations
301  {
302  _points[c] = Point( x, x, z); _weights[c++] = wt;
303  _points[c] = Point( x, -x, z); _weights[c++] = wt;
304  _points[c] = Point(-x, -x, z); _weights[c++] = wt;
305  _points[c] = Point(-x, x, z); _weights[c++] = wt;
306  _points[c] = Point( x, z, -x); _weights[c++] = wt;
307  _points[c] = Point( x, -x, -z); _weights[c++] = wt;
308  _points[c] = Point( x, -z, x); _weights[c++] = wt;
309  _points[c] = Point( z, x, -x); _weights[c++] = wt;
310  _points[c] = Point(-x, x, -z); _weights[c++] = wt;
311  _points[c] = Point(-z, x, x); _weights[c++] = wt;
312  _points[c] = Point( x, -z, -x); _weights[c++] = wt;
313  _points[c] = Point(-z, -x, -x); _weights[c++] = wt;
314  _points[c] = Point(-x, z, -x); _weights[c++] = wt;
315  _points[c] = Point( x, x, -z); _weights[c++] = wt;
316  _points[c] = Point(-x, -x, -z); _weights[c++] = wt;
317  _points[c] = Point( x, z, x); _weights[c++] = wt;
318  _points[c] = Point( z, -x, x); _weights[c++] = wt;
319  _points[c] = Point(-x, -z, x); _weights[c++] = wt;
320  _points[c] = Point(-x, z, x); _weights[c++] = wt;
321  _points[c] = Point( z, -x, -x); _weights[c++] = wt;
322  _points[c] = Point(-z, -x, x); _weights[c++] = wt;
323  _points[c] = Point(-x, -z, -x); _weights[c++] = wt;
324  _points[c] = Point( z, x, x); _weights[c++] = wt;
325  _points[c] = Point(-z, x, -x); _weights[c++] = wt;
326 
327  break;
328  }
329  case 6: // (x,y,z) 24 permutations
330  {
331  _points[c] = Point( x, y, z); _weights[c++] = wt;
332  _points[c] = Point( y, -x, z); _weights[c++] = wt;
333  _points[c] = Point(-x, -y, z); _weights[c++] = wt;
334  _points[c] = Point(-y, x, z); _weights[c++] = wt;
335  _points[c] = Point( x, z, -y); _weights[c++] = wt;
336  _points[c] = Point( x, -y, -z); _weights[c++] = wt;
337  _points[c] = Point( x, -z, y); _weights[c++] = wt;
338  _points[c] = Point( z, y, -x); _weights[c++] = wt;
339  _points[c] = Point(-x, y, -z); _weights[c++] = wt;
340  _points[c] = Point(-z, y, x); _weights[c++] = wt;
341  _points[c] = Point( y, -z, -x); _weights[c++] = wt;
342  _points[c] = Point(-z, -y, -x); _weights[c++] = wt;
343  _points[c] = Point(-y, z, -x); _weights[c++] = wt;
344  _points[c] = Point( y, x, -z); _weights[c++] = wt;
345  _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
346  _points[c] = Point( y, z, x); _weights[c++] = wt;
347  _points[c] = Point( z, -y, x); _weights[c++] = wt;
348  _points[c] = Point(-y, -z, x); _weights[c++] = wt;
349  _points[c] = Point(-x, z, y); _weights[c++] = wt;
350  _points[c] = Point( z, -x, -y); _weights[c++] = wt;
351  _points[c] = Point(-z, -x, y); _weights[c++] = wt;
352  _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
353  _points[c] = Point( z, x, y); _weights[c++] = wt;
354  _points[c] = Point(-z, x, -y); _weights[c++] = wt;
355 
356  break;
357  }
358  default:
359  libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!");
360  } // end switch(rule_id[i])
361  }
362 }
363 
364 } // namespace libMesh
libMesh::QMonomial::clone
virtual std::unique_ptr< QBase > clone() const override
Definition: quadrature_monomial.C:37
libMesh::QMonomial::wissmann_rule
void wissmann_rule(const Real rule_data[][3], const unsigned int n_pts)
Wissmann published three interesting "partially symmetric" rules for integrating degree 4...
Definition: quadrature_monomial.C:42
libMesh::QuadratureType
QuadratureType
Defines an enum for currently available quadrature rules.
Definition: enum_quadrature_type.h:33
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::QMonomial::type
virtual QuadratureType type() const override
Definition: quadrature_monomial.C:32
libMesh::QBase::_points
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:383
libMesh::QBase::_weights
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:389
libMesh::QMonomial::kim_rule
void kim_rule(const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
Rules from Kim and Song, Comm.
Definition: quadrature_monomial.C:210
libMesh::QMonomial::stroud_rule
void stroud_rule(const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
Stroud&#39;s rules for quads and hexes can have one of several different types of symmetry.
Definition: quadrature_monomial.C:62
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:143
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Generated on Wed Apr 17 2024 11:03:50 for libMesh by   doxygen 1.8.14