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
CubatureTensor(std::vector< Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > > cubatures)
Constructor.
virtual int getNumPoints() const
Returns the number of cubature points.
virtual int getDimension() const
Returns dimension of integration domain.
Defines the base class for cubature (integration) rules in Intrepid.
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...
Defines direct cubature (integration) rules in Intrepid.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).