Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Cell/Intrepid_CellToolsDef.hpp
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 #ifndef INTREPID_CELLTOOLSDEF_HPP
00050 #define INTREPID_CELLTOOLSDEF_HPP
00051 
00052 // disable clang warnings
00053 #ifdef __clang__
00054 #pragma clang system_header
00055 #endif
00056 
00057 
00058 namespace Intrepid {
00059 
00060   template<class Scalar>
00061   const FieldContainer<double>& CellTools<Scalar>::getSubcellParametrization(const int subcellDim,
00062                                                                              const shards::CellTopology& parentCell){
00063     
00064 #ifdef HAVE_INTREPID_DEBUG
00065     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
00066                         ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): the specified cell topology does not have a reference cell.");
00067 
00068     TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
00069                            ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): parametrization defined only for 1 and 2-dimensional subcells.");    
00070 #endif
00071     
00072     // Coefficients of the coordinate functions defining the parametrization maps are stored in 
00073     // rank-3 arrays with dimensions (SC, PCD, COEF) where:
00074     //  - SC    is the subcell count of subcells with the specified dimension in the parent cell
00075     //  - PCD   is Parent Cell Dimension, which gives the number of coordinate functions in the map:
00076     //          PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
00077     //          PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
00078     //  - COEF  is number of coefficients needed to specify a coordinate function:
00079     //          COEFF = 2 for edge parametrizations
00080     //          COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
00081     //          are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
00082     //          3 coefficients are sufficient to store Quad face parameterization maps.
00083     //
00084     // Arrays are sized and filled only when parametrization of a particular subcell is requested
00085     // by setSubcellParametrization.
00086     
00087     // Edge maps for 2D non-standard cells: ShellLine and Beam
00088     static FieldContainer<double> lineEdges;        static int lineEdgesSet = 0;
00089     
00090     // Edge maps for 2D standard cells: Triangle and Quadrilateral
00091     static FieldContainer<double> triEdges;         static int triEdgesSet  = 0;
00092     static FieldContainer<double> quadEdges;        static int quadEdgesSet = 0;
00093 
00094     // Edge maps for 3D non-standard cells: Shell Tri and Quad
00095     static FieldContainer<double> shellTriEdges;    static int shellTriEdgesSet  = 0;
00096     static FieldContainer<double> shellQuadEdges;   static int shellQuadEdgesSet = 0;
00097     
00098     // Edge maps for 3D standard cells:
00099     static FieldContainer<double> tetEdges;         static int tetEdgesSet = 0;
00100     static FieldContainer<double> hexEdges;         static int hexEdgesSet = 0;
00101     static FieldContainer<double> pyrEdges;         static int pyrEdgesSet = 0;
00102     static FieldContainer<double> wedgeEdges;       static int wedgeEdgesSet = 0;
00103 
00104 
00105     // Face maps for 3D non-standard cells: Shell Triangle and Quadrilateral
00106     static FieldContainer<double> shellTriFaces;    static int shellTriFacesSet  = 0;
00107     static FieldContainer<double> shellQuadFaces;   static int shellQuadFacesSet = 0;
00108 
00109     // Face maps for 3D standard cells:
00110     static FieldContainer<double> tetFaces;         static int tetFacesSet = 0;
00111     static FieldContainer<double> hexFaces;         static int hexFacesSet = 0;
00112     static FieldContainer<double> pyrFaces;         static int pyrFacesSet = 0;
00113     static FieldContainer<double> wedgeFaces;       static int wedgeFacesSet = 0;
00114 
00115     // Select subcell parametrization according to its parent cell type
00116     switch(parentCell.getKey() ) {
00117       
00118       // Tet cells
00119       case shards::Tetrahedron<4>::key:
00120       case shards::Tetrahedron<8>::key:
00121       case shards::Tetrahedron<10>::key:
00122       case shards::Tetrahedron<11>::key:
00123         if(subcellDim == 2) {
00124           if(!tetFacesSet){
00125             setSubcellParametrization(tetFaces, subcellDim, parentCell);
00126             tetFacesSet = 1;
00127           }
00128           return tetFaces;
00129         }
00130         else if(subcellDim == 1) {
00131           if(!tetEdgesSet){
00132             setSubcellParametrization(tetEdges, subcellDim, parentCell);
00133             tetEdgesSet = 1;
00134           }
00135           return tetEdges;
00136         }
00137         else{
00138           TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument, 
00139                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Tet parametrizations defined for 1 and 2-subcells only");
00140         }
00141         break;
00142         
00143       // Hex cells
00144       case shards::Hexahedron<8>::key:
00145       case shards::Hexahedron<20>::key:
00146       case shards::Hexahedron<27>::key:
00147         if(subcellDim == 2) {
00148           if(!hexFacesSet){
00149             setSubcellParametrization(hexFaces, subcellDim, parentCell);
00150             hexFacesSet = 1;
00151           }
00152           return hexFaces;
00153         }
00154         else if(subcellDim == 1) {
00155           if(!hexEdgesSet){
00156             setSubcellParametrization(hexEdges, subcellDim, parentCell);
00157             hexEdgesSet = 1;
00158           }
00159           return hexEdges;
00160         }
00161         else{
00162           TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument, 
00163                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Hex parametrizations defined for 1 and 2-subcells only");
00164         }
00165         break;
00166         
00167       // Pyramid cells
00168       case shards::Pyramid<5>::key:
00169       case shards::Pyramid<13>::key:
00170       case shards::Pyramid<14>::key:
00171         if(subcellDim == 2) {
00172           if(!pyrFacesSet){
00173             setSubcellParametrization(pyrFaces, subcellDim, parentCell);
00174             pyrFacesSet = 1;
00175           }
00176           return pyrFaces;
00177         }
00178         else if(subcellDim == 1) {
00179           if(!pyrEdgesSet){
00180             setSubcellParametrization(pyrEdges, subcellDim, parentCell);
00181             pyrEdgesSet = 1;
00182           }
00183           return pyrEdges;
00184         }
00185         else {
00186           TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument, 
00187                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Pyramid parametrizations defined for 1 and 2-subcells only");
00188         }
00189         break;
00190         
00191       // Wedge cells
00192       case shards::Wedge<6>::key:
00193       case shards::Wedge<15>::key:
00194       case shards::Wedge<18>::key:
00195         if(subcellDim == 2) {
00196           if(!wedgeFacesSet){
00197             setSubcellParametrization(wedgeFaces, subcellDim, parentCell);
00198             wedgeFacesSet = 1;
00199           }
00200           return wedgeFaces;
00201         }
00202         else if(subcellDim == 1) {
00203           if(!wedgeEdgesSet){
00204             setSubcellParametrization(wedgeEdges, subcellDim, parentCell);
00205             wedgeEdgesSet = 1;
00206           }
00207           return wedgeEdges;
00208         }
00209         else {
00210           TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument, 
00211                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Wedge parametrization defined for 1 and 2-subcells only");
00212         }
00213         break;
00214       //
00215       // Standard 2D cells have only 1-subcells
00216       //
00217       case shards::Triangle<3>::key:
00218       case shards::Triangle<4>::key:
00219       case shards::Triangle<6>::key:
00220         if(subcellDim == 1) {
00221           if(!triEdgesSet){
00222             setSubcellParametrization(triEdges, subcellDim, parentCell);
00223             triEdgesSet = 1;
00224           }
00225           return triEdges;
00226         }
00227         else{
00228           TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00229                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Triangle parametrizations defined for 1-subcells only");
00230         }
00231         break;
00232                   
00233       case shards::Quadrilateral<4>::key:
00234       case shards::Quadrilateral<8>::key:
00235       case shards::Quadrilateral<9>::key:
00236         if(subcellDim == 1) {
00237           if(!quadEdgesSet){
00238             setSubcellParametrization(quadEdges, subcellDim, parentCell);
00239             quadEdgesSet = 1;
00240           }
00241           return quadEdges;
00242         }
00243         else{
00244           TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00245                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Quad parametrizations defined for 1-subcells only");
00246         }
00247         break;        
00248       //
00249       // Non-standard 3D Shell Tri and Quad cells have 1 and 2-subcells. Because they are 3D cells
00250       // can't reuse edge parametrization arrays for 2D Triangle and Quadrilateral.
00251       //
00252       case shards::ShellTriangle<3>::key:
00253       case shards::ShellTriangle<6>::key:
00254         if(subcellDim == 2) {
00255           if(!shellTriFacesSet){
00256             setSubcellParametrization(shellTriFaces, subcellDim, parentCell);
00257             shellTriFacesSet = 1;
00258           }
00259           return shellTriFaces;
00260         }
00261         else if(subcellDim == 1) {
00262           if(!shellTriEdgesSet){
00263             setSubcellParametrization(shellTriEdges, subcellDim, parentCell);
00264             shellTriEdgesSet = 1;
00265           }
00266           return shellTriEdges;
00267         }
00268         else if( subcellDim != 1 || subcellDim != 2){
00269           TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00270                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Shell Triangle parametrizations defined for 1 and 2-subcells only");
00271         }
00272         break;
00273         
00274       case shards::ShellQuadrilateral<4>::key:
00275       case shards::ShellQuadrilateral<8>::key:
00276       case shards::ShellQuadrilateral<9>::key:
00277         if(subcellDim == 2) {
00278           if(!shellQuadFacesSet){
00279             setSubcellParametrization(shellQuadFaces, subcellDim, parentCell);
00280             shellQuadFacesSet = 1;
00281           }
00282           return shellQuadFaces;
00283         }
00284         else if(subcellDim == 1) {
00285           if(!shellQuadEdgesSet){
00286             setSubcellParametrization(shellQuadEdges, subcellDim, parentCell);
00287             shellQuadEdgesSet = 1;
00288           }
00289           return shellQuadEdges;
00290         }
00291         else if( subcellDim != 1 || subcellDim != 2){
00292           TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00293                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Shell Quad parametrizations defined for 1 and 2-subcells only");
00294         }
00295         break;
00296         
00297       //
00298       // Non-standard 2D cells: Shell Lines and Beams have 1-subcells
00299       //
00300       case shards::ShellLine<2>::key:
00301       case shards::ShellLine<3>::key:
00302       case shards::Beam<2>::key:
00303       case shards::Beam<3>::key:
00304         if(subcellDim == 1) {
00305           if(!lineEdgesSet){
00306             setSubcellParametrization(lineEdges, subcellDim, parentCell);
00307             lineEdgesSet = 1;
00308           }
00309           return lineEdges;
00310         }
00311         else{
00312           TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00313                               ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): shell line/beam parametrizations defined for 1-subcells only");
00314         }
00315         break;        
00316 
00317       default:
00318         TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00319                             ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): invalid cell topology.");
00320     }//cell key
00321     // To disable compiler warning, should never be reached
00322     return lineEdges;
00323   }
00324   
00325   
00326   
00327   template<class Scalar>
00328   void CellTools<Scalar>::setSubcellParametrization(FieldContainer<double>&     subcellParametrization,
00329                                                     const int                   subcellDim,
00330                                                     const shards::CellTopology& parentCell)
00331   {
00332 #ifdef HAVE_INTREPID_DEBUG
00333     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
00334                         ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): the specified cell topology does not have a reference cell.");
00335 
00336     TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
00337                         ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): parametrization defined only for 1 and 2-dimensional subcells.");    
00338 #endif
00339                         // subcellParametrization is rank-3 FieldContainer with dimensions (SC, PCD, COEF) where:
00340     //  - SC    is the subcell count of subcells with the specified dimension in the parent cell
00341     //  - PCD   is Parent Cell Dimension, which gives the number of coordinate functions in the map
00342     //          PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
00343     //          PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
00344     //  - COEF  is number of coefficients needed to specify a coordinate function:
00345     //          COEFF = 2 for edge parametrizations
00346     //          COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
00347     //          are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
00348     //          3 coefficients are sufficient to store Quad face parameterization maps.
00349     //  
00350     // Edge parametrization maps [-1,1] to edge defined by (v0, v1)
00351     // Face parametrization maps [-1,1]^2 to quadrilateral face (v0, v1, v2, v3), or
00352     // standard 2-simplex  {(0,0),(1,0),(0,1)} to traingle face (v0, v1, v2).
00353     // This defines orientation-preserving parametrizations with respect to reference edge and
00354     // face orientations induced by their vertex order. 
00355 
00356     // get subcellParametrization dimensions: (sc, pcd, coeff)
00357     unsigned sc    = parentCell.getSubcellCount(subcellDim);
00358     unsigned pcd   = parentCell.getDimension();   
00359     unsigned coeff = (subcellDim == 1) ? 2 : 3;
00360     
00361     
00362     // Resize container
00363     subcellParametrization.resize(sc, pcd, coeff);
00364     
00365     // Edge parametrizations of 2D and 3D cells (shell lines and beams are 2D cells with edges)
00366     if(subcellDim == 1){      
00367       for(unsigned subcellOrd = 0; subcellOrd < sc; subcellOrd++){
00368         
00369         int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
00370         int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
00371         
00372         // vertexK[0] = x_k; vertexK[1] = y_k; vertexK[2] = z_k; z_k = 0 for 2D cells
00373         // Note that ShellLine and Beam are 2D cells!
00374         const double* v0 = getReferenceVertex(parentCell, v0ord);
00375         const double* v1 = getReferenceVertex(parentCell, v1ord);
00376         
00377         // x(t) = (x0 + x1)/2 + t*(x1 - x0)/2 
00378         subcellParametrization(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
00379         subcellParametrization(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
00380         
00381         // y(t) = (y0 + y1)/2 + t*(y1 - y0)/2 
00382         subcellParametrization(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
00383         subcellParametrization(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
00384         
00385         if( pcd == 3 ) {
00386           // z(t) = (z0 + z1)/2 + t*(z1 - z0)/2 
00387           subcellParametrization(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
00388           subcellParametrization(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
00389         }
00390       }// for loop over 1-subcells
00391     }
00392       
00393       // Face parametrizations of 3D cells: (shell Tri and Quad are 3D cells with faces)
00394       // A 3D cell can have both Tri and Quad faces, but because they are affine images of the
00395       // parametrization domain, 3 coefficients are enough to store them in both cases.
00396       else if(subcellDim == 2) {
00397         for(unsigned subcellOrd = 0; subcellOrd < sc; subcellOrd++){
00398           
00399           switch(parentCell.getKey(subcellDim,subcellOrd)){
00400             
00401             // Admissible triangular faces for 3D cells in Shards:
00402             case shards::Triangle<3>::key:
00403             case shards::Triangle<4>::key:
00404             case shards::Triangle<6>::key: 
00405               {
00406                 int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
00407                 int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
00408                 int v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
00409                 const double* v0 = getReferenceVertex(parentCell, v0ord);
00410                 const double* v1 = getReferenceVertex(parentCell, v1ord);
00411                 const double* v2 = getReferenceVertex(parentCell, v2ord);
00412                 
00413                 // x(u,v) = x0 + (x1 - x0)*u + (x2 - x0)*v
00414                 subcellParametrization(subcellOrd, 0, 0) = v0[0];
00415                 subcellParametrization(subcellOrd, 0, 1) = v1[0] - v0[0];
00416                 subcellParametrization(subcellOrd, 0, 2) = v2[0] - v0[0];
00417                   
00418                 // y(u,v) = y0 + (y1 - y0)*u + (y2 - y0)*v
00419                 subcellParametrization(subcellOrd, 1, 0) = v0[1];
00420                 subcellParametrization(subcellOrd, 1, 1) = v1[1] - v0[1];
00421                 subcellParametrization(subcellOrd, 1, 2) = v2[1] - v0[1];
00422                 
00423                 // z(u,v) = z0 + (z1 - z0)*u + (z2 - z0)*v
00424                 subcellParametrization(subcellOrd, 2, 0) = v0[2];
00425                 subcellParametrization(subcellOrd, 2, 1) = v1[2] - v0[2];
00426                 subcellParametrization(subcellOrd, 2, 2) = v2[2] - v0[2];
00427               }
00428               break;
00429               
00430             // Admissible quadrilateral faces for 3D cells in Shards:
00431             case shards::Quadrilateral<4>::key:
00432             case shards::Quadrilateral<8>::key:
00433             case shards::Quadrilateral<9>::key:
00434               {
00435                 int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
00436                 int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
00437                 int v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
00438                 int v3ord = parentCell.getNodeMap(subcellDim, subcellOrd, 3);
00439                 const double* v0 = getReferenceVertex(parentCell, v0ord);
00440                 const double* v1 = getReferenceVertex(parentCell, v1ord);
00441                 const double* v2 = getReferenceVertex(parentCell, v2ord);
00442                 const double* v3 = getReferenceVertex(parentCell, v3ord);
00443                 
00444                 // x(u,v) = (x0+x1+x2+x3)/4+u*(-x0+x1+x2-x3)/4+v*(-x0-x1+x2+x3)/4+uv*(0=x0-x1+x2-x3)/4 
00445                 subcellParametrization(subcellOrd, 0, 0) = ( v0[0] + v1[0] + v2[0] + v3[0])/4.0;
00446                 subcellParametrization(subcellOrd, 0, 1) = (-v0[0] + v1[0] + v2[0] - v3[0])/4.0;
00447                 subcellParametrization(subcellOrd, 0, 2) = (-v0[0] - v1[0] + v2[0] + v3[0])/4.0;
00448                 
00449                 // y(u,v) = (y0+y1+y2+y3)/4+u*(-y0+y1+y2-y3)/4+v*(-y0-y1+y2+y3)/4+uv*(0=y0-y1+y2-y3)/4 
00450                 subcellParametrization(subcellOrd, 1, 0) = ( v0[1] + v1[1] + v2[1] + v3[1])/4.0;
00451                 subcellParametrization(subcellOrd, 1, 1) = (-v0[1] + v1[1] + v2[1] - v3[1])/4.0;
00452                 subcellParametrization(subcellOrd, 1, 2) = (-v0[1] - v1[1] + v2[1] + v3[1])/4.0;
00453                 
00454                 // z(u,v) = (z0+z1+z2+z3)/4+u*(-z0+z1+z2-z3)/4+v*(-z0-z1+z2+z3)/4+uv*(0=z0-z1+z2-z3)/4 
00455                 subcellParametrization(subcellOrd, 2, 0) = ( v0[2] + v1[2] + v2[2] + v3[2])/4.0;
00456                 subcellParametrization(subcellOrd, 2, 1) = (-v0[2] + v1[2] + v2[2] - v3[2])/4.0;
00457                 subcellParametrization(subcellOrd, 2, 2) = (-v0[2] - v1[2] + v2[2] + v3[2])/4.0;
00458               }                
00459               break;
00460             default:
00461               TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00462                                   ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): parametrization not defined for the specified face topology.");
00463              
00464           }// switch face topology key
00465         }// for subcellOrd
00466       }
00467   }
00468   
00469   
00470   
00471   template<class Scalar>
00472   const double* CellTools<Scalar>::getReferenceVertex(const shards::CellTopology& cell,
00473                                                       const int                   vertexOrd){
00474     
00475 #ifdef HAVE_INTREPID_DEBUG
00476     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(cell) ), std::invalid_argument, 
00477                         ">>> ERROR (Intrepid::CellTools::getReferenceVertex): the specified cell topology does not have a reference cell.");
00478     
00479     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= vertexOrd) && (vertexOrd < (int)cell.getVertexCount() ) ), std::invalid_argument,
00480                         ">>> ERROR (Intrepid::CellTools::getReferenceVertex): invalid node ordinal for the specified cell topology. ");
00481 #endif
00482     
00483     // Simply call getReferenceNode with the base topology of the cell
00484     return getReferenceNode(cell.getBaseCellTopologyData(), vertexOrd);
00485   }
00486     
00487   
00488   
00489   template<class Scalar>
00490   template<class ArraySubcellVert>
00491   void CellTools<Scalar>::getReferenceSubcellVertices(ArraySubcellVert &          subcellVertices,
00492                                                       const int                   subcellDim,
00493                                                       const int                   subcellOrd,
00494                                                       const shards::CellTopology& parentCell){
00495 #ifdef HAVE_INTREPID_DEBUG
00496     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
00497                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): the specified cell topology does not have a reference cell.");
00498 
00499     // subcellDim can equal the cell dimension because the cell itself is a valid subcell! In this case
00500     // the method will return all cell cellWorkset.
00501     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellDim) && (subcellDim <= (int)parentCell.getDimension()) ), std::invalid_argument,
00502                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): subcell dimension out of range.");
00503     
00504     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
00505                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): subcell ordinal out of range.");
00506         
00507     // Verify subcellVertices rank and dimensions
00508     {
00509       std::string errmsg = ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices):";
00510       TEUCHOS_TEST_FOR_EXCEPTION( !( requireRankRange(errmsg, subcellVertices, 2, 2) ), std::invalid_argument, errmsg);
00511       
00512       int subcVertexCount = parentCell.getVertexCount(subcellDim, subcellOrd);
00513       int spaceDim = parentCell.getDimension();
00514         
00515       TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellVertices, 0,  subcVertexCount, subcVertexCount) ),
00516                           std::invalid_argument, errmsg);
00517       
00518       TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellVertices, 1,  spaceDim, spaceDim) ),
00519                           std::invalid_argument, errmsg);
00520     }
00521 #endif 
00522     
00523     // Simply call getReferenceNodes with the base topology
00524     getReferenceSubcellNodes(subcellVertices, subcellDim, subcellOrd, parentCell.getBaseCellTopologyData() );
00525   }  
00526 
00527   
00528   
00529   template<class Scalar>
00530   const double* CellTools<Scalar>::getReferenceNode(const shards::CellTopology& cell,
00531                                                     const int                   nodeOrd){
00532     
00533 #ifdef HAVE_INTREPID_DEBUG
00534     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(cell) ), std::invalid_argument, 
00535                         ">>> ERROR (Intrepid::CellTools::getReferenceNode): the specified cell topology does not have a reference cell.");
00536 
00537     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= nodeOrd) && (nodeOrd < (int)cell.getNodeCount() ) ), std::invalid_argument,
00538                         ">>> ERROR (Intrepid::CellTools::getReferenceNode): invalid node ordinal for the specified cell topology. ");
00539 #endif
00540     
00541     // Cartesian coordinates of supported reference cell cellWorkset, padded to three-dimensions.
00542     // Node order follows cell topology definition in Shards
00543     static const double line[2][3] ={
00544       {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0} 
00545     };
00546     static const double line_3[3][3] = {
00547       {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0},     
00548       // Extension node: edge midpoint
00549       { 0.0, 0.0, 0.0}
00550     };
00551     
00552     
00553     // Triangle topologies
00554     static const double triangle[3][3] = {
00555       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0} 
00556     };
00557     static const double triangle_4[4][3] = {
00558       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
00559       // Extension node: cell center
00560       { 1/3, 1/3, 0.0}
00561     };
00562     static const double triangle_6[6][3] = {
00563       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
00564       // Extension cellWorkset: 3 edge midpoints
00565       { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
00566     };
00567     
00568     
00569     // Quadrilateral topologies
00570     static const double quadrilateral[4][3] = {
00571       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}
00572     };
00573     static const double quadrilateral_8[8][3] = {
00574       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
00575       // Extension cellWorkset: 4 edge midpoints
00576       { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}
00577     };
00578     static const double quadrilateral_9[9][3] = {
00579       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
00580       // Extension cellWorkset: 4 edge midpoints + 1 cell center
00581       { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}, { 0.0, 0.0, 0.0}
00582     };
00583     
00584     
00585     // Tetrahedron topologies
00586     static const double tetrahedron[4][3] = {
00587       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
00588     };
00589     static const double tetrahedron_8[8][3] = {
00590       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
00591       // Extension cellWorkset: 4 face centers (do not follow natural face order - see the cell topology!)
00592       { 1/3, 0.0, 1/3}, { 1/3, 1/3, 1/3}, { 1/3, 1/3, 0.0}, { 0.0, 1/3, 1/3} 
00593     };
00594     static const double tetrahedron_10[10][3] = {
00595       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
00596       // Extension cellWorkset: 6 edge midpoints
00597       { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
00598     };
00599 
00600     static const double tetrahedron_11[10][3] = {
00601       { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
00602       // Extension cellWorkset: 6 edge midpoints
00603       { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
00604     };
00605 
00606     
00607     // Hexahedron topologies
00608     static const double hexahedron[8][3] = {
00609       {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
00610       {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0}
00611     };
00612     static const double hexahedron_20[20][3] = {
00613       {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
00614       {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
00615       // Extension cellWorkset: 12 edge midpoints (do not follow natural edge order - see cell topology!)
00616       { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0}, 
00617       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
00618       { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0}
00619     };
00620     static const double hexahedron_27[27][3] = {
00621       {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
00622       {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
00623       // Extension cellWorkset: 12 edge midpoints + 1 cell center + 6 face centers  (do not follow natural subcell order!)
00624       { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0}, 
00625       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
00626       { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0},
00627       { 0.0, 0.0, 0.0},
00628       { 0.0, 0.0,-1.0}, { 0.0, 0.0, 1.0}, {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, {0.0,-1.0, 0.0}, {0.0, 1.0, 0.0} 
00629     };
00630     
00631     
00632     // Pyramid topologies
00633     static const double pyramid[5][3] = {
00634       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
00635     };
00636     static const double pyramid_13[13][3] = {
00637       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
00638       // Extension cellWorkset: 8 edge midpoints 
00639       { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
00640       {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}   
00641     };
00642     static const double pyramid_14[14][3] = {
00643       {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
00644       // Extension cellWorkset: 8 edge midpoints + quadrilateral face midpoint 
00645       { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
00646       {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}, { 0.0, 0.0, 0.0}  
00647     };
00648     
00649     
00650     // Wedge topologies
00651     static const double wedge[6][3] = {
00652       { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0} 
00653     };
00654     static const double wedge_15[15][3] = {
00655       { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
00656       // Extension cellWorkset: 9 edge midpoints (do not follow natural edge order - see cell topology!)
00657       { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
00658       { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0}
00659     };
00660     static const double wedge_18[18][3] = {
00661       { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
00662       // Extension cellWorkset: 9 edge midpoints + 3 quad face centers (do not follow natural subcell order - see cell topology!)
00663       { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
00664       { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0},
00665       { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
00666     };
00667     
00668     
00669     switch(cell.getKey() ) {
00670       
00671       // Base line topologies
00672       case shards::Line<2>::key:
00673       case shards::ShellLine<2>::key:
00674       case shards::Beam<2>::key:
00675         return line[nodeOrd];
00676         break;
00677         
00678       // Extended line topologies
00679       case shards::Line<3>::key:
00680       case shards::ShellLine<3>::key:
00681       case shards::Beam<3>::key:
00682         return line_3[nodeOrd];
00683         break;
00684         
00685         
00686       // Base triangle topologies
00687       case shards::Triangle<3>::key:
00688       case shards::ShellTriangle<3>::key:
00689         return triangle[nodeOrd];
00690         break;
00691         
00692       // Extened Triangle topologies
00693       case shards::Triangle<4>::key:
00694         return triangle_4[nodeOrd];
00695         break;
00696       case shards::Triangle<6>::key:
00697       case shards::ShellTriangle<6>::key:
00698         return triangle_6[nodeOrd];
00699         break;
00700         
00701         
00702       // Base Quadrilateral topologies  
00703       case shards::Quadrilateral<4>::key:
00704       case shards::ShellQuadrilateral<4>::key:
00705         return quadrilateral[nodeOrd];
00706         break;
00707         
00708       // Extended Quadrilateral topologies
00709       case shards::Quadrilateral<8>::key:
00710       case shards::ShellQuadrilateral<8>::key:
00711         return quadrilateral_8[nodeOrd];
00712         break;
00713       case shards::Quadrilateral<9>::key:
00714       case shards::ShellQuadrilateral<9>::key:
00715         return quadrilateral_9[nodeOrd];
00716         break;
00717         
00718         
00719       // Base Tetrahedron topology
00720       case shards::Tetrahedron<4>::key:
00721         return tetrahedron[nodeOrd];
00722         break;
00723         
00724       // Extended Tetrahedron topologies
00725       case shards::Tetrahedron<8>::key:
00726         return tetrahedron_8[nodeOrd];
00727         break;
00728       case shards::Tetrahedron<10>::key:
00729         return tetrahedron_10[nodeOrd];
00730         break;
00731       case shards::Tetrahedron<11>::key:
00732         return tetrahedron_11[nodeOrd];
00733         break;
00734 
00735         
00736       // Base Hexahedron topology
00737       case shards::Hexahedron<8>::key:
00738         return hexahedron[nodeOrd];
00739         break;
00740         
00741       // Extended Hexahedron topologies
00742       case shards::Hexahedron<20>::key:
00743         return hexahedron_20[nodeOrd];
00744         break;
00745       case shards::Hexahedron<27>::key:
00746         return hexahedron_27[nodeOrd];
00747         break;
00748 
00749         
00750       // Base Pyramid topology  
00751       case shards::Pyramid<5>::key:
00752         return pyramid[nodeOrd];
00753         break;
00754         
00755       // Extended pyramid topologies
00756       case shards::Pyramid<13>::key:
00757         return pyramid_13[nodeOrd];
00758         break;
00759      case shards::Pyramid<14>::key:
00760         return pyramid_14[nodeOrd];
00761         break;
00762       
00763         
00764       // Base Wedge topology
00765       case shards::Wedge<6>::key:
00766         return wedge[nodeOrd];
00767         break;
00768         
00769       // Extended Wedge topologies
00770       case shards::Wedge<15>::key:
00771         return wedge_15[nodeOrd];
00772         break;
00773       case shards::Wedge<18>::key:
00774         return wedge_18[nodeOrd];
00775         break;
00776         
00777       default:
00778         TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
00779                             ">>> ERROR (Intrepid::CellTools::getReferenceNode): invalid cell topology.");
00780     }
00781     // To disable compiler warning, should never be reached
00782     return line[0];
00783   }
00784   
00785   
00786   
00787   template<class Scalar>
00788   template<class ArraySubcellNode>
00789   void CellTools<Scalar>::getReferenceSubcellNodes(ArraySubcellNode &          subcellNodes,
00790                                                    const int                   subcellDim,
00791                                                    const int                   subcellOrd,
00792                                                    const shards::CellTopology& parentCell){
00793 #ifdef HAVE_INTREPID_DEBUG
00794     TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
00795                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): the specified cell topology does not have a reference cell.");
00796     
00797     // subcellDim can equal the cell dimension because the cell itself is a valid subcell! In this case
00798     // the method will return all cell cellWorkset.
00799     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellDim) && (subcellDim <= (int)parentCell.getDimension()) ), std::invalid_argument,
00800                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): subcell dimension out of range.");
00801     
00802     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
00803                         ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): subcell ordinal out of range.");
00804     
00805     // Verify subcellNodes rank and dimensions
00806     {
00807       std::string errmsg = ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes):";
00808       TEUCHOS_TEST_FOR_EXCEPTION( !( requireRankRange(errmsg, subcellNodes, 2, 2) ), std::invalid_argument, errmsg);
00809       
00810       int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
00811       int spaceDim = parentCell.getDimension();
00812       
00813       TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellNodes, 0,  subcNodeCount, subcNodeCount) ),
00814                           std::invalid_argument, errmsg);
00815       
00816       TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellNodes, 1,  spaceDim, spaceDim) ),
00817                           std::invalid_argument, errmsg);
00818     }
00819 #endif 
00820     
00821     // Find how many cellWorkset does the specified subcell have.
00822     int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
00823     
00824     // Loop over subcell cellWorkset
00825     for(int subcNodeOrd = 0; subcNodeOrd < subcNodeCount; subcNodeOrd++){
00826       
00827       // Get the node number relative to the parent reference cell
00828       int cellNodeOrd = parentCell.getNodeMap(subcellDim, subcellOrd, subcNodeOrd);
00829             
00830       // Loop over node's Cartesian coordinates
00831       for(int dim = 0; dim < (int)parentCell.getDimension(); dim++){
00832         subcellNodes(subcNodeOrd, dim) = CellTools::getReferenceNode(parentCell, cellNodeOrd)[dim];
00833       }
00834     }
00835   }  
00836   
00837   
00838   
00839   template<class Scalar>
00840   int CellTools<Scalar>::hasReferenceCell(const shards::CellTopology& cell) {
00841     
00842     switch(cell.getKey() ) {
00843       case shards::Line<2>::key:
00844       case shards::Line<3>::key:
00845       case shards::ShellLine<2>::key:
00846       case shards::ShellLine<3>::key:
00847       case shards::Beam<2>::key:
00848       case shards::Beam<3>::key:
00849         
00850       case shards::Triangle<3>::key:
00851       case shards::Triangle<4>::key:
00852       case shards::Triangle<6>::key:
00853       case shards::ShellTriangle<3>::key:
00854       case shards::ShellTriangle<6>::key:
00855         
00856       case shards::Quadrilateral<4>::key:
00857       case shards::Quadrilateral<8>::key:
00858       case shards::Quadrilateral<9>::key:
00859       case shards::ShellQuadrilateral<4>::key:
00860       case shards::ShellQuadrilateral<8>::key:
00861       case shards::ShellQuadrilateral<9>::key:
00862         
00863       case shards::Tetrahedron<4>::key:
00864       case shards::Tetrahedron<8>::key:
00865       case shards::Tetrahedron<10>::key:
00866       case shards::Tetrahedron<11>::key:
00867         
00868       case shards::Hexahedron<8>::key:
00869       case shards::Hexahedron<20>::key:
00870       case shards::Hexahedron<27>::key:
00871         
00872       case shards::Pyramid<5>::key:
00873       case shards::Pyramid<13>::key:
00874       case shards::Pyramid<14>::key:
00875         
00876       case shards::Wedge<6>::key:
00877       case shards::Wedge<15>::key:
00878       case shards::Wedge<18>::key:
00879         return 1;
00880         break;
00881         
00882       default:
00883         return 0;
00884     }
00885     return 0;
00886   }
00887   
00888   //============================================================================================//
00889   //                                                                                            //
00890   //                     Jacobian, inverse Jacobian and Jacobian determinant                    //
00891   //                                                                                            //
00892   //============================================================================================//  
00893   
00894   template<class Scalar>
00895   template<class ArrayJac, class ArrayPoint, class ArrayCell>
00896   void CellTools<Scalar>::setJacobian(ArrayJac &                   jacobian,
00897                                       const ArrayPoint &           points,
00898                                       const ArrayCell  &           cellWorkset,
00899                                       const shards::CellTopology & cellTopo,
00900                                       const int &                  whichCell) 
00901   {
00902     INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell,  cellTopo) );
00903     
00904     int spaceDim  = (int)cellTopo.getDimension();
00905     int numCells  = cellWorkset.dimension(0);
00906     //points can be rank-2 (P,D), or rank-3 (C,P,D)
00907     int numPoints = (points.rank() == 2) ? points.dimension(0) : points.dimension(1);
00908     
00909     // Jacobian is computed using gradients of an appropriate H(grad) basis function: define RCP to the base class
00910     Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
00911     
00912     // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
00913     switch( cellTopo.getKey() ){
00914       
00915       // Standard Base topologies (number of cellWorkset = number of vertices)
00916       case shards::Line<2>::key:
00917         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00918         break;
00919         
00920       case shards::Triangle<3>::key:
00921         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00922         break;
00923         
00924       case shards::Quadrilateral<4>::key:
00925         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00926         break;
00927         
00928       case shards::Tetrahedron<4>::key:
00929         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00930         break;
00931         
00932       case shards::Hexahedron<8>::key:
00933         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00934         break;
00935         
00936       case shards::Wedge<6>::key:
00937         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00938         break;
00939 
00940       case shards::Pyramid<5>::key:
00941             HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
00942             break;
00943         
00944       // Standard Extended topologies
00945       case shards::Triangle<6>::key:    
00946         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00947         break;
00948       case shards::Quadrilateral<9>::key:
00949         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00950         break;
00951         
00952       case shards::Tetrahedron<10>::key:
00953         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00954         break;
00955 
00956       case shards::Tetrahedron<11>::key:
00957         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
00958         break;
00959 
00960       case shards::Hexahedron<20>::key:
00961         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
00962         break;
00963         
00964       case shards::Hexahedron<27>::key:
00965         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00966         break;
00967 
00968       case shards::Wedge<15>::key:
00969         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
00970         break;
00971         
00972       case shards::Wedge<18>::key:
00973         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00974         break;
00975 
00976       case shards::Pyramid<13>::key:
00977         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
00978         break;
00979         
00980         // These extended topologies are not used for mapping purposes
00981       case shards::Quadrilateral<8>::key:
00982         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
00983                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
00984         break;
00985         
00986         // Base and Extended Line, Beam and Shell topologies  
00987       case shards::Line<3>::key:
00988       case shards::Beam<2>::key:
00989       case shards::Beam<3>::key:
00990       case shards::ShellLine<2>::key:
00991       case shards::ShellLine<3>::key:
00992       case shards::ShellTriangle<3>::key:
00993       case shards::ShellTriangle<6>::key:
00994       case shards::ShellQuadrilateral<4>::key:
00995       case shards::ShellQuadrilateral<8>::key:
00996       case shards::ShellQuadrilateral<9>::key:
00997         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
00998                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
00999         break;
01000       default:
01001         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01002                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");        
01003     }// switch  
01004     
01005     // Temp (F,P,D) array for the values of basis functions gradients at the reference points
01006     int basisCardinality = HGRAD_Basis -> getCardinality();
01007     FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
01008     
01009     // Initialize jacobian
01010     for(int i = 0; i < jacobian.size(); i++){
01011       jacobian[i] = 0.0;
01012     }
01013         
01014     // Handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of points arrays.
01015     switch(points.rank()) {
01016       
01017       // refPoints is (P,D): a single or multiple cell jacobians computed for a single set of ref. points
01018       case 2:
01019         {
01020           // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
01021           FieldContainer<Scalar> tempPoints( points.dimension(0), points.dimension(1) );
01022           // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01023           for(int pt = 0; pt < points.dimension(0); pt++){
01024             for(int dm = 0; dm < points.dimension(1) ; dm++){
01025               tempPoints(pt, dm) = points(pt, dm);
01026             }//dm
01027           }//pt
01028           HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
01029           
01030           // The outer loops select the multi-index of the Jacobian entry: cell, point, row, col
01031           // If whichCell = -1, all jacobians are computed, otherwise a single cell jacobian is computed
01032           int cellLoop = (whichCell == -1) ? numCells : 1 ;
01033           
01034           if(whichCell == -1) {
01035             for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01036               for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01037                 for(int row = 0; row < spaceDim; row++){
01038                   for(int col = 0; col < spaceDim; col++){
01039                     
01040                     // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
01041                     for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01042                       jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01043                     } // bfOrd
01044                   } // col
01045                 } // row
01046               } // pointOrd
01047             } // cellOrd
01048           }
01049           else {
01050             for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01051               for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01052                 for(int row = 0; row < spaceDim; row++){
01053                   for(int col = 0; col < spaceDim; col++){
01054                   
01055                     // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
01056                     for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01057                       jacobian(pointOrd, row, col) += cellWorkset(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01058                     } // bfOrd
01059                   } // col
01060                 } // row
01061               } // pointOrd
01062             } // cellOrd
01063           } // if whichcell
01064         }// case 2
01065         break;
01066         
01067         // points is (C,P,D): multiple jacobians computed at multiple point sets, one jacobian per cell  
01068       case 3:
01069         {
01070           // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
01071           FieldContainer<Scalar> tempPoints( points.dimension(1), points.dimension(2) );
01072           
01073           for(int cellOrd = 0; cellOrd < numCells; cellOrd++) {
01074             
01075             // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01076             for(int pt = 0; pt < points.dimension(1); pt++){
01077               for(int dm = 0; dm < points.dimension(2) ; dm++){
01078                 tempPoints(pt, dm) = points(cellOrd, pt, dm);
01079               }//dm
01080             }//pt
01081             
01082             // Compute gradients of basis functions at this set of ref. points
01083             HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
01084             
01085             // Compute jacobians for the point set corresponding to the current cellordinal
01086             for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01087               for(int row = 0; row < spaceDim; row++){
01088                 for(int col = 0; col < spaceDim; col++){
01089                   
01090                   // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same
01091                   for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01092                     jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01093                   } // bfOrd
01094                 } // col
01095               } // row
01096             } // pointOrd
01097           }//cellOrd
01098         }// case 3
01099         
01100         break;
01101         
01102       default:
01103         TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 2) && (points.rank() == 3) ), std::invalid_argument,
01104                             ">>> ERROR (Intrepid::CellTools::setJacobian): rank 2 or 3 required for points array. ");        
01105     }//switch
01106   }
01107   
01108 
01109 
01110 template<class Scalar>
01111 template<class ArrayJacInv, class ArrayJac>
01112 void CellTools<Scalar>::setJacobianInv(ArrayJacInv &     jacobianInv,
01113                                        const ArrayJac &  jacobian) 
01114 {
01115   INTREPID_VALIDATE( validateArguments_setJacobianInv(jacobianInv, jacobian) );
01116 
01117   RealSpaceTools<Scalar>::inverse(jacobianInv, jacobian);
01118 }
01119 
01120 
01121 
01122 template<class Scalar>
01123 template<class ArrayJacDet, class ArrayJac>
01124 void CellTools<Scalar>::setJacobianDet(ArrayJacDet &     jacobianDet,
01125                                        const ArrayJac &  jacobian)
01126 {
01127   INTREPID_VALIDATE( validateArguments_setJacobianDetArgs(jacobianDet, jacobian) );
01128 
01129   RealSpaceTools<Scalar>::det(jacobianDet, jacobian);
01130 }
01131 
01132 //============================================================================================//
01133 //                                                                                            //
01134 //                      Reference-to-physical frame mapping and its inverse                   //
01135 //                                                                                            //
01136 //============================================================================================//
01137 
01138 template<class Scalar>
01139 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
01140 void CellTools<Scalar>::mapToPhysicalFrame(ArrayPhysPoint      &        physPoints,
01141                                            const ArrayRefPoint &        refPoints,
01142                                            const ArrayCell     &        cellWorkset,
01143                                            const shards::CellTopology & cellTopo,
01144                                            const int &                  whichCell)
01145 {
01146   INTREPID_VALIDATE(validateArguments_mapToPhysicalFrame( physPoints, refPoints, cellWorkset, cellTopo, whichCell) );
01147   
01148   int spaceDim  = (int)cellTopo.getDimension();
01149   int numCells  = cellWorkset.dimension(0);
01150   //points can be rank-2 (P,D), or rank-3 (C,P,D)
01151   int numPoints = (refPoints.rank() == 2) ? refPoints.dimension(0) : refPoints.dimension(1);
01152     
01153   // Mapping is computed using an appropriate H(grad) basis function: define RCP to the base class
01154   Teuchos::RCP<Basis<Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
01155   
01156   // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
01157   switch( cellTopo.getKey() ){
01158     
01159     // Standard Base topologies (number of cellWorkset = number of vertices)
01160     case shards::Line<2>::key:
01161       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01162       break;
01163       
01164     case shards::Triangle<3>::key:
01165       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01166       break;
01167       
01168     case shards::Quadrilateral<4>::key:
01169       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01170       break;
01171       
01172     case shards::Tetrahedron<4>::key:
01173       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01174       break;
01175       
01176     case shards::Hexahedron<8>::key:
01177       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01178       break;
01179       
01180     case shards::Wedge<6>::key:
01181       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01182       break;
01183       
01184     case shards::Pyramid<5>::key:
01185           HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01186           break;
01187 
01188     // Standard Extended topologies
01189     case shards::Triangle<6>::key:    
01190       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01191       break;
01192       
01193     case shards::Quadrilateral<9>::key:
01194       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01195       break;
01196       
01197     case shards::Tetrahedron<10>::key:
01198       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01199       break;
01200 
01201     case shards::Tetrahedron<11>::key:
01202       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
01203       break;
01204 
01205     case shards::Hexahedron<20>::key:
01206       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
01207       break;
01208 
01209     case shards::Hexahedron<27>::key:
01210       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01211       break;
01212 
01213     case shards::Wedge<15>::key:
01214       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
01215       break;
01216       
01217     case shards::Wedge<18>::key:
01218       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01219       break;
01220 
01221     case shards::Pyramid<13>::key:
01222       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
01223       break;
01224       
01225     // These extended topologies are not used for mapping purposes
01226     case shards::Quadrilateral<8>::key:
01227       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01228                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
01229       break;
01230       
01231     // Base and Extended Line, Beam and Shell topologies  
01232     case shards::Line<3>::key:
01233     case shards::Beam<2>::key:
01234     case shards::Beam<3>::key:
01235     case shards::ShellLine<2>::key:
01236     case shards::ShellLine<3>::key:
01237     case shards::ShellTriangle<3>::key:
01238     case shards::ShellTriangle<6>::key:
01239     case shards::ShellQuadrilateral<4>::key:
01240     case shards::ShellQuadrilateral<8>::key:
01241     case shards::ShellQuadrilateral<9>::key:
01242       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01243                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
01244       break;
01245     default:
01246       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01247                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported.");        
01248   }// switch  
01249 
01250   // Temp (F,P) array for the values of nodal basis functions at the reference points
01251   int basisCardinality = HGRAD_Basis -> getCardinality();
01252   FieldContainer<Scalar> basisVals(basisCardinality, numPoints);
01253   
01254   // Initialize physPoints
01255   for(int i = 0; i < physPoints.size(); i++){
01256     physPoints[i] = 0.0;
01257   }
01258   
01259   // handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of refPoints
01260   switch(refPoints.rank()) {
01261     
01262     // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
01263     case 2:
01264       {
01265         // getValues requires rank-2 (P,D) input array, but refPoints cannot be passed directly as argument because they are a user type
01266         FieldContainer<Scalar> tempPoints( refPoints.dimension(0), refPoints.dimension(1) );
01267         // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01268         for(int pt = 0; pt < refPoints.dimension(0); pt++){
01269           for(int dm = 0; dm < refPoints.dimension(1) ; dm++){
01270             tempPoints(pt, dm) = refPoints(pt, dm);
01271           }//dm
01272         }//pt
01273         HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
01274         
01275         // If whichCell = -1, ref pt. set is mapped to all cells, otherwise, the set is mapped to one cell only
01276         int cellLoop = (whichCell == -1) ? numCells : 1 ;
01277         
01278         // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
01279         for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01280           for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01281             for(int dim = 0; dim < spaceDim; dim++){
01282               for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01283                 
01284                 if(whichCell == -1){
01285                   physPoints(cellOrd, pointOrd, dim) += cellWorkset(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01286                 }
01287                 else{
01288                   physPoints(pointOrd, dim) += cellWorkset(whichCell, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01289                 }
01290               } // bfOrd
01291             }// dim
01292           }// pointOrd
01293         }//cellOrd
01294       }// case 2
01295       break;
01296       
01297     // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.  
01298     case 3:
01299       {
01300         // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
01301         FieldContainer<Scalar> tempPoints( refPoints.dimension(1), refPoints.dimension(2) );
01302         
01303         // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
01304         for(int cellOrd = 0; cellOrd < numCells; cellOrd++) {
01305           
01306           // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01307           for(int pt = 0; pt < refPoints.dimension(1); pt++){
01308             for(int dm = 0; dm < refPoints.dimension(2) ; dm++){
01309               tempPoints(pt, dm) = refPoints(cellOrd, pt, dm);
01310             }//dm
01311           }//pt
01312           
01313           // Compute basis values for this set of ref. points
01314           HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
01315           
01316           for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01317             for(int dim = 0; dim < spaceDim; dim++){
01318               for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01319                 
01320                 physPoints(cellOrd, pointOrd, dim) += cellWorkset(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01321                 
01322               } // bfOrd
01323             }// dim
01324           }// pointOrd
01325         }//cellOrd        
01326       }// case 3
01327       break;
01328       
01329     default:
01330       TEUCHOS_TEST_FOR_EXCEPTION( !( (refPoints.rank() == 2) && (refPoints.rank() == 3) ), std::invalid_argument,
01331                              ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): rank 2 or 3 required for refPoints array. ");
01332   }
01333 }
01334 
01335 
01336 
01337 template<class Scalar>
01338 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
01339 void CellTools<Scalar>::mapToReferenceFrame(ArrayRefPoint        &        refPoints,
01340                                             const ArrayPhysPoint &        physPoints,
01341                                             const ArrayCell      &        cellWorkset,
01342                                             const shards::CellTopology &  cellTopo,
01343                                             const int &                   whichCell)
01344 {
01345   INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell) );
01346   
01347   int spaceDim  = (int)cellTopo.getDimension();
01348   int numPoints;
01349   int numCells;
01350   
01351   // Define initial guesses to be  the Cell centers of the reference cell topology
01352   FieldContainer<Scalar> cellCenter(spaceDim);
01353   switch( cellTopo.getKey() ){
01354     // Standard Base topologies (number of cellWorkset = number of vertices)
01355     case shards::Line<2>::key:
01356       cellCenter(0) = 0.0;    break;
01357       
01358     case shards::Triangle<3>::key:
01359     case shards::Triangle<6>::key:    
01360       cellCenter(0) = 1./3.;    cellCenter(1) = 1./3.;  break;
01361       
01362     case shards::Quadrilateral<4>::key:
01363     case shards::Quadrilateral<9>::key:
01364       cellCenter(0) = 0.0;      cellCenter(1) = 0.0;    break;
01365       
01366     case shards::Tetrahedron<4>::key:
01367     case shards::Tetrahedron<10>::key:
01368     case shards::Tetrahedron<11>::key:
01369       cellCenter(0) = 1./6.;    cellCenter(1) =  1./6.;    cellCenter(2) =  1./6.;  break;
01370       
01371     case shards::Hexahedron<8>::key:
01372     case shards::Hexahedron<20>::key:
01373     case shards::Hexahedron<27>::key:
01374       cellCenter(0) = 0.0;      cellCenter(1) =  0.0;       cellCenter(2) =  0.0;   break;
01375       
01376     case shards::Wedge<6>::key:
01377     case shards::Wedge<15>::key:
01378     case shards::Wedge<18>::key:
01379       cellCenter(0) = 1./3.;    cellCenter(1) =  1./3.;     cellCenter(2) = 0.0;    break;
01380 
01381     case shards::Pyramid<5>::key:
01382     case shards::Pyramid<13>::key:
01383       cellCenter(0) = 0.;       cellCenter(1) = 0.;         cellCenter(2) = 0.25;    break;
01384       
01385       // These extended topologies are not used for mapping purposes
01386     case shards::Quadrilateral<8>::key:
01387       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01388                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
01389       break;
01390       
01391       // Base and Extended Line, Beam and Shell topologies  
01392     case shards::Line<3>::key:
01393     case shards::Beam<2>::key:
01394     case shards::Beam<3>::key:
01395     case shards::ShellLine<2>::key:
01396     case shards::ShellLine<3>::key:
01397     case shards::ShellTriangle<3>::key:
01398     case shards::ShellTriangle<6>::key:
01399     case shards::ShellQuadrilateral<4>::key:
01400     case shards::ShellQuadrilateral<8>::key:
01401     case shards::ShellQuadrilateral<9>::key:
01402       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01403                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
01404       break;
01405     default:
01406       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01407                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported.");        
01408   }// switch key 
01409   
01410   // Resize initial guess depending on the rank of the physical points array
01411   FieldContainer<Scalar> initGuess;
01412   
01413   // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) initial guess.
01414   if(whichCell == -1){
01415     numPoints = physPoints.dimension(1);
01416     numCells = cellWorkset.dimension(0);
01417     initGuess.resize(numCells, numPoints, spaceDim);
01418     // Set initial guess:
01419     for(int c = 0; c < numCells; c++){
01420       for(int p = 0; p < numPoints; p++){
01421         for(int d = 0; d < spaceDim; d++){
01422           initGuess(c, p, d) = cellCenter(d);
01423         }// d
01424       }// p
01425     }// c
01426   }
01427   // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) initial guess.
01428   else {
01429     numPoints = physPoints.dimension(0);
01430     initGuess.resize(numPoints, spaceDim);
01431     // Set initial guess:
01432     for(int p = 0; p < numPoints; p++){
01433       for(int d = 0; d < spaceDim; d++){
01434         initGuess(p, d) = cellCenter(d);
01435       }// d
01436     }// p
01437   }
01438   
01439   // Call method with initial guess
01440   mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell);  
01441 }
01442   
01443   
01444 
01445 template<class Scalar>
01446 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
01447 void CellTools<Scalar>::mapToReferenceFrameInitGuess(ArrayRefPoint        &        refPoints,
01448                                                      const ArrayInitGuess &        initGuess,
01449                                                      const ArrayPhysPoint &        physPoints,
01450                                                      const ArrayCell      &        cellWorkset,
01451                                                      const shards::CellTopology &  cellTopo,
01452                                                      const int &                   whichCell)
01453 {
01454   INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell) );
01455 
01456   int spaceDim  = (int)cellTopo.getDimension();
01457   int numPoints;
01458   int numCells=0;
01459   
01460   // Temp arrays for Newton iterates and Jacobians. Resize according to rank of ref. point array
01461   FieldContainer<Scalar> xOld;
01462   FieldContainer<Scalar> xTem;  
01463   FieldContainer<Scalar> jacobian;
01464   FieldContainer<Scalar> jacobInv;
01465   FieldContainer<Scalar> error; 
01466   FieldContainer<Scalar> cellCenter(spaceDim);
01467   
01468   // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) temp arrays and (C,P,D,D) Jacobians.
01469   if(whichCell == -1){
01470     numPoints = physPoints.dimension(1);
01471     numCells = cellWorkset.dimension(0);
01472     xOld.resize(numCells, numPoints, spaceDim);
01473     xTem.resize(numCells, numPoints, spaceDim);  
01474     jacobian.resize(numCells,numPoints, spaceDim, spaceDim);
01475     jacobInv.resize(numCells,numPoints, spaceDim, spaceDim);
01476     error.resize(numCells,numPoints); 
01477     // Set initial guess to xOld
01478     for(int c = 0; c < numCells; c++){
01479       for(int p = 0; p < numPoints; p++){
01480         for(int d = 0; d < spaceDim; d++){
01481           xOld(c, p, d) = initGuess(c, p, d);
01482         }// d
01483       }// p
01484     }// c
01485   }
01486   // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) temp arrays and (P,D,D) Jacobians.
01487   else {
01488     numPoints = physPoints.dimension(0);
01489     xOld.resize(numPoints, spaceDim);
01490     xTem.resize(numPoints, spaceDim);  
01491     jacobian.resize(numPoints, spaceDim, spaceDim);
01492     jacobInv.resize(numPoints, spaceDim, spaceDim);
01493     error.resize(numPoints); 
01494     // Set initial guess to xOld
01495     for(int p = 0; p < numPoints; p++){
01496       for(int d = 0; d < spaceDim; d++){
01497         xOld(p, d) = initGuess(p, d);
01498       }// d
01499     }// p
01500   }
01501   
01502   
01503   // Newton method to solve the equation F(refPoints) - physPoints = 0:
01504   // refPoints = xOld - DF^{-1}(xOld)*(F(xOld) - physPoints) = xOld + DF^{-1}(xOld)*(physPoints - F(xOld))
01505   for(int iter = 0; iter < INTREPID_MAX_NEWTON; ++iter) {
01506     
01507     // Jacobians at the old iterates and their inverses. 
01508     setJacobian(jacobian, xOld, cellWorkset, cellTopo, whichCell);
01509     setJacobianInv(jacobInv, jacobian);
01510         
01511     // The Newton step.
01512     mapToPhysicalFrame( xTem, xOld, cellWorkset, cellTopo, whichCell );      // xTem <- F(xOld)
01513     RealSpaceTools<Scalar>::subtract( xTem, physPoints, xTem );        // xTem <- physPoints - F(xOld)
01514     RealSpaceTools<Scalar>::matvec( refPoints, jacobInv, xTem);        // refPoints <- DF^{-1}( physPoints - F(xOld) )
01515     RealSpaceTools<Scalar>::add( refPoints, xOld );                    // refPoints <- DF^{-1}( physPoints - F(xOld) ) + xOld
01516     
01517     // l2 error (Euclidean distance) between old and new iterates: |xOld - xNew|
01518     RealSpaceTools<Scalar>::subtract( xTem, xOld, refPoints );
01519     RealSpaceTools<Scalar>::vectorNorm( error, xTem, NORM_TWO );
01520     
01521     // Average L2 error for a multiple sets of physical points: error is rank-2 (C,P) array 
01522     double totalError;
01523     if(whichCell == -1) {
01524       FieldContainer<Scalar> cellWiseError(numCells);
01525       // error(C,P) -> cellWiseError(P)
01526       RealSpaceTools<Scalar>::vectorNorm( cellWiseError, error, NORM_ONE );
01527       totalError = RealSpaceTools<Scalar>::vectorNorm( cellWiseError, NORM_ONE );
01528     }
01529     //Average L2 error for a single set of physical points: error is rank-1 (P) array
01530     else{
01531       totalError = RealSpaceTools<Scalar>::vectorNorm( error, NORM_ONE );   
01532       totalError = totalError;
01533     }
01534     
01535     // Stopping criterion:
01536     if (totalError < INTREPID_TOL) {
01537       break;
01538     } 
01539     else if ( iter > INTREPID_MAX_NEWTON) {
01540       INTREPID_VALIDATE(std::cout << " Intrepid::CellTools::mapToReferenceFrameInitGuess failed to converge to desired tolerance within " 
01541                       << INTREPID_MAX_NEWTON  << " iterations\n" );
01542       break;
01543     }
01544     
01545     // initialize next Newton step
01546     xOld = refPoints;
01547   } // for(iter)
01548 }
01549 
01550 
01551 
01552 template<class Scalar>
01553 template<class ArraySubcellPoint, class ArrayParamPoint>
01554 void CellTools<Scalar>::mapToReferenceSubcell(ArraySubcellPoint     &       refSubcellPoints,
01555                                               const ArrayParamPoint &       paramPoints,
01556                                               const int                     subcellDim,
01557                                               const int                     subcellOrd,
01558                                               const shards::CellTopology &  parentCell){
01559   
01560   int cellDim = parentCell.getDimension();
01561   int numPts  = paramPoints.dimension(0);
01562 
01563 #ifdef HAVE_INTREPID_DEBUG
01564   TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
01565                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
01566   
01567   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
01568                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-dimensional subcells.");  
01569   
01570   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
01571                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
01572   
01573   // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
01574   std::string errmsg = ">>> ERROR (Intrepid::mapToReferenceSubcell):";
01575   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, refSubcellPoints, 2,2), std::invalid_argument, errmsg);
01576   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, refSubcellPoints, 1, cellDim, cellDim), std::invalid_argument, errmsg);
01577                     
01578   // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
01579   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, paramPoints, 2,2), std::invalid_argument, errmsg);
01580   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, paramPoints, 1, subcellDim, subcellDim), std::invalid_argument, errmsg);    
01581   
01582   // cross check: refSubcellPoints and paramPoints: dimension 0 must match
01583   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refSubcellPoints, 0,  paramPoints, 0), std::invalid_argument, errmsg);      
01584 #endif
01585   
01586   
01587   // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
01588   const FieldContainer<double>& subcellMap = getSubcellParametrization(subcellDim, parentCell);
01589 
01590   // Apply the parametrization map to every point in parameter domain
01591   if(subcellDim == 2) {
01592     for(int pt = 0; pt < numPts; pt++){
01593       double u = paramPoints(pt,0);
01594       double v = paramPoints(pt,1);
01595       
01596       // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
01597       for(int  dim = 0; dim < cellDim; dim++){
01598         refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + \
01599                                     subcellMap(subcellOrd, dim, 1)*u + \
01600                                     subcellMap(subcellOrd, dim, 2)*v;
01601       }
01602     }
01603   }
01604   else if(subcellDim == 1) {    
01605     for(int pt = 0; pt < numPts; pt++){
01606       for(int dim = 0; dim < cellDim; dim++) {
01607         refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + subcellMap(subcellOrd, dim, 1)*paramPoints(pt,0);
01608       }
01609     }
01610   }
01611   else{
01612     TEUCHOS_TEST_FOR_EXCEPTION( !( (subcellDim == 1) || (subcellDim == 2) ), std::invalid_argument, 
01613                         ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
01614   }
01615 }
01616 
01617 
01618 
01619 template<class Scalar>
01620 template<class ArrayEdgeTangent>
01621 void CellTools<Scalar>::getReferenceEdgeTangent(ArrayEdgeTangent &            refEdgeTangent,
01622                                                 const int &                   edgeOrd,
01623                                                 const shards::CellTopology &  parentCell){
01624   
01625   int spaceDim  = parentCell.getDimension();
01626   
01627 #ifdef HAVE_INTREPID_DEBUG
01628   
01629   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01630                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): two or three-dimensional parent cell required");
01631   
01632   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= edgeOrd) && (edgeOrd < (int)parentCell.getSubcellCount(1) ) ), std::invalid_argument,
01633                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): edge ordinal out of bounds");  
01634   
01635   TEUCHOS_TEST_FOR_EXCEPTION( !( refEdgeTangent.size() == spaceDim ), std::invalid_argument,
01636                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): output array size is required to match space dimension");  
01637 #endif
01638   
01639   // Edge parametrizations are computed in setSubcellParametrization and stored in rank-3 array 
01640   // (subcOrd, coordinate, coefficient)
01641   const FieldContainer<double>& edgeMap = getSubcellParametrization(1, parentCell);
01642   
01643   // All ref. edge maps have affine coordinate functions: f_dim(u) = C_0(dim) + C_1(dim)*u, 
01644   //                                     => edge Tangent: -> C_1(*)
01645   refEdgeTangent(0) = edgeMap(edgeOrd, 0, 1);
01646   refEdgeTangent(1) = edgeMap(edgeOrd, 1, 1);
01647   
01648   // Skip last coordinate for 2D parent cells
01649   if(spaceDim == 3) {
01650     refEdgeTangent(2) = edgeMap(edgeOrd, 2, 1);  
01651   }
01652 }
01653 
01654 
01655 
01656 template<class Scalar>
01657 template<class ArrayFaceTangentU, class ArrayFaceTangentV>
01658 void CellTools<Scalar>::getReferenceFaceTangents(ArrayFaceTangentU &           uTan,
01659                                                  ArrayFaceTangentV &           vTan,
01660                                                  const int &                   faceOrd,
01661                                                  const shards::CellTopology &  parentCell){
01662   
01663 #ifdef HAVE_INTREPID_DEBUG
01664   int spaceDim  = parentCell.getDimension();
01665   TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument, 
01666                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): three-dimensional parent cell required");  
01667   
01668   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
01669                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): face ordinal out of bounds");  
01670 
01671   TEUCHOS_TEST_FOR_EXCEPTION( !( (uTan.rank() == 1)  && (vTan.rank() == 1) ), std::invalid_argument,  
01672                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): rank = 1 required for output arrays"); 
01673   
01674   TEUCHOS_TEST_FOR_EXCEPTION( !( uTan.dimension(0) == spaceDim ), std::invalid_argument,
01675                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");  
01676 
01677   TEUCHOS_TEST_FOR_EXCEPTION( !( vTan.dimension(0) == spaceDim ), std::invalid_argument,
01678                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");  
01679 #endif
01680   
01681   // Face parametrizations are computed in setSubcellParametrization and stored in rank-3 array 
01682   // (subcOrd, coordinate, coefficient): retrieve this array
01683   const FieldContainer<double>& faceMap = getSubcellParametrization(2, parentCell);
01684   
01685   /*  All ref. face maps have affine coordinate functions:  f_dim(u,v) = C_0(dim) + C_1(dim)*u + C_2(dim)*v
01686    *                           `   => Tangent vectors are:  uTan -> C_1(*);    vTan -> C_2(*)
01687    */
01688     // set uTan -> C_1(*)
01689     uTan(0) = faceMap(faceOrd, 0, 1);
01690     uTan(1) = faceMap(faceOrd, 1, 1);
01691     uTan(2) = faceMap(faceOrd, 2, 1);
01692     
01693      // set vTan -> C_2(*)
01694     vTan(0) = faceMap(faceOrd, 0, 2);
01695     vTan(1) = faceMap(faceOrd, 1, 2);
01696     vTan(2) = faceMap(faceOrd, 2, 2);
01697 }
01698 
01699 
01700 
01701 template<class Scalar>
01702 template<class ArraySideNormal>
01703 void CellTools<Scalar>::getReferenceSideNormal(ArraySideNormal &             refSideNormal,
01704                                                const int &                   sideOrd,
01705                                                const shards::CellTopology &  parentCell){
01706   int spaceDim  = parentCell.getDimension();
01707 #ifdef HAVE_INTREPID_DEBUG
01708   
01709   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01710                       ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): two or three-dimensional parent cell required");
01711   
01712   // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
01713   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= sideOrd) && (sideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
01714                       ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): side ordinal out of bounds");    
01715 #endif  
01716   
01717   if(spaceDim == 2){
01718     
01719     // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
01720     getReferenceEdgeTangent(refSideNormal, sideOrd, parentCell);
01721     
01722     // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
01723     Scalar temp = refSideNormal(0);
01724     refSideNormal(0) = refSideNormal(1);
01725     refSideNormal(1) = -temp;
01726   }
01727   else{
01728     // 3D parent cell: side = 2D subcell (face), call the face normal method.
01729     getReferenceFaceNormal(refSideNormal, sideOrd, parentCell);
01730   }
01731 }
01732   
01733 
01734 
01735 template<class Scalar>
01736 template<class ArrayFaceNormal>
01737 void CellTools<Scalar>::getReferenceFaceNormal(ArrayFaceNormal &             refFaceNormal,
01738                                                const int &                   faceOrd,
01739                                                const shards::CellTopology &  parentCell){
01740   int spaceDim  = parentCell.getDimension();
01741   
01742 #ifdef HAVE_INTREPID_DEBUG
01743   
01744   TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument, 
01745                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): three-dimensional parent cell required");  
01746   
01747   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
01748                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): face ordinal out of bounds");  
01749   
01750   TEUCHOS_TEST_FOR_EXCEPTION( !( refFaceNormal.rank() == 1 ), std::invalid_argument,  
01751                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): rank = 1 required for output array"); 
01752     
01753   TEUCHOS_TEST_FOR_EXCEPTION( !( refFaceNormal.dimension(0) == spaceDim ), std::invalid_argument,
01754                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): dim0 (spatial dim) must match that of parent cell");  
01755 #endif
01756 
01757   // Reference face normal = vector product of reference face tangents. Allocate temp FC storage:
01758   FieldContainer<Scalar> uTan(spaceDim);
01759   FieldContainer<Scalar> vTan(spaceDim);
01760   getReferenceFaceTangents(uTan, vTan, faceOrd, parentCell);
01761   
01762   // Compute the vector product of the reference face tangents:
01763   RealSpaceTools<Scalar>::vecprod(refFaceNormal, uTan, vTan);
01764 }
01765 
01766 
01767 
01768 template<class Scalar>
01769 template<class ArrayEdgeTangent, class ArrayJac>
01770 void CellTools<Scalar>::getPhysicalEdgeTangents(ArrayEdgeTangent &            edgeTangents,
01771                                                 const ArrayJac &              worksetJacobians,
01772                                                 const int &                   worksetEdgeOrd,
01773                                                 const shards::CellTopology &  parentCell){
01774   int worksetSize = worksetJacobians.dimension(0);
01775   int edgePtCount = worksetJacobians.dimension(1); 
01776   int pCellDim    = parentCell.getDimension();
01777   
01778 #ifdef HAVE_INTREPID_DEBUG
01779   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents):";
01780   
01781   TEUCHOS_TEST_FOR_EXCEPTION( !( (pCellDim == 3) || (pCellDim == 2) ), std::invalid_argument, 
01782                       ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required");  
01783   
01784   // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
01785   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, edgeTangents, 3,3), std::invalid_argument, errmsg);
01786   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, edgeTangents, 2, 2,3), std::invalid_argument, errmsg);
01787  
01788   // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
01789   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01790   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 2,3), std::invalid_argument, errmsg);
01791   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 2,3), std::invalid_argument, errmsg);
01792   
01793   // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
01794   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, edgeTangents, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);      
01795   
01796 #endif
01797   
01798   // Temp storage for constant reference edge tangent: rank-1 (D) arrays
01799   FieldContainer<double> refEdgeTan(pCellDim);
01800   getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
01801   
01802   // Loop over workset faces and edge points
01803   for(int pCell = 0; pCell < worksetSize; pCell++){
01804     for(int pt = 0; pt < edgePtCount; pt++){
01805       
01806       // Apply parent cell Jacobian to ref. edge tangent
01807       for(int i = 0; i < pCellDim; i++){
01808         edgeTangents(pCell, pt, i) = 0.0;
01809         for(int j = 0; j < pCellDim; j++){
01810           edgeTangents(pCell, pt, i) +=  worksetJacobians(pCell, pt, i, j)*refEdgeTan(j);
01811         }// for j
01812       }// for i
01813     }// for pt
01814   }// for pCell
01815 }
01816 
01817 
01818 
01819 template<class Scalar>
01820 template<class ArrayFaceTangentU, class ArrayFaceTangentV, class ArrayJac>
01821 void CellTools<Scalar>::getPhysicalFaceTangents(ArrayFaceTangentU &           faceTanU,
01822                                                 ArrayFaceTangentV &           faceTanV,
01823                                                 const ArrayJac &              worksetJacobians,
01824                                                 const int &                   worksetFaceOrd,
01825                                                 const shards::CellTopology &  parentCell){
01826   int worksetSize = worksetJacobians.dimension(0);
01827   int facePtCount = worksetJacobians.dimension(1); 
01828   int pCellDim    = parentCell.getDimension();
01829   
01830 #ifdef HAVE_INTREPID_DEBUG
01831   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents):";
01832 
01833   TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument, 
01834                       ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");  
01835   
01836   // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
01837   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanU, 3,3), std::invalid_argument, errmsg);
01838   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanU, 2, 3,3), std::invalid_argument, errmsg);
01839 
01840   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanV, 3,3), std::invalid_argument, errmsg);
01841   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanV, 2, 3,3), std::invalid_argument, errmsg);
01842 
01843   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU,  faceTanV), std::invalid_argument, errmsg);      
01844 
01845   // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
01846   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01847   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
01848   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
01849 
01850   // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
01851   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);      
01852 
01853 #endif
01854     
01855   // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
01856   FieldContainer<double> refFaceTanU(pCellDim);
01857   FieldContainer<double> refFaceTanV(pCellDim);
01858   getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
01859 
01860   // Loop over workset faces and face points
01861   for(int pCell = 0; pCell < worksetSize; pCell++){
01862     for(int pt = 0; pt < facePtCount; pt++){
01863       
01864       // Apply parent cell Jacobian to ref. face tangents
01865       for(int dim = 0; dim < pCellDim; dim++){
01866         faceTanU(pCell, pt, dim) = 0.0;
01867         faceTanV(pCell, pt, dim) = 0.0;
01868         
01869         // Unroll loops: parent cell dimension can only be 3
01870         faceTanU(pCell, pt, dim) = \
01871           worksetJacobians(pCell, pt, dim, 0)*refFaceTanU(0) + \
01872           worksetJacobians(pCell, pt, dim, 1)*refFaceTanU(1) + \
01873           worksetJacobians(pCell, pt, dim, 2)*refFaceTanU(2);
01874         faceTanV(pCell, pt, dim) = \
01875           worksetJacobians(pCell, pt, dim, 0)*refFaceTanV(0) + \
01876           worksetJacobians(pCell, pt, dim, 1)*refFaceTanV(1) + \
01877           worksetJacobians(pCell, pt, dim, 2)*refFaceTanV(2);
01878       }// for dim
01879     }// for pt
01880   }// for pCell
01881 }
01882 
01883 
01884 template<class Scalar>
01885 template<class ArraySideNormal, class ArrayJac>
01886 void CellTools<Scalar>::getPhysicalSideNormals(ArraySideNormal &             sideNormals,
01887                                                const ArrayJac &              worksetJacobians,
01888                                                const int &                   worksetSideOrd,
01889                                                const shards::CellTopology &  parentCell){
01890   int worksetSize = worksetJacobians.dimension(0);
01891   int sidePtCount = worksetJacobians.dimension(1);   
01892   int spaceDim  = parentCell.getDimension();
01893   
01894 #ifdef HAVE_INTREPID_DEBUG
01895   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01896                       ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
01897   
01898   // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
01899   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= worksetSideOrd) && (worksetSideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
01900                       ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): side ordinal out of bounds");  
01901 #endif  
01902   
01903   if(spaceDim == 2){
01904 
01905     // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
01906     getPhysicalEdgeTangents(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
01907     
01908     // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
01909     for(int cell = 0; cell < worksetSize; cell++){
01910       for(int pt = 0; pt < sidePtCount; pt++){
01911         Scalar temp = sideNormals(cell, pt, 0);
01912         sideNormals(cell, pt, 0) = sideNormals(cell, pt, 1);
01913         sideNormals(cell, pt, 1) = -temp;
01914       }// for pt
01915     }// for cell
01916   }
01917   else{
01918     // 3D parent cell: side = 2D subcell (face), call the face normal method.
01919     getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
01920   }
01921 }
01922   
01923   
01924   
01925 template<class Scalar>
01926 template<class ArrayFaceNormal, class ArrayJac>
01927 void CellTools<Scalar>::getPhysicalFaceNormals(ArrayFaceNormal &             faceNormals,
01928                                                const ArrayJac &              worksetJacobians,
01929                                                const int &                   worksetFaceOrd,
01930                                                const shards::CellTopology &  parentCell){
01931   int worksetSize = worksetJacobians.dimension(0);
01932   int facePtCount = worksetJacobians.dimension(1); 
01933   int pCellDim    = parentCell.getDimension();
01934   
01935 #ifdef HAVE_INTREPID_DEBUG
01936   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals):";
01937   
01938   TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument, 
01939                       ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required");  
01940   
01941   // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
01942   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceNormals, 3,3), std::invalid_argument, errmsg);
01943   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceNormals, 2, 3,3), std::invalid_argument, errmsg);
01944   
01945   // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
01946   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01947   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
01948   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
01949   
01950   // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
01951   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceNormals, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);        
01952 #endif
01953   
01954   // Temp storage for physical face tangents: rank-3 (C,P,D) arrays
01955   FieldContainer<Scalar> faceTanU(worksetSize, facePtCount, pCellDim);
01956   FieldContainer<Scalar> faceTanV(worksetSize, facePtCount, pCellDim);
01957   getPhysicalFaceTangents(faceTanU, faceTanV, worksetJacobians, worksetFaceOrd, parentCell);
01958   
01959   // Compute the vector product of the physical face tangents:
01960   RealSpaceTools<Scalar>::vecprod(faceNormals, faceTanU, faceTanV);
01961   
01962   
01963 }
01964 
01965 //============================================================================================//
01966 //                                                                                            //
01967 //                                        Inclusion tests                                     //
01968 //                                                                                            //
01969 //============================================================================================//
01970 
01971 
01972 template<class Scalar>
01973 int CellTools<Scalar>::checkPointInclusion(const Scalar*                 point,
01974                                            const int                     pointDim,
01975                                            const shards::CellTopology &  cellTopo,
01976                                            const double &                threshold) {
01977 #ifdef HAVE_INTREPID_DEBUG
01978   TEUCHOS_TEST_FOR_EXCEPTION( !(pointDim == (int)cellTopo.getDimension() ), std::invalid_argument,
01979                       ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Point and cell dimensions do not match. ");
01980 #endif
01981   int testResult = 1;
01982   
01983   // Using these values in the tests effectievly inflates the reference element to a larger one
01984   double minus_one = -1.0 - threshold;
01985   double plus_one  =  1.0 + threshold;
01986   double minus_zero = - threshold;
01987   
01988   // A cell with extended topology has the same reference cell as a cell with base topology. 
01989   // => testing for inclusion in a reference Triangle<> and a reference Triangle<6> relies on 
01990   // on the same set of inequalities. To eliminate unnecessary cases we switch on the base topology
01991   unsigned key = cellTopo.getBaseCellTopologyData() -> key ;
01992   switch( key ) {
01993     
01994     case shards::Line<>::key :
01995       if( !(minus_one <= point[0] && point[0] <= plus_one))  testResult = 0;
01996       break;
01997       
01998     case shards::Triangle<>::key : {
01999       Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1.0 );
02000       if( distance > threshold ) testResult = 0;
02001       break;
02002     }
02003       
02004     case shards::Quadrilateral<>::key :
02005       if(!( (minus_one <= point[0] && point[0] <= plus_one) && \
02006             (minus_one <= point[1] && point[1] <= plus_one) ) ) testResult = 0;   
02007       break;
02008       
02009     case shards::Tetrahedron<>::key : {
02010       Scalar distance = std::max(  std::max(-point[0],-point[1]), \
02011                                    std::max(-point[2], point[0] + point[1] + point[2] - 1)  );
02012       if( distance > threshold ) testResult = 0;
02013       break;
02014     }
02015       
02016     case shards::Hexahedron<>::key :
02017       if(!((minus_one <= point[0] && point[0] <= plus_one) && \
02018            (minus_one <= point[1] && point[1] <= plus_one) && \
02019            (minus_one <= point[2] && point[2] <= plus_one)))  \
02020              testResult = 0;
02021       break;
02022       
02023     // The base of the reference prism is the same as the reference triangle => apply triangle test
02024     // to X and Y coordinates and test whether Z is in [-1,1]
02025     case shards::Wedge<>::key : {
02026       Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1 );
02027       if( distance > threshold  || \
02028           point[2] < minus_one || point[2] > plus_one) \
02029             testResult = 0;
02030       break;
02031     }
02032 
02033     // The base of the reference pyramid is the same as the reference quad cell => a horizontal plane
02034     // through a point P(x,y,z) intersects the pyramid at a quadrilateral that equals the base quad 
02035     // scaled by (1-z) => P(x,y,z) is inside the pyramid <=> (x,y) is in [-1+z,1-z]^2 && 0 <= Z <= 1 
02036     case shards::Pyramid<>::key : {
02037       Scalar left  = minus_one + point[2];
02038       Scalar right = plus_one  - point[2];
02039       if(!( (left       <= point[0] && point[0] <= right) && \
02040             (left       <= point[1] && point[1] <= right) && 
02041             (minus_zero <= point[2] && point[2] <= plus_one) ) )  \
02042              testResult = 0;  
02043       break;
02044     }
02045       
02046     default:
02047       TEUCHOS_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) ||
02048                              (key == shards::Triangle<>::key)  ||
02049                              (key == shards::Quadrilateral<>::key) ||
02050                              (key == shards::Tetrahedron<>::key)  ||
02051                              (key == shards::Hexahedron<>::key)  ||
02052                              (key == shards::Wedge<>::key)  ||
02053                              (key == shards::Pyramid<>::key) ),
02054                           std::invalid_argument,
02055                           ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Invalid cell topology. ");
02056   }
02057   return testResult;
02058 }
02059 
02060 
02061 
02062 template<class Scalar>
02063 template<class ArrayPoint>
02064 int CellTools<Scalar>::checkPointsetInclusion(const ArrayPoint&             points,
02065                                               const shards::CellTopology &  cellTopo, 
02066                                               const double &                threshold) {
02067   
02068   int rank = points.rank();  
02069   
02070 #ifdef HAVE_INTREPID_DEBUG
02071   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= points.rank() ) && (points.rank() <= 3) ), std::invalid_argument,
02072                       ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): rank-1, 2 or 3 required for input points array. ");
02073 
02074   // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
02075   TEUCHOS_TEST_FOR_EXCEPTION( !( points.dimension(rank - 1) == (int)cellTopo.getDimension() ), std::invalid_argument,
02076                       ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): Point and cell dimensions do not match. ");
02077 #endif
02078   
02079   // create temp output array depending on the rank of the input array 
02080   FieldContainer<int> inRefCell;
02081   switch(rank) {
02082     case 1: inRefCell.resize(1); break;
02083     case 2: inRefCell.resize( points.dimension(0) ); break;
02084     case 3: inRefCell.resize( points.dimension(0), points.dimension(1) ); break;
02085   }
02086 
02087   // Call the inclusion method which returns inclusion results for all points
02088   checkPointwiseInclusion(inRefCell, points, cellTopo, threshold);
02089   
02090   // Check if any points were outside, break when finding the first one
02091   int allInside = 1;
02092   for(int i = 0; i < inRefCell.size(); i++ ){
02093     if (inRefCell[i] == 0) {
02094       allInside = 0;
02095       break;
02096     }
02097   }
02098    return allInside;
02099 }
02100 
02101 
02102 
02103 template<class Scalar>
02104 template<class ArrayIncl, class ArrayPoint>
02105 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl &                   inRefCell,
02106                                                 const ArrayPoint &            points,
02107                                                 const shards::CellTopology &  cellTopo, 
02108                                                 const double &                threshold) {
02109   int apRank   = points.rank();
02110   
02111 #ifdef HAVE_INTREPID_DEBUG
02112   
02113   // Verify that points and inRefCell have correct ranks and dimensions
02114   std::string errmsg = ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion):";
02115   if(points.rank() == 1) {
02116     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 1 ), std::invalid_argument, 
02117                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires rank-1 output array.");  
02118     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.dimension(0) == 1), std::invalid_argument,
02119                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires dim0 = 1 for output array.");  
02120   }
02121   else if(points.rank() == 2){
02122     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 1 ), std::invalid_argument, 
02123                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-2 input array requires rank-1 output array.");  
02124     // dimension 0 of the arrays must match
02125     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0,  points, 0), std::invalid_argument, errmsg);
02126   }
02127   else if (points.rank() == 3) {
02128     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 2 ), std::invalid_argument, 
02129                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-3 input array requires rank-2 output array.");  
02130     // dimensions 0 and 1 of the arrays must match
02131     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0,1,  points, 0,1), std::invalid_argument, errmsg);
02132   }
02133   else{
02134     TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 1) || (points.rank() == 2) || (points.rank() == 3) ), std::invalid_argument,
02135                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");      
02136   }    
02137   
02138   // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
02139   TEUCHOS_TEST_FOR_EXCEPTION( !( points.dimension(apRank - 1) == (int)cellTopo.getDimension() ), std::invalid_argument,
02140                       ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Point and cell dimensions do not match. ");
02141   
02142 #endif
02143   
02144   // Initializations
02145   int dim0     = 1;
02146   int dim1     = 1;
02147   int pointDim = 0;
02148   switch(apRank) {
02149     case 1:
02150       pointDim = points.dimension(0);
02151       break;
02152     case 2:
02153       dim1     = points.dimension(0);
02154       pointDim = points.dimension(1);
02155       break;
02156     case 3:
02157       dim0     = points.dimension(0);
02158       dim1     = points.dimension(1);
02159       pointDim = points.dimension(2);
02160       break;
02161     default:
02162       TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= points.rank() ) && (points.rank() <= 3) ), std::invalid_argument,
02163                           ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");      
02164   }// switch
02165   
02166   
02167   // This method can handle up to rank-3 input arrays. The spatial dim must be the last dimension. 
02168   // The method uses [] accessor because array rank is determined at runtime and the appropriate
02169   // (i,j,..,k) accessor is not known. Use of [] requires the following offsets:
02170   //    for input array  = i0*dim1*pointDim + i1*dim1  (computed in 2 pieces: inPtr0 and inPtr1, resp)
02171   //    for output array = i0*dim1                     (computed in one piece: outPtr0)
02172   int inPtr0  = 0;
02173   int inPtr1  = 0;
02174   int outPtr0 = 0;
02175   Scalar point[3] = {0.0, 0.0, 0.0};
02176   
02177   for(int i0 = 0; i0 < dim0; i0++){
02178     outPtr0 = i0*dim1;
02179     inPtr0  = outPtr0*pointDim;
02180     
02181     for(int i1 = 0; i1 < dim1; i1++) {
02182       inPtr1 = inPtr0 + i1*pointDim;      
02183       point[0] = points[inPtr1];
02184       if(pointDim > 1) {
02185         point[1] = points[inPtr1 + 1];
02186         if(pointDim > 2) {
02187           point[2] = points[inPtr1 + 2];
02188           if(pointDim > 3) {
02189             TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= pointDim) && (pointDim <= 3)), std::invalid_argument, 
02190                                 ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Input array specifies invalid point dimension ");      
02191           }
02192         }
02193       } //if(pointDim > 1)
02194       inRefCell[outPtr0 + i1] = checkPointInclusion(point, pointDim, cellTopo, threshold);
02195     } // for (i1)
02196   } // for(i2)
02197 
02198 }  
02199 
02200 
02201 template<class Scalar>
02202 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
02203 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl &                   inCell,
02204                                                 const ArrayPoint &            points,
02205                                                 const ArrayCell &             cellWorkset,
02206                                                 const shards::CellTopology &  cell,
02207                                                 const int &                   whichCell, 
02208                                                 const double &                threshold)
02209 {
02210   INTREPID_VALIDATE( validateArguments_checkPointwiseInclusion(inCell, points, cellWorkset, whichCell, cell) );
02211   
02212   // For cell topologies with reference cells this test maps the points back to the reference cell
02213   // and uses the method for reference cells
02214   unsigned baseKey = cell.getBaseCellTopologyData() -> key;
02215   
02216   switch(baseKey){
02217     
02218     case shards::Line<>::key :
02219     case shards::Triangle<>::key:
02220     case shards::Quadrilateral<>::key :
02221     case shards::Tetrahedron<>::key :
02222     case shards::Hexahedron<>::key :
02223     case shards::Wedge<>::key :
02224     case shards::Pyramid<>::key :
02225       {
02226         FieldContainer<Scalar> refPoints;
02227         
02228         if(points.rank() == 2){
02229           refPoints.resize(points.dimension(0), points.dimension(1) );
02230           mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
02231           checkPointwiseInclusion(inCell, refPoints, cell, threshold );
02232         }
02233         else if(points.rank() == 3){
02234           refPoints.resize(points.dimension(0), points.dimension(1), points.dimension(2) );
02235           mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
02236           checkPointwiseInclusion(inCell, refPoints, cell, threshold );          
02237         }
02238         break;
02239       }
02240     default: 
02241       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
02242                           ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): cell topology not supported");
02243   }// switch
02244   
02245 }
02246 
02247 
02248 //============================================================================================//
02249 //                                                                                            //
02250 //                  Validation of input/output arguments for CellTools methods                //
02251 //                                                                                            //
02252 //============================================================================================//
02253 
02254 template<class Scalar>
02255 template<class ArrayJac, class ArrayPoint, class ArrayCell>
02256 void CellTools<Scalar>::validateArguments_setJacobian(const ArrayJac    &          jacobian,
02257                                                       const ArrayPoint  &          points,
02258                                                       const ArrayCell   &          cellWorkset,
02259                                                       const int &                  whichCell,
02260                                                       const shards::CellTopology & cellTopo){
02261   
02262   // Validate cellWorkset array
02263   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02264                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for cellWorkset array");
02265   
02266   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02267                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) >= 1 required for cellWorkset array");
02268   
02269   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02270                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02271   
02272   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02273                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02274     
02275   // validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02276   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02277                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): whichCell = -1 or a valid cell ordinal is required.");
02278   
02279   
02280   // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
02281   // If rank-2: admissible jacobians: rank-3 (P,D,D) or rank-4 (C,P,D,D); admissible whichCell: -1 (default) or cell ordinal.
02282   if(points.rank() == 2) {
02283     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(0) <= 0), std::invalid_argument,
02284                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) >= 1 required for points array ");
02285     
02286     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(1) != (int)cellTopo.getDimension() ), std::invalid_argument,
02287                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of points array does not match cell dimension");
02288     
02289     // Validate the output array for the Jacobian: if whichCell == -1 all Jacobians are computed, rank-4 (C,P,D,D) required
02290     if(whichCell == -1) {
02291       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.rank() != 4), std::invalid_argument, 
02292                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 4 required for jacobian array");
02293       
02294       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != cellWorkset.dimension(0)), std::invalid_argument,
02295                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of cellWorkset array");
02296       
02297       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(0)), std::invalid_argument,
02298                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 0 of points array");
02299 
02300       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(2) != points.dimension(1)), std::invalid_argument,
02301                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 1 of points array");
02302       
02303       TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(2) == jacobian.dimension(3) ), std::invalid_argument,
02304                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
02305       
02306       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(3) ) && (jacobian.dimension(3) < 4) ), std::invalid_argument,
02307                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
02308     }     
02309     // A single cell Jacobian is computed when whichCell != -1 (whichCell has been already validated), rank-3 (P,D,D) required
02310     else {
02311       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.rank() != 3), std::invalid_argument, 
02312                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for jacobian array");
02313       
02314       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != points.dimension(0)), std::invalid_argument,
02315                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) of jacobian array must equal dim 0 of points array");
02316 
02317       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(1)), std::invalid_argument,
02318                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of jacobian array must equal dim 1 of points array");
02319       
02320       TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(1) == jacobian.dimension(2) ), std::invalid_argument,
02321                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 = dim 2 (same spatial dimensions) required for jacobian array. ");
02322       
02323       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(1) ) && (jacobian.dimension(1) < 4) ), std::invalid_argument,
02324                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 and dim 2 (spatial dimensions) must be between 1 and 3. ");
02325     }
02326   }
02327   // Point array is rank-3 (C,P,D): requires whichCell = -1 and rank-4 (C,P,D,D) jacobians
02328   else if(points.rank() ==3){
02329     std::string errmsg  = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian):";
02330     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02331                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of points array must equal dim 0 of cellWorkset array");
02332 
02333     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(1) <= 0), std::invalid_argument,
02334                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) >= 1 required for points array ");
02335     
02336     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02337                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of points array does not match cell dimension");
02338     
02339     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument,
02340                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): default value whichCell=-1 required for rank-3 input points");
02341     
02342     // rank-4 (C,P,D,D) jacobian required for rank-3 (C,P,D) input points
02343     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, jacobian,  4, 4), std::invalid_argument,errmsg);
02344     
02345     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != points.dimension(0)), std::invalid_argument,
02346                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of points array");
02347     
02348     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(1)), std::invalid_argument,
02349                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 1 of points array");
02350   
02351     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(2) != points.dimension(2)), std::invalid_argument,
02352                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 2 of points array");
02353     
02354     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(2) == jacobian.dimension(3) ), std::invalid_argument,
02355                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
02356     
02357     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(3) ) && (jacobian.dimension(3) < 4) ), std::invalid_argument,
02358                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
02359   }
02360   else {
02361     TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 2) && (points.rank() ==3) ), std::invalid_argument,
02362                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 2 or 3 required for points array");
02363   }  
02364 }
02365 
02366 
02367 
02368 template<class Scalar>
02369 template<class ArrayJacInv, class ArrayJac>
02370 void CellTools<Scalar>::validateArguments_setJacobianInv(const ArrayJacInv & jacobianInv,
02371                                                          const ArrayJac &    jacobian)
02372 {
02373   // Validate input jacobian array: admissible ranks & dimensions are: 
02374   // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
02375   int jacobRank = jacobian.rank();
02376   TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
02377                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
02378   
02379   // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
02380   TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
02381                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
02382   
02383   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
02384                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
02385   
02386   // Validate output jacobianInv array: must have the same rank and dimensions as the input array.
02387   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv):";
02388 
02389   TEUCHOS_TEST_FOR_EXCEPTION( !(requireRankMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
02390   
02391   TEUCHOS_TEST_FOR_EXCEPTION( !(requireDimensionMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
02392 }
02393 
02394 
02395 
02396 template<class Scalar>
02397 template<class ArrayJacDet, class ArrayJac>
02398 void CellTools<Scalar>::validateArguments_setJacobianDetArgs(const ArrayJacDet &  jacobianDet,
02399                                                              const ArrayJac    &  jacobian)
02400 {
02401   // Validate input jacobian array: admissible ranks & dimensions are: 
02402   // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
02403   int jacobRank = jacobian.rank();
02404   TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
02405                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
02406   
02407   // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
02408   TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
02409                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
02410   
02411   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
02412                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
02413 
02414   
02415   // Validate output jacobianDet array: must be rank-2 with dimensions (C,P) if jacobian was rank-4:
02416   if(jacobRank == 4){
02417     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.rank() == 2), std::invalid_argument,
02418                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 2 required for jacobianDet if jacobian is rank-4. ");
02419     
02420     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(0) == jacobian.dimension(0) ), std::invalid_argument,
02421                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of cells) of jacobianDet array must equal dim 0 of jacobian array. ");
02422     
02423     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(1) == jacobian.dimension(1) ), std::invalid_argument,
02424                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 1 (number of points) of jacobianDet array must equal dim 1 of jacobian array.");  
02425   }
02426   
02427   // must be rank-1 with dimension (P) if jacobian was rank-3
02428   else {
02429     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.rank() == 1), std::invalid_argument,
02430                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 1 required for jacobianDet if jacobian is rank-3. ");
02431     
02432     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(0) == jacobian.dimension(0) ), std::invalid_argument,
02433                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of points) of jacobianDet array must equal dim 0 of jacobian array.");  
02434   }
02435 }
02436 
02437 
02438 
02439 template<class Scalar>
02440 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
02441 void CellTools<Scalar>::validateArguments_mapToPhysicalFrame(const ArrayPhysPoint &        physPoints,
02442                                                              const ArrayRefPoint  &        refPoints,
02443                                                              const ArrayCell      &        cellWorkset,
02444                                                              const shards::CellTopology &  cellTopo,
02445                                                              const int&                    whichCell)
02446 {
02447   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame):";
02448   
02449   // Validate cellWorkset array
02450   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02451                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for cellWorkset array");
02452   
02453   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02454                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
02455   
02456   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02457                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02458   
02459   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02460                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02461   
02462     
02463   // validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02464   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02465                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): whichCell = -1 or a valid cell ordinal is required.");
02466   
02467   // Validate refPoints array: can be rank-2 (P,D) or rank-3 (C,P,D) array
02468   // If rank-2: admissible output array is (P,D) or (C,P,D); admissible whichCell: -1 (default) or cell ordinal
02469   if(refPoints.rank() == 2) {
02470     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(0) <= 0), std::invalid_argument,
02471                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) >= 1 required for refPoints array ");
02472     
02473     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(1) != (int)cellTopo.getDimension() ), std::invalid_argument,
02474                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) of refPoints array does not match cell dimension");
02475 
02476     // Validate output array: whichCell = -1 requires rank-3 array with dimensions (C,P,D)  
02477     if(whichCell == -1) {
02478       TEUCHOS_TEST_FOR_EXCEPTION( ( (physPoints.rank() != 3) && (whichCell == -1) ), std::invalid_argument,
02479                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for physPoints array for the default whichCell value");
02480       
02481       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != cellWorkset.dimension(0)), std::invalid_argument,
02482                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of physPoints array must equal dim 0 of cellWorkset array");
02483       
02484       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != refPoints.dimension(0)), std::invalid_argument,
02485                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) of physPoints array must equal dim 0 of refPoints array"); 
02486       
02487       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(2) != (int)cellTopo.getDimension()), std::invalid_argument,
02488                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) does not match cell dimension ");  
02489     }
02490     // 0 <= whichCell < num cells requires rank-2 (P,D) arrays for both refPoints and physPoints
02491     else{
02492       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.rank() != 2), std::invalid_argument,
02493                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 required for physPoints array");
02494       
02495       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != refPoints.dimension(0)), std::invalid_argument,
02496                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) of physPoints array must equal dim 0 of refPoints array"); 
02497       
02498       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != (int)cellTopo.getDimension()), std::invalid_argument,
02499                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) does not match cell dimension ");      
02500     }
02501   }
02502   // refPoints is (C,P,D): requires physPoints to be (C,P,D) and whichCell=-1  (because all cell mappings are applied)
02503   else if(refPoints.rank() == 3) {
02504     
02505     // 1. validate refPoints dimensions and rank
02506     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02507                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of refPoints and cellWorkset arraya are required to match ");
02508 
02509     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(1) <= 0), std::invalid_argument,
02510                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) >= 1 required for refPoints array ");
02511     
02512     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02513                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of refPoints array does not match cell dimension");
02514     
02515     // 2. whichCell  must be -1
02516     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument, 
02517                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): default value is required for rank-3 refPoints array");
02518 
02519     // 3.  physPoints must match rank and dimensions of refPoints
02520     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg );
02521     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg);
02522   }
02523   // if rank is not 2 or 3 throw exception
02524   else {
02525     TEUCHOS_TEST_FOR_EXCEPTION( !( (refPoints.rank() == 2) || (refPoints.rank() == 3) ), std::invalid_argument,
02526                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 or 3 required for refPoints array");
02527   }
02528 }
02529 
02530 
02531 
02532 template<class Scalar>
02533 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
02534 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint  &        refPoints,
02535                                                               const ArrayPhysPoint &        physPoints,
02536                                                               const ArrayCell      &        cellWorkset,
02537                                                               const shards::CellTopology &  cellTopo,
02538                                                               const int&                    whichCell)
02539 {
02540   std::string errmsg  = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02541   std::string errmsg1 = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02542   
02543   // Validate cellWorkset array
02544   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02545                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): rank = 3 required for cellWorkset array");
02546   
02547   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02548                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
02549   
02550   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02551                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02552   
02553   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02554                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02555     
02556   // Validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02557   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02558                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): whichCell = -1 or a valid cell ordinal is required.");
02559   
02560   // Admissible ranks and dimensions of refPoints and physPoints depend on whichCell value:
02561   // default is to map multiple sets of points to multiple sets of points. (C,P,D) arrays required
02562   int validRank;
02563   if(whichCell == -1) {
02564     validRank = 3;
02565     errmsg1 += " default value of whichCell requires rank-3 arrays:";
02566   }
02567   // whichCell is valid cell ordinal => we map single set of pts to a single set of pts. (P,D) arrays required
02568   else{
02569     errmsg1 += " rank-2 arrays required when whichCell is valid cell ordinal";
02570     validRank = 2;
02571   }
02572   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg1, refPoints,  validRank,validRank), std::invalid_argument, errmsg1);
02573   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg1, physPoints, refPoints),           std::invalid_argument, errmsg1);
02574   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg1, refPoints, physPoints),      std::invalid_argument, errmsg1);
02575 }
02576 
02577 
02578 
02579 template<class Scalar>
02580 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
02581 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint  &        refPoints,
02582                                                               const ArrayInitGuess &        initGuess,
02583                                                               const ArrayPhysPoint &        physPoints,
02584                                                               const ArrayCell      &        cellWorkset,
02585                                                               const shards::CellTopology &  cellTopo,
02586                                                               const int&                    whichCell)
02587 {
02588   // Call the method that validates arguments with the default initial guess selection
02589   validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell);
02590   
02591   // Then check initGuess: its rank and dimensions must match those of physPoints.
02592   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02593   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, initGuess, physPoints), std::invalid_argument, errmsg);  
02594 }
02595 
02596 
02597 template<class Scalar>
02598 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
02599 void CellTools<Scalar>::validateArguments_checkPointwiseInclusion(ArrayIncl &                   inCell,
02600                                                                   const ArrayPoint &            physPoints,
02601                                                                   const ArrayCell &             cellWorkset,
02602                                                                   const int &                   whichCell,
02603                                                                   const shards::CellTopology &  cell)
02604 {
02605   // Validate cellWorkset array
02606   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02607                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 3 required for cellWorkset array");
02608   
02609   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02610                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) >= 1 required for cellWorkset array");
02611   
02612   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cell.getSubcellCount(0) ), std::invalid_argument,
02613                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02614   
02615   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cell.getDimension() ), std::invalid_argument,
02616                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02617   
02618   
02619   // Validate whichCell It can be either -1 (default value) or a valid cell ordinal.
02620   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02621                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 or a valid cell ordinal is required.");
02622   
02623   
02624   // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
02625   // If rank-2: admissible inCell is rank-1 (P); admissible whichCell is valid cell ordinal but not -1.
02626   if(physPoints.rank() == 2) {
02627     
02628     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell == -1), std::invalid_argument,
02629                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = a valid cell ordinal is required with rank-2 input array.");
02630 
02631     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) <= 0), std::invalid_argument,
02632                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) >= 1 required for physPoints array ");
02633     
02634     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != (int)cell.getDimension() ), std::invalid_argument,
02635                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (spatial dimension) of physPoints array does not match cell dimension");
02636     
02637     // Validate inCell
02638     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.rank() != 1), std::invalid_argument, 
02639                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 1 required for inCell array");
02640     
02641     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(0) != physPoints.dimension(0)), std::invalid_argument,
02642                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) of inCell array must equal dim 0 of physPoints array");
02643   }
02644   // If rank-3: admissible inCell is rank-2 (C,P); admissible whichCell = -1.
02645   else if (physPoints.rank() == 3){
02646     
02647     TEUCHOS_TEST_FOR_EXCEPTION( !(whichCell == -1), std::invalid_argument,
02648                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 is required with rank-3 input array.");
02649     
02650     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02651                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells)  of physPoints array must equal dim 0 of cellWorkset array ");
02652 
02653     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) <= 0), std::invalid_argument,
02654                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) >= 1 required for physPoints array ");
02655     
02656     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(2) != (int)cell.getDimension() ), std::invalid_argument,
02657                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of physPoints array does not match cell dimension");
02658     
02659     // Validate inCell
02660     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.rank() != 2), std::invalid_argument, 
02661                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 required for inCell array");
02662     
02663     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(0) != physPoints.dimension(0)), std::invalid_argument,
02664                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) of inCell array must equal dim 0 of physPoints array");    
02665 
02666     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(1) != physPoints.dimension(1)), std::invalid_argument,
02667                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) of inCell array must equal dim 1 of physPoints array");    
02668   }
02669   else {
02670     TEUCHOS_TEST_FOR_EXCEPTION( !( (physPoints.rank() == 2) && (physPoints.rank() ==3) ), std::invalid_argument,
02671                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 or 3 required for points array");
02672   }
02673 }
02674 
02675 
02676 
02677 //============================================================================================//
02678 //                                                                                            //
02679 //                                           Debug                                            //
02680 //                                                                                            //
02681 //============================================================================================//
02682 
02683 
02684 template<class Scalar>
02685 void CellTools<Scalar>::printSubcellVertices(const int subcellDim,
02686                                              const int subcellOrd,
02687                                              const shards::CellTopology & parentCell){
02688   
02689   // Get number of vertices for the specified subcell and parent cell dimension
02690   int subcVertexCount = parentCell.getVertexCount(subcellDim, subcellOrd);
02691   int cellDim         = parentCell.getDimension();
02692   
02693   // Allocate space for the subcell vertex coordinates
02694   FieldContainer<double> subcellVertices(subcVertexCount, cellDim);
02695   
02696   // Retrieve the vertex coordinates
02697   getReferenceSubcellVertices(subcellVertices,
02698                               subcellDim,
02699                               subcellOrd,
02700                               parentCell);
02701   
02702   // Print the vertices
02703   std::cout 
02704     << " Subcell " << std::setw(2) << subcellOrd 
02705     <<  " is " << parentCell.getName(subcellDim, subcellOrd) << " with vertices = {";
02706   
02707   // Loop over subcell vertices
02708   for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
02709     std::cout<< "(";
02710     
02711     // Loop over vertex Cartesian coordinates
02712     for(int dim = 0; dim < (int)parentCell.getDimension(); dim++){
02713       std::cout << subcellVertices(subcVertOrd, dim);
02714       if(dim < (int)parentCell.getDimension()-1 ) { std::cout << ","; }
02715     }
02716     std::cout<< ")";
02717     if(subcVertOrd < subcVertexCount - 1) { std::cout << ", "; }
02718   }
02719   std::cout << "}\n";
02720 }
02721   
02722 
02723 template<class Scalar>
02724 template<class ArrayCell>
02725 void CellTools<Scalar>::printWorksetSubcell(const ArrayCell &             cellWorkset,
02726                                             const shards::CellTopology &  parentCell,
02727                                             const int&                    pCellOrd,
02728                                             const int&                    subcellDim,
02729                                             const int&                    subcellOrd,
02730                                             const int&                    fieldWidth){
02731   
02732   // Get the ordinals, relative to reference cell, of subcell cellWorkset
02733   int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
02734   int pCellDim      = parentCell.getDimension();
02735   std::vector<int> subcNodeOrdinals(subcNodeCount);
02736   
02737   for(int i = 0; i < subcNodeCount; i++){
02738     subcNodeOrdinals[i] = parentCell.getNodeMap(subcellDim, subcellOrd, i);
02739   }
02740   
02741   // Loop over parent cells and print subcell cellWorkset
02742   
02743   std::cout 
02744     << " Subcell " << subcellOrd << " on parent cell " << pCellOrd << " is " 
02745     << parentCell.getName(subcellDim, subcellOrd) << " with node(s) \n ({";
02746   
02747   for(int i = 0; i < subcNodeCount; i++){
02748     
02749     // print Cartesian coordinates of the node
02750     for(int dim = 0; dim < pCellDim; dim++){
02751       std::cout
02752       << std::setw(fieldWidth) << std::right << cellWorkset(pCellOrd, subcNodeOrdinals[i], dim); 
02753       if(dim < pCellDim - 1){ std::cout << ","; }
02754     }
02755     std::cout << "}";
02756     if(i < subcNodeCount - 1){ std::cout <<", {"; }
02757   }
02758   std::cout << ")\n\n";
02759 }
02760 
02761 
02762 
02763 } // namespace Intrepid
02764 #endif