52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 #include "Teuchos_Array.hpp"
64 using namespace Intrepid;
66 #define INTREPID_TEST_COMMAND( S ) \
71 catch (std::logic_error err) { \
72 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
73 *outStream << err.what() << '\n'; \
74 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
78 template<
class Scalar>
79 Scalar evalQuad(
int order,
int power, Scalar x[], Scalar w[]) {
84 Q = w[mid]*powl(x[mid],power);
86 for (
int i=0; i<mid; i++) {
87 Q += w[i]*powl(x[i],power)+w[order-i-1]*powl(x[order-i-1],power);
100 template<
class Scalar>
101 Scalar factorial2 (
int n) {
108 value *= (Scalar)n_copy;
114 template<
class Scalar>
115 Scalar chebyshev1(
int power) {
116 Scalar bot, exact, top;
119 for (
int i=2;i<=power;i+=2) {
120 top *= (Scalar)(i-1);
123 exact = M_PI*top/bot;
131 template<
class Scalar>
132 Scalar chebyshev2(
int power) {
133 Scalar bot, exact, top;
136 for (
int i=2;i<=power;i+=2) {
137 top *= (Scalar)(i-1);
140 bot *= (Scalar)(power+2);
141 exact = M_PI*top/bot;
149 int main(
int argc,
char *argv[]) {
150 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
154 int iprint = argc - 1;
155 Teuchos::RCP<std::ostream> outStream;
156 Teuchos::oblackholestream bhs;
158 outStream = Teuchos::rcp(&std::cout,
false);
160 outStream = Teuchos::rcp(&bhs,
false);
163 Teuchos::oblackholestream oldFormatState;
164 oldFormatState.copyfmt(std::cout);
167 <<
"===============================================================================\n" \
169 <<
"| Unit Test (IntrepidBurkardtRules) |\n" \
171 <<
"| 1) the Burkardt rule tests |\n" \
173 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
174 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
176 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
177 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
179 <<
"===============================================================================\n";
185 long double reltol = 1e-8;
186 long double analyticInt = 0, testInt = 0;
190 *outStream <<
"Gauss-Legendre Cubature \n";
191 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = 1\n";
192 for (
int i = 1; i<=maxOrder; i++) {
193 Teuchos::Array<long double> nodes(i), weights(i);
194 IntrepidBurkardtRules::legendre_compute(i,nodes.getRawPtr(),weights.getRawPtr());
195 for (
int j=0; j<=2*i-1; j++) {
199 analyticInt = 2.0/((
long double)j+1.0);
200 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
201 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
202 long double absdiff = std::fabs(analyticInt - testInt);
203 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
204 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
205 << std::scientific << std::setprecision(16) << testInt <<
" "
206 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
207 <<
" " << abstol <<
"\n";
208 if (absdiff > abstol) {
210 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
216 *outStream <<
"Clenshaw-Curtis Cubature \n";
217 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = 1\n";
218 for (
int i = 1; i<=maxOrder; i++) {
219 Teuchos::Array<long double> nodes(i), weights(i);
220 IntrepidBurkardtRules::clenshaw_curtis_compute(i,nodes.getRawPtr(),weights.getRawPtr());
221 for (
int j=0; j<i; j++) {
225 analyticInt = 2.0/((
long double)j+1.0);
226 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
227 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
228 long double absdiff = std::fabs(analyticInt - testInt);
229 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
230 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
231 << std::scientific << std::setprecision(16) << testInt <<
" "
232 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
233 <<
" " << abstol <<
"\n";
234 if (absdiff > abstol) {
236 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
242 *outStream <<
"Fejer Type 2 Cubature \n";
243 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = 1\n";
244 for (
int i = 1; i<=maxOrder; i++) {
245 Teuchos::Array<long double> nodes(i), weights(i);
246 IntrepidBurkardtRules::fejer2_compute(i,nodes.getRawPtr(),weights.getRawPtr());
247 for (
int j=0; j<i; j++) {
251 analyticInt = 2.0/((
long double)j+1.0);
252 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
253 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
254 long double absdiff = std::fabs(analyticInt - testInt);
255 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
256 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
257 << std::scientific << std::setprecision(16) << testInt <<
" "
258 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
259 <<
" " << abstol <<
"\n";
260 if (absdiff > abstol) {
262 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
268 *outStream <<
"Gauss-Patterson Cubature \n";
269 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = 1\n";
270 for (
int l = 1; l<=7; l++) {
271 int i = (int)pow(2.0,(
double)l+1.0)-1;
272 Teuchos::Array<long double> nodes(i), weights(i);
273 IntrepidBurkardtRules::patterson_lookup(i,nodes.getRawPtr(),weights.getRawPtr());
274 for (
int j=0; j<=(1.5*i+0.5); j++) {
278 analyticInt = 2.0/((
long double)j+1.0);
279 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
280 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
281 long double absdiff = std::fabs(analyticInt - testInt);
282 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
283 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
284 << std::scientific << std::setprecision(16) << testInt <<
" "
285 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
286 <<
" " << abstol <<
"\n";
287 if (absdiff > abstol) {
289 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
295 *outStream <<
"Gauss-Chebyshev Type 1 Cubature \n";
296 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = 1/sqrt(1-x^2)\n";
297 for (
int i = 1; i<=maxOrder; i++) {
298 Teuchos::Array<long double> nodes(i), weights(i);
299 IntrepidBurkardtRules::chebyshev1_compute(i,nodes.getRawPtr(),weights.getRawPtr());
300 for (
int j=0; j<=2*i-1; j++) {
301 analyticInt = chebyshev1<long double>(j);
302 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
303 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
304 long double absdiff = std::fabs(analyticInt - testInt);
305 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
306 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
307 << std::scientific << std::setprecision(16) << testInt <<
" "
308 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
309 <<
" " << abstol <<
"\n";
310 if (absdiff > abstol) {
312 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
318 *outStream <<
"Gauss-Chebyshev Type 2 Cubature \n";
319 *outStream <<
"Integrates functions on [-1,1] weighted by w(x) = sqrt(1-x^2)\n";
320 for (
int i = 1; i<=maxOrder; i++) {
321 Teuchos::Array<long double> nodes(i), weights(i);
322 IntrepidBurkardtRules::chebyshev2_compute(i,nodes.getRawPtr(),weights.getRawPtr());
323 for (
int j=0; j<=2*i-1; j++) {
324 analyticInt = chebyshev2<long double>(j);
325 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
326 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
327 long double absdiff = std::fabs(analyticInt - testInt);
328 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
329 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
330 << std::scientific << std::setprecision(16) << testInt <<
" "
331 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
332 <<
" " << abstol <<
"\n";
333 if (absdiff > abstol) {
335 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
341 *outStream <<
"Gauss-Laguerre Cubature \n";
342 *outStream <<
"Integrates functions on [0,oo) weighted by w(x) = exp(-x)\n";
343 for (
int i = 1; i<=maxOrder; i++) {
344 Teuchos::Array<long double> nodes(i), weights(i);
345 IntrepidBurkardtRules::laguerre_compute(i,nodes.getRawPtr(),weights.getRawPtr());
346 for (
int j=0; j<=2*i-1; j++) {
347 analyticInt = tgammal((
long double)(j+1));
348 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
349 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
350 long double absdiff = std::fabs(analyticInt - testInt);
351 *outStream <<
"Cubature order " << std::setw(2) << std::left << i <<
" integrating "
352 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
353 << std::scientific << std::setprecision(16) << testInt <<
" "
354 << analyticInt <<
" " << std::setprecision(4) << absdiff <<
" " <<
"<?"
355 <<
" " << abstol <<
"\n";
356 if (absdiff > abstol) {
358 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
366 *outStream <<
"Gauss-Hermite Cubature \n";
367 *outStream <<
"Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n";
368 for (
int i = 1; i<=maxOrder; i++) {
369 Teuchos::Array<long double> nodes(i), weights(i);
370 IntrepidBurkardtRules::hermite_compute(i,nodes.getRawPtr(),
371 weights.getRawPtr());
372 for (
int j=0; j<=2*i-1; j++) {
376 analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(
long double)j/2.0);
377 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
378 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
379 long double absdiff = std::fabs(analyticInt - testInt);
380 *outStream <<
"Cubature order " << std::setw(2) << std::left << i
382 <<
"x^" << std::setw(2) << std::left << j <<
":"
384 << std::scientific << std::setprecision(16) << testInt
386 << analyticInt <<
" " << std::setprecision(4)
387 << absdiff <<
" " <<
"<?"
388 <<
" " << abstol <<
"\n";
389 if (absdiff > abstol) {
391 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
399 *outStream <<
"Hermite-Genz-Keister Cubature \n";
400 *outStream <<
"Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n";
401 int order[4] = {1,3, 9,19};
402 int max[4] = {1,5,15,29};
403 for (
int l = 0; l<4; l++) {
406 Teuchos::Array<long double> nodes(i), weights(i);
407 IntrepidBurkardtRules::hermite_genz_keister_lookup(i,nodes.getRawPtr(),
408 weights.getRawPtr());
409 for (
int j=0; j<=m; j++) {
413 analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(
long double)j/2.0);
415 analyticInt /= sqrt(M_PI);
416 testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
417 long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
418 long double absdiff = std::fabs(analyticInt - testInt);
419 *outStream <<
"Cubature order " << std::setw(2) << std::left << i
421 <<
"x^" << std::setw(2) << std::left << j <<
":" <<
" "
422 << std::scientific << std::setprecision(16) << testInt
424 << analyticInt <<
" " << std::setprecision(4) << absdiff
426 <<
" " << abstol <<
"\n";
427 if (absdiff > abstol) {
429 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
435 catch (std::logic_error err) {
436 *outStream << err.what() <<
"\n";
442 std::cout <<
"End Result: TEST FAILED\n";
444 std::cout <<
"End Result: TEST PASSED\n";
447 std::cout.copyfmt(oldFormatState);
Header file for integration rules provided by John Burkardt. <>