quadrature.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 // C++ includes
20 
21 // Local includes
22 #include "libmesh/quadrature.h"
23 
24 namespace libMesh
25 {
26 
27 void QBase::init(const ElemType t,
28  unsigned int p)
29 {
30  // check to see if we have already
31  // done the work for this quadrature rule
32  if (t == _type && p == _p_level)
33  return;
34  else
35  {
36  _type = t;
37  _p_level = p;
38  }
39 
40 
41 
42  switch(_dim)
43  {
44  case 0:
45  this->init_0D(_type,_p_level);
46 
47  return;
48 
49  case 1:
50  this->init_1D(_type,_p_level);
51 
52  return;
53 
54  case 2:
55  this->init_2D(_type,_p_level);
56 
57  return;
58 
59  case 3:
60  this->init_3D(_type,_p_level);
61 
62  return;
63 
64  default:
65  libmesh_error();
66  }
67 }
68 
69 
70 
72  unsigned int)
73 {
74  _points.resize(1);
75  _weights.resize(1);
76  _points[0] = Point(0.);
77  _weights[0] = 1.0;
78 }
79 
80 
81 
82 void QBase::scale(std::pair<Real, Real> old_range,
83  std::pair<Real, Real> new_range)
84 {
85  // Make sure we are in 1D
86  libmesh_assert_equal_to (_dim, 1);
87 
88  // Make sure that we have sane ranges
89  libmesh_assert_greater (new_range.second, new_range.first);
90  libmesh_assert_greater (old_range.second, old_range.first);
91 
92  // Make sure there are some points
93  libmesh_assert_greater (_points.size(), 0);
94 
95  // We're mapping from old_range -> new_range
96  for (unsigned int i=0; i<_points.size(); i++)
97  {
98  _points[i](0) =
99  (_points[i](0) - old_range.first) *
100  (new_range.second - new_range.first) /
101  (old_range.second - old_range.first) +
102  new_range.first;
103  }
104 
105  // Compute the scale factor and scale the weights
106  const Real scfact = (new_range.second - new_range.first) /
107  (old_range.second - old_range.first);
108 
109  for (unsigned int i=0; i<_points.size(); i++)
110  _weights[i] *= scfact;
111 }
112 
113 
114 
115 
116 void QBase::tensor_product_quad(const QBase& q1D)
117 {
118 
119  const unsigned int np = q1D.n_points();
120 
121  _points.resize(np * np);
122 
123  _weights.resize(np * np);
124 
125  unsigned int q=0;
126 
127  for (unsigned int j=0; j<np; j++)
128  for (unsigned int i=0; i<np; i++)
129  {
130  _points[q](0) = q1D.qp(i)(0);
131  _points[q](1) = q1D.qp(j)(0);
132 
133  _weights[q] = q1D.w(i)*q1D.w(j);
134 
135  q++;
136  }
137 }
138 
139 
140 
141 
142 
143 void QBase::tensor_product_hex(const QBase& q1D)
144 {
145  const unsigned int np = q1D.n_points();
146 
147  _points.resize(np * np * np);
148 
149  _weights.resize(np * np * np);
150 
151  unsigned int q=0;
152 
153  for (unsigned int k=0; k<np; k++)
154  for (unsigned int j=0; j<np; j++)
155  for (unsigned int i=0; i<np; i++)
156  {
157  _points[q](0) = q1D.qp(i)(0);
158  _points[q](1) = q1D.qp(j)(0);
159  _points[q](2) = q1D.qp(k)(0);
160 
161  _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
162 
163  q++;
164  }
165 }
166 
167 
168 
169 
170 void QBase::tensor_product_prism(const QBase& q1D, const QBase& q2D)
171 {
172  const unsigned int n_points1D = q1D.n_points();
173  const unsigned int n_points2D = q2D.n_points();
174 
175  _points.resize (n_points1D * n_points2D);
176  _weights.resize (n_points1D * n_points2D);
177 
178  unsigned int q=0;
179 
180  for (unsigned int j=0; j<n_points1D; j++)
181  for (unsigned int i=0; i<n_points2D; i++)
182  {
183  _points[q](0) = q2D.qp(i)(0);
184  _points[q](1) = q2D.qp(i)(1);
185  _points[q](2) = q1D.qp(j)(0);
186 
187  _weights[q] = q2D.w(i) * q1D.w(j);
188 
189  q++;
190  }
191 
192 }
193 
194 
195 
196 
197 std::ostream& operator << (std::ostream& os, const QBase& q)
198 {
199  q.print_info(os);
200  return os;
201 }
202 
203 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo