Intrepid
Intrepid_CubaturePolylibDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49namespace Intrepid {
50
51template <class Scalar, class ArrayPoint, class ArrayWeight>
52CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::CubaturePolylib(int degree, EIntrepidPLPoly poly_type, Scalar alpha, Scalar beta) {
53 TEUCHOS_TEST_FOR_EXCEPTION((degree < 0),
54 std::out_of_range,
55 ">>> ERROR (CubaturePolylib): No cubature rule implemented for the desired polynomial degree.");
56 degree_ = degree;
57 dimension_ = 1;
58 poly_type_ = poly_type;
59 alpha_ = alpha;
60 beta_ = beta;
61} // end constructor
62
63
64
65template <class Scalar, class ArrayPoint, class ArrayWeight>
67 return cubature_name_;
68} // end getName
69
70
71
72template <class Scalar, class ArrayPoint, class ArrayWeight>
74 return dimension_;
75} // end dimension
76
77
78
79template <class Scalar, class ArrayPoint, class ArrayWeight>
81 int np = 0;
82 switch (poly_type_) {
83 case PL_GAUSS:
84 np = (degree_+(int)2)/(int)2;
85 break;
86 case PL_GAUSS_RADAU_LEFT:
87 case PL_GAUSS_RADAU_RIGHT:
88 if (degree_ == 0)
89 np = 2;
90 else
91 np = (degree_+(int)3)/(int)2;
92 break;
93 case PL_GAUSS_LOBATTO:
94 np = (degree_+(int)4)/(int)2;
95 break;
96 default:
97 TEUCHOS_TEST_FOR_EXCEPTION((1),
98 std::invalid_argument,
99 ">>> ERROR (CubaturePolylib): Unknown point type argument.");
100 }
101 return np;
102} // end getNumPoints
103
104
105
106template <class Scalar, class ArrayPoint, class ArrayWeight>
107void CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int> & accuracy) const {
108 accuracy.assign(1, degree_);
109} // end getAccuracy
110
111
112
113template <class Scalar, class ArrayPoint, class ArrayWeight>
114const char* CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::cubature_name_ = "INTREPID_CUBATURE_POLYLIB";
115
116
117
118template <class Scalar, class ArrayPoint, class ArrayWeight>
119void CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint & cubPoints, ArrayWeight & cubWeights) const {
120 int numCubPoints = getNumPoints();
121 int cellDim = getDimension();
122 // check size of cubPoints and cubWeights
123 TEUCHOS_TEST_FOR_EXCEPTION( ( ( (int)cubPoints.size() < numCubPoints*cellDim ) || ( (int)cubWeights.size() < numCubPoints ) ),
124 std::out_of_range,
125 ">>> ERROR (CubatureDirect): Insufficient space allocated for cubature points or weights.");
126
127 // temporary storage
128 FieldContainer<Scalar> z(numCubPoints);
129 FieldContainer<Scalar> w(numCubPoints);
130
131 // run Polylib routines
132 switch (poly_type_) {
133 case PL_GAUSS:
134 IntrepidPolylib::zwgj(&z[0], &w[0], numCubPoints, alpha_, beta_);
135 break;
136 case PL_GAUSS_RADAU_LEFT:
137 IntrepidPolylib::zwgrjm(&z[0], &w[0], numCubPoints, alpha_, beta_);
138 break;
139 case PL_GAUSS_RADAU_RIGHT:
140 IntrepidPolylib::zwgrjp(&z[0], &w[0], numCubPoints, alpha_, beta_);
141 break;
142 case PL_GAUSS_LOBATTO:
143 IntrepidPolylib::zwglj(&z[0], &w[0], numCubPoints, alpha_, beta_);
144 break;
145 default:
146 TEUCHOS_TEST_FOR_EXCEPTION((1),
147 std::invalid_argument,
148 ">>> ERROR (CubaturePolylib): Unknown point type argument.");
149 }
150
151 // fill input arrays
152 for (int pointId = 0; pointId < numCubPoints; pointId++) {
153 for (int dim = 0; dim < cellDim; dim++) {
154 cubPoints(pointId,dim) = z[pointId];
155 }
156 cubWeights(pointId) = w[pointId];
157 }
158} // end getCubature
159
160template <class Scalar, class ArrayPoint, class ArrayWeight>
162 ArrayWeight & cubWeights,
163 ArrayPoint& cellCoords) const
164{
165 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
166 ">>> ERROR (CubaturePolylib): Cubature defined in reference space calling method for physical space cubature.");
167}
168
169
170} // end namespace Intrepid
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
CubaturePolylib(int degree=0, EIntrepidPLPoly pt_type=PL_GAUSS, Scalar alpha=0.0, Scalar beta=0.0)
Constructor.
int getNumPoints() const
Returns the number of cubature points.
const char * getName() const
Returns cubature name.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual int getDimension() const
Returns dimension of integration domain.
static void zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
static void zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
static void zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
static void zwgj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Jacobi zeros and weights.