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
55namespace Intrepid{
56
57/************************************************************************
58** Class Definition for class SGPoint
59** Function: Helper Class with cosntruction of Sparse Grid
60*************************************************************************/
61template<class Scalar, int D>
62class SGPoint{
63public:
64 Scalar coords[D];
65
66 SGPoint();
67 SGPoint(Scalar p[D]);
68 bool operator==(const SGPoint<Scalar, D> & right) const;
69 bool operator<(const SGPoint<Scalar, D> & right) const;
70 bool operator>(const SGPoint<Scalar, D> & right) const;
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************************************************************************/
78template<class Scalar, int D, class ArrayPoint=FieldContainer<Scalar>, class ArrayWeight=ArrayPoint>
79class SGNodes{
80public:
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***************************************************************************/
92template<class Scalar, int D>
94{
95 for(int i = 0; i < D; i++)
96 {
97 coords[i] = 0;
98 }
99}
100
101template<class Scalar, int D>
102SGPoint<Scalar, D>::SGPoint(Scalar p[D])
103{
104 for(int i = 0; i < D; i++)
105 {
106 coords[i] = p[i];
107 }
108}
109
110template<class Scalar, int D>
111bool SGPoint<Scalar, D>::operator==(const SGPoint<Scalar, D> & right) const
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
124template<class Scalar, int D>
125bool SGPoint<Scalar, D>::operator<(const SGPoint<Scalar, D> & right) const
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
138template<class Scalar, int D>
139bool SGPoint<Scalar, D>::operator>(const SGPoint<Scalar, D> & right) const
140{
141 if(this < right || this == right)
142 return false;
143
144 return true;
145}
146
147template<class Scalar, int D>
148std::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
162template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
163bool 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
229template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
230void 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/*
244template< class Scalar, int D>
245void 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.