quadrature_gauss.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 // See the files:
25 
26 // quadrature_gauss_1D.C
27 // quadrature_gauss_2D.C
28 // quadrature_gauss_3D.C
29 
30 // for implementation of specific element types.
31 
32 
33 void QGauss::keast_rule(const Real rule_data[][4],
34  const unsigned int n_pts)
35 {
36  // Like the Dunavant rule, the input data should have 4 columns. These columns
37  // have the following format and implied permutations (w=weight).
38  // {a, 0, 0, w} = 1-permutation (a,a,a)
39  // {a, b, 0, w} = 4-permutation (a,b,b), (b,a,b), (b,b,a), (b,b,b)
40  // {a, 0, b, w} = 6-permutation (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a)
41  // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a)
42  // (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a)
43 
44  // Always insert into the points & weights vector relative to the offset
45  unsigned int offset=0;
46 
47 
48  for (unsigned int p=0; p<n_pts; ++p)
49  {
50 
51  // There must always be a non-zero entry to start the row
52  libmesh_assert_not_equal_to (rule_data[p][0], static_cast<Real>(0.0));
53 
54  // A zero weight may imply you did not set up the raw data correctly
55  libmesh_assert_not_equal_to (rule_data[p][3], static_cast<Real>(0.0));
56 
57  // What kind of point is this?
58  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
59  // Two non-zero entries in first 3 cols ? 3-perm point = 3
60  // Three non-zero entries ? 6-perm point = 6
61  unsigned int pointtype=1;
62 
63  if (rule_data[p][1] != static_cast<Real>(0.0))
64  {
65  if (rule_data[p][2] != static_cast<Real>(0.0))
66  pointtype = 12;
67  else
68  pointtype = 4;
69  }
70  else
71  {
72  // The second entry is zero. What about the third?
73  if (rule_data[p][2] != static_cast<Real>(0.0))
74  pointtype = 6;
75  }
76 
77 
78  switch (pointtype)
79  {
80  case 1:
81  {
82  // Be sure we have enough space to insert this point
83  libmesh_assert_less (offset + 0, _points.size());
84 
85  const Real a = rule_data[p][0];
86 
87  // The point has only a single permutation (the centroid!)
88  _points[offset + 0] = Point(a,a,a);
89 
90  // The weight is always the last entry in the row.
91  _weights[offset + 0] = rule_data[p][3];
92 
93  offset += pointtype;
94  break;
95  }
96 
97  case 4:
98  {
99  // Be sure we have enough space to insert these points
100  libmesh_assert_less (offset + 3, _points.size());
101 
102  const Real a = rule_data[p][0];
103  const Real b = rule_data[p][1];
104  const Real wt = rule_data[p][3];
105 
106  // Here it's understood the second entry is to be used twice, and
107  // thus there are three possible permutations.
108  _points[offset + 0] = Point(a,b,b);
109  _points[offset + 1] = Point(b,a,b);
110  _points[offset + 2] = Point(b,b,a);
111  _points[offset + 3] = Point(b,b,b);
112 
113  for (unsigned int j=0; j<pointtype; ++j)
114  _weights[offset + j] = wt;
115 
116  offset += pointtype;
117  break;
118  }
119 
120  case 6:
121  {
122  // Be sure we have enough space to insert these points
123  libmesh_assert_less (offset + 5, _points.size());
124 
125  const Real a = rule_data[p][0];
126  const Real b = rule_data[p][2];
127  const Real wt = rule_data[p][3];
128 
129  // Three individual entries with six permutations.
130  _points[offset + 0] = Point(a,a,b);
131  _points[offset + 1] = Point(a,b,b);
132  _points[offset + 2] = Point(b,b,a);
133  _points[offset + 3] = Point(b,a,b);
134  _points[offset + 4] = Point(b,a,a);
135  _points[offset + 5] = Point(a,b,a);
136 
137  for (unsigned int j=0; j<pointtype; ++j)
138  _weights[offset + j] = wt;
139 
140  offset += pointtype;
141  break;
142  }
143 
144 
145  case 12:
146  {
147  // Be sure we have enough space to insert these points
148  libmesh_assert_less (offset + 11, _points.size());
149 
150  const Real a = rule_data[p][0];
151  const Real b = rule_data[p][1];
152  const Real c = rule_data[p][2];
153  const Real wt = rule_data[p][3];
154 
155  // Three individual entries with six permutations.
156  _points[offset + 0] = Point(a,a,b); _points[offset + 6] = Point(a,b,c);
157  _points[offset + 1] = Point(a,a,c); _points[offset + 7] = Point(a,c,b);
158  _points[offset + 2] = Point(b,a,a); _points[offset + 8] = Point(b,a,c);
159  _points[offset + 3] = Point(c,a,a); _points[offset + 9] = Point(b,c,a);
160  _points[offset + 4] = Point(a,b,a); _points[offset + 10] = Point(c,a,b);
161  _points[offset + 5] = Point(a,c,a); _points[offset + 11] = Point(c,b,a);
162 
163  for (unsigned int j=0; j<pointtype; ++j)
164  _weights[offset + j] = wt;
165 
166  offset += pointtype;
167  break;
168  }
169 
170  default:
171  {
172  libMesh::err << "Don't know what to do with this many permutation points!" << std::endl;
173  libmesh_error();
174  }
175  }
176 
177  }
178 
179 }
180 
181 
182 // A number of different rules for triangles can be described by
183 // permutations of the following types of points:
184 // I: "1"-permutation, (1/3,1/3) (single point only)
185 // II: 3-permutation, (a,a,1-2a)
186 // III: 6-permutation, (a,b,1-a-b)
187 // The weights for a given set of permutations are all the same.
188 void QGauss::dunavant_rule2(const Real* wts,
189  const Real* a,
190  const Real* b,
191  const unsigned int* permutation_ids,
192  unsigned int n_wts)
193 {
194  // Figure out how many total points by summing up the entries
195  // in the permutation_ids array, and resize the _points and _weights
196  // vectors appropriately.
197  unsigned int total_pts = 0;
198  for (unsigned int p=0; p<n_wts; ++p)
199  total_pts += permutation_ids[p];
200 
201  // Resize point and weight vectors appropriately.
202  _points.resize(total_pts);
203  _weights.resize(total_pts);
204 
205  // Always insert into the points & weights vector relative to the offset
206  unsigned int offset=0;
207 
208  for (unsigned int p=0; p<n_wts; ++p)
209  {
210  switch (permutation_ids[p])
211  {
212  case 1:
213  {
214  // The point has only a single permutation (the centroid!)
215  // So we don't even need to look in the a or b arrays.
216  _points [offset + 0] = Point(1.0L/3.0L, 1.0L/3.0L);
217  _weights[offset + 0] = wts[p];
218 
219  offset += 1;
220  break;
221  }
222 
223 
224  case 3:
225  {
226  // For this type of rule, don't need to look in the b array.
227  _points[offset + 0] = Point(a[p], a[p]); // (a,a)
228  _points[offset + 1] = Point(a[p], 1.L-2.L*a[p]); // (a,1-2a)
229  _points[offset + 2] = Point(1.L-2.L*a[p], a[p]); // (1-2a,a)
230 
231  for (unsigned int j=0; j<3; ++j)
232  _weights[offset + j] = wts[p];
233 
234  offset += 3;
235  break;
236  }
237 
238  case 6:
239  {
240  // This type of point uses all 3 arrays...
241  _points[offset + 0] = Point(a[p], b[p]);
242  _points[offset + 1] = Point(b[p], a[p]);
243  _points[offset + 2] = Point(a[p], 1.L-a[p]-b[p]);
244  _points[offset + 3] = Point(1.L-a[p]-b[p], a[p]);
245  _points[offset + 4] = Point(b[p], 1.L-a[p]-b[p]);
246  _points[offset + 5] = Point(1.L-a[p]-b[p], b[p]);
247 
248  for (unsigned int j=0; j<6; ++j)
249  _weights[offset + j] = wts[p];
250 
251  offset += 6;
252  break;
253  }
254 
255  default:
256  {
257  libMesh::err << "Unknown permutation id: " << permutation_ids[p] << "!" << std::endl;
258  libmesh_error();
259  }
260  }
261  }
262 
263 }
264 
265 
266 void QGauss::dunavant_rule(const Real rule_data[][4],
267  const unsigned int n_pts)
268 {
269  // The input data array has 4 columns. The first 3 are the permutation points.
270  // The last column is the weights for a given set of permutation points. A zero
271  // in two of the first 3 columns implies the point is a 1-permutation (centroid).
272  // A zero in one of the first 3 columns implies the point is a 3-permutation.
273  // Otherwise each point is assumed to be a 6-permutation.
274 
275  // Always insert into the points & weights vector relative to the offset
276  unsigned int offset=0;
277 
278 
279  for (unsigned int p=0; p<n_pts; ++p)
280  {
281 
282  // There must always be a non-zero entry to start the row
283  libmesh_assert_not_equal_to ( rule_data[p][0], static_cast<Real>(0.0) );
284 
285  // A zero weight may imply you did not set up the raw data correctly
286  libmesh_assert_not_equal_to ( rule_data[p][3], static_cast<Real>(0.0) );
287 
288  // What kind of point is this?
289  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
290  // Two non-zero entries in first 3 cols ? 3-perm point = 3
291  // Three non-zero entries ? 6-perm point = 6
292  unsigned int pointtype=1;
293 
294  if (rule_data[p][1] != static_cast<Real>(0.0))
295  {
296  if (rule_data[p][2] != static_cast<Real>(0.0))
297  pointtype = 6;
298  else
299  pointtype = 3;
300  }
301 
302  switch (pointtype)
303  {
304  case 1:
305  {
306  // Be sure we have enough space to insert this point
307  libmesh_assert_less (offset + 0, _points.size());
308 
309  // The point has only a single permutation (the centroid!)
310  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][0]);
311 
312  // The weight is always the last entry in the row.
313  _weights[offset + 0] = rule_data[p][3];
314 
315  offset += 1;
316  break;
317  }
318 
319  case 3:
320  {
321  // Be sure we have enough space to insert these points
322  libmesh_assert_less (offset + 2, _points.size());
323 
324  // Here it's understood the second entry is to be used twice, and
325  // thus there are three possible permutations.
326  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
327  _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]);
328  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]);
329 
330  for (unsigned int j=0; j<3; ++j)
331  _weights[offset + j] = rule_data[p][3];
332 
333  offset += 3;
334  break;
335  }
336 
337  case 6:
338  {
339  // Be sure we have enough space to insert these points
340  libmesh_assert_less (offset + 5, _points.size());
341 
342  // Three individual entries with six permutations.
343  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
344  _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]);
345  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]);
346  _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]);
347  _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]);
348  _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]);
349 
350  for (unsigned int j=0; j<6; ++j)
351  _weights[offset + j] = rule_data[p][3];
352 
353  offset += 6;
354  break;
355  }
356 
357  default:
358  {
359  libMesh::err << "Don't know what to do with this many permutation points!" << std::endl;
360  libmesh_error();
361  }
362  }
363  }
364 }
365 
366 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo