Intrepid
Intrepid_CubatureSparseHelper.hpp
1 /*
2 // @HEADER
3 // ************************************************************************
4 //
5 // Intrepid Package
6 // Copyright (2007) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
39 // Denis Ridzal (dridzal@sandia.gov), or
40 // Kara Peterson (kjpeter@sandia.gov)
41 //
42 // ************************************************************************
43 // @HEADER
44 */
45 
46 #ifndef INTREPID_CUBATURE_SPARSE_HELPER_HPP
47 #define INTREPID_CUBATURE_SPARSE_HELPER_HPP
48 
49 #include "Intrepid_ConfigDefs.hpp"
50 #include "Intrepid_Types.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include "Teuchos_Array.hpp"
54 
55 namespace Intrepid{
56 
57 /************************************************************************
58 ** Class Definition for class SGPoint
59 ** Function: Helper Class with cosntruction of Sparse Grid
60 *************************************************************************/
61 template<class Scalar, int D>
62 class SGPoint{
63 public:
64  Scalar coords[D];
65 
66  SGPoint();
67  SGPoint(Scalar p[D]);
68  bool const operator==(const SGPoint<Scalar, D> & right);
69  bool const operator<(const SGPoint<Scalar, D> & right);
70  bool const operator>(const SGPoint<Scalar, D> & right);
71  //friend ostream & operator<<(ostream & o, const SGPoint<D> & p);
72 };
73 
74 /************************************************************************
75 ** Class Definition for class SGNodes
76 ** function: Helper Class with constrution of Sparse Grid
77 ************************************************************************/
78 template<class Scalar, int D, class ArrayPoint=FieldContainer<Scalar>, class ArrayWeight=ArrayPoint>
79 class SGNodes{
80 public:
81  Teuchos::Array< SGPoint<Scalar, D> > nodes;
82  Teuchos::Array< Scalar > weights;
83  bool addNode(Scalar new_node[D], Scalar weight);
84  void copyToArrays(ArrayPoint & cubPoints, ArrayWeight & cubWeights) const;
85  //void copyToTeuchos(Teuchos::Array< Scalar* > & cubPoints, Teuchos::Array<Scalar> & cubWeights) const;
86  int size() const {return nodes.size();}
87 };
88 
89 /**************************************************************************
90 ** Function Definitions for Class SGPoint
91 ***************************************************************************/
92 template<class Scalar, int D>
94 {
95  for(int i = 0; i < D; i++)
96  {
97  coords[i] = 0;
98  }
99 }
100 
101 template<class Scalar, int D>
102 SGPoint<Scalar, D>::SGPoint(Scalar p[D])
103 {
104  for(int i = 0; i < D; i++)
105  {
106  coords[i] = p[i];
107  }
108 }
109 
110 template<class Scalar, int D>
111 bool const SGPoint<Scalar, D>::operator==(const SGPoint<Scalar, D> & right)
112 {
113  bool equal = true;
114 
115  for(int i = 0; i < D; i++)
116  {
117  if(coords[i] != right.coords[i])
118  return false;
119  }
120 
121  return equal;
122 }
123 
124 template<class Scalar, int D>
125 bool const SGPoint<Scalar, D>::operator<(const SGPoint<Scalar, D> & right)
126 {
127  for(int i = 0; i < D; i++)
128  {
129  if(coords[i] < right.coords[i])
130  return true;
131  else if(coords[i] > right.coords[i])
132  return false;
133  }
134 
135  return false;
136 }
137 
138 template<class Scalar, int D>
139 bool const SGPoint<Scalar, D>::operator>(const SGPoint<Scalar, D> & right)
140 {
141  if(this < right || this == right)
142  return false;
143 
144  return true;
145 }
146 
147 template<class Scalar, int D>
148 std::ostream & operator<<(std::ostream & o, SGPoint<Scalar, D> & p)
149 {
150  o << "(";
151  for(int i = 0; i<D;i++)
152  o<< p.coords[i] << " ";
153  o << ")";
154  return o;
155 }
156 
157 
158 /**************************************************************************
159 ** Function Definitions for Class SGNodes
160 ***************************************************************************/
161 
162 template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
163 bool SGNodes<Scalar,D,ArrayPoint,ArrayWeight>::addNode(Scalar new_node[D], Scalar weight)
164 {
165  SGPoint<Scalar, D> new_point(new_node);
166  bool new_and_added = true;
167 
168  if(nodes.size() == 0)
169  {
170  nodes.push_back(new_point);
171  weights.push_back(weight);
172  }
173  else
174  {
175  int left = -1;
176  int right = (int)nodes.size();
177  int mid_node = (int)ceil(nodes.size()/2.0)-1;
178 
179  bool iterate_continue = true;
180 
181  while(iterate_continue)
182  {
183  if(new_point == nodes[mid_node]){
184  weights[mid_node] += weight;
185  iterate_continue = false;
186  new_and_added = false;
187  }
188  else if(new_point < nodes[mid_node]){
189  if(right - left <= 2)
190  {
191  //insert the node into the vector to the left of mid_node
192  nodes.insert(nodes.begin()+mid_node, new_point);
193  weights.insert(weights.begin()+mid_node,weight);
194  iterate_continue = false;
195  }
196  else
197  {
198  right = mid_node;
199  mid_node += (int)ceil((left-mid_node)/2.0);
200  }
201  }
202  else{ //new_point > nodes[mid_node];
203 
204  if(mid_node == (int)nodes.size()-1)
205  {
206  nodes.push_back(new_point);
207  weights.push_back(weight);
208  iterate_continue = false;
209  }
210  else if(right - left <= 2)
211  {
212  //insert the node into the vector to the right of mid_node
213  nodes.insert(nodes.begin()+mid_node+1, new_point);
214  weights.insert(weights.begin()+mid_node+1,weight);
215  iterate_continue = false;
216  }
217  else
218  {
219  left = mid_node;
220  mid_node += (int)ceil((right-mid_node)/2.0);
221  }
222  }
223  }
224  }
225 
226  return new_and_added;
227 }
228 
229 template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
230 void SGNodes<Scalar,D,ArrayPoint,ArrayWeight>::copyToArrays(ArrayPoint & cubPoints, ArrayWeight & cubWeights) const
231 {
232  int numPoints = size();
233 
234  for(int i = 0; i < numPoints; i++)
235  {
236  for (int j=0; j<D; j++) {
237  cubPoints(i,j) = nodes[i].coords[j];
238  }
239  cubWeights(i) = weights[i];
240  }
241 }
242 
243 /*
244 template< class Scalar, int D>
245 void SGNodes<Scalar,D>::copyToTeuchos(Teuchos::Array< Scalar* > & cubPoints, Teuchos::Array<Scalar> & cubWeights) const
246 {
247  int numPoints = size();
248 
249  Scalar tempPoint[D];
250  cubPoints.assign(numPoints,tempPoint);
251  cubWeights.assign(numPoints, 0.0);
252 
253  for(int i = 0; i < numPoints; i++)
254  {
255  cubPoints[i] = nodes[i].coords;
256  cubWeights[i] = weights[i];
257  }
258 }
259 */
260 
261 }
262 #endif
Header file for utility class to provide multidimensional containers.
Contains definitions of custom data types in Intrepid.