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 
51 #include "Intrepid_CubatureDirectLineGauss.hpp"
52 #include "Intrepid_CubatureDirectLineGaussJacobi20.hpp"
53 #include "Intrepid_CubatureDirectTriDefault.hpp"
54 #include "Intrepid_CubatureDirectTetDefault.hpp"
55 #include "Intrepid_CubatureTensor.hpp"
56 #include "Intrepid_CubatureTensorPyr.hpp"
57 #include "Shards_CellTopology.hpp"
58 #include "Teuchos_oblackholestream.hpp"
59 #include "Teuchos_RCP.hpp"
60 #include "Teuchos_GlobalMPISession.hpp"
61 
62 using 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 */
79 double 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 
151 int 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<> >());
405  for (int deg=0; deg<=std::min(INTREPID_CUBATURE_LINE_GAUSS_MAX,INTREPID_CUBATURE_LINE_GAUSSJACOBI20_MAX); deg++) {
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 }
Intrepid::CubatureDirect::getAccuracy
virtual void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
Definition: Intrepid_CubatureDirectDef.hpp:104
Intrepid_CubatureDirectLineGauss.hpp
Header file for the Intrepid::CubatureDirectLineGauss class.
Intrepid::CubatureDirectLineGauss
Defines Gauss integration rules on a line.
Definition: Intrepid_CubatureDirectLineGauss.hpp:71
Intrepid::CubatureDirectTriDefault::getName
const char * getName() const
Returns cubature name.
Definition: Intrepid_CubatureDirectTriDefaultDef.hpp:77
Intrepid::CubatureDirectTetDefault
Defines direct integration rules on a tetrahedron.
Definition: Intrepid_CubatureDirectTetDefault.hpp:71
INTREPID_CUBATURE_TRI_DEFAULT_MAX
#define INTREPID_CUBATURE_TRI_DEFAULT_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct triangle rule of the ...
Definition: Intrepid_CubatureDirectTriDefault.hpp:62
Intrepid_CubatureTensor.hpp
Header file for the Intrepid::CubatureTensor class.
INTREPID_CUBATURE_LINE_GAUSS_MAX
#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...
Definition: Intrepid_CubatureDirectLineGauss.hpp:62
Intrepid_CubatureTensorPyr.hpp
Header file for the Intrepid::CubatureTensorPyr class.
INTREPID_CUBATURE_LINE_GAUSSJACOBI20_MAX
#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...
Definition: Intrepid_CubatureDirectLineGaussJacobi20.hpp:62
INTREPID_CUBATURE_TET_DEFAULT_MAX
#define INTREPID_CUBATURE_TET_DEFAULT_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct tetrahedron rule of t...
Definition: Intrepid_CubatureDirectTetDefault.hpp:62
Intrepid_CubatureDirectTetDefault.hpp
Header file for the Intrepid::CubatureDirectTetDefault class.
Intrepid::CubatureDirectLineGauss::getName
const char * getName() const
Returns cubature name.
Definition: Intrepid_CubatureDirectLineGaussDef.hpp:77
Intrepid_CubatureDirectTriDefault.hpp
Header file for the Intrepid::CubatureDirectTriDefault class.
Intrepid_CubatureDirectLineGaussJacobi20.hpp
Header file for the Intrepid::CubatureDirectLineGaussJacobi20 class.
Intrepid::CubatureTensorPyr
Defines tensor-product cubature (integration) rules in Intrepid.
Definition: Intrepid_CubatureTensorPyr.hpp:60
Intrepid::CubatureDirectTriDefault
Defines direct integration rules on a triangle.
Definition: Intrepid_CubatureDirectTriDefault.hpp:71
Intrepid::CubatureDirectTetDefault::getName
const char * getName() const
Returns cubature name.
Definition: Intrepid_CubatureDirectTetDefaultDef.hpp:77
Intrepid::Cubature
Defines the base class for cubature (integration) rules in Intrepid.
Definition: Intrepid_Cubature.hpp:70
Intrepid::CubatureTensor
Defines tensor-product cubature (integration) rules in Intrepid.
Definition: Intrepid_CubatureTensor.hpp:64
Intrepid::CubatureDirect
Defines direct cubature (integration) rules in Intrepid.
Definition: Intrepid_CubatureDirect.hpp:76
Intrepid::CubatureDirectLineGaussJacobi20
Defines GaussJacobi20 integration rules on a line.
Definition: Intrepid_CubatureDirectLineGaussJacobi20.hpp:71
Generated by   doxygen 1.8.5