57#include "Shards_CellTopology.hpp"
58#include "Teuchos_oblackholestream.hpp"
59#include "Teuchos_RCP.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
62using namespace Intrepid;
64#define INTREPID_TEST_COMMAND( S ) \
69 catch (const std::logic_error & err) { \
70 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
71 *outStream << err.what() << '\n'; \
72 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
79double computeRefVolume(shards::CellTopology & cellTopology,
int cubDegree) {
80 Teuchos::RCP< Cubature<double> > myCub;
83 switch (cellTopology.getBaseCellTopologyData()->key) {
85 case shards::Line<>::key:
86 myCub = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
88 case shards::Triangle<>::key:
89 myCub = Teuchos::rcp(
new CubatureDirectTriDefault<double>(cubDegree));
91 case shards::Tetrahedron<>::key:
92 myCub = Teuchos::rcp(
new CubatureDirectTetDefault<double>(cubDegree));
94 case shards::Quadrilateral<>::key: {
95 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
96 lineCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
98 myCub = Teuchos::rcp(
new CubatureTensor<double>(lineCubs));
101 case shards::Hexahedron<>::key: {
102 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
103 std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(2);
104 lineCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
106 miscCubs[0] = Teuchos::rcp(
new CubatureTensor<double>(lineCubs));
108 myCub = Teuchos::rcp(
new CubatureTensor<double>(miscCubs));
111 case shards::Wedge<>::key: {
112 Teuchos::RCP<CubatureDirect<double> > triCub = Teuchos::rcp(
new CubatureDirectTriDefault<double>(cubDegree));
113 Teuchos::RCP<CubatureDirect<double> > lineCub = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
114 myCub = Teuchos::rcp(
new CubatureTensor<double>(triCub,lineCub));
117 case shards::Pyramid<>::key: {
118 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(3);
119 lineCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
120 lineCubs[1] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(cubDegree));
121 lineCubs[2] = Teuchos::rcp(
new CubatureDirectLineGaussJacobi20<double>(cubDegree));
122 myCub = Teuchos::rcp(
new CubatureTensorPyr<double>(lineCubs));
127 TEUCHOS_TEST_FOR_EXCEPTION( ( (cellTopology.getBaseCellTopologyData()->key != shards::Line<>::key) ||
128 (cellTopology.getBaseCellTopologyData()->key != shards::Triangle<>::key) ||
129 (cellTopology.getBaseCellTopologyData()->key != shards::Tetrahedron<>::key) ||
130 (cellTopology.getBaseCellTopologyData()->key != shards::Quadrilateral<>::key) ||
131 (cellTopology.getBaseCellTopologyData()->key != shards::Hexahedron<>::key) ||
132 (cellTopology.getBaseCellTopologyData()->key != shards::Wedge<>::key) ||
133 (cellTopology.getBaseCellTopologyData()->key != shards::Pyramid<>::key) ),
134 std::invalid_argument,
135 ">>> ERROR (Unit Test -- Cubature -- Volume): Invalid cell type.");
138 int numCubPoints = myCub->getNumPoints();
140 FieldContainer<double> cubPoints( numCubPoints, myCub->getDimension() );
141 FieldContainer<double> cubWeights( numCubPoints );
142 myCub->getCubature(cubPoints, cubWeights);
144 for (
int i=0; i<numCubPoints; i++)
145 vol += cubWeights[i];
151int main(
int argc,
char *argv[]) {
153 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
156 int iprint = argc - 1;
157 Teuchos::RCP<std::ostream> outStream;
158 Teuchos::oblackholestream bhs;
160 outStream = Teuchos::rcp(&std::cout,
false);
162 outStream = Teuchos::rcp(&bhs,
false);
165 Teuchos::oblackholestream oldFormatState;
166 oldFormatState.copyfmt(std::cout);
169 <<
"===============================================================================\n" \
171 <<
"| Unit Test (CubatureDirect,CubatureTensor) |\n" \
173 <<
"| 1) Computing volumes of reference cells |\n" \
175 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
176 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
178 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
179 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
181 <<
"===============================================================================\n"\
182 <<
"| TEST 1: construction and basic functionality |\n"\
183 <<
"===============================================================================\n";
187 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
188 int endThrowNumber = beginThrowNumber + 7;
192 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(-1) );
194 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub;
195 std::string testName =
"INTREPID_CUBATURE_LINE_GAUSS";
196 std::string lineCubName = lineCub.getName();
197 *outStream <<
"\nComparing strings: " << testName <<
" and " << lineCubName <<
"\n\n";
198 TEUCHOS_TEST_FOR_EXCEPTION( (testName != lineCubName), std::logic_error,
"Name mismatch!" ) );
199 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub;
200 std::vector<int> accuracy;
201 lineCub.getAccuracy(accuracy);
202 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error,
"Check member getAccuracy!" ) );
203 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(55);
204 TEUCHOS_TEST_FOR_EXCEPTION( (lineCub.getNumPoints() != 28), std::logic_error,
"Check member getNumPoints!" ) );
205 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub;
206 TEUCHOS_TEST_FOR_EXCEPTION( (lineCub.getDimension() != 1),
208 "Check member dimension!" ) );
210 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(-1) );
212 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub;
213 std::string testName =
"INTREPID_CUBATURE_TRI_DEFAULT";
214 std::string triCubName = triCub.getName();
215 *outStream <<
"\nComparing strings: " << testName <<
" and " << triCubName <<
"\n\n";
216 TEUCHOS_TEST_FOR_EXCEPTION( (testName != triCubName), std::logic_error,
"Name mismatch!" ) );
217 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub;
218 std::vector<int> accuracy;
219 triCub.getAccuracy(accuracy);
220 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error,
"Check member getAccuracy!" ) );
221 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(17);
222 TEUCHOS_TEST_FOR_EXCEPTION( (triCub.getNumPoints() != 61), std::logic_error,
"Check member getNumPoints!" ) );
223 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub;
224 TEUCHOS_TEST_FOR_EXCEPTION( (triCub.getDimension() != 2),
226 "Check member dimension!" ) );
228 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(-1) );
230 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub;
231 std::string testName =
"INTREPID_CUBATURE_TET_DEFAULT";
232 std::string tetCubName = tetCub.getName();
233 *outStream <<
"\nComparing strings: " << testName <<
" and " << tetCubName <<
"\n\n";
234 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
235 TEUCHOS_TEST_FOR_EXCEPTION( (testName != tetCubName), std::logic_error,
"Name mismatch!" ) );
236 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub;
237 std::vector<int> accuracy;
238 tetCub.getAccuracy(accuracy);
239 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error,
"Check member getAccuracy!" ) );
240 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(17);
241 TEUCHOS_TEST_FOR_EXCEPTION( (tetCub.getNumPoints() != 495), std::logic_error,
"Check member getNumPoints!" ) );
242 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub;
243 TEUCHOS_TEST_FOR_EXCEPTION( (tetCub.getDimension() != 3),
245 "Check member getCellTopology!" ) );
248 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(0);
249 CubatureTensor<double> quadCub(lineCubs) );
250 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3);
251 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
252 lineCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(3));
253 lineCubs[1] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(16));
254 miscCubs[0] = Teuchos::rcp(
new CubatureTensor<double>(lineCubs));
255 miscCubs[1] = Teuchos::rcp(
new CubatureDirectTriDefault<double>);
256 miscCubs[2] = Teuchos::rcp(
new CubatureDirectTetDefault<double>(19));
257 CubatureTensor<double> tensorCub(miscCubs);
258 std::vector<int> a(4); a[0]=3; a[1]=16; a[2]=0; a[3]=19;
259 std::vector<int> atest(4);
260 tensorCub.getAccuracy(atest);
261 TEUCHOS_TEST_FOR_EXCEPTION( (a != atest), std::logic_error,
"Check member getAccuracy!" ) );
263 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
264 lineCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(15));
265 lineCubs[1] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(11));
266 CubatureTensor<double> tensorCub(lineCubs);
267 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getNumPoints() != 48), std::logic_error,
"Check member getNumPoints!" ) );
268 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3);
269 miscCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>);
270 miscCubs[1] = Teuchos::rcp(
new CubatureDirectTriDefault<double>);
271 miscCubs[2] = Teuchos::rcp(
new CubatureDirectTetDefault<double>);
272 CubatureTensor<double> tensorCub(miscCubs);
273 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 6), std::logic_error,
"Check member dimension!" ) );
274 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3);
275 miscCubs[0] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(3));
276 miscCubs[1] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(7));
277 miscCubs[2] = Teuchos::rcp(
new CubatureDirectLineGauss<double>(5));
278 CubatureTensor<double> tensorCub(miscCubs);
279 FieldContainer<double> points(tensorCub.getNumPoints(), tensorCub.getDimension());
280 FieldContainer<double> weights(tensorCub.getNumPoints());
281 tensorCub.getCubature(points, weights)
285 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > lineCub = Teuchos::rcp(
new CubatureDirectLineGauss<double>(15));
286 Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(
new CubatureDirectTriDefault<double>(12));
287 CubatureTensor<double> tensorCub(lineCub, triCub);
288 std::vector<int> a(2); a[0] = 15; a[1] = 12;
289 std::vector<int> atest(2);
290 tensorCub.getAccuracy(atest);
291 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 3) || (a != atest),
293 "Check constructormembers dimension and getAccuracy!" ) );
294 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > lineCub = Teuchos::rcp(
new CubatureDirectLineGauss<double>(15));
295 Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(
new CubatureDirectTriDefault<double>(12));
296 CubatureTensor<double> tensorCub(triCub, lineCub, triCub);
297 std::vector<int> a(3); a[0] = 12; a[1] = 15; a[2] = 12;
298 std::vector<int> atest(3);
299 tensorCub.getAccuracy(atest);
300 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 5) || (a != atest),
302 "Check constructor and members dimension and getAccuracy!" ) );
304 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(
new CubatureDirectTriDefault<double>(12));
305 CubatureTensor<double> tensorCub(triCub, 5);
306 std::vector<int> a(5); a[0] = 12; a[1] = 12; a[2] = 12; a[3] = 12; a[4] = 12;
307 std::vector<int> atest(5);
308 tensorCub.getAccuracy(atest);
309 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 10) || (a != atest),
311 "Check constructor and members dimension and getAccuracy!" ) );
312 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) {
316 catch (
const std::logic_error & err) {
317 *outStream << err.what() <<
"\n";
322 <<
"===============================================================================\n"\
323 <<
"| TEST 2: volume computations |\n"\
324 <<
"===============================================================================\n";
326 double testVol = 0.0;
327 double tol = 100.0 * INTREPID_TOL;
330 double volumeList[] = {0.0, 2.0, 1.0/2.0, 4.0, 1.0/6.0, 8.0, 1.0, 4.0/3.0, 32.0};
332 *outStream <<
"\nReference cell volumes:\n\n";
335 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
337 testVol = computeRefVolume(line, deg);
338 *outStream << std::setw(30) <<
"Line volume --> " << std::setw(10) << std::scientific << testVol <<
339 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[1]) <<
"\n";
340 if (std::abs(testVol - volumeList[1]) > tol) {
342 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
346 *outStream <<
"\n\n";
347 shards::CellTopology tri(shards::getCellTopologyData< shards::Triangle<> >());
349 testVol = computeRefVolume(tri, deg);
350 *outStream << std::setw(30) <<
"Triangle volume --> " << std::setw(10) << std::scientific << testVol <<
351 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[2]) <<
"\n";
352 if (std::abs(testVol - volumeList[2]) > tol) {
354 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
358 *outStream <<
"\n\n";
359 shards::CellTopology quad(shards::getCellTopologyData< shards::Quadrilateral<> >());
361 testVol = computeRefVolume(quad, deg);
362 *outStream << std::setw(30) <<
"Quadrilateral volume --> " << std::setw(10) << std::scientific << testVol <<
363 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[3]) <<
"\n";
364 if (std::abs(testVol - volumeList[3]) > tol) {
366 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
370 *outStream <<
"\n\n";
371 shards::CellTopology tet(shards::getCellTopologyData< shards::Tetrahedron<> >());
373 testVol = computeRefVolume(tet, deg);
374 *outStream << std::setw(30) <<
"Tetrahedron volume --> " << std::setw(10) << std::scientific << testVol <<
375 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[4]) <<
"\n";
376 if (std::abs(testVol - volumeList[4]) > tol) {
378 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
381 *outStream <<
"\n\n";
382 shards::CellTopology hex(shards::getCellTopologyData< shards::Hexahedron<> >());
384 testVol = computeRefVolume(hex, deg);
385 *outStream << std::setw(30) <<
"Hexahedron volume --> " << std::setw(10) << std::scientific << testVol <<
386 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[5]) <<
"\n";
387 if (std::abs(testVol - volumeList[5]) > tol) {
389 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
392 *outStream <<
"\n\n";
393 shards::CellTopology wedge(shards::getCellTopologyData< shards::Wedge<> >());
395 testVol = computeRefVolume(wedge, deg);
396 *outStream << std::setw(30) <<
"Wedge volume --> " << std::setw(10) << std::scientific << testVol <<
397 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[6]) <<
"\n";
398 if (std::abs(testVol - volumeList[6]) > tol) {
400 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
403 *outStream <<
"\n\n";
404 shards::CellTopology pyr(shards::getCellTopologyData< shards::Pyramid<> >());
406 testVol = computeRefVolume(pyr, deg);
407 *outStream << std::setw(30) <<
"Pyramid volume --> " << std::setw(10) << std::scientific << testVol <<
408 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[7]) <<
"\n";
409 if (std::abs(testVol - volumeList[7]) > tol) {
411 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
414 *outStream <<
"\n\n";
415 for (
int deg=0; deg<=20; deg++) {
416 Teuchos::RCP<CubatureDirectLineGauss<double> > lineCub = Teuchos::rcp(
new CubatureDirectLineGauss<double>(deg));
417 CubatureTensor<double> hypercubeCub(lineCub, 5);
418 int numCubPoints = hypercubeCub.getNumPoints();
419 FieldContainer<double> cubPoints( numCubPoints, hypercubeCub.getDimension() );
420 FieldContainer<double> cubWeights( numCubPoints );
421 hypercubeCub.getCubature(cubPoints, cubWeights);
423 for (
int i=0; i<numCubPoints; i++)
424 testVol += cubWeights[i];
425 *outStream << std::setw(30) <<
"5-D Hypercube volume --> " << std::setw(10) << std::scientific << testVol <<
426 std::setw(10) <<
"diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[8]) <<
"\n";
427 if (std::abs(testVol - volumeList[8])/std::abs(testVol) > tol) {
429 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
433 catch (
const std::logic_error & err) {
434 *outStream << err.what() <<
"\n";
440 std::cout <<
"End Result: TEST FAILED\n";
442 std::cout <<
"End Result: TEST PASSED\n";
445 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubatureDirectLineGaussJacobi20 class.
#define INTREPID_CUBATURE_LINE_GAUSSJACOBI20_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct line rule of the Gaus...
Header file for the Intrepid::CubatureDirectLineGauss class.
#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...
Header file for the Intrepid::CubatureDirectTetDefault class.
#define INTREPID_CUBATURE_TET_DEFAULT_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct tetrahedron rule of t...
Header file for the Intrepid::CubatureDirectTriDefault class.
#define INTREPID_CUBATURE_TRI_DEFAULT_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct triangle rule of the ...
Header file for the Intrepid::CubatureTensorPyr class.
Header file for the Intrepid::CubatureTensor class.