Intrepid
test_11.cpp
Go to the documentation of this file.
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 
44 
51 //#include "Intrepid_CubatureLineSorted.hpp"
53 #include "Intrepid_Utils.hpp"
54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_GlobalMPISession.hpp"
57 
58 using namespace Intrepid;
59 
60 /*
61  Computes integrals of monomials over a given reference cell.
62 */
63 long double evalQuad(int order, int power, EIntrepidBurkardt rule) {
64 
65  CubatureLineSorted<long double> lineCub(rule,order,false);
66  int size = lineCub.getNumPoints();
67  FieldContainer<long double> cubPoints(size);
68  FieldContainer<long double> cubWeights(size);
69  lineCub.getCubature(cubPoints,cubWeights);
70 
71  int mid = size/2;
72  long double Q = 0.0;
73  if (size%2)
74  Q = cubWeights(mid)*powl(cubPoints(mid),power);
75 
76  for (int i=0; i<mid; i++) {
77  Q += cubWeights(i)*powl(cubPoints(i),power)+
78  cubWeights(size-i-1)*powl(cubPoints(size-i-1),power);
79  }
80  return Q;
81 }
82 
83 long double evalInt(int power, EIntrepidBurkardt rule) {
84  double I = 0;
85 
86  if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2||
87  rule==BURK_LEGENDRE||rule==BURK_PATTERSON ||
88  rule==BURK_TRAPEZOIDAL) {
89  if (power%2)
90  I = 0.0;
91  else
92  I = 2.0/((long double)power+1.0);
93  }
94  else if (rule==BURK_LAGUERRE) {
95  I = tgammal((long double)(power+1));
96  }
97  else if (rule==BURK_CHEBYSHEV1) {
98  long double bot, top;
99  if (!(power%2)) {
100  top = 1; bot = 1;
101  for (int i=2;i<=power;i+=2) {
102  top *= (long double)(i-1);
103  bot *= (long double)i;
104  }
105  I = M_PI*top/bot;
106  }
107  else {
108  I = 0.0;
109  }
110  }
111  else if (rule==BURK_CHEBYSHEV2) {
112  long double bot, top;
113  if (!(power%2)) {
114  top = 1; bot = 1;
115  for (int i=2;i<=power;i+=2) {
116  top *= (long double)(i-1);
117  bot *= (long double)i;
118  }
119  bot *= (long double)(power+2);
120  I = M_PI*top/bot;
121  }
122  else {
123  I = 0.0;
124  }
125  }
126  else if (rule==BURK_HERMITE||rule==BURK_GENZKEISTER) {
127  if (power%2) {
128  I = 0.0;
129  }
130  else {
131  long double value = 1.0;
132  if ((power-1)>=1) {
133  int n_copy = power-1;
134  while (1<n_copy) {
135  value *= (long double)n_copy;
136  n_copy -= 2;
137  }
138  }
139  I = value*sqrt(M_PI)/powl(2.0,(long double)power/2.0);
140  }
141  }
142  return I;
143 }
144 
145 int main(int argc, char *argv[]) {
146 
147  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
148 
149  // This little trick lets us print to std::cout only if
150  // a (dummy) command-line argument is provided.
151  int iprint = argc - 1;
152  Teuchos::RCP<std::ostream> outStream;
153  Teuchos::oblackholestream bhs; // outputs nothing
154  if (iprint > 0)
155  outStream = Teuchos::rcp(&std::cout, false);
156  else
157  outStream = Teuchos::rcp(&bhs, false);
158 
159  // Save the format state of the original std::cout.
160  Teuchos::oblackholestream oldFormatState;
161  oldFormatState.copyfmt(std::cout);
162 
163  *outStream \
164  << "===============================================================================\n" \
165  << "| |\n" \
166  << "| Unit Test (CubatureLineSorted) |\n" \
167  << "| |\n" \
168  << "| 1) Computing integrals of monomials in 1D |\n" \
169  << "| |\n" \
170  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
171  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
172  << "| |\n" \
173  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
174  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
175  << "| |\n" \
176  << "===============================================================================\n"\
177  << "| TEST 11: integrals of monomials in 1D |\n"\
178  << "===============================================================================\n";
179 
180  // internal variables:
181  int errorFlag = 0;
182  long double reltol = 1.0e+05*INTREPID_TOL;
183  int maxDeg = 0;
184  long double analyticInt = 0;
185  long double testInt = 0;
186  int maxOrder = 30;
187 
188  *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
189  // compute and compare integrals
190  try {
191  for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1; rule <= BURK_LAGUERRE; rule++) {
192  *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
193  // compute integrals
194  if (rule==BURK_HERMITE)
195  maxOrder = 10;
196  else if (rule==BURK_TRAPEZOIDAL)
197  maxOrder = 2;
198  else
199  maxOrder = 30;
200 
201  if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER) {
202  for (int i=1; i <= maxOrder; i++) {
203  if ( rule==BURK_CHEBYSHEV1 ||
204  rule==BURK_CHEBYSHEV2 ||
205  rule==BURK_LEGENDRE ||
206  rule==BURK_LAGUERRE ||
207  rule==BURK_HERMITE )
208  maxDeg = 2*i-1;
209  else if ( rule==BURK_CLENSHAWCURTIS ||
210  rule==BURK_FEJER2 ||
211  rule==BURK_TRAPEZOIDAL )
212  maxDeg = i-1;
213 
214  for (int j=0; j <= maxDeg; j++) {
215  analyticInt = evalInt(j,rule);
216  testInt = evalQuad(i,j,rule);
217 
218  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
219  long double absdiff = std::fabs(analyticInt - testInt);
220  *outStream << "Cubature order " << std::setw(2) << std::left
221  << i << " integrating "
222  << "x^" << std::setw(2) << std::left << j << ":"
223  << " "
224  << std::scientific << std::setprecision(16) << testInt
225  << " " << analyticInt
226  << " " << std::setprecision(4) << absdiff << " "
227  << "<?" << " " << abstol
228  << "\n";
229  if (absdiff > abstol) {
230  errorFlag++;
231  *outStream << std::right << std::setw(104)
232  << "^^^^---FAILURE!\n";
233  }
234  } // end for j
235  *outStream << "\n";
236  } // end for i
237  }
238  else if (rule==BURK_PATTERSON) {
239  for (int i=0; i < 8; i++) {
240  int l = (int)std::pow(2.0,(double)i+1.0)-1;
241  if (i==0)
242  maxDeg = 1;
243  else
244  maxDeg = (int)(1.5*(double)l+0.5);
245  for (int j=0; j <= maxDeg; j++) {
246  analyticInt = evalInt(j,rule);
247  testInt = evalQuad(l,j,rule);
248 
249  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
250  long double absdiff = std::fabs(analyticInt - testInt);
251  *outStream << "Cubature order " << std::setw(2) << std::left
252  << l << " integrating "
253  << "x^" << std::setw(2) << std::left << j << ":"
254  << " "
255  << std::scientific << std::setprecision(16) << testInt
256  << " " << analyticInt
257  << " " << std::setprecision(4) << absdiff << " "
258  << "<?" << " " << abstol
259  << "\n";
260  if (absdiff > abstol) {
261  errorFlag++;
262  *outStream << std::right << std::setw(104)
263  << "^^^^---FAILURE!\n";
264  }
265  } // end for j
266  *outStream << "\n";
267  } // end for i
268  }
269  else if (rule==BURK_GENZKEISTER) {
270  reltol *= 1.0e+02;
271  int o_ghk[4] = {1,3, 9,19};
272  int p_ghk[4] = {1,5,15,29};
273  for (int i=0; i < 4; i++) {
274  int l = o_ghk[i];
275  maxDeg = p_ghk[i];
276  for (int j=0; j <= maxDeg; j++) {
277  analyticInt = evalInt(j,rule);
278  testInt = evalQuad(l,j,rule);
279 
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
283  << l << " integrating "
284  << "x^" << std::setw(2) << std::left << j << ":"
285  << " "
286  << std::scientific << std::setprecision(16) << testInt
287  << " " << analyticInt
288  << " " << std::setprecision(4) << absdiff << " "
289  << "<?" << " " << abstol
290  << "\n";
291  if (absdiff > abstol) {
292  errorFlag++;
293  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
294  }
295  } // end for j
296  *outStream << "\n";
297  } // end for i
298  }
299  } // end for rule
300  }
301  catch (std::logic_error err) {
302  *outStream << err.what() << "\n";
303  errorFlag = -1;
304  };
305 
306 
307  if (errorFlag != 0)
308  std::cout << "End Result: TEST FAILED\n";
309  else
310  std::cout << "End Result: TEST PASSED\n";
311 
312  // reset format state of std::cout
313  std::cout.copyfmt(oldFormatState);
314 
315  return errorFlag;
316 }
Utilizes cubature (integration) rules contained in the library sandia_rules (John Burkardt...
Header file for the Intrepid::CubatureLineSorted class.
Intrepid utilities.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...