Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureSparseHelper.hpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //
00005 //                           Intrepid Package
00006 //                 Copyright (2007) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00039 //                    Denis Ridzal  (dridzal@sandia.gov), or
00040 //                    Kara Peterson (kjpeter@sandia.gov)
00041 //
00042 // ************************************************************************
00043 // @HEADER
00044 */
00045 
00046 #ifndef INTREPID_CUBATURE_SPARSE_HELPER_HPP
00047 #define INTREPID_CUBATURE_SPARSE_HELPER_HPP
00048 
00049 #include "Intrepid_ConfigDefs.hpp"
00050 #include "Intrepid_Types.hpp"
00051 #include "Intrepid_FieldContainer.hpp"
00052 #include "Teuchos_Assert.hpp"
00053 #include "Teuchos_Array.hpp"
00054 
00055 namespace Intrepid{
00056 
00057 /************************************************************************
00058 **  Class Definition for class SGPoint
00059 **  Function: Helper Class with cosntruction of Sparse Grid
00060 *************************************************************************/
00061 template<class Scalar, int D>
00062 class SGPoint{
00063 public:
00064   Scalar coords[D];
00065 
00066   SGPoint();
00067   SGPoint(Scalar p[D]);
00068   bool const operator==(const SGPoint<Scalar, D> & right);
00069   bool const operator<(const SGPoint<Scalar, D> & right);
00070   bool const operator>(const SGPoint<Scalar, D> & right);
00071   //friend ostream & operator<<(ostream & o, const SGPoint<D> & p);
00072 };
00073 
00074 /************************************************************************
00075 **  Class Definition for class SGNodes
00076 **  function: Helper Class with constrution of Sparse Grid
00077 ************************************************************************/
00078 template<class Scalar, int D, class ArrayPoint=FieldContainer<Scalar>, class ArrayWeight=ArrayPoint>
00079 class SGNodes{
00080 public:
00081   Teuchos::Array< SGPoint<Scalar, D> > nodes;
00082   Teuchos::Array< Scalar > weights;
00083   bool addNode(Scalar new_node[D], Scalar weight);
00084   void copyToArrays(ArrayPoint & cubPoints, ArrayWeight & cubWeights) const;
00085   //void copyToTeuchos(Teuchos::Array< Scalar* > & cubPoints, Teuchos::Array<Scalar> & cubWeights) const;
00086   int size() const {return nodes.size();}
00087 };
00088 
00089 /**************************************************************************
00090 **  Function Definitions for Class SGPoint
00091 ***************************************************************************/
00092 template<class Scalar, int D>
00093 SGPoint<Scalar,D>::SGPoint()
00094 {
00095   for(int i = 0; i < D; i++)
00096   {
00097     coords[i] = 0;
00098   }
00099 }
00100 
00101 template<class Scalar, int D>
00102 SGPoint<Scalar, D>::SGPoint(Scalar p[D])
00103 {
00104   for(int i = 0; i < D; i++)
00105   {
00106     coords[i] = p[i];
00107   }
00108 }
00109 
00110 template<class Scalar, int D>
00111 bool const SGPoint<Scalar, D>::operator==(const SGPoint<Scalar, D> & right)
00112 {
00113   bool equal = true;
00114 
00115   for(int i = 0; i < D; i++)
00116   {
00117     if(coords[i] != right.coords[i])
00118       return false;
00119   }
00120 
00121   return equal;
00122 }
00123 
00124 template<class Scalar, int D>
00125 bool const SGPoint<Scalar, D>::operator<(const SGPoint<Scalar, D> & right)
00126 {
00127   for(int i = 0; i < D; i++)
00128   {
00129     if(coords[i] < right.coords[i])
00130       return true;
00131     else if(coords[i] > right.coords[i])
00132       return false;
00133   }
00134   
00135   return false;
00136 }
00137 
00138 template<class Scalar, int D>
00139 bool const SGPoint<Scalar, D>::operator>(const SGPoint<Scalar, D> & right)
00140 {
00141   if(this < right || this == right)
00142     return false;
00143 
00144   return true;
00145 }
00146 
00147 template<class Scalar, int D>
00148 std::ostream & operator<<(std::ostream & o, SGPoint<Scalar, D> & p)
00149 {
00150   o << "(";
00151   for(int i = 0; i<D;i++)
00152     o<< p.coords[i] << " ";
00153   o << ")";
00154   return o;
00155 }
00156 
00157 
00158 /**************************************************************************
00159 **  Function Definitions for Class SGNodes
00160 ***************************************************************************/
00161 
00162 template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
00163 bool SGNodes<Scalar,D,ArrayPoint,ArrayWeight>::addNode(Scalar new_node[D], Scalar weight)
00164 {
00165   SGPoint<Scalar, D> new_point(new_node);
00166   bool new_and_added = true;
00167 
00168   if(nodes.size() == 0)
00169   {
00170     nodes.push_back(new_point);
00171     weights.push_back(weight);
00172   }
00173   else
00174   {   
00175     int left = -1;
00176     int right = (int)nodes.size();
00177     int mid_node = (int)ceil(nodes.size()/2.0)-1;
00178 
00179     bool iterate_continue = true;
00180 
00181     while(iterate_continue)
00182     {
00183       if(new_point == nodes[mid_node]){
00184         weights[mid_node] += weight;
00185         iterate_continue = false;
00186         new_and_added = false;  
00187       }
00188       else if(new_point < nodes[mid_node]){
00189         if(right - left <= 2)
00190         {
00191           //insert the node into the vector to the left of mid_node
00192           nodes.insert(nodes.begin()+mid_node, new_point);
00193           weights.insert(weights.begin()+mid_node,weight);
00194           iterate_continue = false;
00195         }
00196         else 
00197         {
00198           right = mid_node;
00199           mid_node += (int)ceil((left-mid_node)/2.0);
00200         }
00201       }
00202       else{ //new_point > nodes[mid_node];
00203 
00204         if(mid_node == (int)nodes.size()-1)
00205         {
00206           nodes.push_back(new_point);
00207           weights.push_back(weight);
00208           iterate_continue = false;
00209         }
00210         else if(right - left <= 2)
00211         {
00212           //insert the node into the vector to the right of mid_node
00213           nodes.insert(nodes.begin()+mid_node+1, new_point);
00214           weights.insert(weights.begin()+mid_node+1,weight);
00215           iterate_continue = false;
00216         }
00217         else 
00218         {
00219           left = mid_node;
00220           mid_node += (int)ceil((right-mid_node)/2.0);
00221         }
00222       }
00223     }
00224   }
00225 
00226   return new_and_added;
00227 }
00228 
00229 template<class Scalar, int D, class ArrayPoint, class ArrayWeight>
00230 void SGNodes<Scalar,D,ArrayPoint,ArrayWeight>::copyToArrays(ArrayPoint & cubPoints, ArrayWeight & cubWeights) const
00231 {
00232   int numPoints = size();
00233 
00234   for(int i = 0; i < numPoints; i++)
00235   {
00236     for (int j=0; j<D; j++) {
00237       cubPoints(i,j) = nodes[i].coords[j];
00238     }
00239     cubWeights(i) = weights[i];
00240   }
00241 }
00242 
00243 /*
00244 template< class Scalar, int D>
00245 void SGNodes<Scalar,D>::copyToTeuchos(Teuchos::Array< Scalar* > & cubPoints, Teuchos::Array<Scalar> & cubWeights) const
00246 {
00247   int numPoints = size();
00248 
00249   Scalar tempPoint[D];
00250   cubPoints.assign(numPoints,tempPoint);
00251   cubWeights.assign(numPoints, 0.0);
00252 
00253   for(int i = 0; i < numPoints; i++)
00254   {
00255     cubPoints[i] = nodes[i].coords;
00256     cubWeights[i] = weights[i];
00257   }
00258 }
00259 */
00260 
00261 }
00262 #endif