Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureSparse.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 #ifndef INTREPID_CUBATURE_SPARSE_HPP
00050 #define INTREPID_CUBATURE_SPARSE_HPP
00051 
00052 #include "Intrepid_ConfigDefs.hpp"
00053 #include "Intrepid_Cubature.hpp"
00054 #include "Intrepid_CubatureDirectLineGauss.hpp"
00055 #include "Intrepid_CubatureSparseHelper.hpp"
00056 #include "Teuchos_Assert.hpp"
00057 
00058 
00063 #define INTREPID_CUBATURE_SPARSE2D_GAUSS_MAX 59
00064 
00069 #define INTREPID_CUBATURE_SPARSE3D_GAUSS_MAX 57
00070 
00071 
00072 namespace Intrepid{
00073 
00074 template<class Scalar, int dimension_, class ArrayPoint = FieldContainer<Scalar>, class ArrayWeight = ArrayPoint>
00075 class CubatureSparse : public Intrepid::Cubature<Scalar,ArrayPoint,ArrayWeight> {
00076   private:
00077 
00078   int level_;
00079 
00080   int numPoints_;
00081 
00082   const int degree_;
00083 
00084   
00085   public:
00086 
00087   ~CubatureSparse() {}
00088 
00089 
00090   CubatureSparse(const int degree);
00091 
00098   virtual void getCubature(ArrayPoint  & cubPoints,
00099                            ArrayWeight & cubWeights) const;
00100 
00103   virtual int getNumPoints() const;
00104 
00107   virtual int getDimension() const;
00108 
00112   virtual void getAccuracy(std::vector<int> & accuracy) const;
00113 
00114 }; // end class CubatureSparse 
00115 
00116 // helper functions
00117 
00118 template<class Scalar, int DIM>
00119 void iterateThroughDimensions(int level,
00120                               int dims_left,
00121                               SGNodes<Scalar,DIM> & cubPointsND,
00122                               Teuchos::Array<Scalar> & partial_node,
00123                               Scalar partial_weight);
00124 
00125 inline int factorial(int num)
00126 {
00127   int answer = 1;
00128   if(num >= 1)
00129   {
00130     while(num > 0)
00131     {
00132       answer = answer*num;
00133       num--;
00134     }
00135   }
00136   else if(num == 0)
00137     answer = 1;
00138   else
00139     answer = -1;
00140 
00141   return answer;
00142 }
00143 
00144 inline double combination(int top, int bot)
00145 {
00146   double answer = factorial(top)/(factorial(bot) * factorial(top-bot));
00147   return answer;
00148 }
00149 
00150 inline int iterateThroughDimensionsForNumCalc(int dims_left,
00151                                               int level,
00152                                               int levels_left,
00153                                               int level_so_far,
00154                                               Teuchos::Array<int> & nodes,
00155                                               int product,
00156                                               bool no_uni_quad)
00157 {
00158   int numNodes = 0;
00159   for(int j = 1; j <= levels_left; j++)
00160   {
00161     bool temp_bool = no_uni_quad;
00162     int temp_knots = nodes[j-1]*product;
00163     int temp_lsf = level_so_far + j;
00164 
00165     if(j==1)
00166       temp_bool = false;
00167 
00168     if(dims_left == 1)
00169     {
00170       if(temp_lsf < level && temp_bool == true)
00171         numNodes += 0;
00172       else
00173       {
00174         numNodes += temp_knots;
00175       }
00176     }
00177     else
00178     {
00179       numNodes += iterateThroughDimensionsForNumCalc(dims_left-1,level, levels_left-j+1, temp_lsf, nodes, temp_knots, temp_bool);
00180     }
00181   }
00182   return numNodes;
00183 }
00184 
00185 inline int calculateNumPoints(int dim, int level)
00186 {
00187   //int* uninum = new int[level];
00188   Teuchos::Array<int> uninum(level);
00189   uninum[0] = 1;
00190   for(int i = 1; i <= level-1; i++)
00191   {
00192     uninum[i] = 2*i;
00193   }
00194 
00195   int numOfNodes = iterateThroughDimensionsForNumCalc(dim, level, level, 0, uninum, 1, true);
00196   return numOfNodes;
00197 }
00198 
00199 
00200 } // end namespace Intrepid
00201 
00202 
00203 // include templated definitions
00204 #include <Intrepid_CubatureSparseDef.hpp>
00205 
00206 #endif