quadrature_grid_3D.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 
19 
20 // Local includes
22 
23 namespace libMesh
24 {
25 
26 
27 void QGrid::init_3D(const ElemType type_in,
28  unsigned int)
29 {
30 #if LIBMESH_DIM == 3
31 
32  //-----------------------------------------------------------------------
33  // 3D quadrature rules
34 
35  // We ignore p - the grid rule is just for experimentation
36  switch (type_in)
37  {
38  //---------------------------------------------
39  // Hex quadrature rules
40  case HEX8:
41  case HEX20:
42  case HEX27:
43  {
44  // We compute the 3D quadrature rule as a tensor
45  // product of the 1D quadrature rule.
46  QGrid q1D(1,_order);
47  q1D.init(EDGE2);
48 
49  tensor_product_hex( q1D );
50 
51  return;
52  }
53 
54 
55 
56  //---------------------------------------------
57  // Tetrahedral quadrature rules
58  case TET4:
59  case TET10:
60  {
61  _points.resize((_order+1)*(_order+2)*(_order+3)/6);
62  _weights.resize((_order+1)*(_order+2)*(_order+3)/6);
63 
64  unsigned int pt = 0;
65  for (int i = 0; i != _order + 1; ++i)
66  {
67  for (int j = 0; j != _order + 1 - i; ++j)
68  {
69  for (int k = 0; k != _order + 1 - i - j; ++k)
70  {
71  _points[pt](0) = (double)i / (double)_order;
72  _points[pt](1) = (double)j / (double)_order;
73  _points[pt](2) = (double)k / (double)_order;
74  _weights[pt] = 1.0 / (double)(_order+1) /
75  (double)(_order+2) / (double)(_order+3);
76  pt++;
77  }
78  }
79  }
80  return;
81  }
82 
83 
84  // Prism quadrature rules
85  case PRISM6:
86  case PRISM15:
87  case PRISM18:
88  {
89  // We compute the 3D quadrature rule as a tensor
90  // product of the 1D quadrature rule and a 2D
91  // triangle quadrature rule
92 
93  QGrid q1D(1,_order);
94  QGrid q2D(2,_order);
95 
96  // Initialize
97  q1D.init(EDGE2);
98  q2D.init(TRI3);
99 
100  tensor_product_prism(q1D, q2D);
101 
102  return;
103  }
104 
105 
106 
107  //---------------------------------------------
108  // Pyramid
109  case PYRAMID5:
110  case PYRAMID14:
111  {
112  _points.resize((_order+1)*(_order+2)*(_order+3)/6);
113  _weights.resize((_order+1)*(_order+2)*(_order+3)/6);
114 
115  unsigned int pt = 0;
116  for (int k = 0; k != _order + 1; ++k)
117  {
118  for (int i = 0; i != _order + 1 - k; ++i)
119  {
120  for (int j = 0; j != _order + 1 - k; ++j)
121  {
122  _points[pt](0) = 2.0 * (double)i / (double)_order
123  - 1.0 + (double)k / (double)_order;
124  _points[pt](1) = 2.0 * (double)j / (double)_order
125  - 1.0 + (double)k / (double)_order;
126  _points[pt](2) = (double)k / (double)_order;
127  _weights[pt] = 1.0 / (double)(_order+1) /
128  (double)(_order+2) / (double)(_order+3);
129  pt++;
130  }
131  }
132  }
133  return;
134  }
135 
136 
137 
138  //---------------------------------------------
139  // Unsupported type
140  default:
141  {
142  libMesh::err << "ERROR: Unsupported type: " << type_in << std::endl;
143  libmesh_error();
144  }
145  }
146 
147  libmesh_error();
148 
149  return;
150 
151 #endif
152 }
153 
154 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo