|
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 00051 #include "Intrepid_BurkardtRules.hpp" 00052 #include "Teuchos_oblackholestream.hpp" 00053 #include "Teuchos_RCP.hpp" 00054 #include "Teuchos_ScalarTraits.hpp" 00055 #include "Teuchos_GlobalMPISession.hpp" 00056 #include "Teuchos_Array.hpp" 00057 00058 //# include <cstdlib> 00059 //# include <iostream> 00060 //# include <cmath> 00061 //# include <iomanip> 00062 00063 using namespace std; 00064 using namespace Intrepid; 00065 00066 #define INTREPID_TEST_COMMAND( S ) \ 00067 { \ 00068 try { \ 00069 S ; \ 00070 } \ 00071 catch (std::logic_error err) { \ 00072 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00073 *outStream << err.what() << '\n'; \ 00074 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00075 }; \ 00076 } 00077 00078 template<class Scalar> 00079 Scalar evalQuad(int order, int power, Scalar x[], Scalar w[]) { 00080 00081 int mid = order/2; 00082 Scalar Q = 0.0; 00083 if (order%2) 00084 Q = w[mid]*powl(x[mid],power); 00085 00086 for (int i=0; i<mid; i++) { 00087 Q += w[i]*powl(x[i],power)+w[order-i-1]*powl(x[order-i-1],power); 00088 } 00089 00090 return Q; 00091 /* 00092 Scalar Q = 0.0; 00093 for (int i=0; i<order; i++) { 00094 Q += w[i]*powl(x[i],power); 00095 } 00096 return Q; 00097 */ 00098 } 00099 00100 template<class Scalar> 00101 Scalar factorial2 (int n) { 00102 Scalar value = 1.0; 00103 if (n<1) 00104 return value; 00105 00106 int n_copy = n; 00107 while (1<n_copy) { 00108 value *= (Scalar)n_copy; 00109 n_copy -= 2; 00110 } 00111 return value; 00112 } 00113 00114 template<class Scalar> 00115 Scalar chebyshev1(int power) { 00116 Scalar bot, exact, top; 00117 if (!(power%2)) { 00118 top = 1; bot = 1; 00119 for (int i=2;i<=power;i+=2) { 00120 top *= (Scalar)(i-1); 00121 bot *= (Scalar)i; 00122 } 00123 exact = M_PI*top/bot; 00124 } 00125 else { 00126 exact = 0.0; 00127 } 00128 return exact; 00129 } 00130 00131 template<class Scalar> 00132 Scalar chebyshev2(int power) { 00133 Scalar bot, exact, top; 00134 if (!(power%2)) { 00135 top = 1; bot = 1; 00136 for (int i=2;i<=power;i+=2) { 00137 top *= (Scalar)(i-1); 00138 bot *= (Scalar)i; 00139 } 00140 bot *= (Scalar)(power+2); 00141 exact = M_PI*top/bot; 00142 } 00143 else { 00144 exact = 0.0; 00145 } 00146 return exact; 00147 } 00148 00149 int main(int argc, char *argv[]) { 00150 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00151 00152 // This little trick lets us print to std::cout only if 00153 // a (dummy) command-line argument is provided. 00154 int iprint = argc - 1; 00155 Teuchos::RCP<std::ostream> outStream; 00156 Teuchos::oblackholestream bhs; // outputs nothing 00157 if (iprint > 0) 00158 outStream = Teuchos::rcp(&std::cout, false); 00159 else 00160 outStream = Teuchos::rcp(&bhs, false); 00161 00162 // Save the format state of the original std::cout. 00163 Teuchos::oblackholestream oldFormatState; 00164 oldFormatState.copyfmt(std::cout); 00165 00166 *outStream \ 00167 << "===============================================================================\n" \ 00168 << "| |\n" \ 00169 << "| Unit Test (IntrepidBurkardtRules) |\n" \ 00170 << "| |\n" \ 00171 << "| 1) the Burkardt rule tests |\n" \ 00172 << "| |\n" \ 00173 << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \ 00174 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00175 << "| |\n" \ 00176 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00177 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00178 << "| |\n" \ 00179 << "===============================================================================\n"; 00180 00181 00182 int errorFlag = 0; 00183 00184 int maxOrder = 30; 00185 long double reltol = 1e-8; 00186 long double analyticInt = 0, testInt = 0; 00187 // compute and compare integrals 00188 try { 00189 00190 *outStream << "Gauss-Legendre Cubature \n"; 00191 *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n"; 00192 for (int i = 1; i<=maxOrder; i++) { 00193 Teuchos::Array<long double> nodes(i), weights(i); 00194 IntrepidBurkardtRules::legendre_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00195 for (int j=0; j<=2*i-1; j++) { 00196 if (j%2) 00197 analyticInt = 0.0; 00198 else 00199 analyticInt = 2.0/((long double)j+1.0); 00200 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00201 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00202 long double absdiff = std::fabs(analyticInt - testInt); 00203 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00204 << "x^" << std::setw(2) << std::left << j << ":" << " " 00205 << std::scientific << std::setprecision(16) << testInt << " " 00206 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00207 << " " << abstol << "\n"; 00208 if (absdiff > abstol) { 00209 errorFlag++; 00210 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00211 } 00212 } 00213 } 00214 *outStream << "\n"; 00215 00216 *outStream << "Clenshaw-Curtis Cubature \n"; 00217 *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n"; 00218 for (int i = 1; i<=maxOrder; i++) { 00219 Teuchos::Array<long double> nodes(i), weights(i); 00220 IntrepidBurkardtRules::clenshaw_curtis_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00221 for (int j=0; j<i; j++) { 00222 if (j%2) 00223 analyticInt = 0.0; 00224 else 00225 analyticInt = 2.0/((long double)j+1.0); 00226 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00227 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00228 long double absdiff = std::fabs(analyticInt - testInt); 00229 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00230 << "x^" << std::setw(2) << std::left << j << ":" << " " 00231 << std::scientific << std::setprecision(16) << testInt << " " 00232 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00233 << " " << abstol << "\n"; 00234 if (absdiff > abstol) { 00235 errorFlag++; 00236 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00237 } 00238 } 00239 } 00240 *outStream << "\n"; 00241 00242 *outStream << "Fejer Type 2 Cubature \n"; 00243 *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n"; 00244 for (int i = 1; i<=maxOrder; i++) { 00245 Teuchos::Array<long double> nodes(i), weights(i); 00246 IntrepidBurkardtRules::fejer2_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00247 for (int j=0; j<i; j++) { 00248 if (j%2) 00249 analyticInt = 0.0; 00250 else 00251 analyticInt = 2.0/((long double)j+1.0); 00252 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00253 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00254 long double absdiff = std::fabs(analyticInt - testInt); 00255 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00256 << "x^" << std::setw(2) << std::left << j << ":" << " " 00257 << std::scientific << std::setprecision(16) << testInt << " " 00258 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00259 << " " << abstol << "\n"; 00260 if (absdiff > abstol) { 00261 errorFlag++; 00262 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00263 } 00264 } 00265 } 00266 *outStream << "\n"; 00267 00268 *outStream << "Gauss-Patterson Cubature \n"; 00269 *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n"; 00270 for (int l = 1; l<=7; l++) { 00271 int i = (int)pow(2.0,(double)l+1.0)-1; 00272 Teuchos::Array<long double> nodes(i), weights(i); 00273 IntrepidBurkardtRules::patterson_lookup(i,nodes.getRawPtr(),weights.getRawPtr()); 00274 for (int j=0; j<=(1.5*i+0.5); j++) { 00275 if (j%2) 00276 analyticInt = 0.0; 00277 else 00278 analyticInt = 2.0/((long double)j+1.0); 00279 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00280 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00281 long double absdiff = std::fabs(analyticInt - testInt); 00282 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00283 << "x^" << std::setw(2) << std::left << j << ":" << " " 00284 << std::scientific << std::setprecision(16) << testInt << " " 00285 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00286 << " " << abstol << "\n"; 00287 if (absdiff > abstol) { 00288 errorFlag++; 00289 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00290 } 00291 } 00292 } 00293 *outStream << "\n"; 00294 00295 *outStream << "Gauss-Chebyshev Type 1 Cubature \n"; 00296 *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1/sqrt(1-x^2)\n"; 00297 for (int i = 1; i<=maxOrder; i++) { 00298 Teuchos::Array<long double> nodes(i), weights(i); 00299 IntrepidBurkardtRules::chebyshev1_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00300 for (int j=0; j<=2*i-1; j++) { 00301 analyticInt = chebyshev1<long double>(j); 00302 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00303 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00304 long double absdiff = std::fabs(analyticInt - testInt); 00305 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00306 << "x^" << std::setw(2) << std::left << j << ":" << " " 00307 << std::scientific << std::setprecision(16) << testInt << " " 00308 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00309 << " " << abstol << "\n"; 00310 if (absdiff > abstol) { 00311 errorFlag++; 00312 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00313 } 00314 } 00315 } 00316 *outStream << "\n"; 00317 00318 *outStream << "Gauss-Chebyshev Type 2 Cubature \n"; 00319 *outStream << "Integrates functions on [-1,1] weighted by w(x) = sqrt(1-x^2)\n"; 00320 for (int i = 1; i<=maxOrder; i++) { 00321 Teuchos::Array<long double> nodes(i), weights(i); 00322 IntrepidBurkardtRules::chebyshev2_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00323 for (int j=0; j<=2*i-1; j++) { 00324 analyticInt = chebyshev2<long double>(j); 00325 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00326 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00327 long double absdiff = std::fabs(analyticInt - testInt); 00328 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00329 << "x^" << std::setw(2) << std::left << j << ":" << " " 00330 << std::scientific << std::setprecision(16) << testInt << " " 00331 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00332 << " " << abstol << "\n"; 00333 if (absdiff > abstol) { 00334 errorFlag++; 00335 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00336 } 00337 } 00338 } 00339 *outStream << "\n"; 00340 00341 *outStream << "Gauss-Laguerre Cubature \n"; 00342 *outStream << "Integrates functions on [0,oo) weighted by w(x) = exp(-x)\n"; 00343 for (int i = 1; i<=maxOrder; i++) { 00344 Teuchos::Array<long double> nodes(i), weights(i); 00345 IntrepidBurkardtRules::laguerre_compute(i,nodes.getRawPtr(),weights.getRawPtr()); 00346 for (int j=0; j<=2*i-1; j++) { 00347 analyticInt = tgammal((long double)(j+1)); 00348 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00349 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00350 long double absdiff = std::fabs(analyticInt - testInt); 00351 *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating " 00352 << "x^" << std::setw(2) << std::left << j << ":" << " " 00353 << std::scientific << std::setprecision(16) << testInt << " " 00354 << analyticInt << " " << std::setprecision(4) << absdiff << " " << "<?" 00355 << " " << abstol << "\n"; 00356 if (absdiff > abstol) { 00357 errorFlag++; 00358 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00359 } 00360 } 00361 } 00362 *outStream << "\n"; 00363 00364 maxOrder = 10; 00365 00366 *outStream << "Gauss-Hermite Cubature \n"; 00367 *outStream << "Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n"; 00368 for (int i = 1; i<=maxOrder; i++) { 00369 Teuchos::Array<long double> nodes(i), weights(i); 00370 IntrepidBurkardtRules::hermite_compute(i,nodes.getRawPtr(), 00371 weights.getRawPtr()); 00372 for (int j=0; j<=2*i-1; j++) { 00373 if (j%2) 00374 analyticInt = 0.0; 00375 else 00376 analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(long double)j/2.0); 00377 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00378 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00379 long double absdiff = std::fabs(analyticInt - testInt); 00380 *outStream << "Cubature order " << std::setw(2) << std::left << i 00381 << " integrating " 00382 << "x^" << std::setw(2) << std::left << j << ":" 00383 << " " 00384 << std::scientific << std::setprecision(16) << testInt 00385 << " " 00386 << analyticInt << " " << std::setprecision(4) 00387 << absdiff << " " << "<?" 00388 << " " << abstol << "\n"; 00389 if (absdiff > abstol) { 00390 errorFlag++; 00391 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00392 } 00393 } 00394 } 00395 *outStream << "\n"; 00396 00397 reltol = 1e-6; 00398 00399 *outStream << "Hermite-Genz-Keister Cubature \n"; 00400 *outStream << "Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n"; 00401 int order[4] = {1,3, 9,19}; 00402 int max[4] = {1,5,15,29}; 00403 for (int l = 0; l<4; l++) { 00404 int i = order[l]; 00405 int m = max[l]; 00406 Teuchos::Array<long double> nodes(i), weights(i); 00407 IntrepidBurkardtRules::hermite_genz_keister_lookup(i,nodes.getRawPtr(), 00408 weights.getRawPtr()); 00409 for (int j=0; j<=m; j++) { 00410 if (j%2) 00411 analyticInt = 0.0; 00412 else 00413 analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(long double)j/2.0); 00414 if (i>=36) 00415 analyticInt /= sqrt(M_PI); 00416 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr()); 00417 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) ); 00418 long double absdiff = std::fabs(analyticInt - testInt); 00419 *outStream << "Cubature order " << std::setw(2) << std::left << i 00420 << " integrating " 00421 << "x^" << std::setw(2) << std::left << j << ":" << " " 00422 << std::scientific << std::setprecision(16) << testInt 00423 << " " 00424 << analyticInt << " " << std::setprecision(4) << absdiff 00425 << " " << "<?" 00426 << " " << abstol << "\n"; 00427 if (absdiff > abstol) { 00428 errorFlag++; 00429 *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n"; 00430 } 00431 } 00432 } 00433 *outStream << "\n"; 00434 } 00435 catch (std::logic_error err) { 00436 *outStream << err.what() << "\n"; 00437 errorFlag = -1; 00438 }; 00439 00440 00441 if (errorFlag != 0) 00442 std::cout << "End Result: TEST FAILED\n"; 00443 else 00444 std::cout << "End Result: TEST PASSED\n"; 00445 00446 // reset format state of std::cout 00447 std::cout.copyfmt(oldFormatState); 00448 00449 return errorFlag; 00450 } 00451
1.7.6.1