53 #include "Teuchos_oblackholestream.hpp"
54 #include "Teuchos_RCP.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
57 using namespace Intrepid;
70 polydeg[0] = xDeg; polydeg[1] = yDeg; polydeg[2] = zDeg;
72 val *= std::pow(p(i),polydeg[i]);
81 double computeIntegral(
int cubDegree,
int xDeg,
int yDeg) {
86 int cubDim = myCub.getDimension();
87 int numCubPoints = myCub.getNumPoints();
93 myCub.getCubature(cubPoints, cubWeights);
95 for (
int i=0; i<numCubPoints; i++) {
96 for (
int j=0; j<cubDim; j++) {
97 point(j) = cubPoints(i,j);
99 val += computeMonomial(point, xDeg, yDeg)*cubWeights(i);
106 int main(
int argc,
char *argv[]) {
108 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
112 int iprint = argc - 1;
113 Teuchos::RCP<std::ostream> outStream;
114 Teuchos::oblackholestream bhs;
116 outStream = Teuchos::rcp(&std::cout,
false);
118 outStream = Teuchos::rcp(&bhs,
false);
121 Teuchos::oblackholestream oldFormatState;
122 oldFormatState.copyfmt(std::cout);
125 <<
"===============================================================================\n" \
127 <<
"| Unit Test (CubatureSparse) |\n" \
129 <<
"| 1) Computing integrals of monomials on reference cells in 2D |\n" \
131 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
132 <<
"| Denis Ridzal (dridzal@sandia.gov) or |\n" \
133 <<
"| Matthew Keegan (mskeega@sandia.gov) |\n" \
135 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
136 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
138 <<
"===============================================================================\n"\
139 <<
"| TEST 1: integrals of monomials in 2D for Sparse Grid Construction |\n"\
140 <<
"===============================================================================\n";
146 Teuchos::Array< Teuchos::Array<double> > testInt;
147 Teuchos::Array< Teuchos::Array<double> > analyticInt;
148 Teuchos::Array<double> tmparray(1);
149 double reltol = 1.0e+03 * INTREPID_TOL;
152 int numPoly = (maxDeg+1)*(maxDeg+2)/2;
153 int numAnalytic = (maxOffset+1)*(maxOffset+2)/2;
154 testInt.assign(numPoly, tmparray);
155 analyticInt.assign(numAnalytic, tmparray);
158 std::string basedir =
"./data";
159 std::stringstream namestream;
160 std::string filename;
161 namestream << basedir <<
"/QUAD_integrals" <<
".dat";
162 namestream >> filename;
166 *outStream <<
"\nIntegrals of monomials:\n";
168 std::ifstream filecompare(&filename[0]);
170 for (
int cubDeg=0; cubDeg <= maxDeg; cubDeg++) {
172 testInt[cubDeg].resize((cubDeg+1)*(cubDeg+2)/2);
173 for (
int xDeg=0; xDeg <= cubDeg; xDeg++) {
174 for (
int yDeg=0; yDeg <= cubDeg-xDeg; yDeg++) {
175 testInt[cubDeg][polyCt] = computeIntegral(cubDeg, xDeg, yDeg);
182 if (filecompare.is_open()) {
183 getAnalytic(analyticInt, filecompare);
188 for (
int cubDeg=0; cubDeg <= maxDeg; cubDeg++) {
191 for (
int xDeg=0; xDeg <= cubDeg; xDeg++) {
192 for (
int yDeg=0; yDeg <= cubDeg-xDeg; yDeg++) {
193 double abstol = ( analyticInt[polyCt+offset][0] == 0.0 ? reltol : std::fabs(reltol*analyticInt[polyCt+offset][0]) );
194 double absdiff = std::fabs(analyticInt[polyCt+offset][0] - testInt[cubDeg][polyCt]);
195 *outStream <<
"Cubature order " << std::setw(2) << std::left << cubDeg <<
" integrating "
196 <<
"x^" << std::setw(2) << std::left << xDeg <<
" * y^" << std::setw(2) << yDeg <<
":" <<
" "
197 << std::scientific << std::setprecision(16) << testInt[cubDeg][polyCt] <<
" " << analyticInt[polyCt+offset][0] <<
" "
198 << std::setprecision(4) << absdiff <<
" " <<
"<?" <<
" " << abstol <<
"\n";
199 if (absdiff > abstol) {
201 *outStream << std::right << std::setw(111) <<
"^^^^---FAILURE!\n";
205 offset = offset + maxOffset - cubDeg;
211 catch (std::logic_error err) {
212 *outStream << err.what() <<
"\n";
218 std::cout <<
"End Result: TEST FAILED\n";
220 std::cout <<
"End Result: TEST PASSED\n";
223 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubatureSparse class.
int dimension(const int whichDim) const
Returns the specified dimension.
#define INTREPID_CUBATURE_LINE_GAUSS_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct line rule of the Gaus...