Intrepid
Intrepid_CubatureSparseDef.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 
50 namespace Intrepid {
51 
52 /**************************************************************************
53 ** Function Definitions for Class CubatureSparse
54 ***************************************************************************/
55 
56 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
57 CubatureSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::CubatureSparse(const int degree) :
58  degree_(degree) {
59 
60  if(dimension_ == 2)
61  {
62  if(degree == 1)
63  level_ = 1;
64  else if(degree <= 3)
65  level_ = 2;
66  else if(degree <= 7)
67  level_ = 3;
68  else if(degree <= 11)
69  level_ = 4;
70  else if(degree <= 15)
71  level_ = 5;
72  else if(degree <= 19)
73  level_ = 6;
74  else if(degree <= 23)
75  level_ = 7;
76  else if(degree <= 27)
77  level_ = 8;
78  else if(degree <= 31)
79  level_ = 9;
80  else if(degree <= 35)
81  level_ = 10;
82  else if(degree <= 39)
83  level_ = 11;
84  else if(degree <= 43)
85  level_ = 12;
86  else if(degree <= 47)
87  level_ = 13;
88  else if(degree <= 51)
89  level_ = 14;
90  else if(degree <= 55)
91  level_ = 15;
92  else if(degree <= 59)
93  level_ = 16;
94  }
95  else if(dimension_ == 3)
96  {
97  if(degree == 1)
98  level_ = 1;
99  else if(degree <= 3)
100  level_ = 2;
101  else if(degree <= 5)
102  level_ = 3;
103  else if(degree <= 9)
104  level_ = 4;
105  else if(degree <= 13)
106  level_ = 5;
107  else if(degree <= 17)
108  level_ = 6;
109  else if(degree <= 21)
110  level_ = 7;
111  else if(degree <= 25)
112  level_ = 8;
113  else if(degree <= 29)
114  level_ = 9;
115  else if(degree <= 33)
116  level_ = 10;
117  else if(degree <= 37)
118  level_ = 11;
119  else if(degree <= 41)
120  level_ = 12;
121  else if(degree <= 45)
122  level_ = 13;
123  else if(degree <= 49)
124  level_ = 14;
125  else if(degree <= 53)
126  level_ = 15;
127  else if(degree <= 57)
128  level_ = 16;
129  }
130 
131  numPoints_ = calculateNumPoints(dimension_,level_);
132 }
133 
134 
135 
136 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
138  ArrayWeight & cubWeights) const{
139  Teuchos::Array<Scalar> dummy_point(1);
140  dummy_point[0] = 0.0;
141  Scalar dummy_weight = 1.0;
143 
144  iterateThroughDimensions(level_, dimension_, grid, dummy_point, dummy_weight);
145 
146  grid.copyToArrays(cubPoints, cubWeights);
147 } // end getCubature
148 
149 template<class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
151  ArrayWeight& cubWeights,
152  ArrayPoint& cellCoords) const
153 {
154  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
155  ">>> ERROR (CubatureSparse): Cubature defined in reference space calling method for physical space cubature.");
156 }
157 
158 
159 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
161  return numPoints_;
162 } // end getNumPoints
163 
164 
165 
166 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
168  return dimension_;
169 } // end dimension
170 
171 
172 
173 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
175  accuracy.assign(1, degree_);
176 } //end getAccuracy
177 
178 
179 
180 /************************************************************************
181 ** Function Definition for iterateThroughDimensions()
182 ** and its helper functions
183 *************************************************************************/
184 
185 template< class Scalar, int DIM>
186 void iterateThroughDimensions(int level,
187  int dims_left,
188  SGNodes<Scalar,DIM> & cubPointsND,
189  Teuchos::Array<Scalar> & partial_node,
190  Scalar partial_weight)
191 {
192  int l = level;
193  int d = DIM;
194  int add_on = d - dims_left;
195  int start = dims_left > 1 ? 1 : (int)std::max(1, l);
196  int end = l + add_on;
197 
198  for(int k_i = start; k_i <= end; k_i++)
199  {
200  /*******************
201  ** Slow-Gauss
202  ********************/
203  int order1D = 2*k_i-1;
204 
205  /*******************
206  ** Fast-Gauss
207  ********************/
208  //int order1D = (int)pow(2,k_i) - 1;
209 
210  int cubDegree1D = 2*order1D - 1;
211  CubatureDirectLineGauss<Scalar> Cub1D(cubDegree1D);
212  FieldContainer<Scalar> cubPoints1D(order1D, 1);
213  FieldContainer<Scalar> cubWeights1D(order1D);
214 
215  Cub1D.getCubature(cubPoints1D, cubWeights1D);
216 
217  for(int node1D = 0; node1D < order1D; node1D++)
218  {
219  Teuchos::Array<Scalar> node(d-dims_left+1);
220  node[d-dims_left] = cubPoints1D(node1D,0);
221  for(int m = 0; m < d-dims_left; m++)
222  node[m] = partial_node[m];
223 
224  Scalar weight = cubWeights1D(node1D)*partial_weight;
225 
226  if(dims_left != 1)
227  {
228  iterateThroughDimensions(l - k_i, dims_left-1, cubPointsND, node, weight);
229  }
230  else
231  {
232  weight = pow(-1.0, end - k_i)*combination(d-1, k_i - l)*weight;
233  cubPointsND.addNode(&node[0], weight);
234  }
235  }
236  }
237 }
238 
239 } // end namespace Intrepid
Defines Gauss integration rules on a line.
virtual void getAccuracy(std::vector< int > &accuracy) const
Returns algebraic accuracy (e.g. max. degree of polynomial that is integrated exactly).
virtual int getNumPoints() const
Returns the number of cubature points.
virtual int getDimension() const
Returns dimension of the integration domain.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).