|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00035 #include "Intrepid_PointTools.hpp" 00036 #include "Intrepid_FieldContainer.hpp" 00037 #include "Teuchos_oblackholestream.hpp" 00038 #include "Teuchos_GlobalMPISession.hpp" 00039 #include "Shards_CellTopology.hpp" 00040 00041 using namespace std; 00042 using namespace Intrepid; 00043 00044 00045 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00046 { \ 00047 ++nException; \ 00048 try { \ 00049 S ; \ 00050 } \ 00051 catch (std::logic_error err) { \ 00052 ++throwCounter; \ 00053 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00054 *outStream << err.what() << '\n'; \ 00055 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00056 }; \ 00057 } 00058 00059 00060 int main(int argc, char *argv[]) { 00061 00062 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00063 00064 // This little trick lets us print to std::cout only if 00065 // a (dummy) command-line argument is provided. 00066 int iprint = argc - 1; 00067 Teuchos::RCP<std::ostream> outStream; 00068 Teuchos::oblackholestream bhs; // outputs nothing 00069 if (iprint > 0) 00070 outStream = Teuchos::rcp(&std::cout, false); 00071 else 00072 outStream = Teuchos::rcp(&bhs, false); 00073 00074 // Save the format state of the original std::cout. 00075 Teuchos::oblackholestream oldFormatState; 00076 oldFormatState.copyfmt(std::cout); 00077 00078 *outStream \ 00079 << "===============================================================================\n" \ 00080 << "| |\n" \ 00081 << "| Unit Test (PointTools) |\n" \ 00082 << "| |\n" \ 00083 << "| 1) Construction of equispaced and warped lattices on simplices |\n" \ 00084 << "| |\n" \ 00085 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00086 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \ 00087 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \ 00088 << "| |\n" \ 00089 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00090 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00091 << "| |\n" \ 00092 << "===============================================================================\n"; 00093 00094 00095 00096 int errorFlag = 0; 00097 00098 shards::CellTopology myLine_2( shards::getCellTopologyData< shards::Line<2> >() ); 00099 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() ); 00100 shards::CellTopology myTet_4( shards::getCellTopologyData< shards::Tetrahedron<4> >() ); 00101 00102 00103 *outStream \ 00104 << "\n" 00105 << "===============================================================================\n"\ 00106 << "| TEST 1: size of lattices |\n"\ 00107 << "===============================================================================\n"; 00108 00109 try{ 00110 // first try the lattices with offset = 0. This is a spot-check 00111 00112 if (PointTools::getLatticeSize( myLine_2 , 4 , 0 ) != 5) { 00113 errorFlag++; 00114 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00115 *outStream << " size of 4th order lattice on a line with no offset: " << PointTools::getLatticeSize( myLine_2 , 4 , 0 ) << "\n"; 00116 *outStream << " should be 5\n"; 00117 } 00118 00119 00120 if (PointTools::getLatticeSize( myTri_3 , 3 , 0 ) != 10) { 00121 errorFlag++; 00122 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00123 *outStream << " size of 3rd order lattice on a line with no offset: " << PointTools::getLatticeSize( myTri_3 , 3 , 0 ) << "\n"; 00124 *outStream << " should be 10\n"; 00125 } 00126 00127 00128 if (PointTools::getLatticeSize( myTet_4 , 3 , 0 ) != 20) { 00129 errorFlag++; 00130 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00131 *outStream << " size of 3rd order lattice on a tet with no offset: " << PointTools::getLatticeSize( myTet_4 , 3 , 0 ) << "\n"; 00132 *outStream << " should be 20\n"; 00133 } 00134 00135 00136 // check with the offset = 1 00137 if (PointTools::getLatticeSize( myLine_2 , 3 , 1 ) != 2) { 00138 errorFlag++; 00139 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00140 *outStream << " size of 3rd order lattice on a line with 1 offset: " << PointTools::getLatticeSize( myLine_2 , 3 , 1 ) << "\n"; 00141 *outStream << " should be 2\n"; 00142 } 00143 00144 if (PointTools::getLatticeSize( myTri_3 , 4 , 1 ) != 3) { 00145 errorFlag++; 00146 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00147 *outStream << " size of 4th order lattice on a triangle with 1 offset: " << PointTools::getLatticeSize( myTri_3 , 4 , 1 ) << "\n"; 00148 *outStream << " should be 3\n"; 00149 } 00150 00151 if (PointTools::getLatticeSize( myTet_4 , 5 , 1 ) != 4) { 00152 errorFlag++; 00153 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00154 *outStream << " size of 5th order lattice on a tetrahedron with 1 offset: " << PointTools::getLatticeSize( myTet_4 , 5 , 1 ) << "\n"; 00155 *outStream << " should be 4\n"; 00156 } 00157 00158 } 00159 catch (std::exception &err) { 00160 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00161 *outStream << err.what() << '\n'; 00162 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00163 errorFlag = -1000; 00164 }; 00165 00166 // Now verify that we throw an exception on some of the non-supported cell types. 00167 00168 *outStream \ 00169 << "\n" 00170 << "===============================================================================\n"\ 00171 << "| TEST 2: check for unsupported cell types \n"\ 00172 << "===============================================================================\n"; 00173 try{ 00174 try { 00175 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Quadrilateral<4> >() , 3 , 0 ); 00176 } 00177 catch (std::invalid_argument err) { 00178 *outStream << err.what() << "\n"; 00179 } 00180 00181 try { 00182 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Hexahedron<8> >() , 3 , 0 ); 00183 } 00184 catch (std::invalid_argument err) { 00185 *outStream << err.what() << "\n"; 00186 } 00187 } 00188 catch (std::exception &err) { 00189 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00190 *outStream << err.what() << '\n'; 00191 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00192 errorFlag = -1000; 00193 }; 00194 00195 // Check for malformed point arrays 00196 00197 00198 *outStream \ 00199 << "\n" 00200 << "===============================================================================\n"\ 00201 << "| TEST 2: malformed point arrays \n"\ 00202 << "===============================================================================\n"; 00203 try{ 00204 // line: not enough points allocated 00205 try { 00206 FieldContainer<double> pts(3,1); 00207 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 ); 00208 } 00209 catch (std::invalid_argument err) { 00210 *outStream << err.what() << "\n"; 00211 } 00212 // line: wrong dimension for points 00213 try { 00214 FieldContainer<double> pts(6,2); 00215 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 ); 00216 } 00217 catch (std::invalid_argument err) { 00218 *outStream << err.what() << "\n"; 00219 } 00220 // triangle: too many points allocated 00221 try { 00222 FieldContainer<double> pts(4,2); 00223 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 1 ); 00224 } 00225 catch (std::invalid_argument err) { 00226 *outStream << err.what() << "\n"; 00227 } 00228 // triangle: wrong dimension for points 00229 try { 00230 FieldContainer<double> pts(6,1); 00231 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 0 ); 00232 } 00233 catch (std::invalid_argument err) { 00234 *outStream << err.what() << "\n"; 00235 } 00236 // tetrahedron: not enough points allocated 00237 try { 00238 FieldContainer<double> pts(4,2); 00239 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 2 , 0 ); 00240 } 00241 catch (std::invalid_argument err) { 00242 *outStream << err.what() << "\n"; 00243 } 00244 // tetrahedron: wrong dimension for points 00245 try { 00246 FieldContainer<double> pts(4,2); 00247 PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 1 , 0 ); 00248 } 00249 catch (std::invalid_argument err) { 00250 *outStream << err.what() << "\n"; 00251 } 00252 00253 } 00254 catch (std::exception &err) { 00255 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00256 *outStream << err.what() << '\n'; 00257 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00258 errorFlag = -1000; 00259 }; 00260 00261 00262 *outStream \ 00263 << "\n" 00264 << "===============================================================================\n"\ 00265 << "| TEST 3: check values of triangular lattice compared to Warburton's code \n"\ 00266 << "===============================================================================\n"; 00267 try { 00268 // triangle case 00269 const int order = 4; 00270 const int offset = 0; 00271 int numPts = PointTools::getLatticeSize( myTri_3 , order , offset ); 00272 int ptDim = 2; 00273 FieldContainer<double> warpBlendPts( numPts , ptDim ); 00274 PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts , 00275 myTri_3 , 00276 order , 00277 offset , 00278 POINTTYPE_WARPBLEND ); 00279 FieldContainer<double> verts( 1, 3 , 2 ); 00280 verts(0,0,0) = -1.0; 00281 verts(0,0,1) = -1.0/sqrt(3.0); 00282 verts(0,1,0) = 1.0; 00283 verts(0,1,1) = -1.0/sqrt(3.0); 00284 verts(0,2,0) = 0.0; 00285 verts(0,2,1) = 2.0/sqrt(3.0); 00286 00287 // holds points on the equilateral triangle 00288 FieldContainer<double> warpBlendMappedPts( numPts , ptDim ); 00289 00290 CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts , 00291 warpBlendPts , 00292 verts , 00293 myTri_3 , 00294 0 ); 00295 00296 // Values from MATLAB code 00297 double points[] = { -1.000000000000000 , -0.577350269189626 , 00298 -0.654653670707977 , -0.577350269189626 , 00299 -0.000000000000000 , -0.577350269189626 , 00300 0.654653670707977 , -0.577350269189626 , 00301 1.000000000000000 , -0.577350269189626 , 00302 -0.827326835353989 , -0.278271574919028 , 00303 -0.327375261332958 , -0.189010195256608 , 00304 0.327375261332958 , -0.189010195256608 , 00305 0.827326835353989, -0.278271574919028, 00306 -0.500000000000000, 0.288675134594813, 00307 0.000000000000000, 0.378020390513215, 00308 0.500000000000000, 0.288675134594813, 00309 -0.172673164646011, 0.855621844108654, 00310 0.172673164646011, 0.855621844108654, 00311 0, 1.154700538379252 }; 00312 00313 // compare 00314 for (int i=0;i<numPts;i++) { 00315 for (int j=0;j<2;j++) { 00316 int l = 2*i+j; 00317 if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) { 00318 errorFlag++; 00319 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00320 00321 // Output the multi-index of the value where the error is: 00322 *outStream << " At multi-index { "; 00323 *outStream << i << " ";*outStream << j << " "; 00324 *outStream << "} computed value: " << warpBlendMappedPts(i,j) 00325 << " but correct value: " << points[l] << "\n"; 00326 } 00327 } 00328 } 00329 00330 00331 } 00332 catch (std::exception &err) { 00333 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00334 *outStream << err.what() << '\n'; 00335 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00336 errorFlag = -1000; 00337 }; 00338 00339 00340 *outStream \ 00341 << "\n" 00342 << "===============================================================================\n"\ 00343 << "| TEST 4: check values of tetrahedral lattice compared to Warburton's code \n"\ 00344 << "===============================================================================\n"; 00345 try { 00346 // triangle case 00347 const int order = 6; 00348 const int offset = 0; 00349 int numPts = PointTools::getLatticeSize( myTet_4 , order , offset ); 00350 int ptDim = 3; 00351 FieldContainer<double> warpBlendPts( numPts , ptDim ); 00352 PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts , 00353 myTet_4 , 00354 order , 00355 offset , 00356 POINTTYPE_WARPBLEND ); 00357 00358 FieldContainer<double> verts(1,4,3); 00359 verts(0,0,0) = -1.0; 00360 verts(0,0,1) = -1.0/sqrt(3.0); 00361 verts(0,0,2) = -1.0/sqrt(6.0); 00362 verts(0,1,0) = 1.0; 00363 verts(0,1,1) = -1.0/sqrt(3.0); 00364 verts(0,1,2) = -1.0/sqrt(6.0); 00365 verts(0,2,0) = 0.0; 00366 verts(0,2,1) = 2.0 / sqrt(3.0); 00367 verts(0,2,2) = -1.0/sqrt(6.0); 00368 verts(0,3,0) = 0.0; 00369 verts(0,3,1) = 0.0; 00370 verts(0,3,2) = 3.0 / sqrt(6.0); 00371 00372 00373 // points on the equilateral tet 00374 FieldContainer<double> warpBlendMappedPts( numPts , ptDim ); 00375 00376 CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts , 00377 warpBlendPts , 00378 verts , 00379 myTet_4 , 00380 0 ); 00381 00382 // Values from MATLAB code 00383 double points[] = { -1.000000000000000, -0.577350269189626, -0.408248290463863, 00384 -0.830223896278567, -0.577350269189626, -0.408248290463863, 00385 -0.468848793470714, -0.577350269189626, -0.408248290463863, 00386 -0.000000000000000, -0.577350269189626, -0.408248290463863, 00387 0.468848793470714, -0.577350269189626, -0.408248290463863, 00388 0.830223896278567, -0.577350269189626, -0.408248290463863, 00389 1.000000000000000, -0.577350269189626, -0.408248290463863, 00390 -0.915111948139283, -0.430319850411323, -0.408248290463863, 00391 -0.660434383303427, -0.381301968982318, -0.408248290463863, 00392 -0.239932664820086, -0.368405260495326, -0.408248290463863, 00393 0.239932664820086, -0.368405260495326, -0.408248290463863, 00394 0.660434383303426, -0.381301968982318, -0.408248290463863, 00395 0.915111948139283, -0.430319850411323, -0.408248290463863, 00396 -0.734424396735357, -0.117359831084509, -0.408248290463863, 00397 -0.439014646886819, -0.023585152684228, -0.408248290463863, 00398 -0.000000000000000, -0.000000000000000, -0.408248290463863, 00399 0.439014646886819, -0.023585152684228, -0.408248290463863, 00400 0.734424396735357, -0.117359831084509, -0.408248290463863, 00401 -0.500000000000000, 0.288675134594813, -0.408248290463863, 00402 -0.199081982066733, 0.391990413179555, -0.408248290463863, 00403 0.199081982066733, 0.391990413179555, -0.408248290463863, 00404 0.500000000000000, 0.288675134594813, -0.408248290463863, 00405 -0.265575603264643, 0.694710100274135, -0.408248290463863, 00406 -0.000000000000000, 0.762603937964635, -0.408248290463863, 00407 0.265575603264643, 0.694710100274135, -0.408248290463863, 00408 -0.084888051860716, 1.007670119600949, -0.408248290463863, 00409 0.084888051860716, 1.007670119600949, -0.408248290463863, 00410 0, 1.154700538379252, -0.408248290463863, 00411 -0.915111948139284, -0.528340129596858, -0.269626682252082, 00412 -0.660434383303427, -0.512000835787190, -0.223412180441618, 00413 -0.239932664820086, -0.507701932958193, -0.211253047073435, 00414 0.239932664820086, -0.507701932958193, -0.211253047073435, 00415 0.660434383303426, -0.512000835787190, -0.223412180441618, 00416 0.915111948139284, -0.528340129596858, -0.269626682252082, 00417 -0.773622922202284, -0.315952535579882, -0.223412180441618, 00418 -0.421605613935553, -0.243414114697549, -0.172119771139157, 00419 -0.000000000000000, -0.224211101329670, -0.158541190167514, 00420 0.421605613935553, -0.243414114697549, -0.172119771139157, 00421 0.773622922202284, -0.315952535579882, -0.223412180441618, 00422 -0.559649103902302, 0.046063183547205, -0.211253047073435, 00423 -0.194172509561981, 0.112105550664835, -0.158541190167514, 00424 0.194172509561981, 0.112105550664835, -0.158541190167514, 00425 0.559649103902302, 0.046063183547205, -0.211253047073435, 00426 -0.319716439082216, 0.461638749410988, -0.211253047073435, 00427 -0.000000000000000, 0.486828229395098, -0.172119771139157, 00428 0.319716439082216, 0.461638749410988, -0.211253047073435, 00429 -0.113188538898858, 0.827953371367071, -0.223412180441618, 00430 0.113188538898858, 0.827953371367071, -0.223412180441618, 00431 -0.000000000000000, 1.056680259193716, -0.269626682252082, 00432 -0.734424396735357, -0.424020123154587, 0.025434853622935, 00433 -0.439014646886819, -0.392761897021160, 0.113846468290170, 00434 -0.000000000000000, -0.384900179459751, 0.136082763487954, 00435 0.439014646886819, -0.392761897021160, 0.113846468290170, 00436 0.734424396735357, -0.424020123154587, 0.025434853622935, 00437 -0.559649103902302, -0.183816888326860, 0.113846468290170, 00438 -0.194172509561981, -0.112105550664835, 0.158541190167514, 00439 0.194172509561981, -0.112105550664835, 0.158541190167514, 00440 0.559649103902302, -0.183816888326860, 0.113846468290170, 00441 -0.333333333333333, 0.192450089729875, 0.136082763487954, 00442 -0.000000000000000, 0.224211101329670, 0.158541190167514, 00443 0.333333333333333, 0.192450089729875, 0.136082763487954, 00444 -0.120634457015483, 0.576578785348020, 0.113846468290170, 00445 0.120634457015482, 0.576578785348020, 0.113846468290170, 00446 -0.000000000000000, 0.848040246309174, 0.025434853622935, 00447 -0.500000000000000, -0.288675134594813, 0.408248290463863, 00448 -0.199081982066733, -0.254236708399899, 0.505654869247127, 00449 0.199081982066733, -0.254236708399899, 0.505654869247127, 00450 0.500000000000000, -0.288675134594813, 0.408248290463863, 00451 -0.319716439082216, -0.045291699705599, 0.505654869247127, 00452 -0.000000000000000, -0.000000000000000, 0.516359313417471, 00453 0.319716439082216, -0.045291699705599, 0.505654869247127, 00454 -0.120634457015483, 0.299528408105498, 0.505654869247127, 00455 0.120634457015483, 0.299528408105498, 0.505654869247127, 00456 -0.000000000000000, 0.577350269189626, 0.408248290463863, 00457 -0.265575603264643, -0.153330146035039, 0.791061727304791, 00458 -0.000000000000000, -0.130698866804872, 0.855072651347100, 00459 0.265575603264643, -0.153330146035039, 0.791061727304791, 00460 -0.113188538898858, 0.065349433402436, 0.855072651347100, 00461 0.113188538898858, 0.065349433402436, 0.855072651347099, 00462 0, 0.306660292070078, 0.791061727304791, 00463 -0.084888051860717, -0.049010139592768, 1.086123263179808, 00464 0.084888051860717, -0.049010139592768, 1.086123263179808, 00465 0.000000000000000, 0.098020279185535, 1.086123263179808, 00466 0, 0, 1.224744871391589 00467 }; 00468 00469 // compare 00470 for (int i=0;i<numPts;i++) { 00471 for (int j=0;j<ptDim;j++) { 00472 int l = ptDim*i+j; 00473 if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) { 00474 errorFlag++; 00475 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00476 00477 // Output the multi-index of the value where the error is: 00478 *outStream << " At multi-index { "; 00479 *outStream << i << " ";*outStream << j << " "; 00480 *outStream << "} computed value: " << warpBlendMappedPts(i,j) 00481 << " but correct value: " << points[l] << "\n"; 00482 } 00483 } 00484 } 00485 00486 00487 } 00488 catch (std::exception &err) { 00489 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00490 *outStream << err.what() << '\n'; 00491 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00492 errorFlag = -1000; 00493 }; 00494 00495 00496 00497 if (errorFlag != 0) 00498 std::cout << "End Result: TEST FAILED\n"; 00499 else 00500 std::cout << "End Result: TEST PASSED\n"; 00501 00502 // reset format state of std::cout 00503 std::cout.copyfmt(oldFormatState); 00504 00505 return errorFlag; 00506 }
1.7.4