Intrepid
http://trilinos.sandia.gov/packages/docs/r10.8/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureTensorSortedDef.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 namespace Intrepid {
00050 
00051 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00052 CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
00053                                                int numPoints, int dimension) {
00054   /*
00055     This constructor initializes the nodes and weights for an Ndim quadrature 
00056     rule and sets the nodes and weights lists to zero.
00057   */
00058   points_.clear(); weights_.clear();
00059   numPoints_ = numPoints;
00060   dimension_ = dimension;
00061 }
00062 
00063 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00064 CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
00065                      CubatureLineSorted<Scalar> & cubLine) {
00066   /*
00067     This constructor takes a 1D rule and casts it as a tensor product rule.
00068   */
00069   dimension_ = 1;
00070   numPoints_ = cubLine.getNumPoints();
00071   degree_.resize(1);
00072   cubLine.getAccuracy(degree_);
00073 
00074   int loc = 0;
00075   std::vector<Scalar> node(1,0.0);
00076   typename std::map<Scalar,int>::iterator it;
00077   points_.clear(); weights_.clear(); 
00078   for (it = cubLine.begin(); it != cubLine.end(); it++) {
00079     node[0] = cubLine.getNode(it);
00080     points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
00081     weights_.push_back(cubLine.getWeight(it->second));
00082     loc++;
00083   }
00084 }
00085 
00086 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00087 CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
00088                  int dimension, std::vector<int> numPoints1D, 
00089                  std::vector<EIntrepidBurkardt> rule1D, bool isNormalized) {
00090   /*
00091     This constructor builds a tensor product rule according to quadInfo myRule.
00092   */  
00093   TEST_FOR_EXCEPTION((dimension!=(int)numPoints1D.size()||
00094                       dimension!=(int)rule1D.size()),std::out_of_range,
00095            ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
00096 
00097   dimension_ = dimension;  
00098   degree_.resize(dimension);
00099   std::vector<int> degree(1,0);
00100   CubatureTensorSorted<Scalar> newRule(0,1);
00101   for (int i=0; i<dimension; i++) {
00102     // Compute 1D rules
00103     CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints1D[i],isNormalized);
00104     rule1.getAccuracy(degree);
00105     degree_[i] = degree[0];
00106     // Build Tensor Rule
00107     newRule = kron_prod<Scalar>(newRule,rule1);
00108   }
00109   numPoints_ = newRule.getNumPoints();
00110   typename std::map<std::vector<Scalar>,int>::iterator it;
00111   points_.clear(); weights_.clear();
00112   int loc = 0;
00113   std::vector<Scalar> node(dimension_,0.0);
00114   for (it=newRule.begin(); it!=newRule.end(); it++) {
00115     node = it->first;
00116     points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
00117     weights_.push_back(newRule.getWeight(node));
00118     loc++;
00119   }
00120 } 
00121 
00122 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00123 CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
00124                         int dimension, std::vector<int> numPoints1D, 
00125                         std::vector<EIntrepidBurkardt> rule1D, 
00126                         std::vector<EIntrepidGrowth> growth1D, 
00127                         bool isNormalized) {
00128   /*
00129     This constructor builds a tensor product rule according to quadInfo myRule.
00130   */  
00131   TEST_FOR_EXCEPTION((dimension!=(int)numPoints1D.size()||
00132                       dimension!=(int)rule1D.size()||
00133                       dimension!=(int)growth1D.size()),std::out_of_range,
00134            ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
00135   dimension_ = dimension;  
00136   degree_.resize(dimension);
00137   std::vector<int> degree(1);
00138   CubatureTensorSorted<Scalar> newRule(0,1);
00139   for (int i=0; i<dimension; i++) {
00140     // Compute 1D rules
00141     int numPoints = growthRule1D(numPoints1D[i],growth1D[i],rule1D[i]);
00142     CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
00143     rule1.getAccuracy(degree);
00144     degree_[i] = degree[0];
00145     // Build Tensor Rule
00146     newRule = kron_prod<Scalar>(newRule,rule1);
00147   }
00148   numPoints_ = newRule.getNumPoints();
00149 
00150   typename std::map<std::vector<Scalar>,int>::iterator it;
00151   points_.clear(); weights_.clear();
00152   int loc = 0;
00153   std::vector<Scalar> node;
00154   for (it=newRule.begin(); it!=newRule.end(); it++) {
00155     node = it->first;
00156     points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
00157     weights_.push_back(newRule.getWeight(node));
00158     loc++;
00159   }
00160 } 
00161 
00162 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00163 CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
00164                          int dimension, int maxNumPoints, 
00165                          std::vector<EIntrepidBurkardt> rule1D, 
00166                          std::vector<EIntrepidGrowth> growth1D, 
00167                          bool isNormalized) {
00168   /*
00169     This constructor builds a tensor product rule according to quadInfo myRule.
00170   */  
00171   TEST_FOR_EXCEPTION((dimension!=(int)rule1D.size()||
00172                       dimension!=(int)growth1D.size()),std::out_of_range,
00173             ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
00174   dimension_ = dimension;
00175   degree_.resize(dimension);
00176   std::vector<int> degree(1);
00177   CubatureTensorSorted<Scalar> newRule(0,1);
00178   for (int i=0; i<dimension; i++) {
00179     // Compute 1D rules   
00180     int numPoints = growthRule1D(maxNumPoints,growth1D[i],rule1D[i]);
00181     CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
00182     rule1.getAccuracy(degree);
00183     degree_[i] = degree[0];
00184     // Build Tensor Rule
00185     newRule = kron_prod<Scalar>(newRule,rule1);
00186   }
00187   numPoints_ = newRule.getNumPoints();
00188  
00189   typename std::map<std::vector<Scalar>,int>::iterator it;
00190   points_.clear(); weights_.clear();
00191   int loc = 0;
00192   std::vector<Scalar> node;
00193   for (it=newRule.begin(); it!=newRule.end(); it++) {
00194     node = it->first;
00195     points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
00196     weights_.push_back(newRule.getWeight(node));
00197     loc++;
00198   }
00199 } 
00200 
00201 /* =========================================================================
00202                      Access Operator - ruleTP
00203    ========================================================================= */
00204 template <class Scalar, class ArrayPoint, class ArrayWeight>
00205 int CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
00206   return numPoints_;
00207 } // end getNumPoints
00208 
00209 template <class Scalar, class ArrayPoint, class ArrayWeight>
00210 void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(
00211                                            std::vector<int> & accuracy) const {
00212   accuracy = degree_;
00213 } // end getAccuracy
00214 
00215 template <class Scalar, class ArrayPoint, class ArrayWeight>
00216 void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(
00217                       ArrayPoint & cubPoints, ArrayWeight & cubWeights) const {
00218 
00219   typename std::map<std::vector<Scalar>,int>::const_iterator it;
00220   for (it=points_.begin(); it!=points_.end();it++) {
00221     for (int j=0; j<dimension_; j++) {
00222       cubPoints(it->second,j)  = it->first[j];
00223     }
00224     cubWeights(it->second) = weights_[it->second];
00225   }
00226 
00227   /*
00228   typename std::map<std::vector<Scalar>,int>::const_iterator it = 
00229     points_.begin();
00230   for (int i=0; i<numPoints_; i++) {
00231     for (int j=0; j<dimension_; j++) {
00232       cubPoints(i,j)  = it->first[j];
00233     }
00234     cubWeights(i) = weights_[it->second];
00235     it++;
00236   }
00237   */
00238 } // end getCubature
00239 
00240 template <class Scalar, class ArrayPoint, class ArrayWeight>
00241 int CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
00242   return dimension_;
00243 } // end getDimension
00244 
00245 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00246 typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::begin() {
00247   return points_.begin();
00248 }
00249 
00250 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00251 typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::end() {
00252   return points_.end();
00253 }
00254 
00255 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00256 void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::insert(
00257                     typename std::map<std::vector<Scalar>,int>::iterator it, 
00258                     std::vector<Scalar> point, 
00259                     Scalar weight) {
00260   points_.insert(it,std::pair<std::vector<Scalar>,int>(point,
00261                                                        (int)points_.size()));
00262   weights_.push_back(weight);
00263   numPoints_++;
00264   return;
00265 }
00266 
00267 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00268 std::vector<Scalar> CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getNode(typename std::map<std::vector<Scalar>,int>::iterator it) {
00269   /*
00270     Access node for ruleTP
00271   */ 
00272   return it->first;
00273 }
00274 
00275 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00276 Scalar CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
00277                                                                    int node) { 
00278   /*
00279     Access weight for ruleTP
00280   */   
00281   return weights_[node];
00282 }
00283 
00284 template <class Scalar, class ArrayPoint, class ArrayWeight> 
00285 Scalar CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
00286                                                    std::vector<Scalar> point) {
00287   //
00288   //  Access weight for ruleTP
00289   //   
00290   return weights_[points_[point]];
00291 }
00292 
00293 template <class Scalar, class ArrayPoint, class ArrayWeight>
00294 void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::update(
00295                              Scalar alpha2, 
00296                              CubatureTensorSorted<Scalar> & cubRule2, 
00297                              Scalar alpha1) {
00298   
00299   // Initialize an iterator on std::map<std::vector<Scalar>,Scalar>
00300   typename std::map<std::vector<Scalar>,int>::iterator it;
00301   
00302   // Temporary Container for updated rule
00303   typename std::map<std::vector<Scalar>,int> newPoints;
00304   std::vector<Scalar> newWeights(0,0.0);
00305   std::vector<Scalar> node(dimension_,0.0);
00306    int loc = 0;
00307 
00308   // Intersection of rule1 and rule2
00309   typename std::map<std::vector<Scalar>,int> inter; 
00310   std::set_intersection(points_.begin(),points_.end(),
00311                         cubRule2.begin(),cubRule2.end(),
00312                         inserter(inter,inter.begin()),inter.value_comp());
00313   for (it=inter.begin(); it!=inter.end(); it++) { 
00314     node = it->first;
00315     newWeights.push_back( alpha1*weights_[it->second]
00316                          +alpha2*cubRule2.getWeight(node));
00317     newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));
00318     //points_.erase(node); cubRule2.erase(node);
00319     loc++;    
00320   }
00321   int isize = inter.size(); 
00322 
00323   // Set Difference rule1 \ rule2
00324   int size = points_.size();
00325   if (isize!=size) {
00326     typename std::map<std::vector<Scalar>,int> diff1; 
00327     std::set_difference(points_.begin(),points_.end(),
00328                         cubRule2.begin(),cubRule2.end(),
00329                         inserter(diff1,diff1.begin()),diff1.value_comp());
00330     for (it=diff1.begin(); it!=diff1.end(); it++) {      
00331       node = it->first;
00332       newWeights.push_back(alpha1*weights_[it->second]);
00333       newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));    
00334       loc++;
00335     }  
00336   }
00337 
00338   // Set Difference rule2 \ rule1
00339   size = cubRule2.getNumPoints();
00340   if (isize!=size) {
00341     typename std::map<std::vector<Scalar>,int> diff2; 
00342     std::set_difference(cubRule2.begin(),cubRule2.end(),
00343                         points_.begin(),points_.end(),
00344                         inserter(diff2,diff2.begin()),diff2.value_comp());
00345     for (it=diff2.begin(); it!=diff2.end(); it++) {      
00346       node = it->first;
00347       newWeights.push_back(alpha2*cubRule2.getWeight(it->second));
00348       newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));  
00349       loc++;
00350     }        
00351   }  
00352  
00353   points_.clear();  points_.insert(newPoints.begin(),newPoints.end());
00354   weights_.clear(); weights_.assign(newWeights.begin(),newWeights.end());
00355   numPoints_ = (int)points_.size(); 
00356 }
00357 
00358 template <class Scalar, class ArrayPoint, class ArrayWeight>
00359 void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::normalize(){
00360   Scalar sum = 0.0;
00361   
00362   typename std::vector<Scalar>::iterator it;
00363   for (it=weights_.begin(); it!=weights_.end(); it++)
00364     sum += *it;
00365 
00366   for (it=weights_.begin(); it!=weights_.end(); it++)
00367     *it /= sum;
00368 }
00369 
00370 
00371 /* ======================================================================
00372                              Kronecker Products
00373    ====================================================================== */ 
00374 template <class Scalar>
00375 CubatureTensorSorted<Scalar> kron_prod(CubatureTensorSorted<Scalar> & rule1, 
00376                                        CubatureLineSorted<Scalar> & rule2) {
00377   /* 
00378     Compute the Kronecker Product of a Tensor Product rule and a 1D rule.
00379   */
00380   int s1   = rule1.getNumPoints();
00381   // int s2   = rule2.getNumPoints();
00382   int Ndim = rule1.getDimension();
00383 
00384   if (s1==0) {
00385     CubatureTensorSorted<Scalar> TPrule(rule2);
00386     return TPrule;
00387   }
00388   else {
00389     // Initialize Arrays Containing Updated Nodes and Weights
00390     CubatureTensorSorted<Scalar> TPrule(0,Ndim+1); 
00391 
00392     Scalar weight = 0.0;
00393     Scalar node2  = 0.0;
00394 
00395     // Perform Kronecker Products
00396     // Compute Kronecker Product of Nodes
00397     typename std::map<std::vector<Scalar>,int>::iterator it = TPrule.begin();
00398     typename std::map<std::vector<Scalar>,int>::iterator it_i;
00399     typename std::map<Scalar,int>::iterator it_j;
00400     for (it_i=rule1.begin(); it_i!=rule1.end(); it_i++) {
00401       for (it_j=rule2.begin(); it_j!=rule2.end(); it_j++) {
00402         std::vector<Scalar> node = rule1.getNode(it_i);
00403         node2   = rule2.getNode(it_j);
00404         weight  = rule1.getWeight(node)*rule2.getWeight(node2);
00405         node.push_back(node2);
00406         TPrule.insert(it,node,weight);
00407         it = TPrule.end();
00408       }
00409     }
00410     return TPrule;
00411   }
00412 }
00413 } // end Intrepid namespace