|
Intrepid
|
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 00044 00051 #include "Intrepid_AdaptiveSparseGrid.hpp" 00052 //#include "Intrepid_CubatureLineSorted.hpp" 00053 #include "Intrepid_Utils.hpp" 00054 #include "Teuchos_oblackholestream.hpp" 00055 #include "Teuchos_RCP.hpp" 00056 #include "Teuchos_RefCountPtr.hpp" 00057 #include "Teuchos_GlobalMPISession.hpp" 00058 00059 using namespace Intrepid; 00060 00061 00062 template<class Scalar> 00063 class StdVector { 00064 private: 00065 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_; 00066 00067 public: 00068 00069 StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec ) 00070 : std_vec_(std_vec) {} 00071 00072 Teuchos::RefCountPtr<StdVector<Scalar> > Create() const { 00073 return Teuchos::rcp( new StdVector<Scalar>( 00074 Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0)))); 00075 } 00076 00077 void Update( StdVector<Scalar> & s ) { 00078 int dimension = (int)(std_vec_->size()); 00079 for (int i=0; i<dimension; i++) 00080 (*std_vec_)[i] += s[i]; 00081 } 00082 00083 void Update( Scalar alpha, StdVector<Scalar> & s ) { 00084 int dimension = (int)(std_vec_->size()); 00085 for (int i=0; i<dimension; i++) 00086 (*std_vec_)[i] += alpha*s[i]; 00087 } 00088 00089 Scalar operator[](int i) { 00090 return (*std_vec_)[i]; 00091 } 00092 00093 void clear() { 00094 std_vec_->clear(); 00095 } 00096 00097 void resize(int n, Scalar p) { 00098 std_vec_->resize(n,p); 00099 } 00100 00101 int size() { 00102 return (int)std_vec_->size(); 00103 } 00104 00105 void Set( Scalar alpha ) { 00106 int dimension = (int)(std_vec_->size()); 00107 for (int i=0; i<dimension; i++) 00108 (*std_vec_)[i] = alpha; 00109 } 00110 }; 00111 00112 template<class Scalar, class UserVector> 00113 class ASGdata : 00114 public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> { 00115 public: 00116 ~ASGdata() {} 00117 00118 ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D, 00119 std::vector<EIntrepidGrowth> growth1D, int maxLevel, 00120 bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>( 00121 dimension,rule1D,growth1D,maxLevel,isNormalized) {} 00122 00123 void eval_integrand(UserVector & output, std::vector<Scalar> & input) { 00124 output.clear(); output.resize(1,powl(input[0]+input[1],(long double)6.0)); 00125 } 00126 Scalar error_indicator(UserVector & input) { 00127 int dimension = (int)input.size(); 00128 Scalar norm2 = 0.0; 00129 for (int i=0; i<dimension; i++) 00130 norm2 += input[i]*input[i]; 00131 00132 Scalar ID = AdaptiveSparseGridInterface<Scalar,UserVector>:: 00133 getInitialDiff(); 00134 norm2 = std::sqrt(norm2)/ID; 00135 return norm2; 00136 } 00137 }; 00138 00139 long double adaptSG(StdVector<long double> & iv, 00140 std::multimap<long double,std::vector<int> > & activeIndex, 00141 std::set<std::vector<int> > & oldIndex, 00142 AdaptiveSparseGridInterface<long double,StdVector<long double> > & problem_data, 00143 CubatureTensorSorted<long double> & cubRule, 00144 long double TOL) { 00145 00146 // Construct a Container for the adapted rule 00147 int dimension = problem_data.getDimension(); 00148 std::vector<int> index(dimension,1); 00149 00150 // Initialize global error indicator 00151 long double eta = 1.0; 00152 00153 // Initialize the Active index set 00154 activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index)); 00155 00156 // Perform Adaptation 00157 while (eta > TOL) { 00158 eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid( 00159 activeIndex,oldIndex, 00160 iv,cubRule, 00161 eta,problem_data); 00162 } 00163 cubRule.normalize(); 00164 return eta; 00165 } 00166 00167 long double evalQuad(CubatureTensorSorted<long double> & lineCub) { 00168 00169 int size = lineCub.getNumPoints(); 00170 int dimension = lineCub.getDimension(); 00171 FieldContainer<long double> cubPoints(size,dimension); 00172 FieldContainer<long double> cubWeights(size); 00173 lineCub.getCubature(cubPoints,cubWeights); 00174 00175 long double Q = 0.0; 00176 for (int k=0; k<size; k++) 00177 Q += cubWeights(k)*powl(cubPoints(k,0)+cubPoints(k,1),(long double)6.0); 00178 00179 return Q; 00180 } 00181 00182 int main(int argc, char *argv[]) { 00183 00184 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00185 00186 // This little trick lets us print to std::cout only if 00187 // a (dummy) command-line argument is provided. 00188 int iprint = argc - 1; 00189 Teuchos::RCP<std::ostream> outStream; 00190 Teuchos::oblackholestream bhs; // outputs nothing 00191 if (iprint > 0) 00192 outStream = Teuchos::rcp(&std::cout, false); 00193 else 00194 outStream = Teuchos::rcp(&bhs, false); 00195 00196 // Save the format state of the original std::cout. 00197 Teuchos::oblackholestream oldFormatState; 00198 oldFormatState.copyfmt(std::cout); 00199 00200 *outStream \ 00201 << "===============================================================================\n" \ 00202 << "| |\n" \ 00203 << "| Unit Test (AdaptiveSparseGrid) |\n" \ 00204 << "| |\n" \ 00205 << "| 1) Integrate a sum of Gaussians in 2D and compare index sets. |\n" \ 00206 << "| |\n" \ 00207 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \ 00208 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00209 << "| |\n" \ 00210 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00211 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00212 << "| |\n" \ 00213 << "===============================================================================\n"\ 00214 << "| TEST 23: Compare index sets for different instances of refine grid |\n"\ 00215 << "===============================================================================\n"; 00216 00217 00218 // internal variables: 00219 int errorFlag = 0; 00220 long double TOL = INTREPID_TOL; 00221 int dimension = 2; 00222 int maxLevel = 4; 00223 bool isNormalized = false; 00224 00225 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS); 00226 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP); 00227 00228 ASGdata<long double,StdVector<long double> > problem_data(dimension,rule1D, 00229 growth1D,maxLevel, 00230 isNormalized); 00231 Teuchos::RCP<std::vector<long double> > integralValue = 00232 Teuchos::rcp(new std::vector<long double>(1,0.0)); 00233 StdVector<long double> sol(integralValue); sol.Set(0.0); 00234 problem_data.init(sol); 00235 00236 try { 00237 00238 // Initialize the index sets 00239 std::multimap<long double,std::vector<int> > activeIndex1; 00240 std::set<std::vector<int> > oldIndex1; 00241 std::vector<int> index(dimension,1); 00242 CubatureTensorSorted<long double> adaptedRule(dimension,index,rule1D, 00243 growth1D,isNormalized); 00244 adaptSG(sol,activeIndex1,oldIndex1,problem_data,adaptedRule,TOL); 00245 long double Q1 = sol[0]; 00246 00247 CubatureTensorSorted<long double> fullRule(0,dimension); 00248 AdaptiveSparseGrid<long double,StdVector<long double> >::buildSparseGrid( 00249 fullRule,dimension, 00250 maxLevel,rule1D, 00251 growth1D,isNormalized); 00252 long double Q2 = evalQuad(fullRule); 00253 fullRule.normalize(); 00254 00255 long double diff = fabs(Q1-Q2); 00256 00257 *outStream << "Q1 = " << Q1 << " Q2 = " << Q2 00258 << " |Q1-Q2| = " << diff << "\n"; 00259 00260 int size1 = adaptedRule.getNumPoints(); 00261 FieldContainer<long double> aPoints(size1,dimension); 00262 FieldContainer<long double> aWeights(size1); 00263 adaptedRule.getCubature(aPoints,aWeights); 00264 00265 *outStream << "\n\nAdapted Rule Nodes and Weights\n"; 00266 for (int i=0; i<size1; i++) 00267 *outStream << aPoints(i,0) << "\t" << aPoints(i,1) 00268 << "\t" << aWeights(i) << "\n"; 00269 00270 int size2 = fullRule.getNumPoints(); 00271 FieldContainer<long double> fPoints(size2,dimension); 00272 FieldContainer<long double> fWeights(size2); 00273 fullRule.getCubature(fPoints,fWeights); 00274 00275 *outStream << "\n\nFull Rule Nodes and Weights\n"; 00276 for (int i=0; i<size2; i++) 00277 *outStream << fPoints(i,0) << "\t" << fPoints(i,1) 00278 << "\t" << fWeights(i) << "\n"; 00279 00280 *outStream << "\n\nSize of adapted rule = " << size1 00281 << " Size of full rule = " << size2 << "\n"; 00282 if (diff > TOL*fabs(Q2)||size1!=size2) { 00283 errorFlag++; 00284 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n"; 00285 } 00286 else { 00287 long double sum1 = 0.0, sum2 = 0.0; 00288 for (int i=0; i<size1; i++) { 00289 //diff = fabs(fWeights(i)-aWeights(i)); 00290 sum1 += fWeights(i); 00291 sum2 += aWeights(i); 00292 } 00293 *outStream << "Check if weights are normalized:" 00294 << " Adapted Rule Sum = " << sum2 00295 << " Full Rule Sum = " << sum1 << "\n"; 00296 if (fabs(sum1-1.0) > TOL || fabs(sum2-1.0) > TOL) { 00297 errorFlag++; 00298 *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n"; 00299 } 00300 } 00301 } 00302 catch (std::logic_error err) { 00303 *outStream << err.what() << "\n"; 00304 errorFlag = -1; 00305 }; 00306 00307 if (errorFlag != 0) 00308 std::cout << "End Result: TEST FAILED\n"; 00309 else 00310 std::cout << "End Result: TEST PASSED\n"; 00311 00312 // reset format state of std::cout 00313 std::cout.copyfmt(oldFormatState); 00314 00315 return errorFlag; 00316 }
1.7.6.1