quadrature_grid_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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  const unsigned int np = (_order+1)*(_order+2)*(_order+3)/6;
62  // Master tet has 1x1 triangle base, height 1, so volume = 1/6
63  const Real weight = Real(1)/Real(6)/np;
64  const Real dx = Real(1)/(_order+1);
65  _points.resize(np);
66  _weights.resize(np);
67 
68  unsigned int pt = 0;
69  for (int i = 0; i != _order + 1; ++i)
70  {
71  for (int j = 0; j != _order + 1 - i; ++j)
72  {
73  for (int k = 0; k != _order + 1 - i - j; ++k)
74  {
75  _points[pt](0) = (i+0.5)*dx;
76  _points[pt](1) = (j+0.5)*dx;
77  _points[pt](2) = (k+0.5)*dx;
78  _weights[pt] = weight;
79  pt++;
80  }
81  }
82  }
83  return;
84  }
85 
86 
87  // Prism quadrature rules
88  case PRISM6:
89  case PRISM15:
90  case PRISM18:
91  {
92  // We compute the 3D quadrature rule as a tensor
93  // product of the 1D quadrature rule and a 2D
94  // triangle quadrature rule
95 
96  QGrid q1D(1,_order);
97  QGrid q2D(2,_order);
98 
99  // Initialize
100  q1D.init(EDGE2);
101  q2D.init(TRI3);
102 
103  tensor_product_prism(q1D, q2D);
104 
105  return;
106  }
107 
108 
109 
110  //---------------------------------------------
111  // Pyramid
112  case PYRAMID5:
113  case PYRAMID13:
114  case PYRAMID14:
115  {
116  const unsigned int np = (_order+1)*(_order+2)*(_order+3)/6;
117  _points.resize(np);
118  _weights.resize(np);
119  // Master pyramid has 2x2 base, height 1, so volume = 4/3
120  const Real weight = Real(4)/Real(3)/np;
121  const Real dx = Real(2)/(_order+1);
122  const Real dz = Real(1)/(_order+1);
123 
124  unsigned int pt = 0;
125  for (int k = 0; k != _order + 1; ++k)
126  {
127  for (int i = 0; i != _order + 1 - k; ++i)
128  {
129  for (int j = 0; j != _order + 1 - k; ++j)
130  {
131  _points[pt](0) = (i+0.5)*dx-1.0 +
132  (k+0.5)*dz;
133  _points[pt](1) = (j+0.5)*dx-1.0 +
134  (k+0.5)*dz;
135  _points[pt](2) = (k+0.5)*dz;
136  _weights[pt] = weight;
137  pt++;
138  }
139  }
140  }
141  return;
142  }
143 
144 
145 
146  //---------------------------------------------
147  // Unsupported type
148  default:
149  libmesh_error_msg("ERROR: Unsupported type: " << type_in);
150  }
151 #endif
152 }
153 
154 } // namespace libMesh
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28
std::vector< Point > _points
Definition: quadrature.h:349
std::vector< Real > _weights
Definition: quadrature.h:355
void tensor_product_prism(const QBase &q1D, const QBase &q2D)
Definition: quadrature.C:181
virtual void init_3D(const ElemType _type=INVALID_ELEM, unsigned int p_level=0) override
void tensor_product_hex(const QBase &q1D)
Definition: quadrature.C:154
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:233
Implements grid-based quadrature rules suitable for non-smooth functions.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real