1#include "Teuchos_Array.hpp"
2#include "Teuchos_RCP.hpp"
3#include "Teuchos_oblackholestream.hpp"
4#include "Teuchos_GlobalMPISession.hpp"
18#define INTREPID_TEST_COMMAND( S ) \
23 catch (const std::logic_error & err) { \
24 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
25 *outStream << err.what() << '\n'; \
26 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
31int main(
int argc ,
char **argv )
33 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
36 int iprint = argc - 1;
38 Teuchos::RCP<std::ostream> outStream;
39 Teuchos::oblackholestream bhs;
42 outStream = Teuchos::rcp(&std::cout,
false);
44 outStream = Teuchos::rcp(&bhs,
false);
47 Teuchos::oblackholestream oldFormatState;
48 oldFormatState.copyfmt(std::cout);
52 <<
"===============================================================================\n" \
54 <<
"| Unit Test TensorProductSpace Tools |\n" \
56 <<
"| Tests sum-factored polynomial evaluation and integration |\n" \
58 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
59 <<
"| Denis Ridzal (dridzal@sandia.gov) or |\n" \
60 <<
"| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
62 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
63 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
65 <<
"===============================================================================\n";
71 Array<RCP<TensorBasis<double,FieldContainer<double> > > > basesByDim(4);
72 Array<RCP<FieldContainer<double> > > cubPtsByDim(4);
79 cpl.getCubature( cubPoints, cubWeights );
87 for (
int j=0;j<cpl.getNumPoints();j++)
89 for (
int i=0;i<cpl.getNumPoints();i++)
91 int index = j*cpl.getNumPoints() + i;
92 quadPts(index,0) = cubPoints(i,0);
93 quadPts(index,1) = cubPoints(j,0);
98 for (
int k=0;k<cpl.getNumPoints();k++)
100 for (
int j=0;j<cpl.getNumPoints();j++)
102 for (
int i=0;i<cpl.getNumPoints();i++)
104 int index = k* cpl.getNumPoints() * cpl.getNumPoints() + j*cpl.getNumPoints() + i;
105 cubPts(index,0) = cubPoints(i,0);
106 cubPts(index,1) = cubPoints(j,0);
107 cubPts(index,2) = cubPoints(k,0);
112 cubPtsByDim[2] = Teuchos::rcp( &quadPts ,
false );
113 cubPtsByDim[3] = Teuchos::rcp( &cubPts ,
false );
117 Array<Array<RCP<Basis<double,FieldContainer<double> > > > > &bases = basesByDim[space_dim]->getBases();
123 Array<RCP<FieldContainer<double> > > pts( space_dim );
124 pts[0] = Teuchos::rcp( &cubPoints,
false );
125 for (
int i=1;i<space_dim;i++)
130 Array<RCP<FieldContainer<double> > > wts(space_dim);
131 wts[0] = Teuchos::rcp( &cubWeights ,
false );
132 for (
int i=1;i<space_dim;i++)
138 cpl.getNumPoints() );
140 cpl.getNumPoints() );
142 cpl.getNumPoints(), 1 );
144 cpl.getNumPoints(), 1 );
146 bases[0][0]->getValues( Phix , cubPoints, Intrepid::OPERATOR_VALUE );
147 bases[0][1]->getValues( Phiy , cubPoints, Intrepid::OPERATOR_VALUE );
148 bases[0][0]->getValues( DPhix , cubPoints, Intrepid::OPERATOR_D1 );
149 bases[0][1]->getValues( DPhiy , cubPoints, Intrepid::OPERATOR_D1 );
151 Array<RCP<FieldContainer<double> > > basisVals(2);
152 basisVals[0] = Teuchos::rcp( &Phix ,
false );
153 basisVals[1] = Teuchos::rcp( &Phiy ,
false );
155 Array<RCP<FieldContainer<double> > > basisDVals(2);
156 basisDVals[0] = Teuchos::rcp( &DPhix ,
false );
157 basisDVals[1] = Teuchos::rcp( &DPhiy ,
false );
175 basesByDim[space_dim]->getCardinality());
177 basesByDim[space_dim]->getCardinality(),
181 basesByDim[space_dim]->getValues( fullVals ,
183 Intrepid::OPERATOR_VALUE );
184 basesByDim[space_dim]->getValues( fullGrads ,
186 Intrepid::OPERATOR_GRAD );
188 for (
int i=0;i<fullVals.dimension(1);i++)
190 if (std::abs( fullVals(0,i) - vals(0,0,i) ) > Intrepid::INTREPID_TOL )
193 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
196 *outStream <<
" Evaluating first bf At multi-index { ";
198 *outStream <<
"} brute force value: " << fullVals(0,i)
199 <<
" but tensor-product value: " << vals(0,i) <<
"\n";
200 *outStream <<
"Difference: " << std::abs( fullVals(0,i) - vals(0,i) ) <<
"\n";
204 for (
int i=0;i<fullGrads.dimension(1);i++)
206 for (
int j=0;j<fullGrads.dimension(2);j++)
208 if (std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) > Intrepid::INTREPID_TOL )
211 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
214 *outStream <<
" Evaluating first bf At multi-index { ";
215 *outStream << i <<
" " << j;
216 *outStream <<
"} brute force value: " << fullGrads(0,i,j)
217 <<
" but tensor-product value: " << grads(0,0,i,j) <<
"\n";
218 *outStream <<
"Difference: " << std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) <<
"\n";
229 for (
int i=0;i<basesByDim[2]->getCardinality();i++)
231 momentsNaive(0,i) = 0.0;
232 for (
int qpty=0;qpty<cubPoints.dimension(0);qpty++)
234 for (
int qptx=0;qptx<cubPoints.dimension(0);qptx++)
236 momentsNaive(0,i) += cubWeights(qpty) * cubWeights(qptx) *
237 vals( 0, 0, qpty*cubPoints.dimension(0)+qptx )
238 * fullVals(i,qpty*cubPoints.dimension(0)+qptx);
248 for (
int j=0;j<momentsClever.dimension(0);j++)
253 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
256 *outStream <<
" At multi-index { ";
257 *outStream <<
" " << j;
258 *outStream <<
"} brute force value: " << momentsNaive(0,j)
259 <<
" but sum-factored value: " << momentsClever(0,0,j) <<
"\n";
260 *outStream <<
"Difference: " << std::abs( momentsNaive(0,j) - momentsClever(0,0,j) ) <<
"\n";
266 std::cout <<
"End Result: TEST FAILED\n";
270 std::cout <<
"End Result: TEST PASSED\n";
274 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
Contains definitions of custom data types in Intrepid.
static const double INTREPID_TOL
General purpose tolerance in, e.g., internal Newton's method to invert ref to phys maps.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
An abstract base class that defines interface for bases that are tensor products of simpler bases.