53#include "Teuchos_oblackholestream.hpp"
54#include "Teuchos_RCP.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
57#define INTREPID_CUBATURE_LINE_MAX 61
59using namespace Intrepid;
71 polydeg[0] = xDeg; polydeg[1] = yDeg; polydeg[2] = zDeg;
73 val *= std::pow(p(i),polydeg[i]);
82double computeIntegral(
int cubDegree,
int polyDegree, EIntrepidPLPoly poly_type) {
87 int cubDim = lineCub.getDimension();
89 int numCubPoints = lineCub.getNumPoints();
95 lineCub.getCubature(cubPoints, cubWeights);
97 for (
int i=0; i<numCubPoints; i++) {
98 for (
int j=0; j<cubDim; j++) {
99 point(j) = cubPoints(i,j);
101 val += computeMonomial(point, polyDegree)*cubWeights(i);
108int main(
int argc,
char *argv[]) {
110 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
114 int iprint = argc - 1;
115 Teuchos::RCP<std::ostream> outStream;
116 Teuchos::oblackholestream bhs;
118 outStream = Teuchos::rcp(&std::cout,
false);
120 outStream = Teuchos::rcp(&bhs,
false);
123 Teuchos::oblackholestream oldFormatState;
124 oldFormatState.copyfmt(std::cout);
127 <<
"===============================================================================\n" \
129 <<
"| Unit Test (CubaturePolylib) |\n" \
131 <<
"| 1) Computing integrals of monomials on reference cells in 1D |\n" \
133 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
134 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
136 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
137 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
139 <<
"===============================================================================\n"\
140 <<
"| TEST 1: integrals of monomials in 1D |\n"\
141 <<
"===============================================================================\n";
145 Teuchos::Array< Teuchos::Array<double> > testInt;
146 Teuchos::Array< Teuchos::Array<double> > analyticInt;
147 Teuchos::Array<double> tmparray(1);
148 double reltol = 1.0e+03 * INTREPID_TOL;
149 testInt.assign(INTREPID_CUBATURE_LINE_MAX+1, tmparray);
150 analyticInt.assign(INTREPID_CUBATURE_LINE_MAX+1, tmparray);
153 std::string basedir =
"./data";
154 std::stringstream namestream;
155 std::string filename;
156 namestream << basedir <<
"/EDGE_integrals" <<
".dat";
157 namestream >> filename;
158 std::ifstream filecompare(&filename[0]);
160 *outStream <<
"\nIntegrals of monomials on a reference line (edge):\n";
164 for (EIntrepidPLPoly poly_type=PL_GAUSS; poly_type <= PL_GAUSS_LOBATTO; poly_type++) {
166 for (
int cubDeg=0; cubDeg <= INTREPID_CUBATURE_LINE_MAX; cubDeg++) {
167 testInt[cubDeg].resize(cubDeg+1);
168 for (
int polyDeg=0; polyDeg <= cubDeg; polyDeg++) {
169 testInt[cubDeg][polyDeg] = computeIntegral(cubDeg, polyDeg, poly_type);
173 if (filecompare.is_open()) {
179 for (
int cubDeg=0; cubDeg <= INTREPID_CUBATURE_LINE_MAX; cubDeg++) {
180 for (
int polyDeg=0; polyDeg <= cubDeg; polyDeg++) {
181 double abstol = ( analyticInt[polyDeg][0] == 0.0 ? reltol : std::fabs(reltol*analyticInt[polyDeg][0]) );
182 double absdiff = std::fabs(analyticInt[polyDeg][0] - testInt[cubDeg][polyDeg]);
183 *outStream <<
"Cubature order " << std::setw(2) << std::left << cubDeg <<
" integrating "
184 <<
"x^" << std::setw(2) << std::left << polyDeg <<
":" <<
" "
185 << std::scientific << std::setprecision(16) << testInt[cubDeg][polyDeg] <<
" " << analyticInt[polyDeg][0] <<
" "
186 << std::setprecision(4) << absdiff <<
" " <<
"<?" <<
" " << abstol <<
"\n";
187 if (absdiff > abstol) {
189 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
196 catch (
const std::logic_error & err) {
197 *outStream << err.what() <<
"\n";
203 std::cout <<
"End Result: TEST FAILED\n";
205 std::cout <<
"End Result: TEST PASSED\n";
208 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubaturePolylib class.
void getAnalytic(Teuchos::Array< Teuchos::Array< Scalar > > &testMat, std::ifstream &inputFile, TypeOfExactData analyticDataType=INTREPID_UTILS_FRACTION)
Loads analytic values stored in a file into the matrix testMat, where the output matrix is an array o...
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int dimension(const int whichDim) const
Returns the specified dimension.