54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_GlobalMPISession.hpp"
58 using namespace Intrepid;
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) {
69 int size = lineCub.getNumPoints();
72 lineCub.getCubature(cubPoints,cubWeights);
80 int l1 = growthRule1D(order[0],growth[0],rule[0]);
81 int l2 = growthRule1D(order[1],growth[1],rule[1]);
85 for (
int i=0; i<l1; i++) {
86 locMid = i*l1+mid2; Qi = 0.0;
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;
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]);
101 Qi *= powl(cubPoints(locMid,0),power[0]);
128 long double evalInt(
int dimension, std::vector<int> power,
129 std::vector<EIntrepidBurkardt> rule) {
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) {
139 I *= 2.0/((
long double)power[i]+1.0);
141 else if (rule[i]==BURK_LAGUERRE) {
142 I *= tgammal((
long double)(power[i]+1));
144 else if (rule[i]==BURK_CHEBYSHEV1) {
145 long double bot, top;
148 for (
int j=2;j<=power[i];j+=2) {
149 top *= (
long double)(j-1);
150 bot *= (
long double)j;
158 else if (rule[i]==BURK_CHEBYSHEV2) {
159 long double bot, top;
162 for (
int j=2;j<=power[i];j+=2) {
163 top *= (
long double)(j-1);
164 bot *= (
long double)j;
166 bot *= (
long double)(power[i]+2);
173 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
178 long double value = 1.0;
179 if ((power[i]-1)>=1) {
180 int n_copy = power[i]-1;
182 value *= (
long double)n_copy;
186 I *= value*sqrt(M_PI)/powl(2.0,(
long double)power[i]/2.0);
193 int main(
int argc,
char *argv[]) {
195 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
199 int iprint = argc - 1;
200 Teuchos::RCP<std::ostream> outStream;
201 Teuchos::oblackholestream bhs;
203 outStream = Teuchos::rcp(&std::cout,
false);
205 outStream = Teuchos::rcp(&bhs,
false);
208 Teuchos::oblackholestream oldFormatState;
209 oldFormatState.copyfmt(std::cout);
212 <<
"===============================================================================\n" \
214 <<
"| Unit Test (CubatureTensorSorted) |\n" \
216 <<
"| 1) Computing integrals of monomials in 2D |\n" \
218 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
219 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
221 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
222 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
224 <<
"===============================================================================\n"\
225 <<
"| TEST 13: integrals of monomials in 2D - Anisotropic with growth rules |\n"\
226 <<
"===============================================================================\n";
232 long double reltol = 1.0e+06*INTREPID_TOL;
235 long double analyticInt = 0;
236 long double testInt = 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);
244 *outStream <<
"\nIntegrals of monomials on a reference line (edge):\n";
247 for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1;rule<=BURK_LAGUERRE;rule++) {
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 ) {
263 else if ( rule==BURK_CLENSHAWCURTIS ||
264 rule==BURK_FEJER2 ) {
269 order[0] = i; order[1] = i;
270 for (
int j=0; j <= maxDegx; j++) {
272 for (
int k=0; k <= maxDegy; k++) {
274 analyticInt = evalInt(dimension, power, rule1);
275 testInt = evalQuad(power,dimension,order,rule1,growth);
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
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) {
291 *outStream << std::right << std::setw(104)
292 <<
"^^^^---FAILURE!\n";
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);
310 maxDegx = (int)(1.5*(
double)l1-0.5);
311 maxDegy = (int)(1.5*(
double)l2-0.5);
314 order[0] = i; order[1] = i;
315 for (
int j=0; j <= maxDegx; j++) {
317 for (
int k=0; k <= maxDegy; k++) {
319 analyticInt = evalInt(dimension, power, rule1);
320 testInt = evalQuad(power,dimension,order,rule1,growth);
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) {
336 *outStream << std::right << std::setw(104)
337 <<
"^^^^---FAILURE!\n";
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);
355 maxDegx = (int)(1.5*(
double)l1-0.5);
356 maxDegy = (int)(1.5*(
double)l2-0.5);
359 order[0] = i; order[1] = i;
360 for (
int j=0; j <= maxDegx; j++) {
362 for (
int k=0; k <= maxDegy; k++) {
364 analyticInt = evalInt(dimension, power, rule1);
365 testInt = evalQuad(power,dimension,order,rule1,growth);
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
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) {
381 *outStream << std::right << std::setw(104)
382 <<
"^^^^---FAILURE!\n";
392 catch (std::logic_error err) {
393 *outStream << err.what() <<
"\n";
399 std::cout <<
"End Result: TEST FAILED\n";
401 std::cout <<
"End Result: TEST PASSED\n";
404 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...