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 
49 namespace Intrepid {
50 
51 template <class Scalar, class ArrayPoint, class ArrayWeight>
52 CubaturePolylib<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 
65 template <class Scalar, class ArrayPoint, class ArrayWeight>
67  return cubature_name_;
68 } // end getName
69 
70 
71 
72 template <class Scalar, class ArrayPoint, class ArrayWeight>
74  return dimension_;
75 } // end dimension
76 
77 
78 
79 template <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 
106 template <class Scalar, class ArrayPoint, class ArrayWeight>
107 void CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int> & accuracy) const {
108  accuracy.assign(1, degree_);
109 } // end getAccuracy
110 
111 
112 
113 template <class Scalar, class ArrayPoint, class ArrayWeight>
114 const char* CubaturePolylib<Scalar,ArrayPoint,ArrayWeight>::cubature_name_ = "INTREPID_CUBATURE_POLYLIB";
115 
116 
117 
118 template <class Scalar, class ArrayPoint, class ArrayWeight>
119 void 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 
160 template <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
171 
172 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
173 #ifdef __GNUC__
174 #warning "The Intrepid package is deprecated"
175 #endif
176 #endif
177 
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
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.
virtual int getDimension() const
Returns dimension of integration domain.
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.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
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.
int getNumPoints() const
Returns the number of cubature points.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
CubaturePolylib(int degree=0, EIntrepidPLPoly pt_type=PL_GAUSS, Scalar alpha=0.0, Scalar beta=0.0)
Constructor.
const char * getName() const
Returns cubature name.
static void zwgj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Jacobi zeros and weights.