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
50namespace Intrepid {
51
52/**************************************************************************
53** Function Definitions for Class CubatureGenSparse
54***************************************************************************/
55
56template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
57CubatureGenSparse<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
220template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
222 ArrayWeight & cubWeights) const{
223 grid.copyToArrays(cubPoints, cubWeights);
224} // end getCubature
225
226template<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
236template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
238 return numPoints_;
239} // end getNumPoints
240
241
242
243template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
245 return dimension_;
246} // end dimension
247
248
249
250template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
252 accuracy.assign(1, degree_);
253} //end getAccuracy
254
255
256} // end namespace Intrepid