Intrepid
test_22.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 
53 #include "Intrepid_Utils.hpp"
54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.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, int dimension, int order,
65  std::vector<EIntrepidBurkardt> rule,
66  std::vector<EIntrepidGrowth> growth) {
67 
68  CubatureTensorSorted<long double> lineCub(0,dimension);
70  lineCub,dimension,
71  order,rule,
72  growth,false);
73 
74  int size = lineCub.getNumPoints();
75  FieldContainer<long double> cubPoints(size,dimension);
76  FieldContainer<long double> cubWeights(size);
77  lineCub.getCubature(cubPoints,cubWeights);
78 
79  int mid = size/2;
80  long double Q = 0.0;
81  if (size%2) {
82  Q = cubWeights(mid);
83  for (int i=0; i<dimension; i++) {
84  Q *= powl(cubPoints(mid,i),(long double)power[i]);
85  }
86  }
87 
88  for (int i=0; i<mid; i++) {
89  long double value1 = cubWeights(i);
90  long double value2 = cubWeights(size-i-1);
91  for (int j=0; j<dimension; j++) {
92  value1 *= powl(cubPoints(i,j),(long double)power[j]);
93  value2 *= powl(cubPoints(size-i-1,j),(long double)power[j]);
94  }
95  Q += value1+value2;
96  }
97  return Q;
98 }
99 
100 long double evalInt(int dimension, std::vector<int> power,
101  std::vector<EIntrepidBurkardt> rule) {
102  long double I = 1.0;
103 
104  for (int i=0; i<dimension; i++) {
105  if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2||
106  rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON ||
107  rule[i]==BURK_TRAPEZOIDAL) {
108  if (power[i]%2)
109  I *= 0.0;
110  else
111  I *= 2.0/((long double)power[i]+1.0);
112  }
113  else if (rule[i]==BURK_LAGUERRE) {
114  I *= tgammal((long double)(power[i]+1));
115  }
116  else if (rule[i]==BURK_CHEBYSHEV1) {
117  long double bot, top;
118  if (!(power[i]%2)) {
119  top = 1; bot = 1;
120  for (int j=2;j<=power[i];j+=2) {
121  top *= (long double)(j-1);
122  bot *= (long double)j;
123  }
124  I *= M_PI*top/bot;
125  }
126  else {
127  I *= 0.0;
128  }
129  }
130  else if (rule[i]==BURK_CHEBYSHEV2) {
131  long double bot, top;
132  if (!(power[i]%2)) {
133  top = 1; bot = 1;
134  for (int j=2;j<=power[i];j+=2) {
135  top *= (long double)(j-1);
136  bot *= (long double)j;
137  }
138  bot *= (long double)(power[i]+2);
139  I *= M_PI*top/bot;
140  }
141  else {
142  I *= 0.0;
143  }
144  }
145  else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
146  if (power[i]%2) {
147  I *= 0.0;
148  }
149  else {
150  long double value = 1.0;
151  if ((power[i]-1)>=1) {
152  int n_copy = power[i]-1;
153  while (1<n_copy) {
154  value *= (long double)n_copy;
155  n_copy -= 2;
156  }
157  }
158  I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0);
159  }
160  }
161  }
162  return I;
163 }
164 
165 int main(int argc, char *argv[]) {
166 
167  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
168 
169  // This little trick lets us print to std::cout only if
170  // a (dummy) command-line argument is provided.
171  int iprint = argc - 1;
172  Teuchos::RCP<std::ostream> outStream;
173  Teuchos::oblackholestream bhs; // outputs nothing
174  if (iprint > 0)
175  outStream = Teuchos::rcp(&std::cout, false);
176  else
177  outStream = Teuchos::rcp(&bhs, false);
178 
179  // Save the format state of the original std::cout.
180  Teuchos::oblackholestream oldFormatState;
181  oldFormatState.copyfmt(std::cout);
182 
183  *outStream \
184  << "===============================================================================\n" \
185  << "| |\n" \
186  << "| Unit Test (CubatureTensorSorted) |\n" \
187  << "| |\n" \
188  << "| 1) Computing integrals of monomials in 2D |\n" \
189  << "| |\n" \
190  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
191  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
192  << "| |\n" \
193  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
194  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
195  << "| |\n" \
196  << "===============================================================================\n"\
197  << "| TEST 22: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\
198  << "===============================================================================\n";
199 
200 
201  // internal variables:
202  int dimension = 2;
203  int errorFlag = 0;
204  long double reltol = 1.0e+02*INTREPID_TOL;
205  int maxDeg = 0;
206  long double analyticInt = 0;
207  long double testInt = 0;
208  int maxOrder = 6;
209  std::vector<int> power(2,0);
210  std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
211  std::vector<EIntrepidGrowth> growth1(2,GROWTH_FULLEXP);
212 
213  *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
214  // compute and compare integrals
215  try {
216  for (int i=0; i<=maxOrder; i++) {
217  maxDeg = i-1;
218  for (int j=0; j <= maxDeg; j++) {
219  power[0] = j;
220  for (int k=0; k <= maxDeg; k++) {
221  power[1] = k;
222  if (j+k < maxDeg) {
223  analyticInt = evalInt(dimension, power, rule1);
224  testInt = evalQuad(power,dimension,maxOrder,rule1,growth1);
225 
226  long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
227  long double absdiff = std::fabs(analyticInt - testInt);
228  *outStream << "Cubature order " << std::setw(2) << std::left << i
229  << " integrating " << "x^" << std::setw(2)
230  << std::left << j << "y^" << std::setw(2) << std::left
231  << k << ":" << " " << std::scientific
232  << std::setprecision(16) << testInt
233  << " " << analyticInt << " "
234  << std::setprecision(4) << absdiff << " "
235  << "<?" << " " << abstol << "\n";
236  if (absdiff > abstol) {
237  errorFlag++;
238  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
239  }
240  }
241  } // end for k
242  *outStream << "\n";
243  } // end for j
244  *outStream << "\n";
245  } // end for i
246  }
247  catch (std::logic_error err) {
248  *outStream << err.what() << "\n";
249  errorFlag = -1;
250  };
251 
252 
253  if (errorFlag != 0)
254  std::cout << "End Result: TEST FAILED\n";
255  else
256  std::cout << "End Result: TEST PASSED\n";
257 
258  // reset format state of std::cout
259  std::cout.copyfmt(oldFormatState);
260 
261  return errorFlag;
262 }
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Header file for the Intrepid::AdaptiveSparseGrid class.
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...