quadrature_simpson_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 
28 
29 
30 void QSimpson::init_3D(const ElemType type_in,
31  unsigned int)
32 {
33 #if LIBMESH_DIM == 3
34 
35  //-----------------------------------------------------------------------
36  // 3D quadrature rules
37  switch (type_in)
38  {
39  //---------------------------------------------
40  // Hex quadrature rules
41  case HEX8:
42  case HEX20:
43  case HEX27:
44  {
45  // We compute the 3D quadrature rule as a tensor
46  // product of the 1D quadrature rule.
47  QSimpson q1D(1);
48  q1D.init(EDGE2);
49 
50  tensor_product_hex( q1D );
51 
52  return;
53  }
54 
55 
56 
57  //---------------------------------------------
58  // Tetrahedral quadrature rules
59  case TET4:
60  case TET10:
61  {
62  // This rule is created by combining 8 subtets
63  // which use the trapezoidal rule. The weights
64  // may seem a bit odd, but they are correct,
65  // and should add up to 1/6, the volume of the
66  // reference tet. The points of this rule are
67  // at the nodal points of the TET10, allowing
68  // you to generate diagonal element stiffness
69  // matrices when using quadratic elements.
70  // It should be able to integrate something
71  // better than linears, but I'm not sure how
72  // high.
73 
74  _points.resize(10);
75  _weights.resize(10);
76 
77  _points[0](0) = 0.; _points[5](0) = .5;
78  _points[0](1) = 0.; _points[5](1) = .5;
79  _points[0](2) = 0.; _points[5](2) = 0.;
80 
81  _points[1](0) = 1.; _points[6](0) = 0.;
82  _points[1](1) = 0.; _points[6](1) = .5;
83  _points[1](2) = 0.; _points[6](2) = 0.;
84 
85  _points[2](0) = 0.; _points[7](0) = 0.;
86  _points[2](1) = 1.; _points[7](1) = 0.;
87  _points[2](2) = 0.; _points[7](2) = .5;
88 
89  _points[3](0) = 0.; _points[8](0) = .5;
90  _points[3](1) = 0.; _points[8](1) = 0.;
91  _points[3](2) = 1.; _points[8](2) = .5;
92 
93  _points[4](0) = .5; _points[9](0) = 0.;
94  _points[4](1) = 0.; _points[9](1) = .5;
95  _points[4](2) = 0.; _points[9](2) = .5;
96 
97 
98  _weights[0] = .0052083333333333333333333333333333333333333333; // 1./192.
99  _weights[1] = _weights[0];
100  _weights[2] = _weights[0];
101  _weights[3] = _weights[0];
102 
103  _weights[4] = .0243055555555555555555555555555555555555555555; // 14./576.
104  _weights[5] = _weights[4];
105  _weights[6] = _weights[4];
106  _weights[7] = _weights[4];
107  _weights[8] = _weights[4];
108  _weights[9] = _weights[4];
109 
110  return;
111  }
112 
113 
114 
115  //---------------------------------------------
116  // Prism quadrature rules
117  case PRISM6:
118  case PRISM15:
119  case PRISM18:
120  {
121  // We compute the 3D quadrature rule as a tensor
122  // product of the 1D quadrature rule and a 2D
123  // triangle quadrature rule
124 
125  QSimpson q1D(1);
126  QSimpson q2D(2);
127 
128  // Initialize
129  q1D.init(EDGE2);
130  q2D.init(TRI3);
131 
132  tensor_product_prism(q1D, q2D);
133 
134  return;
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:07 UTC

Hosted By:
SourceForge.net Logo