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