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 
49 namespace Intrepid {
50 
51 template <class Scalar, class ArrayPoint, class ArrayWeight>
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 
83 template <class Scalar, class ArrayPoint, class ArrayWeight>
85  Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature2) {
86  cubatures_.resize(2);
87  cubatures_[0] = cubature1;
88  cubatures_[1] = cubature2;
89 
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 
101 template <class Scalar, class ArrayPoint, class ArrayWeight>
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];
114  }
115 
116  dimension_ = cubatures_[0]->getDimension() + cubatures_[1]->getDimension() + cubatures_[2]->getDimension();
117 }
118 
119 
120 
121 template <class Scalar, class ArrayPoint, class ArrayWeight>
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 
137 template <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 
197 template<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 
207 template <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 
218 template <class Scalar, class ArrayPoint, class ArrayWeight>
220  return dimension_;
221 } // end dimension
222 
223 
224 
225 template <class Scalar, class ArrayPoint, class ArrayWeight>
226 void CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int> & degree) const {
227  degree = degree_;
228 } // end getAccuracy
229 
230 } // end namespace Intrepid
virtual int getDimension() const
Returns dimension of integration domain.
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 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.
Defines the base class for cubature (integration) rules in Intrepid.
Defines direct cubature (integration) rules in Intrepid.