Intrepid
test_13.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 
52 //#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(std::vector<int> power,
64  int dimension, std::vector<int> order,
65  std::vector<EIntrepidBurkardt> rule,
66  std::vector<EIntrepidGrowth> growth) {
67 
68  CubatureTensorSorted<long double> lineCub(dimension,order,rule,growth,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  // for (int i=0; i<size; i++) {
75  // std::cout << cubPoints(i,0) << " " << cubPoints(i,1) << std::endl;
76  // }
77 
78  long double Q = 0.0;
79  long double Qi = 0.0;
80  int l1 = growthRule1D(order[0],growth[0],rule[0]);
81  int l2 = growthRule1D(order[1],growth[1],rule[1]);
82  int mid2 = l2/2;
83  int locMid = 0;
84  int cnt = 0;
85  for (int i=0; i<l1; i++) {
86  locMid = i*l1+mid2; Qi = 0.0;
87  if (l2%2) {
88  Qi = cubWeights(locMid)*powl(cubPoints(locMid,1),power[1]); cnt++;
89  for (int j=1; j<=mid2; j++) {
90  Qi += cubWeights(locMid-j)*powl(cubPoints(locMid-j,1),power[1])
91  +cubWeights(locMid+j)*powl(cubPoints(locMid+j,1),power[1]); cnt += 2;
92  }
93  }
94  else {
95  for (int j=0; j<mid2; j++) {
96  Qi += cubWeights(locMid-j)*powl(cubPoints(locMid-j,1),power[1])
97  +cubWeights(locMid+j+1)*powl(cubPoints(locMid+j+1,1),power[1]);
98  cnt += 2;
99  }
100  }
101  Qi *= powl(cubPoints(locMid,0),power[0]);
102  Q += Qi;
103  }
104  return Q;
105  /*
106  int mid = size/2;
107  long double Q = 0.0;
108  if (size%2) {
109  Q = cubWeights(mid);
110  for (int i=0; i<dimension; i++) {
111  Q *= powl(cubPoints(mid,i),power[i]);
112  }
113  }
114 
115  for (int i=0; i<mid; i++) {
116  long double value1 = cubWeights(i);
117  long double value2 = cubWeights(size-i-1);
118  for (int j=0; j<dimension; j++) {
119  value1 *= powl(cubPoints(i,j),power[j]);
120  value2 *= powl(cubPoints(size-i-1,j),power[j]);
121  }
122  Q += value1+value2;
123  }
124  return Q;
125  */
126 }
127 
128 long double evalInt(int dimension, std::vector<int> power,
129  std::vector<EIntrepidBurkardt> rule) {
130  long double I = 1.0;
131 
132  for (int i=0; i<dimension; i++) {
133  if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2||
134  rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON ||
135  rule[i]==BURK_TRAPEZOIDAL) {
136  if (power[i]%2)
137  I *= 0.0;
138  else
139  I *= 2.0/((long double)power[i]+1.0);
140  }
141  else if (rule[i]==BURK_LAGUERRE) {
142  I *= tgammal((long double)(power[i]+1));
143  }
144  else if (rule[i]==BURK_CHEBYSHEV1) {
145  long double bot, top;
146  if (!(power[i]%2)) {
147  top = 1; bot = 1;
148  for (int j=2;j<=power[i];j+=2) {
149  top *= (long double)(j-1);
150  bot *= (long double)j;
151  }
152  I *= M_PI*top/bot;
153  }
154  else {
155  I *= 0.0;
156  }
157  }
158  else if (rule[i]==BURK_CHEBYSHEV2) {
159  long double bot, top;
160  if (!(power[i]%2)) {
161  top = 1; bot = 1;
162  for (int j=2;j<=power[i];j+=2) {
163  top *= (long double)(j-1);
164  bot *= (long double)j;
165  }
166  bot *= (long double)(power[i]+2);
167  I *= M_PI*top/bot;
168  }
169  else {
170  I *= 0.0;
171  }
172  }
173  else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
174  if (power[i]%2) {
175  I *= 0.0;
176  }
177  else {
178  long double value = 1.0;
179  if ((power[i]-1)>=1) {
180  int n_copy = power[i]-1;
181  while (1<n_copy) {
182  value *= (long double)n_copy;
183  n_copy -= 2;
184  }
185  }
186  I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0);
187  }
188  }
189  }
190  return I;
191 }
192 
193 int main(int argc, char *argv[]) {
194 
195  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
196 
197  // This little trick lets us print to std::cout only if
198  // a (dummy) command-line argument is provided.
199  int iprint = argc - 1;
200  Teuchos::RCP<std::ostream> outStream;
201  Teuchos::oblackholestream bhs; // outputs nothing
202  if (iprint > 0)
203  outStream = Teuchos::rcp(&std::cout, false);
204  else
205  outStream = Teuchos::rcp(&bhs, false);
206 
207  // Save the format state of the original std::cout.
208  Teuchos::oblackholestream oldFormatState;
209  oldFormatState.copyfmt(std::cout);
210 
211  *outStream \
212  << "===============================================================================\n" \
213  << "| |\n" \
214  << "| Unit Test (CubatureTensorSorted) |\n" \
215  << "| |\n" \
216  << "| 1) Computing integrals of monomials in 2D |\n" \
217  << "| |\n" \
218  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
219  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
220  << "| |\n" \
221  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
222  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
223  << "| |\n" \
224  << "===============================================================================\n"\
225  << "| TEST 13: integrals of monomials in 2D - Anisotropic with growth rules |\n"\
226  << "===============================================================================\n";
227 
228 
229  // internal variables:
230  int dimension = 2;
231  int errorFlag = 0;
232  long double reltol = 1.0e+06*INTREPID_TOL;
233  int maxDegx = 0;
234  int maxDegy = 0;
235  long double analyticInt = 0;
236  long double testInt = 0;
237  int maxOrder = 3;
238  int l1 = 0, l2 = 0;
239  std::vector<int> power(2,0);
240  std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
241  std::vector<int> order(2,0);
242  std::vector<EIntrepidGrowth> growth(2,GROWTH_FULLEXP);
243 
244  *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
245  // compute and compare integrals
246  try {
247  for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1;rule<=BURK_LAGUERRE;rule++) {
248  // compute integrals
249  rule1[0] = rule; rule1[1] = rule;
250  if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER&&rule!=BURK_TRAPEZOIDAL){
251  *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
252  for (int i=1; i <= maxOrder; i++) {
253  l1 = growthRule1D(i,growth[0],rule);
254  l2 = growthRule1D(i,growth[1],rule);
255  if ( rule==BURK_CHEBYSHEV1 ||
256  rule==BURK_CHEBYSHEV2 ||
257  rule==BURK_LEGENDRE ||
258  rule==BURK_LAGUERRE ||
259  rule==BURK_HERMITE ) {
260  maxDegx = 2*l1-1;
261  maxDegy = 2*l2-1;
262  }
263  else if ( rule==BURK_CLENSHAWCURTIS ||
264  rule==BURK_FEJER2 ) {
265  maxDegx = l1-1;
266  maxDegy = l2-1;
267  }
268 
269  order[0] = i; order[1] = i;
270  for (int j=0; j <= maxDegx; j++) {
271  power[0] = j;
272  for (int k=0; k <= maxDegy; k++) {
273  power[1] = k;
274  analyticInt = evalInt(dimension, power, rule1);
275  testInt = evalQuad(power,dimension,order,rule1,growth);
276 
277  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
278  long double absdiff = std::fabs(analyticInt - testInt);
279  *outStream << "Cubature order (" << std::setw(2) << std::left
280  << l1 << ", " << std::setw(2) << std::left << l2
281  << ") integrating "
282  << "x^" << std::setw(2) << std::left << j << "y^"
283  << std::setw(2) << std::left
284  << k << ":" << " " << std::scientific
285  << std::setprecision(16) << testInt
286  << " " << analyticInt << " "
287  << std::setprecision(4) << absdiff << " "
288  << "<?" << " " << abstol << "\n";
289  if (absdiff > abstol) {
290  errorFlag++;
291  *outStream << std::right << std::setw(104)
292  << "^^^^---FAILURE!\n";
293  }
294  } // end for k
295  *outStream << "\n";
296  } // end for j
297  *outStream << "\n";
298  } // end for i
299  }
300  else if (rule==BURK_PATTERSON) {
301  *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
302  for (int i=1; i <= 3; i++) {
303  l1 = growthRule1D(i,growth[0],rule);
304  l2 = growthRule1D(i,growth[1],rule);
305  if (i==0) {
306  maxDegx = 1;
307  maxDegy = 1;
308  }
309  else {
310  maxDegx = (int)(1.5*(double)l1-0.5);
311  maxDegy = (int)(1.5*(double)l2-0.5);
312  }
313 
314  order[0] = i; order[1] = i;
315  for (int j=0; j <= maxDegx; j++) {
316  power[0] = j;
317  for (int k=0; k <= maxDegy; k++) {
318  power[1] = k;
319  analyticInt = evalInt(dimension, power, rule1);
320  testInt = evalQuad(power,dimension,order,rule1,growth);
321 
322  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
323  long double absdiff = std::fabs(analyticInt - testInt);
324  *outStream << "Cubature order (" << std::setw(2) << std::left
325  << l1 << ", " << std::setw(2) << std::left
326  << l2 << ") integrating "
327  << "x^" << std::setw(2) << std::left << j
328  << "y^" << std::setw(2) << std::left
329  << k << ":" << " " << std::scientific
330  << std::setprecision(16) << testInt
331  << " " << analyticInt << " "
332  << std::setprecision(4) << absdiff << " "
333  << "<?" << " " << abstol << "\n";
334  if (absdiff > abstol) {
335  errorFlag++;
336  *outStream << std::right << std::setw(104)
337  << "^^^^---FAILURE!\n";
338  }
339  } // end for k
340  *outStream << "\n";
341  } // end for j
342  *outStream << "\n";
343  } // end for i
344  }
345  else if (rule==BURK_GENZKEISTER) {
346  *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
347  for (int i=1; i <= 3; i++) {
348  l1 = growthRule1D(i,growth[0],rule);
349  l2 = growthRule1D(i,growth[1],rule);
350  if (i==0) {
351  maxDegx = 1;
352  maxDegy = 1;
353  }
354  else {
355  maxDegx = (int)(1.5*(double)l1-0.5);
356  maxDegy = (int)(1.5*(double)l2-0.5);
357  }
358 
359  order[0] = i; order[1] = i;
360  for (int j=0; j <= maxDegx; j++) {
361  power[0] = j;
362  for (int k=0; k <= maxDegy; k++) {
363  power[1] = k;
364  analyticInt = evalInt(dimension, power, rule1);
365  testInt = evalQuad(power,dimension,order,rule1,growth);
366 
367  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
368  long double absdiff = std::fabs(analyticInt - testInt);
369  *outStream << "Cubature order (" << std::setw(2) << std::left
370  << l1 << ", " << std::setw(2) << std::left << l2
371  << ") integrating "
372  << "x^" << std::setw(2) << std::left << j << "y^"
373  << std::setw(2) << std::left
374  << k << ":" << " " << std::scientific
375  << std::setprecision(16) << testInt
376  << " " << analyticInt << " "
377  << std::setprecision(4) << absdiff << " "
378  << "<?" << " " << abstol << "\n";
379  if (absdiff > abstol) {
380  errorFlag++;
381  *outStream << std::right << std::setw(104)
382  << "^^^^---FAILURE!\n";
383  }
384  } // end for k
385  *outStream << "\n";
386  } // end for j
387  *outStream << "\n";
388  } // end for i
389  }
390  } // end for rule
391  }
392  catch (std::logic_error err) {
393  *outStream << err.what() << "\n";
394  errorFlag = -1;
395  };
396 
397 
398  if (errorFlag != 0)
399  std::cout << "End Result: TEST FAILED\n";
400  else
401  std::cout << "End Result: TEST PASSED\n";
402 
403  // reset format state of std::cout
404  std::cout.copyfmt(oldFormatState);
405 
406  return errorFlag;
407 }
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...