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