Intrepid
test_01.cpp
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
51 #include "Intrepid_BurkardtRules.hpp"
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"
57 
58 //# include <cstdlib>
59 //# include <iostream>
60 //# include <cmath>
61 //# include <iomanip>
62 
63 using namespace std;
64 using namespace Intrepid;
65 
66 #define INTREPID_TEST_COMMAND( S ) \
67 { \
68  try { \
69  S ; \
70  } \
71  catch (const std::logic_error & err) { \
72  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
73  *outStream << err.what() << '\n'; \
74  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
75  }; \
76 }
77 
78 template<class Scalar>
79 Scalar evalQuad(int order, int power, Scalar x[], Scalar w[]) {
80 
81  int mid = order/2;
82  Scalar Q = 0.0;
83  if (order%2)
84  Q = w[mid]*powl(x[mid],power);
85 
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);
88  }
89 
90  return Q;
91  /*
92  Scalar Q = 0.0;
93  for (int i=0; i<order; i++) {
94  Q += w[i]*powl(x[i],power);
95  }
96  return Q;
97  */
98 }
99 
100 template<class Scalar>
101 Scalar factorial2 (int n) {
102  Scalar value = 1.0;
103  if (n<1)
104  return value;
105 
106  int n_copy = n;
107  while (1<n_copy) {
108  value *= (Scalar)n_copy;
109  n_copy -= 2;
110  }
111  return value;
112 }
113 
114 template<class Scalar>
115 Scalar chebyshev1(int power) {
116  Scalar bot, exact, top;
117  if (!(power%2)) {
118  top = 1; bot = 1;
119  for (int i=2;i<=power;i+=2) {
120  top *= (Scalar)(i-1);
121  bot *= (Scalar)i;
122  }
123  exact = M_PI*top/bot;
124  }
125  else {
126  exact = 0.0;
127  }
128  return exact;
129 }
130 
131 template<class Scalar>
132 Scalar chebyshev2(int power) {
133  Scalar bot, exact, top;
134  if (!(power%2)) {
135  top = 1; bot = 1;
136  for (int i=2;i<=power;i+=2) {
137  top *= (Scalar)(i-1);
138  bot *= (Scalar)i;
139  }
140  bot *= (Scalar)(power+2);
141  exact = M_PI*top/bot;
142  }
143  else {
144  exact = 0.0;
145  }
146  return exact;
147 }
148 
149 int main(int argc, char *argv[]) {
150  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
151 
152  // This little trick lets us print to std::cout only if
153  // a (dummy) command-line argument is provided.
154  int iprint = argc - 1;
155  Teuchos::RCP<std::ostream> outStream;
156  Teuchos::oblackholestream bhs; // outputs nothing
157  if (iprint > 0)
158  outStream = Teuchos::rcp(&std::cout, false);
159  else
160  outStream = Teuchos::rcp(&bhs, false);
161 
162  // Save the format state of the original std::cout.
163  Teuchos::oblackholestream oldFormatState;
164  oldFormatState.copyfmt(std::cout);
165 
166  *outStream \
167  << "===============================================================================\n" \
168  << "| |\n" \
169  << "| Unit Test (IntrepidBurkardtRules) |\n" \
170  << "| |\n" \
171  << "| 1) the Burkardt rule tests |\n" \
172  << "| |\n" \
173  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
174  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
175  << "| |\n" \
176  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
177  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
178  << "| |\n" \
179  << "===============================================================================\n";
180 
181 
182  int errorFlag = 0;
183 
184  int maxOrder = 30;
185  long double reltol = 1e-8;
186  long double analyticInt = 0, testInt = 0;
187  // compute and compare integrals
188  try {
189 
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++) {
196  if (j%2)
197  analyticInt = 0.0;
198  else
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) {
209  errorFlag++;
210  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
211  }
212  }
213  }
214  *outStream << "\n";
215 
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++) {
222  if (j%2)
223  analyticInt = 0.0;
224  else
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) {
235  errorFlag++;
236  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
237  }
238  }
239  }
240  *outStream << "\n";
241 
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++) {
248  if (j%2)
249  analyticInt = 0.0;
250  else
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) {
261  errorFlag++;
262  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
263  }
264  }
265  }
266  *outStream << "\n";
267 
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++) {
275  if (j%2)
276  analyticInt = 0.0;
277  else
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) {
288  errorFlag++;
289  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
290  }
291  }
292  }
293  *outStream << "\n";
294 
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) {
311  errorFlag++;
312  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
313  }
314  }
315  }
316  *outStream << "\n";
317 
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) {
334  errorFlag++;
335  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
336  }
337  }
338  }
339  *outStream << "\n";
340 
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) {
357  errorFlag++;
358  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
359  }
360  }
361  }
362  *outStream << "\n";
363 
364  maxOrder = 10;
365 
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++) {
373  if (j%2)
374  analyticInt = 0.0;
375  else
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
381  << " integrating "
382  << "x^" << std::setw(2) << std::left << j << ":"
383  << " "
384  << std::scientific << std::setprecision(16) << testInt
385  << " "
386  << analyticInt << " " << std::setprecision(4)
387  << absdiff << " " << "<?"
388  << " " << abstol << "\n";
389  if (absdiff > abstol) {
390  errorFlag++;
391  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
392  }
393  }
394  }
395  *outStream << "\n";
396 
397  reltol = 1e-6;
398 
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++) {
404  int i = order[l];
405  int m = max[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++) {
410  if (j%2)
411  analyticInt = 0.0;
412  else
413  analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(long double)j/2.0);
414  if (i>=36)
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
420  << " integrating "
421  << "x^" << std::setw(2) << std::left << j << ":" << " "
422  << std::scientific << std::setprecision(16) << testInt
423  << " "
424  << analyticInt << " " << std::setprecision(4) << absdiff
425  << " " << "<?"
426  << " " << abstol << "\n";
427  if (absdiff > abstol) {
428  errorFlag++;
429  *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
430  }
431  }
432  }
433  *outStream << "\n";
434  }
435  catch (const std::logic_error & err) {
436  *outStream << err.what() << "\n";
437  errorFlag = -1;
438  };
439 
440 
441  if (errorFlag != 0)
442  std::cout << "End Result: TEST FAILED\n";
443  else
444  std::cout << "End Result: TEST PASSED\n";
445 
446  // reset format state of std::cout
447  std::cout.copyfmt(oldFormatState);
448 
449  return errorFlag;
450 }
451 
Intrepid_BurkardtRules.hpp
Header file for integration rules provided by John Burkardt. &lt;&gt;
Generated by   doxygen 1.8.5