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