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