Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Integration/Intrepid_AdaptiveSparseGridDef.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 UserVector>
00052 bool AdaptiveSparseGrid<Scalar,UserVector>::isAdmissible(
00053                std::vector<int> index, 
00054                int direction, 
00055                std::set<std::vector<int> > inOldIndex,
00056                AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
00057   /*
00058     Check if inOldIndex remains admissible if index is added, i.e.
00059        index-ek in inOldIndex for all k=1,...,dim.
00060   */
00061   int dimension = problem_data.getDimension();
00062   for (int i=0; i<dimension; i++) {
00063     if (index[i]>1 && i!=direction) {
00064       index[i]--;
00065       if (!inOldIndex.count(index)) {
00066         return false;
00067       }
00068       index[i]++;
00069     }
00070   }
00071   return true;
00072 }
00073 
00074 template<class Scalar, class UserVector>
00075 void AdaptiveSparseGrid<Scalar,UserVector>::build_diffRule(
00076                CubatureTensorSorted<Scalar> & outRule, 
00077                std::vector<int> index,
00078                AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
00079                
00080   int numPoints     = 0;
00081   int dimension     = problem_data.getDimension();
00082   bool isNormalized = problem_data.isNormalized();
00083   std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
00084   std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
00085 
00086   for (int i=0; i<dimension; i++) {
00087     // Compute 1D rules
00088     numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
00089     CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
00090  
00091     if (numPoints!=1) { // Compute differential rule
00092       numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
00093       CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
00094       diffRule1.update(-1.0,rule1,1.0);
00095     }
00096     // Build Tensor Rule
00097     outRule = kron_prod<Scalar>(outRule,diffRule1); 
00098   }
00099 }
00100 
00101 template<class Scalar, class UserVector>
00102 void AdaptiveSparseGrid<Scalar,UserVector>::build_diffRule(
00103                CubatureTensorSorted<Scalar> & outRule, 
00104                std::vector<int> index, int dimension,
00105                std::vector<EIntrepidBurkardt> rule1D, 
00106                std::vector<EIntrepidGrowth> growth1D,
00107                bool isNormalized) {
00108                
00109   int numPoints = 0;
00110   for (int i=0; i<dimension; i++) {
00111     // Compute 1D rules
00112     numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
00113     CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
00114  
00115     if (numPoints!=1) { // Differential Rule
00116       numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
00117       CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
00118       diffRule1.update(-1.0,rule1,1.0);
00119     }
00120     // Build Tensor Rule
00121     outRule = kron_prod<Scalar>(outRule,diffRule1); 
00122   }
00123 }
00124  
00125 // Update Index Set - no knowledge of active or old indices
00126 template<class Scalar, class UserVector>
00127 Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
00128                  typename std::multimap<Scalar,std::vector<int> >  & indexSet, 
00129                  UserVector & integralValue,
00130                  AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
00131   
00132   int dimension = problem_data.getDimension();
00133   std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
00134   std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
00135   
00136   // Copy Multimap into a Set for ease of use
00137   typename std::multimap<Scalar,std::vector<int> >::iterator it;
00138   std::set<std::vector<int> > oldSet;
00139   std::set<std::vector<int> >::iterator it1(oldSet.begin());
00140   for (it=indexSet.begin(); it!=indexSet.end(); it++) {
00141     oldSet.insert(it1,it->second);
00142     it1++;
00143   }
00144   indexSet.clear();
00145   
00146   // Find Possible Active Points
00147   int flag = 1;
00148   std::vector<int> index(dimension,0);
00149   typename std::multimap<Scalar,std::vector<int> > activeSet;
00150   for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
00151     index = *it1;
00152     for (int i=0; i<dimension; i++) {
00153       index[i]++;
00154       flag = (int)(!oldSet.count(index));
00155       index[i]--;
00156       if (flag) {
00157         activeSet.insert(std::pair<Scalar,std::vector<int> >(1.0,index));
00158         oldSet.erase(it1);
00159         break;
00160       }
00161     }
00162   }
00163 
00164   // Compute local and global error indicators for active set
00165   typename std::multimap<Scalar,std::vector<int> >::iterator it2;
00166   Scalar eta = 0.0;
00167   Scalar G   = 0.0;
00168   Teuchos::RCP<UserVector> s = integralValue.Create();
00169   for (it2=activeSet.begin(); it2!=activeSet.end(); it2++) {
00170     // Build Differential Quarature Rule
00171     index = it2->second;
00172     CubatureTensorSorted<Scalar> diffRule(0,dimension);
00173     build_diffRule(diffRule,index,problem_data);
00174     
00175     // Apply Rule to function
00176     problem_data.eval_cubature(*s,diffRule);
00177     
00178     // Update local error indicator and index set
00179     G  = problem_data.error_indicator(*s);
00180     activeSet.erase(it2);
00181     activeSet.insert(it2,std::pair<Scalar,std::vector<int> >(G,index));
00182     eta += G;
00183   }
00184 
00185   // Refine Sparse Grid
00186   eta = refine_grid(activeSet,oldSet,integralValue,eta,
00187                     dimension,rule1D,growth1D);
00188 
00189   // Insert New Active and Old Index sets into indexSet
00190   indexSet.insert(activeSet.begin(),activeSet.end());
00191   for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
00192     index = *it1;
00193     indexSet.insert(std::pair<Scalar,std::vector<int> >(-1.0,index));
00194   }
00195     
00196   return eta;
00197 }
00198 
00199 // Update index set and output integral
00200 template<class Scalar, class UserVector> 
00201 Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
00202           typename std::multimap<Scalar,std::vector<int> > & activeIndex, 
00203           std::set<std::vector<int> > & oldIndex, 
00204           UserVector & integralValue,
00205           Scalar globalErrorIndicator,
00206           AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
00207  
00208  TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
00209               ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");  
00210 
00211   int dimension = problem_data.getDimension();
00212   std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
00213   std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
00214 
00215   // Initialize Flags
00216   bool maxLevelFlag     = true;
00217   bool isAdmissibleFlag = true;
00218 
00219   // Initialize Cubature Rule
00220   Teuchos::RCP<UserVector> s = integralValue.Create();
00221   // Initialize iterator at end of inOldIndex
00222   std::set<std::vector<int> >::iterator it1(oldIndex.end()); 
00223 
00224   // Initialize iterator at end of inActiveIndex
00225   typename std::multimap<Scalar,std::vector<int> >::iterator it;
00226 
00227   // Obtain Global Error Indicator as sum of key values of inActiveIndex
00228   Scalar eta = globalErrorIndicator;
00229 
00230   // Select Index to refine
00231   it = activeIndex.end(); it--;        // Decrement to position of final value
00232   Scalar G               = it->first;  // Largest Error Indicator is at End 
00233   eta                   -= G;          // Update global error indicator
00234   std::vector<int> index = it->second; // Get Corresponding index
00235   activeIndex.erase(it);               // Erase Index from active index set
00236   // Insert Index into old index set
00237   oldIndex.insert(it1,index); it1 = oldIndex.end();   
00238 
00239   // Refinement process
00240   for (int k=0; k<dimension; k++) {
00241     index[k]++; // index + ek
00242     // Check Max Level
00243     maxLevelFlag = problem_data.max_level(index);
00244     if (maxLevelFlag) {
00245       // Check Admissibility
00246       isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
00247       if (isAdmissibleFlag) { // If admissible
00248         // Build Differential Quarature Rule
00249         CubatureTensorSorted<Scalar> diffRule(0,dimension);
00250         build_diffRule(diffRule,index,problem_data);    
00251 
00252         // Apply Rule to function
00253         problem_data.eval_cubature(*s,diffRule);
00254         
00255         // Update integral value
00256         integralValue.Update(*s);
00257 
00258         // Update local error indicator and index set
00259         G  = problem_data.error_indicator(*s); 
00260         if (activeIndex.end()!=activeIndex.begin()) 
00261           activeIndex.insert(activeIndex.end()--,
00262                            std::pair<Scalar,std::vector<int> >(G,index));
00263         else
00264           activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
00265 
00266         // Update global error indicators
00267         eta += G;
00268       }
00269     }
00270     else { // Max Level Exceeded 
00271       //std::cout << "Maximum Level Exceeded" << std::endl;
00272     }
00273     index[k]--;
00274   }
00275   return eta;
00276 }
00277 
00278 // Update index set and output integral/sparse grid
00279 template<class Scalar, class UserVector> 
00280 Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
00281         typename std::multimap<Scalar,std::vector<int> > & activeIndex, 
00282         std::set<std::vector<int> > & oldIndex, 
00283         UserVector & integralValue,
00284         CubatureTensorSorted<Scalar> & cubRule,
00285         Scalar globalErrorIndicator,
00286         AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
00287 
00288   TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
00289               ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");  
00290 
00291   int dimension = problem_data.getDimension();
00292   std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
00293   std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
00294 
00295   // Initialize Flags
00296   bool maxLevelFlag     = true;
00297   bool isAdmissibleFlag = true;
00298 
00299   // Initialize Cubature Rule
00300   Teuchos::RCP<UserVector> s = integralValue.Create();
00301 
00302   // Initialize iterator at end of inOldIndex
00303   std::set<std::vector<int> >::iterator it1(oldIndex.end());  
00304 
00305   // Initialize iterator at end of inActiveIndex
00306   typename std::multimap<Scalar,std::vector<int> >::iterator it;
00307 
00308   // Obtain Global Error Indicator as sum of key values of inActiveIndex
00309   Scalar eta = globalErrorIndicator;
00310 
00311   // Select Index to refine
00312   it = activeIndex.end(); it--;        // Decrement to position of final value 
00313   Scalar G               = it->first;  // Largest Error Indicator is at End
00314   eta                   -= G;          // Update global error indicator
00315   std::vector<int> index = it->second; // Get Corresponding index
00316   activeIndex.erase(it);               // Erase Index from active index set
00317   // Insert Index into old index set
00318   oldIndex.insert(it1,index); it1 = oldIndex.end(); 
00319   
00320   // Refinement process
00321   for (int k=0; k<dimension; k++) {
00322     index[k]++; // index + ek
00323     // Check Max Level
00324     maxLevelFlag = problem_data.max_level(index);
00325     if (maxLevelFlag) {
00326       // Check Admissibility
00327       isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
00328       if (isAdmissibleFlag) { // If admissible
00329         // Build Differential Quarature Rule
00330         CubatureTensorSorted<Scalar> diffRule(0,dimension);
00331         build_diffRule(diffRule,index,problem_data);
00332         
00333         // Apply Rule to function
00334         problem_data.eval_cubature(*s,diffRule);
00335         
00336         // Update integral value
00337         integralValue.Update(*s);
00338         
00339         // Update local error indicator and index set
00340         G  = problem_data.error_indicator(*s);  
00341         if (activeIndex.end()!=activeIndex.begin()) 
00342           activeIndex.insert(activeIndex.end()--,
00343                            std::pair<Scalar,std::vector<int> >(G,index));
00344         else
00345           activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
00346                 
00347         // Update global error indicators
00348         eta += G;
00349 
00350         // Update adapted quadrature rule nodes and weights
00351         cubRule.update(1.0,diffRule,1.0);
00352       }
00353     }
00354     else { // Max Level Exceeded 
00355       //std::cout << "Maximum Level Exceeded" << std::endl;
00356     }
00357     index[k]--;
00358   }
00359   return eta;
00360 }
00361 
00362 template<class Scalar, class UserVector>
00363 void AdaptiveSparseGrid<Scalar,UserVector>::buildSparseGrid(
00364               CubatureTensorSorted<Scalar> & output,
00365               int dimension, int maxlevel,
00366               std::vector<EIntrepidBurkardt> rule1D, 
00367               std::vector<EIntrepidGrowth> growth1D,
00368               bool isNormalized) {
00369 
00370   if (dimension == 2) {
00371     std::vector<int> index(dimension,0);
00372     for (int i=0; i<maxlevel; i++) {
00373       for (int j=0; j<maxlevel; j++) {
00374         if(i+j+dimension <= maxlevel+dimension-1) {
00375           index[0] = i+1; index[1] = j+1; 
00376           CubatureTensorSorted<Scalar> diffRule(0,dimension);
00377           build_diffRule(diffRule,index,dimension,rule1D,growth1D,isNormalized);
00378           output.update(1.0,diffRule,1.0);
00379         }
00380       }
00381     } 
00382   }
00383   else if (dimension == 3) {    
00384     std::vector<int> index(dimension,0);
00385     for (int i=0; i<maxlevel; i++) {
00386       for (int j=0; j<maxlevel; j++) {
00387         for (int k=0; k<maxlevel; k++) {
00388           if(i+j+k+dimension <= maxlevel+dimension-1) {
00389             index[0] = i+1; index[1] = j+1; index[2] = k+1;
00390             CubatureTensorSorted<Scalar> diffRule(0,dimension);
00391             build_diffRule(diffRule,index,dimension,rule1D,
00392                           growth1D,isNormalized);
00393             output.update(1.0,diffRule,1.0);
00394           }
00395         }
00396       }
00397     } 
00398   }
00399   else if (dimension == 4) {
00400     std::vector<int> index(dimension,0);
00401     for (int i=0; i<maxlevel; i++) {
00402       for (int j=0; j<maxlevel; j++) {
00403         for (int k=0; k<maxlevel; k++) {
00404           for (int l=0; l<maxlevel; l++) {
00405             if(i+j+k+l+dimension <= maxlevel+dimension-1) {
00406               index[0] = i+1; index[1] = j+1; index[2] = k+1; index[3] = l+1;
00407               CubatureTensorSorted<Scalar> diffRule(0,dimension);
00408               build_diffRule(diffRule,index,dimension,rule1D,
00409                             growth1D,isNormalized);
00410               output.update(1.0,diffRule,1.0);
00411             }
00412           }
00413         }
00414       } 
00415     }
00416   }
00417   else if (dimension == 5) {
00418     std::vector<int> index(dimension,0);
00419     for (int i=0; i<maxlevel; i++) {
00420       for (int j=0; j<maxlevel; j++) {
00421         for (int k=0; k<maxlevel; k++) {
00422           for (int l=0; l<maxlevel; l++) {
00423             for (int m=0; m<maxlevel; m++) {
00424               if(i+j+k+l+m+dimension <= maxlevel+dimension-1) {
00425                 index[0] = i+1; index[1] = j+1; index[2] = k+1; 
00426                 index[3] = l+1; index[4] = m+1;
00427                 CubatureTensorSorted<Scalar> diffRule(0,dimension);
00428                 build_diffRule(diffRule,index,dimension,rule1D,
00429                               growth1D,isNormalized);
00430                 output.update(1.0,diffRule,1.0);
00431               }
00432             }
00433           }
00434         } 
00435       }
00436     }
00437   }
00438   else 
00439     std::cout << "Dimension Must Be Less Than 5\n";
00440 }
00441 
00442 } // end Intrepid namespace