Intrepid
test_01.cpp
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 
49 #include "Intrepid_PointTools.hpp"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
53 #include "Shards_CellTopology.hpp"
54 
55 using namespace std;
56 using namespace Intrepid;
57 
58 
59 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
60 { \
61  ++nException; \
62  try { \
63  S ; \
64  } \
65  catch (std::logic_error err) { \
66  ++throwCounter; \
67  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
68  *outStream << err.what() << '\n'; \
69  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
70  }; \
71 }
72 
73 
74 int main(int argc, char *argv[]) {
75 
76  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77 
78  // This little trick lets us print to std::cout only if
79  // a (dummy) command-line argument is provided.
80  int iprint = argc - 1;
81  Teuchos::RCP<std::ostream> outStream;
82  Teuchos::oblackholestream bhs; // outputs nothing
83  if (iprint > 0)
84  outStream = Teuchos::rcp(&std::cout, false);
85  else
86  outStream = Teuchos::rcp(&bhs, false);
87 
88  // Save the format state of the original std::cout.
89  Teuchos::oblackholestream oldFormatState;
90  oldFormatState.copyfmt(std::cout);
91 
92  *outStream \
93  << "===============================================================================\n" \
94  << "| |\n" \
95  << "| Unit Test (PointTools) |\n" \
96  << "| |\n" \
97  << "| 1) Construction of equispaced and warped lattices on simplices |\n" \
98  << "| |\n" \
99  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
100  << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
101  << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \
102  << "| |\n" \
103  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105  << "| |\n" \
106  << "===============================================================================\n";
107 
108 
109 
110  int errorFlag = 0;
111 
112  shards::CellTopology myLine_2( shards::getCellTopologyData< shards::Line<2> >() );
113  shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
114  shards::CellTopology myTet_4( shards::getCellTopologyData< shards::Tetrahedron<4> >() );
115 
116 
117  *outStream \
118  << "\n"
119  << "===============================================================================\n"\
120  << "| TEST 1: size of lattices |\n"\
121  << "===============================================================================\n";
122 
123  try{
124  // first try the lattices with offset = 0. This is a spot-check
125 
126  if (PointTools::getLatticeSize( myLine_2 , 4 , 0 ) != 5) {
127  errorFlag++;
128  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
129  *outStream << " size of 4th order lattice on a line with no offset: " << PointTools::getLatticeSize( myLine_2 , 4 , 0 ) << "\n";
130  *outStream << " should be 5\n";
131  }
132 
133 
134  if (PointTools::getLatticeSize( myTri_3 , 3 , 0 ) != 10) {
135  errorFlag++;
136  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
137  *outStream << " size of 3rd order lattice on a line with no offset: " << PointTools::getLatticeSize( myTri_3 , 3 , 0 ) << "\n";
138  *outStream << " should be 10\n";
139  }
140 
141 
142  if (PointTools::getLatticeSize( myTet_4 , 3 , 0 ) != 20) {
143  errorFlag++;
144  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
145  *outStream << " size of 3rd order lattice on a tet with no offset: " << PointTools::getLatticeSize( myTet_4 , 3 , 0 ) << "\n";
146  *outStream << " should be 20\n";
147  }
148 
149 
150  // check with the offset = 1
151  if (PointTools::getLatticeSize( myLine_2 , 3 , 1 ) != 2) {
152  errorFlag++;
153  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
154  *outStream << " size of 3rd order lattice on a line with 1 offset: " << PointTools::getLatticeSize( myLine_2 , 3 , 1 ) << "\n";
155  *outStream << " should be 2\n";
156  }
157 
158  if (PointTools::getLatticeSize( myTri_3 , 4 , 1 ) != 3) {
159  errorFlag++;
160  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
161  *outStream << " size of 4th order lattice on a triangle with 1 offset: " << PointTools::getLatticeSize( myTri_3 , 4 , 1 ) << "\n";
162  *outStream << " should be 3\n";
163  }
164 
165  if (PointTools::getLatticeSize( myTet_4 , 5 , 1 ) != 4) {
166  errorFlag++;
167  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
168  *outStream << " size of 5th order lattice on a tetrahedron with 1 offset: " << PointTools::getLatticeSize( myTet_4 , 5 , 1 ) << "\n";
169  *outStream << " should be 4\n";
170  }
171 
172  }
173  catch (std::exception &err) {
174  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
175  *outStream << err.what() << '\n';
176  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
177  errorFlag = -1000;
178  };
179 
180  // Now verify that we throw an exception on some of the non-supported cell types.
181 
182  *outStream \
183  << "\n"
184  << "===============================================================================\n"\
185  << "| TEST 2: check for unsupported cell types \n"\
186  << "===============================================================================\n";
187  try{
188  try {
189  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Quadrilateral<4> >() , 3 , 0 );
190  }
191  catch (std::invalid_argument err) {
192  *outStream << err.what() << "\n";
193  }
194 
195  try {
196  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Hexahedron<8> >() , 3 , 0 );
197  }
198  catch (std::invalid_argument err) {
199  *outStream << err.what() << "\n";
200  }
201  }
202  catch (std::exception &err) {
203  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
204  *outStream << err.what() << '\n';
205  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
206  errorFlag = -1000;
207  };
208 
209  // Check for malformed point arrays
210 
211 
212  *outStream \
213  << "\n"
214  << "===============================================================================\n"\
215  << "| TEST 2: malformed point arrays \n"\
216  << "===============================================================================\n";
217  try{
218  // line: not enough points allocated
219  try {
220  FieldContainer<double> pts(3,1);
221  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
222  }
223  catch (std::invalid_argument err) {
224  *outStream << err.what() << "\n";
225  }
226  // line: wrong dimension for points
227  try {
228  FieldContainer<double> pts(6,2);
229  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
230  }
231  catch (std::invalid_argument err) {
232  *outStream << err.what() << "\n";
233  }
234  // triangle: too many points allocated
235  try {
236  FieldContainer<double> pts(4,2);
237  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 1 );
238  }
239  catch (std::invalid_argument err) {
240  *outStream << err.what() << "\n";
241  }
242  // triangle: wrong dimension for points
243  try {
244  FieldContainer<double> pts(6,1);
245  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 0 );
246  }
247  catch (std::invalid_argument err) {
248  *outStream << err.what() << "\n";
249  }
250  // tetrahedron: not enough points allocated
251  try {
252  FieldContainer<double> pts(4,2);
253  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 2 , 0 );
254  }
255  catch (std::invalid_argument err) {
256  *outStream << err.what() << "\n";
257  }
258  // tetrahedron: wrong dimension for points
259  try {
260  FieldContainer<double> pts(4,2);
261  PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 1 , 0 );
262  }
263  catch (std::invalid_argument err) {
264  *outStream << err.what() << "\n";
265  }
266 
267  }
268  catch (std::exception &err) {
269  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
270  *outStream << err.what() << '\n';
271  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
272  errorFlag = -1000;
273  };
274 
275 
276  *outStream \
277  << "\n"
278  << "===============================================================================\n"\
279  << "| TEST 3: check values of triangular lattice compared to Warburton's code \n"\
280  << "===============================================================================\n";
281  try {
282  // triangle case
283  const int order = 4;
284  const int offset = 0;
285  int numPts = PointTools::getLatticeSize( myTri_3 , order , offset );
286  int ptDim = 2;
287  FieldContainer<double> warpBlendPts( numPts , ptDim );
288  PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
289  myTri_3 ,
290  order ,
291  offset ,
292  POINTTYPE_WARPBLEND );
293  FieldContainer<double> verts( 1, 3 , 2 );
294  verts(0,0,0) = -1.0;
295  verts(0,0,1) = -1.0/sqrt(3.0);
296  verts(0,1,0) = 1.0;
297  verts(0,1,1) = -1.0/sqrt(3.0);
298  verts(0,2,0) = 0.0;
299  verts(0,2,1) = 2.0/sqrt(3.0);
300 
301  // holds points on the equilateral triangle
302  FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
303 
304  CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
305  warpBlendPts ,
306  verts ,
307  myTri_3 ,
308  0 );
309 
310  // Values from MATLAB code
311  double points[] = { -1.000000000000000 , -0.577350269189626 ,
312  -0.654653670707977 , -0.577350269189626 ,
313  -0.000000000000000 , -0.577350269189626 ,
314  0.654653670707977 , -0.577350269189626 ,
315  1.000000000000000 , -0.577350269189626 ,
316  -0.827326835353989 , -0.278271574919028 ,
317  -0.327375261332958 , -0.189010195256608 ,
318  0.327375261332958 , -0.189010195256608 ,
319  0.827326835353989, -0.278271574919028,
320  -0.500000000000000, 0.288675134594813,
321  0.000000000000000, 0.378020390513215,
322  0.500000000000000, 0.288675134594813,
323  -0.172673164646011, 0.855621844108654,
324  0.172673164646011, 0.855621844108654,
325  0, 1.154700538379252 };
326 
327  // compare
328  for (int i=0;i<numPts;i++) {
329  for (int j=0;j<2;j++) {
330  int l = 2*i+j;
331  if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
332  errorFlag++;
333  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
334 
335  // Output the multi-index of the value where the error is:
336  *outStream << " At multi-index { ";
337  *outStream << i << " ";*outStream << j << " ";
338  *outStream << "} computed value: " << warpBlendMappedPts(i,j)
339  << " but correct value: " << points[l] << "\n";
340  }
341  }
342  }
343 
344 
345  }
346  catch (std::exception &err) {
347  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
348  *outStream << err.what() << '\n';
349  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
350  errorFlag = -1000;
351  };
352 
353 
354  *outStream \
355  << "\n"
356  << "===============================================================================\n"\
357  << "| TEST 4: check values of tetrahedral lattice compared to Warburton's code \n"\
358  << "===============================================================================\n";
359  try {
360  // triangle case
361  const int order = 6;
362  const int offset = 0;
363  int numPts = PointTools::getLatticeSize( myTet_4 , order , offset );
364  int ptDim = 3;
365  FieldContainer<double> warpBlendPts( numPts , ptDim );
366  PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
367  myTet_4 ,
368  order ,
369  offset ,
370  POINTTYPE_WARPBLEND );
371 
372  FieldContainer<double> verts(1,4,3);
373  verts(0,0,0) = -1.0;
374  verts(0,0,1) = -1.0/sqrt(3.0);
375  verts(0,0,2) = -1.0/sqrt(6.0);
376  verts(0,1,0) = 1.0;
377  verts(0,1,1) = -1.0/sqrt(3.0);
378  verts(0,1,2) = -1.0/sqrt(6.0);
379  verts(0,2,0) = 0.0;
380  verts(0,2,1) = 2.0 / sqrt(3.0);
381  verts(0,2,2) = -1.0/sqrt(6.0);
382  verts(0,3,0) = 0.0;
383  verts(0,3,1) = 0.0;
384  verts(0,3,2) = 3.0 / sqrt(6.0);
385 
386 
387  // points on the equilateral tet
388  FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
389 
390  CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
391  warpBlendPts ,
392  verts ,
393  myTet_4 ,
394  0 );
395 
396  // Values from MATLAB code
397  double points[] = { -1.000000000000000, -0.577350269189626, -0.408248290463863,
398  -0.830223896278567, -0.577350269189626, -0.408248290463863,
399  -0.468848793470714, -0.577350269189626, -0.408248290463863,
400  -0.000000000000000, -0.577350269189626, -0.408248290463863,
401  0.468848793470714, -0.577350269189626, -0.408248290463863,
402  0.830223896278567, -0.577350269189626, -0.408248290463863,
403  1.000000000000000, -0.577350269189626, -0.408248290463863,
404  -0.915111948139283, -0.430319850411323, -0.408248290463863,
405  -0.660434383303427, -0.381301968982318, -0.408248290463863,
406  -0.239932664820086, -0.368405260495326, -0.408248290463863,
407  0.239932664820086, -0.368405260495326, -0.408248290463863,
408  0.660434383303426, -0.381301968982318, -0.408248290463863,
409  0.915111948139283, -0.430319850411323, -0.408248290463863,
410  -0.734424396735357, -0.117359831084509, -0.408248290463863,
411  -0.439014646886819, -0.023585152684228, -0.408248290463863,
412  -0.000000000000000, -0.000000000000000, -0.408248290463863,
413  0.439014646886819, -0.023585152684228, -0.408248290463863,
414  0.734424396735357, -0.117359831084509, -0.408248290463863,
415  -0.500000000000000, 0.288675134594813, -0.408248290463863,
416  -0.199081982066733, 0.391990413179555, -0.408248290463863,
417  0.199081982066733, 0.391990413179555, -0.408248290463863,
418  0.500000000000000, 0.288675134594813, -0.408248290463863,
419  -0.265575603264643, 0.694710100274135, -0.408248290463863,
420  -0.000000000000000, 0.762603937964635, -0.408248290463863,
421  0.265575603264643, 0.694710100274135, -0.408248290463863,
422  -0.084888051860716, 1.007670119600949, -0.408248290463863,
423  0.084888051860716, 1.007670119600949, -0.408248290463863,
424  0, 1.154700538379252, -0.408248290463863,
425  -0.915111948139284, -0.528340129596858, -0.269626682252082,
426  -0.660434383303427, -0.512000835787190, -0.223412180441618,
427  -0.239932664820086, -0.507701932958193, -0.211253047073435,
428  0.239932664820086, -0.507701932958193, -0.211253047073435,
429  0.660434383303426, -0.512000835787190, -0.223412180441618,
430  0.915111948139284, -0.528340129596858, -0.269626682252082,
431  -0.773622922202284, -0.315952535579882, -0.223412180441618,
432  -0.421605613935553, -0.243414114697549, -0.172119771139157,
433  -0.000000000000000, -0.224211101329670, -0.158541190167514,
434  0.421605613935553, -0.243414114697549, -0.172119771139157,
435  0.773622922202284, -0.315952535579882, -0.223412180441618,
436  -0.559649103902302, 0.046063183547205, -0.211253047073435,
437  -0.194172509561981, 0.112105550664835, -0.158541190167514,
438  0.194172509561981, 0.112105550664835, -0.158541190167514,
439  0.559649103902302, 0.046063183547205, -0.211253047073435,
440  -0.319716439082216, 0.461638749410988, -0.211253047073435,
441  -0.000000000000000, 0.486828229395098, -0.172119771139157,
442  0.319716439082216, 0.461638749410988, -0.211253047073435,
443  -0.113188538898858, 0.827953371367071, -0.223412180441618,
444  0.113188538898858, 0.827953371367071, -0.223412180441618,
445  -0.000000000000000, 1.056680259193716, -0.269626682252082,
446  -0.734424396735357, -0.424020123154587, 0.025434853622935,
447  -0.439014646886819, -0.392761897021160, 0.113846468290170,
448  -0.000000000000000, -0.384900179459751, 0.136082763487954,
449  0.439014646886819, -0.392761897021160, 0.113846468290170,
450  0.734424396735357, -0.424020123154587, 0.025434853622935,
451  -0.559649103902302, -0.183816888326860, 0.113846468290170,
452  -0.194172509561981, -0.112105550664835, 0.158541190167514,
453  0.194172509561981, -0.112105550664835, 0.158541190167514,
454  0.559649103902302, -0.183816888326860, 0.113846468290170,
455  -0.333333333333333, 0.192450089729875, 0.136082763487954,
456  -0.000000000000000, 0.224211101329670, 0.158541190167514,
457  0.333333333333333, 0.192450089729875, 0.136082763487954,
458  -0.120634457015483, 0.576578785348020, 0.113846468290170,
459  0.120634457015482, 0.576578785348020, 0.113846468290170,
460  -0.000000000000000, 0.848040246309174, 0.025434853622935,
461  -0.500000000000000, -0.288675134594813, 0.408248290463863,
462  -0.199081982066733, -0.254236708399899, 0.505654869247127,
463  0.199081982066733, -0.254236708399899, 0.505654869247127,
464  0.500000000000000, -0.288675134594813, 0.408248290463863,
465  -0.319716439082216, -0.045291699705599, 0.505654869247127,
466  -0.000000000000000, -0.000000000000000, 0.516359313417471,
467  0.319716439082216, -0.045291699705599, 0.505654869247127,
468  -0.120634457015483, 0.299528408105498, 0.505654869247127,
469  0.120634457015483, 0.299528408105498, 0.505654869247127,
470  -0.000000000000000, 0.577350269189626, 0.408248290463863,
471  -0.265575603264643, -0.153330146035039, 0.791061727304791,
472  -0.000000000000000, -0.130698866804872, 0.855072651347100,
473  0.265575603264643, -0.153330146035039, 0.791061727304791,
474  -0.113188538898858, 0.065349433402436, 0.855072651347100,
475  0.113188538898858, 0.065349433402436, 0.855072651347099,
476  0, 0.306660292070078, 0.791061727304791,
477  -0.084888051860717, -0.049010139592768, 1.086123263179808,
478  0.084888051860717, -0.049010139592768, 1.086123263179808,
479  0.000000000000000, 0.098020279185535, 1.086123263179808,
480  0, 0, 1.224744871391589
481  };
482 
483  // compare
484  for (int i=0;i<numPts;i++) {
485  for (int j=0;j<ptDim;j++) {
486  int l = ptDim*i+j;
487  if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
488  errorFlag++;
489  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
490 
491  // Output the multi-index of the value where the error is:
492  *outStream << " At multi-index { ";
493  *outStream << i << " ";*outStream << j << " ";
494  *outStream << "} computed value: " << warpBlendMappedPts(i,j)
495  << " but correct value: " << points[l] << "\n";
496  }
497  }
498  }
499 
500 
501  }
502  catch (std::exception &err) {
503  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
504  *outStream << err.what() << '\n';
505  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
506  errorFlag = -1000;
507  };
508 
509 
510 
511  if (errorFlag != 0)
512  std::cout << "End Result: TEST FAILED\n";
513  else
514  std::cout << "End Result: TEST PASSED\n";
515 
516  // reset format state of std::cout
517  std::cout.copyfmt(oldFormatState);
518 
519  return errorFlag;
520 }
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:111
Generated by   doxygen 1.8.5