54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
64 long double evalQuad(std::vector<int> power,
int dimension,
int order,
65 std::vector<EIntrepidBurkardt> rule,
66 std::vector<EIntrepidGrowth> growth) {
74 int size = lineCub.getNumPoints();
77 lineCub.getCubature(cubPoints,cubWeights);
83 for (
int i=0; i<dimension; i++) {
84 Q *= powl(cubPoints(mid,i),(
long double)power[i]);
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]);
100 long double evalInt(
int dimension, std::vector<int> power,
101 std::vector<EIntrepidBurkardt> rule) {
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) {
111 I *= 2.0/((
long double)power[i]+1.0);
113 else if (rule[i]==BURK_LAGUERRE) {
114 I *= tgammal((
long double)(power[i]+1));
116 else if (rule[i]==BURK_CHEBYSHEV1) {
117 long double bot, top;
120 for (
int j=2;j<=power[i];j+=2) {
121 top *= (
long double)(j-1);
122 bot *= (
long double)j;
130 else if (rule[i]==BURK_CHEBYSHEV2) {
131 long double bot, top;
134 for (
int j=2;j<=power[i];j+=2) {
135 top *= (
long double)(j-1);
136 bot *= (
long double)j;
138 bot *= (
long double)(power[i]+2);
145 else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
150 long double value = 1.0;
151 if ((power[i]-1)>=1) {
152 int n_copy = power[i]-1;
154 value *= (
long double)n_copy;
158 I *= value*sqrt(M_PI)/powl(2.0,(
long double)power[i]/2.0);
165 int main(
int argc,
char *argv[]) {
167 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
171 int iprint = argc - 1;
172 Teuchos::RCP<std::ostream> outStream;
173 Teuchos::oblackholestream bhs;
175 outStream = Teuchos::rcp(&std::cout,
false);
177 outStream = Teuchos::rcp(&bhs,
false);
180 Teuchos::oblackholestream oldFormatState;
181 oldFormatState.copyfmt(std::cout);
184 <<
"===============================================================================\n" \
186 <<
"| Unit Test (CubatureTensorSorted) |\n" \
188 <<
"| 1) Computing integrals of monomials in 2D |\n" \
190 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
191 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
193 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
194 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
196 <<
"===============================================================================\n"\
197 <<
"| TEST 22: integrals of monomials in 2D - Anisotropic but no growth rules |\n"\
198 <<
"===============================================================================\n";
204 long double reltol = 1.0e+02*INTREPID_TOL;
206 long double analyticInt = 0;
207 long double testInt = 0;
209 std::vector<int> power(2,0);
210 std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
211 std::vector<EIntrepidGrowth> growth1(2,GROWTH_FULLEXP);
213 *outStream <<
"\nIntegrals of monomials on a reference line (edge):\n";
216 for (
int i=0; i<=maxOrder; i++) {
218 for (
int j=0; j <= maxDeg; j++) {
220 for (
int k=0; k <= maxDeg; k++) {
223 analyticInt = evalInt(dimension, power, rule1);
224 testInt = evalQuad(power,dimension,maxOrder,rule1,growth1);
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) {
238 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
247 catch (std::logic_error err) {
248 *outStream << err.what() <<
"\n";
254 std::cout <<
"End Result: TEST FAILED\n";
256 std::cout <<
"End Result: TEST PASSED\n";
259 std::cout.copyfmt(oldFormatState);
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Header file for the Intrepid::AdaptiveSparseGrid class.
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...