Intrepid
test_01.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44
57#include "Shards_CellTopology.hpp"
58#include "Teuchos_oblackholestream.hpp"
59#include "Teuchos_RCP.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
61
62using namespace Intrepid;
63
64#define INTREPID_TEST_COMMAND( S ) \
65{ \
66 try { \
67 S ; \
68 } \
69 catch (const std::logic_error & err) { \
70 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
71 *outStream << err.what() << '\n'; \
72 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
73 }; \
74}
75
76/*
77 Computes volume of a given reference cell.
78*/
79double computeRefVolume(shards::CellTopology & cellTopology, int cubDegree) {
80 Teuchos::RCP< Cubature<double> > myCub;
81 double vol = 0.0;
82
83 switch (cellTopology.getBaseCellTopologyData()->key) {
84
85 case shards::Line<>::key:
86 myCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree));
87 break;
88 case shards::Triangle<>::key:
89 myCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(cubDegree));
90 break;
91 case shards::Tetrahedron<>::key:
92 myCub = Teuchos::rcp(new CubatureDirectTetDefault<double>(cubDegree));
93 break;
94 case shards::Quadrilateral<>::key: {
95 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2);
96 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree));
97 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+7) % INTREPID_CUBATURE_LINE_GAUSS_MAX));
98 myCub = Teuchos::rcp(new CubatureTensor<double>(lineCubs));
99 }
100 break;
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));
105 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+7) % INTREPID_CUBATURE_LINE_GAUSS_MAX));
106 miscCubs[0] = Teuchos::rcp(new CubatureTensor<double>(lineCubs));
107 miscCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+25) % INTREPID_CUBATURE_LINE_GAUSS_MAX));
108 myCub = Teuchos::rcp(new CubatureTensor<double>(miscCubs));
109 }
110 break;
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));
115 }
116 break;
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));
123 }
124 break;
125
126 default:
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.");
136 } // end switch
137
138 int numCubPoints = myCub->getNumPoints();
139
140 FieldContainer<double> cubPoints( numCubPoints, myCub->getDimension() );
141 FieldContainer<double> cubWeights( numCubPoints );
142 myCub->getCubature(cubPoints, cubWeights);
143
144 for (int i=0; i<numCubPoints; i++)
145 vol += cubWeights[i];
146
147 return vol;
148}
149
150
151int main(int argc, char *argv[]) {
152
153 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
154 // This little trick lets us print to std::cout only if
155 // a (dummy) command-line argument is provided.
156 int iprint = argc - 1;
157 Teuchos::RCP<std::ostream> outStream;
158 Teuchos::oblackholestream bhs; // outputs nothing
159 if (iprint > 0)
160 outStream = Teuchos::rcp(&std::cout, false);
161 else
162 outStream = Teuchos::rcp(&bhs, false);
163
164 // Save the format state of the original std::cout.
165 Teuchos::oblackholestream oldFormatState;
166 oldFormatState.copyfmt(std::cout);
167
168 *outStream \
169 << "===============================================================================\n" \
170 << "| |\n" \
171 << "| Unit Test (CubatureDirect,CubatureTensor) |\n" \
172 << "| |\n" \
173 << "| 1) Computing volumes of reference cells |\n" \
174 << "| |\n" \
175 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
176 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
177 << "| |\n" \
178 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
179 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
180 << "| |\n" \
181 << "===============================================================================\n"\
182 << "| TEST 1: construction and basic functionality |\n"\
183 << "===============================================================================\n";
184
185 int errorFlag = 0;
186
187 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
188 int endThrowNumber = beginThrowNumber + 7;
189
190 try {
191 /* Line cubature. */
192 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(-1) );
193 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(INTREPID_CUBATURE_LINE_GAUSS_MAX+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),
207 std::logic_error,
208 "Check member dimension!" ) );
209 /* Triangle cubature. */
210 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(-1) );
211 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(INTREPID_CUBATURE_TRI_DEFAULT_MAX+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),
225 std::logic_error,
226 "Check member dimension!" ) );
227 /* Tetrahedron cubature. */
228 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(-1) );
229 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(INTREPID_CUBATURE_TET_DEFAULT_MAX+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),
244 std::logic_error,
245 "Check member getCellTopology!" ) );
246
247 /* Tensor cubature. */
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!" ) );
262
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)
282 )
283
284
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),
292 std::logic_error,
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),
301 std::logic_error,
302 "Check constructor and members dimension and getAccuracy!" ) );
303
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),
310 std::logic_error,
311 "Check constructor and members dimension and getAccuracy!" ) );
312 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) {
313 errorFlag = -1000;
314 }
315 }
316 catch (const std::logic_error & err) {
317 *outStream << err.what() << "\n";
318 errorFlag = -1000;
319 };
320
321 *outStream \
322 << "===============================================================================\n"\
323 << "| TEST 2: volume computations |\n"\
324 << "===============================================================================\n";
325
326 double testVol = 0.0;
327 double tol = 100.0 * INTREPID_TOL;
328
329 // list of analytic volume values, listed in the enumerated reference cell order up to CELL_HEXAPRISM
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};
331
332 *outStream << "\nReference cell volumes:\n\n";
333
334 try {
335 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
336 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) {
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) {
341 errorFlag = 1;
342 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
343 }
344 }
345
346 *outStream << "\n\n";
347 shards::CellTopology tri(shards::getCellTopologyData< shards::Triangle<> >());
348 for (int deg=0; deg<=INTREPID_CUBATURE_TRI_DEFAULT_MAX; deg++) {
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) {
353 errorFlag = 1;
354 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
355 }
356 }
357
358 *outStream << "\n\n";
359 shards::CellTopology quad(shards::getCellTopologyData< shards::Quadrilateral<> >());
360 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) {
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) {
365 errorFlag = 1;
366 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
367 }
368 }
369
370 *outStream << "\n\n";
371 shards::CellTopology tet(shards::getCellTopologyData< shards::Tetrahedron<> >());
372 for (int deg=0; deg<=INTREPID_CUBATURE_TET_DEFAULT_MAX; deg++) {
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) {
377 errorFlag = 1;
378 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
379 }
380 }
381 *outStream << "\n\n";
382 shards::CellTopology hex(shards::getCellTopologyData< shards::Hexahedron<> >());
383 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) {
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) {
388 errorFlag = 1;
389 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
390 }
391 }
392 *outStream << "\n\n";
393 shards::CellTopology wedge(shards::getCellTopologyData< shards::Wedge<> >());
394 for (int deg=0; deg<=std::min(INTREPID_CUBATURE_LINE_GAUSS_MAX,INTREPID_CUBATURE_TRI_DEFAULT_MAX); deg++) {
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) {
399 errorFlag = 1;
400 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
401 }
402 }
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) {
410 errorFlag = 1;
411 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
412 }
413 }
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);
422 testVol = 0;
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) {
428 errorFlag = 1;
429 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
430 }
431 }
432 }
433 catch (const std::logic_error & err) {
434 *outStream << err.what() << "\n";
435 errorFlag = -1;
436 };
437
438
439 if (errorFlag != 0)
440 std::cout << "End Result: TEST FAILED\n";
441 else
442 std::cout << "End Result: TEST PASSED\n";
443
444 // reset format state of std::cout
445 std::cout.copyfmt(oldFormatState);
446
447 return errorFlag;
448}
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.