|
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 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00044 00049 #include "Intrepid_CellTools.hpp" 00050 #include "Intrepid_FieldContainer.hpp" 00051 #include "Intrepid_DefaultCubatureFactory.hpp" 00052 #include "Shards_CellTopology.hpp" 00053 00054 #include "Teuchos_oblackholestream.hpp" 00055 #include "Teuchos_RCP.hpp" 00056 #include "Teuchos_GlobalMPISession.hpp" 00057 #include "Teuchos_ScalarTraits.hpp" 00058 00059 using namespace std; 00060 using namespace Intrepid; 00061 using namespace shards; 00062 00063 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00064 { \ 00065 ++nException; \ 00066 try { \ 00067 S ; \ 00068 } \ 00069 catch (std::logic_error err) { \ 00070 ++throwCounter; \ 00071 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00072 *outStream << err.what() << '\n'; \ 00073 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00074 }; \ 00075 } 00076 00077 00078 00079 int main(int argc, char *argv[]) { 00080 00081 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00082 00083 typedef CellTools<double> CellTools; 00084 typedef shards::CellTopology CellTopology; 00085 00086 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided. 00087 int iprint = argc - 1; 00088 00089 Teuchos::RCP<std::ostream> outStream; 00090 Teuchos::oblackholestream bhs; // outputs nothing 00091 00092 if (iprint > 0) 00093 outStream = Teuchos::rcp(&std::cout, false); 00094 else 00095 outStream = Teuchos::rcp(&bhs, false); 00096 00097 // Save the format state of the original std::cout. 00098 Teuchos::oblackholestream oldFormatState; 00099 oldFormatState.copyfmt(std::cout); 00100 00101 *outStream \ 00102 << "===============================================================================\n" \ 00103 << "| |\n" \ 00104 << "| Unit Test CellTools |\n" \ 00105 << "| |\n" \ 00106 << "| 1) Mapping to and from reference cells with base and extended topologies|\n" \ 00107 << "| using default initial guesses when computing the inverse F^{-1} |\n" \ 00108 << "| 2) Repeat all tests from 1) using user-defined initial guess for F^{-1} |\n" \ 00109 << "| 3) Exception testing |\n" \ 00110 << "| |\n" \ 00111 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \ 00112 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \ 00113 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \ 00114 << "| |\n" \ 00115 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00116 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00117 << "| |\n" \ 00118 << "===============================================================================\n"; 00119 00120 int errorFlag = 0; 00121 00122 // Collect all supported cell topologies 00123 std::vector<shards::CellTopology> supportedTopologies; 00124 supportedTopologies.push_back(shards::getCellTopologyData<Triangle<3> >() ); 00125 supportedTopologies.push_back(shards::getCellTopologyData<Triangle<6> >() ); 00126 supportedTopologies.push_back(shards::getCellTopologyData<Quadrilateral<4> >() ); 00127 supportedTopologies.push_back(shards::getCellTopologyData<Quadrilateral<9> >() ); 00128 supportedTopologies.push_back(shards::getCellTopologyData<Tetrahedron<4> >() ); 00129 supportedTopologies.push_back(shards::getCellTopologyData<Tetrahedron<10> >() ); 00130 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<8> >() ); 00131 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<20> >() ); 00132 supportedTopologies.push_back(shards::getCellTopologyData<Hexahedron<27> >() ); 00133 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<6> >() ); 00134 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<15> >() ); 00135 supportedTopologies.push_back(shards::getCellTopologyData<Wedge<18> >() ); 00136 supportedTopologies.push_back(shards::getCellTopologyData<Pyramid<5> >() ); 00137 supportedTopologies.push_back(shards::getCellTopologyData<Pyramid<13> >() ); 00138 00139 // Declare iterator to loop over the cell topologies 00140 std::vector<shards::CellTopology>::iterator topo_iterator; 00141 00142 // Test 1 scope 00143 try{ 00144 00145 *outStream \ 00146 << "\n" 00147 << "===============================================================================\n"\ 00148 << "| Test 1: computing F(x) and F^{-1}(x) using default initial guesses. |\n"\ 00149 << "===============================================================================\n\n"; 00150 /* 00151 * Test summary: 00152 * 00153 * A reference point set is mapped to physical frame and then back to reference frame. 00154 * Test passes if the final set of points matches the first set of points. The cell workset 00155 * is generated by perturbing randomly the cellWorkset of a reference cell with the specified 00156 * cell topology. 00157 * 00158 */ 00159 // Declare arrays for cell workset and point sets. We will have 10 cells in the wset. and 10 pts per pt. set 00160 FieldContainer<double> cellWorkset; // physical cell workset 00161 FieldContainer<double> refPoints; // reference point set(s) 00162 FieldContainer<double> physPoints; // physical point set(s) 00163 FieldContainer<double> controlPoints; // preimages: physical points mapped back to ref. frame 00164 00165 // We will use cubature factory to get some points on the reference cells. Declare necessary arrays 00166 DefaultCubatureFactory<double> cubFactory; 00167 FieldContainer<double> cubPoints; 00168 FieldContainer<double> cubWeights; 00169 00170 // Initialize number of cells in the cell workset 00171 int numCells = 10; 00172 00173 00174 // Loop over cell topologies, make cell workset for each one by perturbing the cellWorkset & test methods 00175 for(topo_iterator = supportedTopologies.begin(); topo_iterator != supportedTopologies.end(); ++topo_iterator){ 00176 00177 // 1. Define a single reference point set using cubature factory with order 4 cubature 00178 Teuchos::RCP<Cubature<double> > cellCubature = cubFactory.create( (*topo_iterator), 4); 00179 int cubDim = cellCubature -> getDimension(); 00180 int numPts = cellCubature -> getNumPoints(); 00181 cubPoints.resize(numPts, cubDim); 00182 cubWeights.resize(numPts); 00183 cellCubature -> getCubature(cubPoints, cubWeights); 00184 00185 // 2. Define a cell workset by perturbing the cellWorkset of the reference cell with the specified topology 00186 // 2.1 Resize dimensions of the rank-3 (C,N,D) cell workset array for the current topology 00187 int numNodes = (*topo_iterator).getNodeCount(); 00188 int cellDim = (*topo_iterator).getDimension(); 00189 cellWorkset.resize(numCells, numNodes, cellDim); 00190 00191 // 2.2 Copy cellWorkset of the reference cell with the same topology to temp rank-2 (N,D) array 00192 FieldContainer<double> refCellNodes(numNodes, cellDim ); 00193 CellTools::getReferenceSubcellNodes(refCellNodes, cellDim, 0, (*topo_iterator) ); 00194 00195 // 2.3 Create randomly perturbed version of the reference cell and save in the cell workset array 00196 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00197 00198 // Move vertices +/-0.125 along their axes. Gives nondegenerate cells for base and extended topologies 00199 for(int nodeOrd = 0; nodeOrd < numNodes; nodeOrd++){ 00200 for(int d = 0; d < cellDim; d++){ 00201 double delta = Teuchos::ScalarTraits<double>::random()/16.0; 00202 cellWorkset(cellOrd, nodeOrd, d) = refCellNodes(nodeOrd, d) + delta; 00203 } // d 00204 }// nodeOrd 00205 }// cellOrd 00206 /* 00207 * 3.1 Test 1: single point set to single physical cell: map ref. point set in rank-2 (P,D) array 00208 * to a physical point set in rank-2 (P,D) array for a specified cell ordinal. Use the cub. 00209 * points array for this test. Resize physPoints and controlPoints to rank-2 (P,D) arrays. 00210 */ 00211 physPoints.resize(numPts, cubDim); 00212 controlPoints.resize(numPts, cubDim); 00213 00214 *outStream 00215 << " Mapping a set of " << numPts << " points to one cell in a workset of " << numCells << " " 00216 << (*topo_iterator).getName() << " cells. \n"; 00217 00218 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00219 00220 // Forward map:: requires cell ordinal 00221 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator), cellOrd); 00222 // Inverse map: requires cell ordinal 00223 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator), cellOrd); 00224 00225 // Points in controlPoints should match the originals in cubPoints up to a tolerance 00226 for(int pt = 0; pt < numPts; pt++){ 00227 for(int d = 0; d < cellDim; d++){ 00228 00229 if( abs( controlPoints(pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00230 errorFlag++; 00231 *outStream 00232 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00233 << " Mapping a single point set to a single physical cell in a workset failed for: \n" 00234 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00235 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00236 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00237 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00238 << " Original value = " << cubPoints(pt, d) << "\n" 00239 << " F^{-1}F(P_d) = " << controlPoints(pt, d) <<"\n"; 00240 } 00241 }// d 00242 }// pt 00243 }// cellOrd 00244 /* 00245 * 3.2 Test 2: single point set to multiple physical cells: map ref. point set in rank-2 (P,D) array 00246 * to a physical point set in rank-3 (C, P,D) array for all cell ordinals. Use the cub. 00247 * points array for this test. Resize physPoints and controlPoints to rank-3 (C,P,D) arrays. 00248 */ 00249 physPoints.clear(); 00250 controlPoints.clear(); 00251 physPoints.resize(numCells, numPts, cubDim); 00252 controlPoints.resize(numCells, numPts, cubDim); 00253 00254 *outStream 00255 << " Mapping a set of " << numPts << " points to all cells in workset of " << numCells << " " 00256 << (*topo_iterator).getName() << " cells. \n"; 00257 00258 // Forward map: do not specify cell ordinal 00259 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator)); 00260 // Inverse map: do not specify cell ordinal 00261 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator)); 00262 00263 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance 00264 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00265 for(int pt = 0; pt < numPts; pt++){ 00266 for(int d = 0; d < cellDim; d++){ 00267 00268 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00269 errorFlag++; 00270 *outStream 00271 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00272 << " Mapping a single point set to all physical cells in a workset failed for: \n" 00273 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00274 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00275 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00276 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00277 << " Original value = " << cubPoints(pt, d) << "\n" 00278 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n"; 00279 } 00280 }// d 00281 }// pt 00282 }// cellOrd 00283 /* 00284 * 3.3 Test 3: multiple point sets to multiple physical cells: map ref. point sets in rank-3 (C,P,D) array 00285 * to physical point sets in rank-3 (C, P,D) array for all cell ordinals. The (C,P,D) array 00286 * with reference point sets is obtained by cloning the cubature point array. 00287 */ 00288 physPoints.clear(); 00289 controlPoints.clear(); 00290 refPoints.resize(numCells, numPts, cubDim); 00291 physPoints.resize(numCells, numPts, cubDim); 00292 controlPoints.resize(numCells, numPts, cubDim); 00293 00294 // Clone cubature points in refPoints: 00295 for(int c = 0; c < numCells; c++){ 00296 for(int pt = 0; pt < numPts; pt++){ 00297 for(int d = 0; d < cellDim; d++){ 00298 refPoints(c, pt, d) = cubPoints(pt, d); 00299 }// d 00300 }// pt 00301 }// c 00302 00303 *outStream 00304 << " Mapping " << numCells << " sets of " << numPts << " points to corresponding cells in workset of " << numCells << " " 00305 << (*topo_iterator).getName() << " cells. \n"; 00306 00307 // Forward map: do not specify cell ordinal 00308 CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator)); 00309 // Inverse map: do not specify cell ordinal 00310 CellTools::mapToReferenceFrame(controlPoints, physPoints, cellWorkset, (*topo_iterator)); 00311 00312 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance 00313 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00314 for(int pt = 0; pt < numPts; pt++){ 00315 for(int d = 0; d < cellDim; d++){ 00316 00317 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00318 errorFlag++; 00319 *outStream 00320 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00321 << " Mapping multiple point sets to corresponding physical cells in a workset failed for: \n" 00322 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00323 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00324 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00325 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00326 << " Original value = " << refPoints(cellOrd, pt, d) << "\n" 00327 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n"; 00328 } 00329 }// d 00330 }// pt 00331 }// cellOrd 00332 }// topo_iterator 00333 }// try using default initial guesses branch 00334 00335 /************************************************************************************************* 00336 * Wrap up test: check if the test broke down unexpectedly due to an exception * 00337 ************************************************************************************************/ 00338 00339 catch (std::logic_error err) { 00340 *outStream << err.what() << "\n"; 00341 errorFlag = -1000; 00342 }; 00343 00344 00345 // Test 2: repeat all of the above using the original points as user-defined initial guess 00346 try{ 00347 00348 *outStream \ 00349 << "\n" 00350 << "===============================================================================\n"\ 00351 << "| Test 2: computing F(x) and F^{-1}(x) using user-defined initial guess. |\n"\ 00352 << "===============================================================================\n\n"; 00353 /* 00354 * Test summary: 00355 * 00356 * Repeats all parts of Test 1 using user-defined initial guesses in the computation of F^{-1}. 00357 * The guesses are simply the exact solutions (the original set of points we started with) 00358 * and so, this tests runs much faster because Newton converges in a single iteration. 00359 * 00360 */ 00361 00362 // Declare arrays for cell workset and point sets. We will have 10 cells in the wset. and 10 pts per pt. set 00363 FieldContainer<double> cellWorkset; // physical cell workset 00364 FieldContainer<double> physPoints; // physical point set(s) 00365 FieldContainer<double> controlPoints; // preimages: physical points mapped back to ref. frame 00366 FieldContainer<double> initialGuess; // User-defined initial guesses for F^{-1} 00367 00368 // We will use cubature factory to get some points on the reference cells. Declare necessary arrays 00369 DefaultCubatureFactory<double> cubFactory; 00370 FieldContainer<double> cubPoints; 00371 FieldContainer<double> cubWeights; 00372 00373 // Initialize number of cells in the cell workset 00374 int numCells = 10; 00375 00376 00377 // Loop over cell topologies, make cell workset for each one by perturbing the cellWorkset & test methods 00378 for(topo_iterator = supportedTopologies.begin(); topo_iterator != supportedTopologies.end(); ++topo_iterator){ 00379 00380 // 1. Define a single reference point set using cubature factory with order 6 cubature 00381 Teuchos::RCP<Cubature<double> > cellCubature = cubFactory.create( (*topo_iterator), 4); 00382 int cubDim = cellCubature -> getDimension(); 00383 int numPts = cellCubature -> getNumPoints(); 00384 cubPoints.resize(numPts, cubDim); 00385 cubWeights.resize(numPts); 00386 cellCubature -> getCubature(cubPoints, cubWeights); 00387 00388 // 2. Define a cell workset by perturbing the cellWorkset of the reference cell with the specified topology 00389 // 2.1 Resize dimensions of the rank-3 (C,N,D) cell workset array for the current topology 00390 int numNodes = (*topo_iterator).getNodeCount(); 00391 int cellDim = (*topo_iterator).getDimension(); 00392 cellWorkset.resize(numCells, numNodes, cellDim); 00393 00394 // 2.2 Copy cellWorkset of the reference cell with the same topology to temp rank-2 (N,D) array 00395 FieldContainer<double> refCellNodes(numNodes, cellDim ); 00396 CellTools::getReferenceSubcellNodes(refCellNodes, cellDim, 0, (*topo_iterator) ); 00397 00398 // 2.3 Create randomly perturbed version of the reference cell and save in the cell workset array 00399 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00400 00401 // Move vertices +/-0.125 along their axes. Gives nondegenerate cells for base and extended topologies 00402 for(int nodeOrd = 0; nodeOrd < numNodes; nodeOrd++){ 00403 for(int d = 0; d < cellDim; d++){ 00404 double delta = Teuchos::ScalarTraits<double>::random()/16.0; 00405 cellWorkset(cellOrd, nodeOrd, d) = refCellNodes(nodeOrd, d) + delta; 00406 } // d 00407 }// nodeOrd 00408 }// cellOrd 00409 /* 00410 * 3.1 Test 1: single point set to single physical cell: map ref. point set in rank-2 (P,D) array 00411 * to a physical point set in rank-2 (P,D) array for a specified cell ordinal. Use the cub. 00412 * points array for this test. Resize physPoints and controlPoints to rank-2 (P,D) arrays. 00413 */ 00414 physPoints.resize(numPts, cubDim); 00415 controlPoints.resize(numPts, cubDim); 00416 00417 *outStream 00418 << " Mapping a set of " << numPts << " points to one cell in a workset of " << numCells << " " 00419 << (*topo_iterator).getName() << " cells. \n"; 00420 00421 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00422 00423 // Forward map:: requires cell ordinal 00424 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator), cellOrd); 00425 // Inverse map: requires cell ordinal. Use cubPoints as initial guess 00426 CellTools::mapToReferenceFrameInitGuess(controlPoints, cubPoints, physPoints, cellWorkset, (*topo_iterator), cellOrd); 00427 00428 // Points in controlPoints should match the originals in cubPoints up to a tolerance 00429 for(int pt = 0; pt < numPts; pt++){ 00430 for(int d = 0; d < cellDim; d++){ 00431 00432 if( abs( controlPoints(pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00433 errorFlag++; 00434 *outStream 00435 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00436 << " Mapping a single point set to a single physical cell in a workset failed for: \n" 00437 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00438 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00439 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00440 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00441 << " Original value = " << cubPoints(pt, d) << "\n" 00442 << " F^{-1}F(P_d) = " << controlPoints(pt, d) <<"\n"; 00443 } 00444 }// d 00445 }// pt 00446 }// cellOrd 00447 /* 00448 * 3.2 Test 2: single point set to multiple physical cells: map ref. point set in rank-2 (P,D) array 00449 * to a physical point set in rank-3 (C, P,D) array for all cell ordinals. Use the cub. 00450 * points array for this test. Resize physPoints and controlPoints to rank-3 (C,P,D) arrays. 00451 */ 00452 physPoints.clear(); 00453 controlPoints.clear(); 00454 physPoints.resize(numCells, numPts, cubDim); 00455 controlPoints.resize(numCells, numPts, cubDim); 00456 00457 // Clone cubature points in initialGuess: 00458 initialGuess.resize(numCells, numPts, cubDim); 00459 for(int c = 0; c < numCells; c++){ 00460 for(int pt = 0; pt < numPts; pt++){ 00461 for(int d = 0; d < cellDim; d++){ 00462 initialGuess(c, pt, d) = cubPoints(pt, d); 00463 }// d 00464 }// pt 00465 }// c 00466 00467 *outStream 00468 << " Mapping a set of " << numPts << " points to all cells in workset of " << numCells << " " 00469 << (*topo_iterator).getName() << " cells. \n"; 00470 00471 // Forward map: do not specify cell ordinal 00472 CellTools::mapToPhysicalFrame(physPoints, cubPoints, cellWorkset, (*topo_iterator)); 00473 // Inverse map: do not specify cell ordinal 00474 CellTools::mapToReferenceFrameInitGuess(controlPoints, initialGuess, physPoints, cellWorkset, (*topo_iterator)); 00475 00476 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance 00477 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00478 for(int pt = 0; pt < numPts; pt++){ 00479 for(int d = 0; d < cellDim; d++){ 00480 00481 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00482 errorFlag++; 00483 *outStream 00484 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00485 << " Mapping a single point set to all physical cells in a workset failed for: \n" 00486 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00487 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00488 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00489 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00490 << " Original value = " << cubPoints(pt, d) << "\n" 00491 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n"; 00492 } 00493 }// d 00494 }// pt 00495 }// cellOrd 00496 /* 00497 * 3.3 Test 3: multiple point sets to multiple physical cells: map ref. point sets in rank-3 (C,P,D) array 00498 * to physical point sets in rank-3 (C, P,D) array for all cell ordinals. The initialGuess 00499 * array from last test is used as the required (C,P,D) array of reference points for the 00500 * forward map and as the user-defined initial guess array for the inverse map 00501 */ 00502 physPoints.clear(); 00503 controlPoints.clear(); 00504 physPoints.resize(numCells, numPts, cubDim); 00505 controlPoints.resize(numCells, numPts, cubDim); 00506 00507 *outStream 00508 << " Mapping " << numCells << " sets of " << numPts << " points to corresponding cells in workset of " << numCells << " " 00509 << (*topo_iterator).getName() << " cells. \n"; 00510 00511 // Forward map: do not specify cell ordinal 00512 CellTools::mapToPhysicalFrame(physPoints, initialGuess, cellWorkset, (*topo_iterator)); 00513 // Inverse map: do not specify cell ordinal 00514 CellTools::mapToReferenceFrameInitGuess(controlPoints, initialGuess, physPoints, cellWorkset, (*topo_iterator)); 00515 00516 // Check: points in controlPoints should match the originals in cubPoints up to a tolerance 00517 for(int cellOrd = 0; cellOrd < numCells; cellOrd++){ 00518 for(int pt = 0; pt < numPts; pt++){ 00519 for(int d = 0; d < cellDim; d++){ 00520 00521 if( abs( controlPoints(cellOrd, pt, d) - cubPoints(pt, d) ) > 100.0*INTREPID_TOL ){ 00522 errorFlag++; 00523 *outStream 00524 << std::setw(70) << "^^^^----FAILURE!" << "\n" 00525 << " Mapping multiple point sets to corresponding physical cells in a workset failed for: \n" 00526 << " Cell Topology = " << (*topo_iterator).getName() << "\n" 00527 << " Physical cell ordinal in workset = " << cellOrd << "\n" 00528 << " Reference point ordinal = " << setprecision(12) << pt << "\n" 00529 << " At reference point coordinate = " << setprecision(12) << d << "\n" 00530 << " Original value = " << initialGuess(cellOrd, pt, d) << "\n" 00531 << " F^{-1}F(P_d) = " << controlPoints(cellOrd, pt, d) <<"\n"; 00532 } 00533 }// d 00534 }// pt 00535 }// cellOrd 00536 } //topo-iterator 00537 }// try user-defined initial guess 00538 00539 /************************************************************************************************* 00540 * Wrap up test: check if the test broke down unexpectedly due to an exception * 00541 ************************************************************************************************/ 00542 00543 catch (std::logic_error err) { 00544 *outStream << err.what() << "\n"; 00545 errorFlag = -1000; 00546 }; 00547 00548 *outStream \ 00549 << "\n" 00550 << "===============================================================================\n"\ 00551 << "| Test 3: Exception testing - only when HAVE_INTREPID_DEBUG is defined. |\n"\ 00552 << "===============================================================================\n\n"; 00553 /* 00554 * Test summary: 00555 * Calls methods of CellTools class with incorrectly configured arguments. This test is run only 00556 * in debug mode because many of the exceptions are checked only in that mode. 00557 * 00558 */ 00559 00560 // Initialize throw counter for exception testing 00561 int nException = 0; 00562 int throwCounter = 0; 00563 00564 try { 00565 00566 #ifdef HAVE_INTREPID_DEBUG 00567 // Some arbitrary dimensions 00568 int C = 10; 00569 int P = 21; 00570 int N; 00571 int D; 00572 int V; 00573 00574 // Array arguments 00575 FieldContainer<double> jacobian; 00576 FieldContainer<double> jacobianInv; 00577 FieldContainer<double> jacobianDet; 00578 FieldContainer<double> points; 00579 FieldContainer<double> cellWorkset; 00580 FieldContainer<double> physPoints; 00581 FieldContainer<double> refPoints; 00582 FieldContainer<double> initGuess; 00583 00584 /*********************************************************************************************** 00585 * Exception tests for setJacobian method * 00586 **********************************************************************************************/ 00587 00588 // Use the second cell topology for these tests (Triangle<6>) 00589 topo_iterator = supportedTopologies.begin() + 1; 00590 D = (*topo_iterator).getDimension(); 00591 N = (*topo_iterator).getNodeCount(); 00592 V = (*topo_iterator).getVertexCount(); 00593 00594 // 1. incorrect jacobian rank 00595 jacobian.resize(C, P, D); 00596 points.resize(P, D); 00597 cellWorkset.resize(C, N, D); 00598 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00599 throwCounter, nException ); 00600 00601 // 2. Incorrect cellWorkset rank 00602 cellWorkset.resize(C, D); 00603 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00604 throwCounter, nException ); 00605 00606 // 3. Incorrect points rank 00607 cellWorkset.resize(C, N, D); 00608 points.resize(D); 00609 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00610 throwCounter, nException ); 00611 00612 // 4. points rank incompatible with whichCell = valid cell ordinal 00613 points.resize(C, P, D); 00614 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator), 0 ), 00615 throwCounter, nException ); 00616 00617 // 5. Non-matching dim 00618 jacobian.resize(C, P, D, D); 00619 points.resize(C, P, D - 1); 00620 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00621 throwCounter, nException ); 00622 00623 // 6. Non-matching dim 00624 jacobian.resize(C, P, D, D); 00625 points.resize(C, P - 1, D); 00626 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00627 throwCounter, nException ); 00628 00629 // 7. Non-matching dim 00630 jacobian.resize(C, P, D, D); 00631 points.resize(C - 1, P, D); 00632 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00633 throwCounter, nException ); 00634 00635 // 8. Non-matching dim 00636 jacobian.resize(C, P, D, D); 00637 points.resize(C, P, D); 00638 cellWorkset.resize(C, N, D - 1); 00639 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00640 throwCounter, nException ); 00641 00642 // 9. Non-matching dim 00643 jacobian.resize(C, P, D, D); 00644 points.resize(C, P, D); 00645 cellWorkset.resize(C - 1, N, D); 00646 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00647 throwCounter, nException ); 00648 00649 // 10. Incompatible ranks 00650 jacobian.resize(C, D, D); 00651 points.resize(C, P, D); 00652 cellWorkset.resize(C, N, D); 00653 INTREPID_TEST_COMMAND( CellTools::setJacobian(jacobian, points, cellWorkset, (*topo_iterator) ), 00654 throwCounter, nException ); 00655 00656 /*********************************************************************************************** 00657 * Exception tests for setJacobianInv method * 00658 **********************************************************************************************/ 00659 00660 // 11. incompatible ranks 00661 jacobian.resize(C, P, D, D); 00662 jacobianInv.resize(P, D, D); 00663 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00664 throwCounter, nException ); 00665 00666 // 12. incorrect ranks 00667 jacobian.resize(D, D); 00668 jacobianInv.resize(D, D); 00669 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00670 throwCounter, nException ); 00671 00672 // 13. nonmatching dimensions 00673 jacobian.resize(C, P, D, D - 1); 00674 jacobianInv.resize(C, P, D, D); 00675 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00676 throwCounter, nException ); 00677 00678 // 14. nonmatching dimensions 00679 jacobian.resize(C, P, D - 1, D); 00680 jacobianInv.resize(C, P, D, D); 00681 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00682 throwCounter, nException ); 00683 00684 // 15. nonmatching dimensions 00685 jacobian.resize(C, P - 1, D, D); 00686 jacobianInv.resize(C, P, D, D); 00687 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00688 throwCounter, nException ); 00689 00690 // 16. nonmatching dimensions 00691 jacobian.resize(C - 1, P, D, D); 00692 jacobianInv.resize(C, P, D, D); 00693 INTREPID_TEST_COMMAND( CellTools::setJacobianInv(jacobianInv, jacobian), 00694 throwCounter, nException ); 00695 00696 /*********************************************************************************************** 00697 * Exception tests for setJacobianDet method * 00698 **********************************************************************************************/ 00699 00700 // 17. Incompatible ranks 00701 jacobian.resize(C, P, D, D); 00702 jacobianDet.resize(C, P, D); 00703 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00704 throwCounter, nException ); 00705 00706 // 18. Incompatible ranks 00707 jacobian.resize(P, D, D); 00708 jacobianDet.resize(C, P); 00709 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00710 throwCounter, nException ); 00711 00712 // 19. Incorrect rank 00713 jacobian.resize(D, D); 00714 jacobianDet.resize(C, P); 00715 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00716 throwCounter, nException ); 00717 00718 // 20. Incorrect rank 00719 jacobian.resize(C, P, D, D); 00720 jacobianDet.resize(C); 00721 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00722 throwCounter, nException ); 00723 00724 // 21. Incorrect dimension 00725 jacobian.resize(C, P, D, D); 00726 jacobianDet.resize(C, P-1); 00727 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00728 throwCounter, nException ); 00729 00730 // 22. Incorrect dimension 00731 jacobian.resize(C - 1, P, D, D); 00732 jacobianDet.resize(C, P); 00733 INTREPID_TEST_COMMAND( CellTools::setJacobianDet(jacobianDet, jacobian), 00734 throwCounter, nException ); 00735 00736 /*********************************************************************************************** 00737 * Exception tests for mapToPhysicalFrame method * 00738 **********************************************************************************************/ 00739 00740 // 23. Incorrect refPoint rank 00741 refPoints.resize(P); 00742 physPoints.resize(P, D); 00743 cellWorkset.resize(C, N, D); 00744 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00745 throwCounter, nException ); 00746 // 24. Incorrect workset rank 00747 cellWorkset.resize(P, D); 00748 refPoints.resize(P, D); 00749 physPoints.resize(C, P, D); 00750 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00751 throwCounter, nException ); 00752 00753 // 25. Incompatible ranks 00754 refPoints.resize(C, P, D); 00755 physPoints.resize(P, D); 00756 cellWorkset.resize(C, N, D); 00757 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00758 throwCounter, nException ); 00759 00760 // 26. Incompatible dimensions 00761 refPoints.resize(C, P, D); 00762 physPoints.resize(C, P, D - 1); 00763 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00764 throwCounter, nException ); 00765 00766 // 27. Incompatible dimensions 00767 refPoints.resize(C, P, D); 00768 physPoints.resize(C, P - 1, D); 00769 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00770 throwCounter, nException ); 00771 00772 // 28. Incompatible dimensions 00773 refPoints.resize(C, P, D); 00774 physPoints.resize(C - 1, P, D); 00775 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator) ), 00776 throwCounter, nException ); 00777 00778 // 29. Incorrect physPoints rank when whichCell is valid cell ordinal 00779 refPoints.resize(P, D); 00780 physPoints.resize(C, P, D); 00781 INTREPID_TEST_COMMAND( CellTools::mapToPhysicalFrame(physPoints, refPoints, cellWorkset, (*topo_iterator), 0 ), 00782 throwCounter, nException ); 00783 00784 /*********************************************************************************************** 00785 * Exception tests for mapToReferenceFrame method (with default initial guesses) * 00786 **********************************************************************************************/ 00787 00788 // 30. incompatible ranks 00789 refPoints.resize(C, P, D); 00790 physPoints.resize(P, D); 00791 cellWorkset.resize(C, N, D); 00792 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator) ), 00793 throwCounter, nException ); 00794 00795 // 31. Incompatible ranks with whichCell = valid cell ordinal 00796 refPoints.resize(C, P, D); 00797 physPoints.resize(C, P, D); 00798 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator), 0 ), 00799 throwCounter, nException ); 00800 00801 // 32. Incompatible ranks with whichCell = -1 (default) 00802 refPoints.resize(P, D); 00803 physPoints.resize(P, D); 00804 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)), 00805 throwCounter, nException ); 00806 00807 // 33. Nonmatching dimensions 00808 refPoints.resize(C, P, D - 1); 00809 physPoints.resize(C, P, D); 00810 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)), 00811 throwCounter, nException ); 00812 00813 // 34. Nonmatching dimensions 00814 refPoints.resize(C, P - 1, D); 00815 physPoints.resize(C, P, D); 00816 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)), 00817 throwCounter, nException ); 00818 00819 // 35. Nonmatching dimensions 00820 refPoints.resize(C - 1, P, D); 00821 physPoints.resize(C, P, D); 00822 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)), 00823 throwCounter, nException ); 00824 00825 // 36. Incorrect rank for cellWorkset 00826 refPoints.resize(C, P, D); 00827 physPoints.resize(C, P, D); 00828 cellWorkset.resize(C, N); 00829 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrame(refPoints, physPoints, cellWorkset, (*topo_iterator)), 00830 throwCounter, nException ); 00831 00832 /*********************************************************************************************** 00833 * Exception tests for mapToReferenceFrameInitGuess method (initial guess is a parameter) * 00834 **********************************************************************************************/ 00835 00836 // 37. Incompatible ranks 00837 refPoints.resize(C, P, D); 00838 physPoints.resize(C, P, D); 00839 initGuess.resize(P, D); 00840 cellWorkset.resize(C, N, D); 00841 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator) ), 00842 throwCounter, nException ); 00843 00844 // 38. Incompatible ranks when whichCell is valid ordinal 00845 refPoints.resize(P, D); 00846 physPoints.resize(P, D); 00847 initGuess.resize(C, P, D); 00848 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator), 0), 00849 throwCounter, nException ); 00850 00851 // 39. Nonmatching dimensions 00852 refPoints.resize(C, P, D); 00853 physPoints.resize(C, P, D); 00854 initGuess.resize(C, P, D - 1); 00855 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)), 00856 throwCounter, nException ); 00857 00858 // 40. Nonmatching dimensions 00859 initGuess.resize(C, P - 1, D); 00860 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)), 00861 throwCounter, nException ); 00862 00863 // 41. Nonmatching dimensions 00864 initGuess.resize(C - 1, P, D); 00865 INTREPID_TEST_COMMAND( CellTools::mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, (*topo_iterator)), 00866 throwCounter, nException ); 00867 00868 /*********************************************************************************************** 00869 * Exception tests for mapToReferenceSubcell method * 00870 **********************************************************************************************/ 00871 00872 FieldContainer<double> refSubcellPoints; 00873 FieldContainer<double> paramPoints; 00874 int subcellDim = 2; 00875 int subcellOrd = 0; 00876 00877 // This should set cell topology to Tetrahedron<10> so that we have real edges and faces. 00878 topo_iterator += 5; 00879 D = (*topo_iterator).getDimension(); 00880 00881 // 42. Incorrect array rank 00882 refSubcellPoints.resize(P,3); 00883 paramPoints.resize(P); 00884 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)), 00885 throwCounter, nException ); 00886 00887 // 43. Incorrect array rank 00888 refSubcellPoints.resize(P); 00889 paramPoints.resize(P, 2); 00890 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)), 00891 throwCounter, nException ); 00892 00893 // 44. Incorrect array dimension for face of 3D cell (should be 3) 00894 refSubcellPoints.resize(P, 2); 00895 paramPoints.resize(P, 2); 00896 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)), 00897 throwCounter, nException ); 00898 00899 // 45. Incorrect array dimension for parametrization domain of a face of 3D cell (should be 2) 00900 refSubcellPoints.resize(P, 3); 00901 paramPoints.resize(P, 3); 00902 INTREPID_TEST_COMMAND( CellTools::mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellDim, subcellOrd, (*topo_iterator)), 00903 throwCounter, nException ); 00904 00905 /*********************************************************************************************** 00906 * Exception tests for getReferenceEdgeTangent method * 00907 **********************************************************************************************/ 00908 00909 FieldContainer<double> refEdgeTangent; 00910 00911 // 46. Incorrect rank 00912 refEdgeTangent.resize(C,P,D); 00913 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 0, (*topo_iterator)), 00914 throwCounter, nException ); 00915 00916 // 47. Incorrect dimension D for Tet<10> cell 00917 refEdgeTangent.resize(2); 00918 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 0, (*topo_iterator)), 00919 throwCounter, nException ); 00920 00921 // 48. Invalid edge ordinal for Tet<10> 00922 refEdgeTangent.resize(C,P,D); 00923 INTREPID_TEST_COMMAND( CellTools::getReferenceEdgeTangent(refEdgeTangent, 10, (*topo_iterator)), 00924 throwCounter, nException ); 00925 00926 /*********************************************************************************************** 00927 * Exception tests for getReferenceFaceTangents method * 00928 **********************************************************************************************/ 00929 00930 FieldContainer<double> refFaceTanU; 00931 FieldContainer<double> refFaceTanV; 00932 00933 // 49. Incorrect rank 00934 refFaceTanU.resize(P, D); 00935 refFaceTanV.resize(D); 00936 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)), 00937 throwCounter, nException ); 00938 00939 // 50. Incorrect rank 00940 refFaceTanU.resize(D); 00941 refFaceTanV.resize(P, D); 00942 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)), 00943 throwCounter, nException ); 00944 00945 // 51. Incorrect dimension for 3D cell 00946 refFaceTanU.resize(D - 1); 00947 refFaceTanV.resize(D); 00948 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)), 00949 throwCounter, nException ); 00950 00951 // 52. Incorrect dimension for 3D cell 00952 refFaceTanU.resize(D); 00953 refFaceTanV.resize(D - 1); 00954 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 0, (*topo_iterator)), 00955 throwCounter, nException ); 00956 00957 // 53. Invalid face ordinal 00958 refFaceTanU.resize(D); 00959 refFaceTanV.resize(D); 00960 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceTangents(refFaceTanU, refFaceTanV, 10, (*topo_iterator)), 00961 throwCounter, nException ); 00962 00963 /*********************************************************************************************** 00964 * Exception tests for getReferenceSide/FaceNormal methods * 00965 **********************************************************************************************/ 00966 00967 FieldContainer<double> refSideNormal; 00968 00969 // 54-55. Incorrect rank 00970 refSideNormal.resize(C,P); 00971 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)), 00972 throwCounter, nException ); 00973 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)), 00974 throwCounter, nException ); 00975 00976 // 56-57. Incorrect dimension for 3D cell 00977 refSideNormal.resize(D - 1); 00978 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)), 00979 throwCounter, nException ); 00980 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)), 00981 throwCounter, nException ); 00982 00983 // 58-59. Invalid side ordinal for Tet<10> 00984 refSideNormal.resize(D); 00985 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 10, (*topo_iterator)), 00986 throwCounter, nException ); 00987 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 10, (*topo_iterator)), 00988 throwCounter, nException ); 00989 00990 // 60. Incorrect dimension for 2D cell: reset topo_iterator to the first cell in supportedTopologies which is Tri<3> 00991 topo_iterator = supportedTopologies.begin(); 00992 D = (*topo_iterator).getDimension(); 00993 refSideNormal.resize(D - 1); 00994 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 0, (*topo_iterator)), 00995 throwCounter, nException ); 00996 00997 // 61. Invalid side ordinal for Tri<3> 00998 refSideNormal.resize(D); 00999 INTREPID_TEST_COMMAND( CellTools::getReferenceSideNormal(refSideNormal, 10, (*topo_iterator)), 01000 throwCounter, nException ); 01001 01002 // 62. Cannot call the "face" method for 2D cells 01003 refSideNormal.resize(D); 01004 INTREPID_TEST_COMMAND( CellTools::getReferenceFaceNormal(refSideNormal, 0, (*topo_iterator)), 01005 throwCounter, nException ); 01006 01007 /*********************************************************************************************** 01008 * Exception tests for checkPoint/Pointset/PointwiseInclusion methods * 01009 **********************************************************************************************/ 01010 points.resize(2,3,3,4); 01011 FieldContainer<int> inCell; 01012 01013 // 63. Point dimension does not match cell topology 01014 double * point = 0; 01015 INTREPID_TEST_COMMAND(CellTools::checkPointInclusion(point, (*topo_iterator).getDimension() + 1, (*topo_iterator) ), 01016 throwCounter, nException ); 01017 01018 // 64. Invalid cell topology 01019 CellTopology pentagon_5(shards::getCellTopologyData<shards::Pentagon<> >() ); 01020 INTREPID_TEST_COMMAND(CellTools::checkPointInclusion(point, pentagon_5.getDimension(), pentagon_5 ), 01021 throwCounter, nException ); 01022 01023 // 65. Incorrect spatial dimension of points 01024 points.resize(10, 10, (*topo_iterator).getDimension() + 1); 01025 INTREPID_TEST_COMMAND(CellTools::checkPointsetInclusion(points, (*topo_iterator) ), 01026 throwCounter, nException ); 01027 01028 // 66. Incorrect rank of input array 01029 points.resize(10,10,10,3); 01030 INTREPID_TEST_COMMAND(CellTools::checkPointsetInclusion(points, (*topo_iterator) ), 01031 throwCounter, nException ); 01032 01033 // 67. Incorrect rank of output array 01034 points.resize(10,10,(*topo_iterator).getDimension() ); 01035 inCell.resize(10); 01036 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01037 throwCounter, nException ); 01038 01039 // 68. Incorrect rank of output array 01040 points.resize(10, (*topo_iterator).getDimension() ); 01041 inCell.resize(10, 10); 01042 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01043 throwCounter, nException ); 01044 01045 // 69. Incorrect rank of output array 01046 points.resize((*topo_iterator).getDimension() ); 01047 inCell.resize(10, 10); 01048 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01049 throwCounter, nException ); 01050 01051 // 70. Incorrect dimension of output array 01052 points.resize(10, 10, (*topo_iterator).getDimension() ); 01053 inCell.resize(10, 9); 01054 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01055 throwCounter, nException ); 01056 01057 // 71. Incorrect dimension of output array 01058 points.resize(10, 10, (*topo_iterator).getDimension() ); 01059 inCell.resize(9, 10); 01060 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01061 throwCounter, nException ); 01062 01063 // 72. Incorrect dimension of output array 01064 points.resize(10, (*topo_iterator).getDimension() ); 01065 inCell.resize(9); 01066 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01067 throwCounter, nException ); 01068 01069 // 73. Incorrect spatial dimension of input array 01070 points.resize(10, 10, (*topo_iterator).getDimension() + 1); 01071 inCell.resize(10, 10); 01072 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01073 throwCounter, nException ); 01074 01075 // 74. Incorrect rank of input array. 01076 points.resize(10,10,10,3); 01077 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, points, (*topo_iterator) ), 01078 throwCounter, nException ); 01079 01080 01081 physPoints.resize(C, P, D); 01082 inCell.resize(C, P); 01083 // 75. Invalid rank of cellWorkset 01084 cellWorkset.resize(C, N, D, D); 01085 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ), 01086 throwCounter, nException ); 01087 01088 // 76. Invalid dimension 1 (node count) of cellWorkset 01089 cellWorkset.resize(C, N + 1, D); 01090 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ), 01091 throwCounter, nException ); 01092 01093 // 77. Invalid dimension 2 (spatial dimension) of cellWorkset 01094 cellWorkset.resize(C, N, D + 1); 01095 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ), 01096 throwCounter, nException ); 01097 01098 // 78. Invalid whichCell value (exceeds cell count in the workset) 01099 cellWorkset.resize(C, N, D); 01100 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), C + 1 ), 01101 throwCounter, nException ); 01102 01103 // 79. Invalid whichCell for rank-3 physPoints (must be -1, here it is valid cell ordinal) 01104 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0 ), 01105 throwCounter, nException ); 01106 01107 // 80. Invalid whichCell for rank-2 physPoints (must be a valid cell ordinal, here it is the default -1) 01108 physPoints.resize(P, D); 01109 inCell.resize(P); 01110 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator) ), 01111 throwCounter, nException ); 01112 01113 // 81. Incompatible ranks of I/O arrays 01114 physPoints.resize(C, P, D); 01115 inCell.resize(P); 01116 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)), 01117 throwCounter, nException ); 01118 01119 // 82. Incompatible ranks of I/O arrays 01120 physPoints.resize(P, D); 01121 inCell.resize(C, P); 01122 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0), 01123 throwCounter, nException ); 01124 01125 // 83. Incompatible dimensions of I/O arrays 01126 physPoints.resize(C, P, D); 01127 inCell.resize(C, P + 1); 01128 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)), 01129 throwCounter, nException ); 01130 01131 // 84. Incompatible dimensions of I/O arrays: rank-3 Input 01132 physPoints.resize(C + 1, P, D); 01133 inCell.resize(C, P); 01134 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator)), 01135 throwCounter, nException ); 01136 01137 // 85. Incompatible dimensions of I/O arrays: rank-2 Input 01138 physPoints.resize(P, D); 01139 inCell.resize(P + 1); 01140 INTREPID_TEST_COMMAND(CellTools::checkPointwiseInclusion(inCell, physPoints, cellWorkset, (*topo_iterator), 0 ), 01141 throwCounter, nException ); 01142 01143 01144 /*********************************************************************************************** 01145 * Exception tests for getReferenceVertex/vertices/Node/Nodes methods * 01146 **********************************************************************************************/ 01147 01148 FieldContainer<double> subcellNodes; 01149 01150 // 86-89. Cell does not have reference cell 01151 INTREPID_TEST_COMMAND(CellTools::getReferenceVertex(pentagon_5, 0), throwCounter, nException); 01152 INTREPID_TEST_COMMAND(CellTools::getReferenceNode(pentagon_5, 0), throwCounter, nException); 01153 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, 0, 0, pentagon_5), throwCounter, nException); 01154 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, 0, 0, pentagon_5), throwCounter, nException); 01155 01156 // Use last cell topology (Wedge<18>) for these tests 01157 topo_iterator = supportedTopologies.end() - 1; 01158 D = (*topo_iterator).getDimension(); 01159 int subcDim = D - 1; 01160 int S = (*topo_iterator).getSubcellCount(subcDim); 01161 V = (*topo_iterator).getVertexCount(subcDim, S - 1); 01162 subcellNodes.resize(V, D); 01163 // 90. subcell ordinal out of range 01164 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S + 1, (*topo_iterator)), 01165 throwCounter, nException); 01166 01167 // 91. subcell dim out of range 01168 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, D + 1, S, (*topo_iterator)), 01169 throwCounter, nException); 01170 01171 // 92. Incorrect rank for subcellNodes 01172 subcellNodes.resize(V, D, D); 01173 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01174 throwCounter, nException); 01175 01176 // 93. Incorrect dimension for subcellNodes 01177 subcellNodes.resize(V - 1, D); 01178 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01179 throwCounter, nException); 01180 01181 // 94. Incorrect dimension for subcellNodes 01182 subcellNodes.resize(V, D - 1); 01183 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellVertices(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01184 throwCounter, nException); 01185 01186 01187 N = (*topo_iterator).getNodeCount(subcDim, S - 1); 01188 subcellNodes.resize(N, D); 01189 // 95. subcell ordinal out of range 01190 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S + 1, (*topo_iterator)), 01191 throwCounter, nException); 01192 01193 // 96. subcell dim out of range 01194 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, D + 1, S, (*topo_iterator)), 01195 throwCounter, nException); 01196 01197 // 97. Incorrect rank for subcellNodes 01198 subcellNodes.resize(N, D, D); 01199 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01200 throwCounter, nException); 01201 01202 // 98. Incorrect dimension for subcellNodes 01203 subcellNodes.resize(N - 1, D); 01204 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01205 throwCounter, nException); 01206 01207 // 99. Incorrect dimension for subcellNodes 01208 subcellNodes.resize(N, D - 1); 01209 INTREPID_TEST_COMMAND(CellTools::getReferenceSubcellNodes(subcellNodes, subcDim, S - 1, (*topo_iterator)), 01210 throwCounter, nException); 01211 01212 #endif 01213 } // try exception testing 01214 01215 /************************************************************************************************* 01216 * Wrap up test: check if the test broke down unexpectedly due to an exception * 01217 ************************************************************************************************/ 01218 01219 catch(std::logic_error err) { 01220 *outStream << err.what() << "\n"; 01221 errorFlag = -1000; 01222 } 01223 01224 // Check if number of thrown exceptions matches the one we expect 01225 if (throwCounter != nException) { 01226 errorFlag++; 01227 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 01228 } 01229 01230 01231 if (errorFlag != 0) 01232 std::cout << "End Result: TEST FAILED\n"; 01233 else 01234 std::cout << "End Result: TEST PASSED\n"; 01235 01236 // reset format state of std::cout 01237 std::cout.copyfmt(oldFormatState); 01238 01239 return errorFlag; 01240 } 01241 01242 01243 01244 01245 01246 01247
1.7.6.1