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