quadrature_monomial.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 
20 
21 namespace libMesh
22 {
23 
24 
25 // See the files:
26 
27 // quadrature_monomial_2D.C
28 // quadrature_monomial_3D.C
29 
30 // for implementation of specific element types.
31 
32 // Constructor
33 QMonomial::QMonomial(const unsigned int d,
34  const Order o) : QBase(d,o)
35 {
36 }
37 
38 
39 // Destructor
41 {
42 }
43 
44 void QMonomial::wissmann_rule(const Real rule_data[][3],
45  const unsigned int n_pts)
46 {
47  for (unsigned int i=0, c=0; i<n_pts; ++i)
48  {
49  _points[c] = Point( rule_data[i][0], rule_data[i][1] );
50  _weights[c++] = rule_data[i][2];
51 
52  // This may be an (x1,x2) -> (-x1,x2) point, in which case
53  // we will also generate the mirror point using the same weight.
54  if (rule_data[i][0] != static_cast<Real>(0.0))
55  {
56  _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
57  _weights[c++] = rule_data[i][2];
58  }
59  }
60 }
61 
62 
63 
64 void QMonomial::stroud_rule(const Real rule_data[][3],
65  const unsigned int* rule_symmetry,
66  const unsigned int n_pts)
67 {
68  for (unsigned int i=0, c=0; i<n_pts; ++i)
69  {
70  const Real
71  x=rule_data[i][0],
72  y=rule_data[i][1],
73  wt=rule_data[i][2];
74 
75  switch(rule_symmetry[i])
76  {
77  case 0: // Single point (no symmetry)
78  {
79  _points[c] = Point( x, y);
80  _weights[c++] = wt;
81 
82  break;
83  }
84  case 1: // Fully-symmetric (x,y)
85  {
86  _points[c] = Point( x, y);
87  _weights[c++] = wt;
88 
89  _points[c] = Point(-x, y);
90  _weights[c++] = wt;
91 
92  _points[c] = Point( x,-y);
93  _weights[c++] = wt;
94 
95  _points[c] = Point(-x,-y);
96  _weights[c++] = wt;
97 
98  _points[c] = Point( y, x);
99  _weights[c++] = wt;
100 
101  _points[c] = Point(-y, x);
102  _weights[c++] = wt;
103 
104  _points[c] = Point( y,-x);
105  _weights[c++] = wt;
106 
107  _points[c] = Point(-y,-x);
108  _weights[c++] = wt;
109 
110  break;
111  }
112  case 2: // Fully-symmetric (x,x)
113  {
114  _points[c] = Point( x, x);
115  _weights[c++] = wt;
116 
117  _points[c] = Point(-x, x);
118  _weights[c++] = wt;
119 
120  _points[c] = Point( x,-x);
121  _weights[c++] = wt;
122 
123  _points[c] = Point(-x,-x);
124  _weights[c++] = wt;
125 
126  break;
127  }
128  case 3: // Fully-symmetric (x,0)
129  {
130  libmesh_assert_equal_to (y, 0.0);
131 
132  _points[c] = Point( x,0.);
133  _weights[c++] = wt;
134 
135  _points[c] = Point(-x,0.);
136  _weights[c++] = wt;
137 
138  _points[c] = Point(0., x);
139  _weights[c++] = wt;
140 
141  _points[c] = Point(0.,-x);
142  _weights[c++] = wt;
143 
144  break;
145  }
146  case 4: // Rotational invariant
147  {
148  _points[c] = Point( x, y);
149  _weights[c++] = wt;
150 
151  _points[c] = Point(-x,-y);
152  _weights[c++] = wt;
153 
154  _points[c] = Point(-y, x);
155  _weights[c++] = wt;
156 
157  _points[c] = Point( y,-x);
158  _weights[c++] = wt;
159 
160  break;
161  }
162  case 5: // Partial symmetry (Wissman's rules)
163  {
164  libmesh_assert_not_equal_to (x, 0.0);
165 
166  _points[c] = Point( x, y);
167  _weights[c++] = wt;
168 
169  _points[c] = Point(-x, y);
170  _weights[c++] = wt;
171 
172  break;
173  }
174  case 6: // Rectangular symmetry
175  {
176  _points[c] = Point( x, y);
177  _weights[c++] = wt;
178 
179  _points[c] = Point(-x, y);
180  _weights[c++] = wt;
181 
182  _points[c] = Point(-x,-y);
183  _weights[c++] = wt;
184 
185  _points[c] = Point( x,-y);
186  _weights[c++] = wt;
187 
188  break;
189  }
190  case 7: // Central symmetry
191  {
192  libmesh_assert_equal_to (x, 0.0);
193  libmesh_assert_not_equal_to (y, 0.0);
194 
195  _points[c] = Point(0., y);
196  _weights[c++] = wt;
197 
198  _points[c] = Point(0.,-y);
199  _weights[c++] = wt;
200 
201  break;
202  }
203  default:
204  {
205  libMesh::err << "Unknown symmetry!" << std::endl;
206  libmesh_error();
207  }
208  } // end switch(rule_symmetry[i])
209  }
210 }
211 
212 
213 
214 
215 void QMonomial::kim_rule(const Real rule_data[][4],
216  const unsigned int* rule_id,
217  const unsigned int n_pts)
218 {
219  for (unsigned int i=0, c=0; i<n_pts; ++i)
220  {
221  const Real
222  x=rule_data[i][0],
223  y=rule_data[i][1],
224  z=rule_data[i][2],
225  wt=rule_data[i][3];
226 
227  switch(rule_id[i])
228  {
229  case 0: // (0,0,0) 1 permutation
230  {
231  _points[c] = Point( x, y, z); _weights[c++] = wt;
232 
233  break;
234  }
235  case 1: // (x,0,0) 6 permutations
236  {
237  _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
238  _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
239  _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
240  _points[c] = Point(0., x, 0.); _weights[c++] = wt;
241  _points[c] = Point(0., 0., -x); _weights[c++] = wt;
242  _points[c] = Point(0., 0., x); _weights[c++] = wt;
243 
244  break;
245  }
246  case 2: // (x,x,0) 12 permutations
247  {
248  _points[c] = Point( x, x, 0.); _weights[c++] = wt;
249  _points[c] = Point( x, -x, 0.); _weights[c++] = wt;
250  _points[c] = Point(-x, -x, 0.); _weights[c++] = wt;
251  _points[c] = Point(-x, x, 0.); _weights[c++] = wt;
252  _points[c] = Point( x, 0., -x); _weights[c++] = wt;
253  _points[c] = Point( x, 0., x); _weights[c++] = wt;
254  _points[c] = Point(0., x, -x); _weights[c++] = wt;
255  _points[c] = Point(0., x, x); _weights[c++] = wt;
256  _points[c] = Point(0., -x, -x); _weights[c++] = wt;
257  _points[c] = Point(-x, 0., -x); _weights[c++] = wt;
258  _points[c] = Point(0., -x, x); _weights[c++] = wt;
259  _points[c] = Point(-x, 0., x); _weights[c++] = wt;
260 
261  break;
262  }
263  case 3: // (x,y,0) 24 permutations
264  {
265  _points[c] = Point( x, y, 0.); _weights[c++] = wt;
266  _points[c] = Point( y, -x, 0.); _weights[c++] = wt;
267  _points[c] = Point(-x, -y, 0.); _weights[c++] = wt;
268  _points[c] = Point(-y, x, 0.); _weights[c++] = wt;
269  _points[c] = Point( x, 0., -y); _weights[c++] = wt;
270  _points[c] = Point( x, -y, 0.); _weights[c++] = wt;
271  _points[c] = Point( x, 0., y); _weights[c++] = wt;
272  _points[c] = Point(0., y, -x); _weights[c++] = wt;
273  _points[c] = Point(-x, y, 0.); _weights[c++] = wt;
274  _points[c] = Point(0., y, x); _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( y, x, 0.); _weights[c++] = wt;
279  _points[c] = Point(-y, -x, 0.); _weights[c++] = wt;
280  _points[c] = Point( y, 0., x); _weights[c++] = wt;
281  _points[c] = Point(0., -y, x); _weights[c++] = wt;
282  _points[c] = Point(-y, 0., x); _weights[c++] = wt;
283  _points[c] = Point(-x, 0., y); _weights[c++] = wt;
284  _points[c] = Point(0., -x, -y); _weights[c++] = wt;
285  _points[c] = Point(0., -x, y); _weights[c++] = wt;
286  _points[c] = Point(-x, 0., -y); _weights[c++] = wt;
287  _points[c] = Point(0., x, y); _weights[c++] = wt;
288  _points[c] = Point(0., x, -y); _weights[c++] = wt;
289 
290  break;
291  }
292  case 4: // (x,x,x) 8 permutations
293  {
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  _points[c] = Point(-x, x, x); _weights[c++] = wt;
298  _points[c] = Point( x, x, -x); _weights[c++] = wt;
299  _points[c] = Point( x, -x, -x); _weights[c++] = wt;
300  _points[c] = Point(-x, x, -x); _weights[c++] = wt;
301  _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
302 
303  break;
304  }
305  case 5: // (x,x,z) 24 permutations
306  {
307  _points[c] = Point( x, x, z); _weights[c++] = wt;
308  _points[c] = Point( x, -x, z); _weights[c++] = wt;
309  _points[c] = Point(-x, -x, z); _weights[c++] = wt;
310  _points[c] = Point(-x, x, z); _weights[c++] = wt;
311  _points[c] = Point( x, z, -x); _weights[c++] = wt;
312  _points[c] = Point( x, -x, -z); _weights[c++] = wt;
313  _points[c] = Point( x, -z, x); _weights[c++] = wt;
314  _points[c] = Point( z, x, -x); _weights[c++] = wt;
315  _points[c] = Point(-x, x, -z); _weights[c++] = wt;
316  _points[c] = Point(-z, x, x); _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, x, -z); _weights[c++] = wt;
321  _points[c] = Point(-x, -x, -z); _weights[c++] = wt;
322  _points[c] = Point( x, z, x); _weights[c++] = wt;
323  _points[c] = Point( z, -x, x); _weights[c++] = wt;
324  _points[c] = Point(-x, -z, x); _weights[c++] = wt;
325  _points[c] = Point(-x, z, x); _weights[c++] = wt;
326  _points[c] = Point( z, -x, -x); _weights[c++] = wt;
327  _points[c] = Point(-z, -x, x); _weights[c++] = wt;
328  _points[c] = Point(-x, -z, -x); _weights[c++] = wt;
329  _points[c] = Point( z, x, x); _weights[c++] = wt;
330  _points[c] = Point(-z, x, -x); _weights[c++] = wt;
331 
332  break;
333  }
334  case 6: // (x,y,z) 24 permutations
335  {
336  _points[c] = Point( x, y, z); _weights[c++] = wt;
337  _points[c] = Point( y, -x, z); _weights[c++] = wt;
338  _points[c] = Point(-x, -y, z); _weights[c++] = wt;
339  _points[c] = Point(-y, x, z); _weights[c++] = wt;
340  _points[c] = Point( x, z, -y); _weights[c++] = wt;
341  _points[c] = Point( x, -y, -z); _weights[c++] = wt;
342  _points[c] = Point( x, -z, y); _weights[c++] = wt;
343  _points[c] = Point( z, y, -x); _weights[c++] = wt;
344  _points[c] = Point(-x, y, -z); _weights[c++] = wt;
345  _points[c] = Point(-z, y, x); _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( y, x, -z); _weights[c++] = wt;
350  _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
351  _points[c] = Point( y, z, x); _weights[c++] = wt;
352  _points[c] = Point( z, -y, x); _weights[c++] = wt;
353  _points[c] = Point(-y, -z, x); _weights[c++] = wt;
354  _points[c] = Point(-x, z, y); _weights[c++] = wt;
355  _points[c] = Point( z, -x, -y); _weights[c++] = wt;
356  _points[c] = Point(-z, -x, y); _weights[c++] = wt;
357  _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
358  _points[c] = Point( z, x, y); _weights[c++] = wt;
359  _points[c] = Point(-z, x, -y); _weights[c++] = wt;
360 
361  break;
362  }
363  default:
364  {
365  libMesh::err << "Unknown rule ID: " << rule_id[i] << "!" << std::endl;
366  libmesh_error();
367  }
368  } // end switch(rule_id[i])
369  }
370 }
371 
372 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo