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"; \
31 int 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 );
164 Intrepid::TensorProductSpaceTools::evaluate<double,FieldContainer<double>,
FieldContainer<double>,FieldContainer<double> >( vals , coeff , basisVals );
166 FieldContainer<double> grads( 1 , 1, pts[0]->size() * pts[1]->size() , 2 );
168 Intrepid::TensorProductSpaceTools::evaluateGradient<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( grads ,
174 FieldContainer<double> fullVals(basesByDim[space_dim]->getCardinality(),
175 basesByDim[space_dim]->getCardinality());
176 FieldContainer<double> fullGrads(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";
228 FieldContainer<double> momentsNaive(1,basesByDim[2]->getCardinality());
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);
243 FieldContainer<double> momentsClever(1,1,basesByDim[space_dim]->getCardinality());
244 Intrepid::TensorProductSpaceTools::moments<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( momentsClever ,
248 for (
int j=0;j<momentsClever.dimension(0);j++)
250 if (std::abs( momentsClever(0,0,j) - momentsNaive(0,j) ) > Intrepid::INTREPID_TOL )
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);
Contains definitions of custom data types in Intrepid.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
An abstract base class that defines interface for bases that are tensor products of simpler bases...
Header file for the Intrepid::CubaturePolylib class.