Intrepid
Intrepid_CubatureGenSparseDef.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 CubatureGenSparse
54 ***************************************************************************/
55 
56 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
57 CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::CubatureGenSparse(const int degree) :
58  degree_(degree) {
59 
60  SGNodes<int, dimension_> list;
61  SGNodes<int,dimension_> bigger_rules;
62 
63  bool continue_making_first_list = true;
64  bool more_bigger_rules = true;
65 
66  int poly_exp[dimension_];
67  int level[dimension_];
68  int temp_big_rule[dimension_];
69 
70  for(int i = 0; i<dimension_; i++){
71  poly_exp[i] = 0;
72  temp_big_rule[i] = 0;
73  }
74 
75  while(continue_making_first_list){
76  for(int i = 0; i < dimension_; i++)
77  {
78  int max_exp = 0;
79  if(i == 0)
80  max_exp = std::max(degree_,1) - Sum(poly_exp,1,dimension_-1);
81  else if(i == dimension_ -1)
82  max_exp = std::max(degree_,1) - Sum(poly_exp,0,dimension_-2);
83  else
84  max_exp = std::max(degree_,1) - Sum(poly_exp,0,dimension_-1) + poly_exp[i];
85 
86  if(poly_exp[i] < max_exp)
87  {
88  poly_exp[i]++;
89  break;
90  }
91  else
92  {
93  if(i == dimension_-1)
94  continue_making_first_list = false;
95  else
96  poly_exp[i] = 0;
97 
98  }
99  }
100 
101  if(continue_making_first_list)
102  {
103  for(int j = 0; j < dimension_;j++)
104  {
105  /*******************
106  ** Slow-Gauss
107  ********************/
108  level[j] = (int)std::ceil((((Scalar)poly_exp[j])+3.0)/4.0);
109  /*******************
110  ** Fast-Gauss
111  ********************/
112  //level[j] = intstd::ceil(std::log(poly_exp[j]+3)/std::log(2) - 1);
113  }
114  list.addNode(level,1);
115 
116  }
117  }
118 
119 
120  while(more_bigger_rules)
121  {
122  bigger_rules.addNode(temp_big_rule,1);
123 
124  for(int i = 0; i < dimension_; i++)
125  {
126  if(temp_big_rule[i] == 0){
127  temp_big_rule[i] = 1;
128  break;
129  }
130  else{
131  if(i == dimension_-1)
132  more_bigger_rules = false;
133  else
134  temp_big_rule[i] = 0;
135  }
136  }
137  }
138 
139  for(int x = 0; x < list.size(); x++){
140  for(int y = 0; y < bigger_rules.size(); y++)
141  {
142  SGPoint<int, dimension_> next_rule;
143  for(int t = 0; t < dimension_; t++)
144  next_rule.coords[t] = list.nodes[x].coords[t] + bigger_rules.nodes[y].coords[t];
145 
146  bool is_in_set = false;
147  for(int z = 0; z < list.size(); z++)
148  {
149  if(next_rule == list.nodes[z]){
150  is_in_set = true;
151  break;
152  }
153  }
154 
155  if(is_in_set)
156  {
157  int big_sum[dimension_];
158  for(int i = 0; i < dimension_; i++)
159  big_sum[i] = bigger_rules.nodes[y].coords[i];
160  Scalar coeff = std::pow(-1.0, Sum(big_sum, 0, dimension_-1));
161 
162  Scalar point[dimension_];
163  int point_record[dimension_];
164 
165  for(int j = 0; j<dimension_; j++)
166  point_record[j] = 1;
167 
168  bool more_points = true;
169 
170  while(more_points)
171  {
172  Scalar weight = 1.0;
173 
174  for(int w = 0; w < dimension_; w++){
175  /*******************
176  ** Slow-Gauss
177  ********************/
178  int order1D = 2*list.nodes[x].coords[w]-1;
179  /*******************
180  ** Fast-Gauss
181  ********************/
182  //int order1D = (int)std::pow(2.0,next_rule.coords[w]) - 1;
183 
184  int cubDegree1D = 2*order1D - 1;
185  CubatureDirectLineGauss<Scalar> Cub1D(cubDegree1D);
186  FieldContainer<Scalar> cubPoints1D(order1D, 1);
187  FieldContainer<Scalar> cubWeights1D(order1D);
188 
189  Cub1D.getCubature(cubPoints1D, cubWeights1D);
190 
191  point[w] = cubPoints1D(point_record[w]-1, 0);
192  weight = weight * cubWeights1D(point_record[w]-1);
193  }
194  weight = weight*coeff;
195  grid.addNode(point, weight);
196 
197  for(int v = 0; v < dimension_; v++)
198  {
199  if(point_record[v] < 2*list.nodes[x].coords[v]-1){
200  (point_record[v])++;
201  break;
202  }
203  else{
204  if(v == dimension_-1)
205  more_points = false;
206  else
207  point_record[v] = 1;
208  }
209  }
210  }
211  }
212  }
213  }
214 
215  numPoints_ = grid.size();
216 }
217 
218 
219 
220 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
222  ArrayWeight & cubWeights) const{
223  grid.copyToArrays(cubPoints, cubWeights);
224 } // end getCubature
225 
226 template<class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
228  ArrayWeight& cubWeights,
229  ArrayPoint& cellCoords) const
230 {
231  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
232  ">>> ERROR (CubatureGenSparse): Cubature defined in reference space calling method for physical space cubature.");
233 }
234 
235 
236 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
238  return numPoints_;
239 } // end getNumPoints
240 
241 
242 
243 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
245  return dimension_;
246 } // end dimension
247 
248 
249 
250 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
252  accuracy.assign(1, degree_);
253 } //end getAccuracy
254 
255 
256 } // end namespace Intrepid
virtual int getNumPoints() const
Returns the number of cubature points.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual void getAccuracy(std::vector< int > &accuracy) const
Returns algebraic accuracy (e.g. max. degree of polynomial that is integrated exactly).
virtual int getDimension() const
Returns dimension of the integration domain.