Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/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 // 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