64int main(
int argc,
char *argv[]) {
66 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
69 int iprint = argc - 1;
71 Teuchos::RCP<std::ostream> outStream;
72 Teuchos::oblackholestream bhs;
75 outStream = Teuchos::rcp(&std::cout,
false);
77 outStream = Teuchos::rcp(&bhs,
false);
80 Teuchos::oblackholestream oldFormatState;
81 oldFormatState.copyfmt(std::cout);
84 <<
"===============================================================================\n" \
86 <<
"| Unit Test OrthogonalBases |\n" \
88 <<
"| 1) Tests orthogonality of triangular orthogonal basis (Dubiner) |\n" \
90 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
91 <<
"| Denis Ridzal (dridzal@sandia.gov) or |\n" \
92 <<
"| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
94 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
95 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
97 <<
"===============================================================================\n";
103 CubatureDirectTriDefault<double,FieldContainer<double> > myCub(20);
104 FieldContainer<double> cubPts( myCub.getNumPoints() , 2 );
105 FieldContainer<double> cubWts( myCub.getNumPoints() );
107 myCub.getCubature( cubPts , cubWts );
111 const int polydim = (deg+1)*(deg+2)/2;
112 FieldContainer<double> basisAtCubPts( polydim , myCub.getNumPoints() );
113 OrthogonalBases::tabulateTriangle<double,FieldContainer<double>,FieldContainer<double> >( cubPts , deg , basisAtCubPts );
116 for (
int i=0;i<polydim;i++) {
117 for (
int j=0;j<polydim;j++) {
119 for (
int k=0;k<myCub.getNumPoints();k++) {
120 cur += cubWts(k) * basisAtCubPts( i , k ) * basisAtCubPts( j , k );
122 if (i != j && fabs( cur ) > INTREPID_TOL) {
125 else if (i == j && fabs( cur ) < INTREPID_TOL ) {
134 std::cout <<
"End Result: TEST FAILED\n";
136 std::cout <<
"End Result: TEST PASSED\n";
139 std::cout.copyfmt(oldFormatState);