Intrepid
Intrepid_CubatureTensorDef.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>
52CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(std::vector< Teuchos::RCP<Cubature<Scalar,ArrayPoint,ArrayWeight> > > cubatures) {
53 unsigned numCubs = cubatures.size();
54 TEUCHOS_TEST_FOR_EXCEPTION( (numCubs < 1),
55 std::out_of_range,
56 ">>> ERROR (CubatureTensor): Input cubature array must be of size 1 or larger.");
57
58 cubatures_ = cubatures;
59
60 unsigned numDegrees = 0;
61 for (unsigned i=0; i<numCubs; i++) {
62 std::vector<int> tmp;
63 cubatures[i]->getAccuracy(tmp);
64 numDegrees += tmp.size();
65 }
66
67 degree_.assign(numDegrees, 0);
68 int count = 0;
69 dimension_ = 0;
70 for (unsigned i=0; i<numCubs; i++) {
71 std::vector<int> tmp;
72 cubatures[i]->getAccuracy(tmp);
73 for (unsigned j=0; j<tmp.size(); j++) {
74 degree_[count] = tmp[j];
75 count++;
76 }
77 dimension_ += cubatures[i]->getDimension();
78 }
79}
80
81
82
83template <class Scalar, class ArrayPoint, class ArrayWeight>
84CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature1,
85 Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature2) {
86 cubatures_.resize(2);
87 cubatures_[0] = cubature1;
88 cubatures_[1] = cubature2;
90 degree_.assign(2, 0);
91 for (unsigned i=0; i<2; i++){
92 std::vector<int> d;
93 cubatures_[i]->getAccuracy(d); degree_[i] = d[0];
94 }
95
96 dimension_ = cubatures_[0]->getDimension() + cubatures_[1]->getDimension();
97}
98
99
100
101template <class Scalar, class ArrayPoint, class ArrayWeight>
102CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature1,
103 Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature2,
104 Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature3) {
105 cubatures_.resize(3);
106 cubatures_[0] = cubature1;
107 cubatures_[1] = cubature2;
108 cubatures_[2] = cubature3;
109
110 degree_.assign(3, 0);
111 for (unsigned i=0; i<3; i++){
112 std::vector<int> d;
113 cubatures_[i]->getAccuracy(d); degree_[i] = d[0];
115
116 dimension_ = cubatures_[0]->getDimension() + cubatures_[1]->getDimension() + cubatures_[2]->getDimension();
117}
118
119
120
121template <class Scalar, class ArrayPoint, class ArrayWeight>
122CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature, int n) {
123 cubatures_.resize(n);
124 for (int i=0; i<n; i++) {
125 cubatures_[i] = cubature;
126 }
127
128 std::vector<int> d;
129 cubatures_[0]->getAccuracy(d);
130 degree_.assign(n,d[0]);
131
132 dimension_ = cubatures_[0]->getDimension()*n;
133}
134
135
136
137template <class Scalar, class ArrayPoint, class ArrayWeight>
139 ArrayWeight & cubWeights) const {
140 int numCubPoints = getNumPoints();
141 int cubDim = getDimension();
142 // check size of cubPoints and cubWeights
143 TEUCHOS_TEST_FOR_EXCEPTION( ( ( (int)cubPoints.size() < numCubPoints*cubDim ) || ( (int)cubWeights.size() < numCubPoints ) ),
144 std::out_of_range,
145 ">>> ERROR (CubatureTensor): Insufficient space allocated for cubature points or weights.");
146
147 unsigned numCubs = cubatures_.size();
148 std::vector<unsigned> numLocPoints(numCubs);
149 std::vector<unsigned> locDim(numCubs);
150 std::vector< FieldContainer<Scalar> > points(numCubs);
151 std::vector< FieldContainer<Scalar> > weights(numCubs);
152
153 // extract required points and weights
154 for (unsigned i=0; i<numCubs; i++) {
155
156 numLocPoints[i] = cubatures_[i]->getNumPoints();
157 locDim[i] = cubatures_[i]->getDimension();
158 points[i].resize(numLocPoints[i], locDim[i]);
159 weights[i].resize(numLocPoints[i]);
160
161 // cubPoints and cubWeights are used here only for temporary data retrieval
162 cubatures_[i]->getCubature(cubPoints, cubWeights);
163 for (unsigned pt=0; pt<numLocPoints[i]; pt++) {
164 for (unsigned d=0; d<locDim[i]; d++) {
165 points[i](pt,d) = cubPoints(pt,d);
166 weights[i](pt) = cubWeights(pt);
167 }
168 }
169
170 }
171
172 // reset all weights to 1.0
173 for (int i=0; i<numCubPoints; i++) {
174 cubWeights(i) = (Scalar)1.0;
175 }
176
177 // fill tensor-product cubature
178 int globDimCounter = 0;
179 int shift = 1;
180 for (unsigned i=0; i<numCubs; i++) {
181
182 for (int j=0; j<numCubPoints; j++) {
183 /* int itmp = ((j*shift) % numCubPoints) + (j / (numCubPoints/shift)); // equivalent, but numerically unstable */
184 int itmp = (j % (numCubPoints/shift))*shift + (j / (numCubPoints/shift));
185 for (unsigned k=0; k<locDim[i]; k++) {
186 cubPoints(itmp , globDimCounter+k) = points[i](j % numLocPoints[i], k);
187 }
188 cubWeights( itmp ) *= weights[i](j % numLocPoints[i]);
189 }
190
191 shift *= numLocPoints[i];
192 globDimCounter += locDim[i];
193 }
194
195} // end getCubature
196
197template<class Scalar, class ArrayPoint, class ArrayWeight>
199 ArrayWeight& cubWeights,
200 ArrayPoint& cellCoords) const
201{
202 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
203 ">>> ERROR (CubatureTensor): Cubature defined in reference space calling method for physical space cubature.");
204}
205
206
207template <class Scalar, class ArrayPoint, class ArrayWeight>
209 unsigned numCubs = cubatures_.size();
210 int numCubPoints = 1;
211 for (unsigned i=0; i<numCubs; i++) {
212 numCubPoints *= cubatures_[i]->getNumPoints();
213 }
214 return numCubPoints;
215} // end getNumPoints
216
217
218template <class Scalar, class ArrayPoint, class ArrayWeight>
220 return dimension_;
221} // end dimension
222
223
224
225template <class Scalar, class ArrayPoint, class ArrayWeight>
227 degree = degree_;
228} // end getAccuracy
229
230} // end namespace Intrepid
virtual int getNumPoints() const
Returns the number of cubature points.
virtual void getAccuracy(std::vector< int > &degree) const
Returns max. degree of polynomials that are integrated exactly. The return vector has the size of the...
virtual int getDimension() const
Returns dimension of integration domain.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
CubatureTensor(std::vector< Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > > cubatures)
Constructor.