Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureGenSparseDef.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 
00050 namespace Intrepid {
00051 
00052 /**************************************************************************
00053 **  Function Definitions for Class CubatureGenSparse
00054 ***************************************************************************/
00055 
00056 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
00057 CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::CubatureGenSparse(const int degree) :
00058     degree_(degree) {
00059 
00060   SGNodes<int, dimension_> list;
00061   SGNodes<int,dimension_> bigger_rules;
00062 
00063   bool continue_making_first_list = true;
00064   bool more_bigger_rules = true;
00065 
00066   int poly_exp[dimension_];
00067   int level[dimension_];
00068   int temp_big_rule[dimension_];
00069   
00070   for(int i = 0; i<dimension_; i++){
00071     poly_exp[i] = 0;
00072     temp_big_rule[i] = 0;
00073   }
00074 
00075   while(continue_making_first_list){
00076     for(int i = 0; i < dimension_; i++)
00077     {
00078       int max_exp = 0;
00079       if(i == 0)
00080         max_exp = std::max(degree_,1) - Sum(poly_exp,1,dimension_-1);
00081       else if(i == dimension_ -1)
00082         max_exp = std::max(degree_,1) - Sum(poly_exp,0,dimension_-2);
00083       else
00084         max_exp = std::max(degree_,1) - Sum(poly_exp,0,dimension_-1) + poly_exp[i];
00085 
00086       if(poly_exp[i] < max_exp)
00087       {
00088         poly_exp[i]++;
00089         break;
00090       }
00091       else
00092       {
00093         if(i == dimension_-1)
00094           continue_making_first_list = false;
00095         else
00096           poly_exp[i] = 0;
00097           
00098       }
00099     }
00100 
00101     if(continue_making_first_list)
00102     {
00103       for(int j = 0; j < dimension_;j++)
00104       {
00105         /*******************
00106         **  Slow-Gauss
00107         ********************/
00108         level[j] = (int)std::ceil((((Scalar)poly_exp[j])+3.0)/4.0);
00109         /*******************
00110         **  Fast-Gauss
00111         ********************/
00112         //level[j] = intstd::ceil(std::log(poly_exp[j]+3)/std::log(2) - 1);
00113       }
00114       list.addNode(level,1);
00115       
00116     }
00117   }
00118 
00119 
00120   while(more_bigger_rules)
00121   {
00122     bigger_rules.addNode(temp_big_rule,1);
00123 
00124     for(int i = 0; i < dimension_; i++)
00125     {
00126       if(temp_big_rule[i] == 0){
00127         temp_big_rule[i] = 1;
00128         break;
00129       }
00130       else{
00131         if(i == dimension_-1)
00132           more_bigger_rules = false;
00133         else
00134           temp_big_rule[i] = 0;
00135       }
00136     }
00137   } 
00138 
00139   for(int x = 0; x < list.size(); x++){
00140     for(int y = 0; y < bigger_rules.size(); y++)
00141     { 
00142       SGPoint<int, dimension_> next_rule;
00143       for(int t = 0; t < dimension_; t++)
00144         next_rule.coords[t] = list.nodes[x].coords[t] + bigger_rules.nodes[y].coords[t];
00145 
00146       bool is_in_set = false;
00147       for(int z = 0; z < list.size(); z++)
00148       {
00149         if(next_rule == list.nodes[z]){
00150           is_in_set = true;
00151           break;
00152         }
00153       }
00154 
00155       if(is_in_set)
00156       {
00157         int big_sum[dimension_];
00158         for(int i = 0; i < dimension_; i++)
00159           big_sum[i] = bigger_rules.nodes[y].coords[i];
00160         Scalar coeff = std::pow(-1.0, Sum(big_sum, 0, dimension_-1));
00161         
00162         Scalar point[dimension_];
00163         int point_record[dimension_];
00164 
00165         for(int j = 0; j<dimension_; j++)
00166           point_record[j] = 1;
00167 
00168         bool more_points = true;
00169 
00170         while(more_points)
00171         {
00172           Scalar weight = 1.0;
00173         
00174           for(int w = 0; w < dimension_; w++){
00175             /*******************
00176             **  Slow-Gauss
00177             ********************/
00178             int order1D = 2*list.nodes[x].coords[w]-1;
00179             /*******************
00180             **  Fast-Gauss
00181             ********************/
00182             //int order1D = (int)std::pow(2.0,next_rule.coords[w]) - 1;
00183 
00184             int cubDegree1D = 2*order1D - 1;
00185             CubatureDirectLineGauss<Scalar> Cub1D(cubDegree1D);
00186             FieldContainer<Scalar> cubPoints1D(order1D, 1);
00187             FieldContainer<Scalar> cubWeights1D(order1D);
00188 
00189             Cub1D.getCubature(cubPoints1D, cubWeights1D);
00190 
00191             point[w] =  cubPoints1D(point_record[w]-1, 0);
00192             weight = weight * cubWeights1D(point_record[w]-1);
00193           }     
00194           weight = weight*coeff;
00195           grid.addNode(point, weight);
00196 
00197           for(int v = 0; v < dimension_; v++)
00198           {
00199             if(point_record[v] < 2*list.nodes[x].coords[v]-1){
00200               (point_record[v])++;
00201               break;
00202             }
00203             else{
00204               if(v == dimension_-1)
00205                 more_points = false;
00206               else
00207                 point_record[v] = 1;
00208             }
00209           }
00210         }
00211       }
00212     }
00213   }
00214 
00215   numPoints_ = grid.size();
00216 }
00217 
00218 
00219 
00220 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
00221 void CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint  & cubPoints,
00222                                                                               ArrayWeight & cubWeights) const{
00223   grid.copyToArrays(cubPoints, cubWeights);
00224 } // end getCubature
00225 
00226 
00227 
00228 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
00229 int CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::getNumPoints() const {
00230   return numPoints_;
00231 } // end getNumPoints
00232 
00233 
00234 
00235 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
00236 int CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::getDimension() const {
00237   return dimension_;
00238 } // end dimension
00239 
00240 
00241 
00242 template <class Scalar, int dimension_, class ArrayPoint, class ArrayWeight>
00243 void CubatureGenSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int> & accuracy) const {
00244   accuracy.assign(1, degree_);
00245 } //end getAccuracy
00246 
00247 
00248 } // end namespace Intrepid