55 #include "Teuchos_oblackholestream.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
64 long double evalQuad(std::vector<int> power,
65 int dimension, std::vector<int> order,
66 std::vector<EIntrepidBurkardt> rule) {
69 int size = lineCub.getNumPoints();
72 lineCub.getCubature(cubPoints,cubWeights);
78 for (
int i=0; i<dimension; i++) {
79 Q *= powl(cubPoints(mid,i),power[i]);
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]);
95 long double evalInt(
int dimension, std::vector<int> power,
96 std::vector<EIntrepidBurkardt> rule) {
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) {
106 I *= 2.0/((
long double)power[i]+1.0);
108 else if (rule[i]==BURK_LAGUERRE) {
109 I *= tgammal((
long double)(power[i]+1));
111 else if (rule[i]==BURK_CHEBYSHEV1) {
112 long double bot, top;
115 for (
int j=2;j<=power[i];j+=2) {
116 top *= (
long double)(j-1);
117 bot *= (
long double)j;
125 else if (rule[i]==BURK_CHEBYSHEV2) {
126 long double bot, top;
129 for (
int j=2;j<=power[i];j+=2) {
130 top *= (
long double)(j-1);
131 bot *= (
long double)j;
133 bot *= (
long double)(power[i]+2);
140 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
145 long double value = 1.0;
146 if ((power[i]-1)>=1) {
147 int n_copy = power[i]-1;
149 value *= (
long double)n_copy;
153 I *= value*sqrt(M_PI)/powl(2.0,(
long double)power[i]/2.0);
160 int main(
int argc,
char *argv[]) {
162 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
166 int iprint = argc - 1;
167 Teuchos::RCP<std::ostream> outStream;
168 Teuchos::oblackholestream bhs;
170 outStream = Teuchos::rcp(&std::cout,
false);
172 outStream = Teuchos::rcp(&bhs,
false);
175 Teuchos::oblackholestream oldFormatState;
176 oldFormatState.copyfmt(std::cout);
179 <<
"===============================================================================\n" \
181 <<
"| Unit Test (CubatureTensorSorted) |\n" \
183 <<
"| 1) Computing integrals of monomials in 2D |\n" \
185 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
186 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
188 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
189 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
191 <<
"===============================================================================\n"\
192 <<
"| TEST 12: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\
193 <<
"===============================================================================\n";
199 long double reltol = 1.0e+05*INTREPID_TOL;
201 long double analyticInt = 0;
202 long double testInt = 0;
204 std::vector<int> power(2,0);
205 std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
206 std::vector<int> order(2,0);
208 *outStream <<
"\nIntegrals of monomials on a reference line (edge):\n";
211 for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1;rule<=BURK_LAGUERRE;rule++) {
212 *outStream <<
"Testing " << EIntrepidBurkardtToString(rule) <<
"\n";
214 if (rule==BURK_HERMITE)
216 else if (rule==BURK_TRAPEZOIDAL)
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 ||
230 else if ( rule==BURK_CLENSHAWCURTIS ||
232 rule==BURK_TRAPEZOIDAL )
235 order[0] = i; order[1] = i;
236 for (
int j=0; j <= maxDeg; j++) {
238 for (
int k=0; k <= maxDeg; k++) {
240 analyticInt = evalInt(dimension, power, rule1);
241 testInt = evalQuad(power,dimension,order,rule1);
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) {
256 *outStream << std::right << std::setw(104)
257 <<
"^^^^---FAILURE!\n";
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;
271 maxDeg = (int)(1.5*(
double)l+0.5);
273 order[0] = l; order[1] = l;
274 for (
int j=0; j <= maxDeg; j++) {
276 for (
int k=0; k <= maxDeg; k++) {
278 analyticInt = evalInt(dimension, power, rule1);
279 testInt = evalQuad(power,dimension,order,rule1);
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) {
294 *outStream << std::right << std::setw(104)
295 <<
"^^^^---FAILURE!\n";
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++) {
310 maxDeg = (int)(1.5*(
double)l+0.5);
312 order[0] = l; order[1] = l;
313 for (
int j=0; j <= maxDeg; j++) {
315 for (
int k=0; k <= maxDeg; k++) {
317 analyticInt = evalInt(dimension, power, rule1);
318 testInt = evalQuad(power,dimension,order,rule1);
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) {
333 *outStream << std::right << std::setw(104)
334 <<
"^^^^---FAILURE!\n";
344 catch (std::logic_error err) {
345 *outStream << err.what() <<
"\n";
351 std::cout <<
"End Result: TEST FAILED\n";
353 std::cout <<
"End Result: TEST PASSED\n";
356 std::cout.copyfmt(oldFormatState);
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...