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<27>::key:
00961         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00962         break;
00963         
00964       case shards::Wedge<18>::key:
00965         HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
00966         break;
00967         
00968         // These extended topologies are not used for mapping purposes
00969       case shards::Quadrilateral<8>::key:
00970       case shards::Hexahedron<20>::key:
00971       case shards::Wedge<15>::key:
00972         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
00973                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
00974         break;
00975         
00976         // Base and Extended Line, Beam and Shell topologies  
00977       case shards::Line<3>::key:
00978       case shards::Beam<2>::key:
00979       case shards::Beam<3>::key:
00980       case shards::ShellLine<2>::key:
00981       case shards::ShellLine<3>::key:
00982       case shards::ShellTriangle<3>::key:
00983       case shards::ShellTriangle<6>::key:
00984       case shards::ShellQuadrilateral<4>::key:
00985       case shards::ShellQuadrilateral<8>::key:
00986       case shards::ShellQuadrilateral<9>::key:
00987         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
00988                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
00989         break;
00990       default:
00991         TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
00992                             ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");        
00993     }// switch  
00994     
00995     // Temp (F,P,D) array for the values of basis functions gradients at the reference points
00996     int basisCardinality = HGRAD_Basis -> getCardinality();
00997     FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
00998     
00999     // Initialize jacobian
01000     for(int i = 0; i < jacobian.size(); i++){
01001       jacobian[i] = 0.0;
01002     }
01003         
01004     // Handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of points arrays.
01005     switch(points.rank()) {
01006       
01007       // refPoints is (P,D): a single or multiple cell jacobians computed for a single set of ref. points
01008       case 2:
01009         {
01010           // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
01011           FieldContainer<Scalar> tempPoints( points.dimension(0), points.dimension(1) );
01012           // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01013           for(int pt = 0; pt < points.dimension(0); pt++){
01014             for(int dm = 0; dm < points.dimension(1) ; dm++){
01015               tempPoints(pt, dm) = points(pt, dm);
01016             }//dm
01017           }//pt
01018           HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
01019           
01020           // The outer loops select the multi-index of the Jacobian entry: cell, point, row, col
01021           // If whichCell = -1, all jacobians are computed, otherwise a single cell jacobian is computed
01022           int cellLoop = (whichCell == -1) ? numCells : 1 ;
01023           
01024           if(whichCell == -1) {
01025             for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01026               for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01027                 for(int row = 0; row < spaceDim; row++){
01028                   for(int col = 0; col < spaceDim; col++){
01029                     
01030                     // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
01031                     for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01032                       jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01033                     } // bfOrd
01034                   } // col
01035                 } // row
01036               } // pointOrd
01037             } // cellOrd
01038           }
01039           else {
01040             for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01041               for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01042                 for(int row = 0; row < spaceDim; row++){
01043                   for(int col = 0; col < spaceDim; col++){
01044                   
01045                     // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
01046                     for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01047                       jacobian(pointOrd, row, col) += cellWorkset(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01048                     } // bfOrd
01049                   } // col
01050                 } // row
01051               } // pointOrd
01052             } // cellOrd
01053           } // if whichcell
01054         }// case 2
01055         break;
01056         
01057         // points is (C,P,D): multiple jacobians computed at multiple point sets, one jacobian per cell  
01058       case 3:
01059         {
01060           // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
01061           FieldContainer<Scalar> tempPoints( points.dimension(1), points.dimension(2) );
01062           
01063           for(int cellOrd = 0; cellOrd < numCells; cellOrd++) {
01064             
01065             // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01066             for(int pt = 0; pt < points.dimension(1); pt++){
01067               for(int dm = 0; dm < points.dimension(2) ; dm++){
01068                 tempPoints(pt, dm) = points(cellOrd, pt, dm);
01069               }//dm
01070             }//pt
01071             
01072             // Compute gradients of basis functions at this set of ref. points
01073             HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
01074             
01075             // Compute jacobians for the point set corresponding to the current cellordinal
01076             for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01077               for(int row = 0; row < spaceDim; row++){
01078                 for(int col = 0; col < spaceDim; col++){
01079                   
01080                   // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same
01081                   for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01082                     jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
01083                   } // bfOrd
01084                 } // col
01085               } // row
01086             } // pointOrd
01087           }//cellOrd
01088         }// case 3
01089         
01090         break;
01091         
01092       default:
01093         TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 2) && (points.rank() == 3) ), std::invalid_argument,
01094                             ">>> ERROR (Intrepid::CellTools::setJacobian): rank 2 or 3 required for points array. ");        
01095     }//switch
01096   }
01097   
01098 
01099 
01100 template<class Scalar>
01101 template<class ArrayJacInv, class ArrayJac>
01102 void CellTools<Scalar>::setJacobianInv(ArrayJacInv &     jacobianInv,
01103                                        const ArrayJac &  jacobian) 
01104 {
01105   INTREPID_VALIDATE( validateArguments_setJacobianInv(jacobianInv, jacobian) );
01106 
01107   RealSpaceTools<Scalar>::inverse(jacobianInv, jacobian);
01108 }
01109 
01110 
01111 
01112 template<class Scalar>
01113 template<class ArrayJacDet, class ArrayJac>
01114 void CellTools<Scalar>::setJacobianDet(ArrayJacDet &     jacobianDet,
01115                                        const ArrayJac &  jacobian)
01116 {
01117   INTREPID_VALIDATE( validateArguments_setJacobianDetArgs(jacobianDet, jacobian) );
01118 
01119   RealSpaceTools<Scalar>::det(jacobianDet, jacobian);
01120 }
01121 
01122 //============================================================================================//
01123 //                                                                                            //
01124 //                      Reference-to-physical frame mapping and its inverse                   //
01125 //                                                                                            //
01126 //============================================================================================//
01127 
01128 template<class Scalar>
01129 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
01130 void CellTools<Scalar>::mapToPhysicalFrame(ArrayPhysPoint      &        physPoints,
01131                                            const ArrayRefPoint &        refPoints,
01132                                            const ArrayCell     &        cellWorkset,
01133                                            const shards::CellTopology & cellTopo,
01134                                            const int &                  whichCell)
01135 {
01136   INTREPID_VALIDATE(validateArguments_mapToPhysicalFrame( physPoints, refPoints, cellWorkset, cellTopo, whichCell) );
01137   
01138   int spaceDim  = (int)cellTopo.getDimension();
01139   int numCells  = cellWorkset.dimension(0);
01140   //points can be rank-2 (P,D), or rank-3 (C,P,D)
01141   int numPoints = (refPoints.rank() == 2) ? refPoints.dimension(0) : refPoints.dimension(1);
01142     
01143   // Mapping is computed using an appropriate H(grad) basis function: define RCP to the base class
01144   Teuchos::RCP<Basis<Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
01145   
01146   // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
01147   switch( cellTopo.getKey() ){
01148     
01149     // Standard Base topologies (number of cellWorkset = number of vertices)
01150     case shards::Line<2>::key:
01151       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01152       break;
01153       
01154     case shards::Triangle<3>::key:
01155       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01156       break;
01157       
01158     case shards::Quadrilateral<4>::key:
01159       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01160       break;
01161       
01162     case shards::Tetrahedron<4>::key:
01163       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01164       break;
01165       
01166     case shards::Hexahedron<8>::key:
01167       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01168       break;
01169       
01170     case shards::Wedge<6>::key:
01171       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01172       break;
01173       
01174     case shards::Pyramid<5>::key:
01175           HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
01176           break;
01177 
01178     // Standard Extended topologies
01179     case shards::Triangle<6>::key:    
01180       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01181       break;
01182       
01183     case shards::Quadrilateral<9>::key:
01184       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01185       break;
01186       
01187     case shards::Tetrahedron<10>::key:
01188       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01189       break;
01190 
01191     case shards::Tetrahedron<11>::key:
01192       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
01193       break;
01194 
01195     case shards::Hexahedron<27>::key:
01196       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01197       break;
01198       
01199     case shards::Wedge<18>::key:
01200       HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
01201       break;
01202       
01203     // These extended topologies are not used for mapping purposes
01204     case shards::Quadrilateral<8>::key:
01205     case shards::Hexahedron<20>::key:
01206     case shards::Wedge<15>::key:
01207       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01208                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
01209       break;
01210       
01211     // Base and Extended Line, Beam and Shell topologies  
01212     case shards::Line<3>::key:
01213     case shards::Beam<2>::key:
01214     case shards::Beam<3>::key:
01215     case shards::ShellLine<2>::key:
01216     case shards::ShellLine<3>::key:
01217     case shards::ShellTriangle<3>::key:
01218     case shards::ShellTriangle<6>::key:
01219     case shards::ShellQuadrilateral<4>::key:
01220     case shards::ShellQuadrilateral<8>::key:
01221     case shards::ShellQuadrilateral<9>::key:
01222       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01223                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
01224       break;
01225     default:
01226       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01227                           ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported.");        
01228   }// switch  
01229 
01230   // Temp (F,P) array for the values of nodal basis functions at the reference points
01231   int basisCardinality = HGRAD_Basis -> getCardinality();
01232   FieldContainer<Scalar> basisVals(basisCardinality, numPoints);
01233   
01234   // Initialize physPoints
01235   for(int i = 0; i < physPoints.size(); i++){
01236     physPoints[i] = 0.0;
01237   }
01238   
01239   // handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of refPoints
01240   switch(refPoints.rank()) {
01241     
01242     // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
01243     case 2:
01244       {
01245         // getValues requires rank-2 (P,D) input array, but refPoints cannot be passed directly as argument because they are a user type
01246         FieldContainer<Scalar> tempPoints( refPoints.dimension(0), refPoints.dimension(1) );
01247         // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01248         for(int pt = 0; pt < refPoints.dimension(0); pt++){
01249           for(int dm = 0; dm < refPoints.dimension(1) ; dm++){
01250             tempPoints(pt, dm) = refPoints(pt, dm);
01251           }//dm
01252         }//pt
01253         HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
01254         
01255         // If whichCell = -1, ref pt. set is mapped to all cells, otherwise, the set is mapped to one cell only
01256         int cellLoop = (whichCell == -1) ? numCells : 1 ;
01257         
01258         // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
01259         for(int cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
01260           for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01261             for(int dim = 0; dim < spaceDim; dim++){
01262               for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01263                 
01264                 if(whichCell == -1){
01265                   physPoints(cellOrd, pointOrd, dim) += cellWorkset(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01266                 }
01267                 else{
01268                   physPoints(pointOrd, dim) += cellWorkset(whichCell, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01269                 }
01270               } // bfOrd
01271             }// dim
01272           }// pointOrd
01273         }//cellOrd
01274       }// case 2
01275       break;
01276       
01277     // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.  
01278     case 3:
01279       {
01280         // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
01281         FieldContainer<Scalar> tempPoints( refPoints.dimension(1), refPoints.dimension(2) );
01282         
01283         // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
01284         for(int cellOrd = 0; cellOrd < numCells; cellOrd++) {
01285           
01286           // Copy point set corresponding to this cell oridinal to the temp (P,D) array
01287           for(int pt = 0; pt < refPoints.dimension(1); pt++){
01288             for(int dm = 0; dm < refPoints.dimension(2) ; dm++){
01289               tempPoints(pt, dm) = refPoints(cellOrd, pt, dm);
01290             }//dm
01291           }//pt
01292           
01293           // Compute basis values for this set of ref. points
01294           HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
01295           
01296           for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
01297             for(int dim = 0; dim < spaceDim; dim++){
01298               for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
01299                 
01300                 physPoints(cellOrd, pointOrd, dim) += cellWorkset(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
01301                 
01302               } // bfOrd
01303             }// dim
01304           }// pointOrd
01305         }//cellOrd        
01306       }// case 3
01307       break;
01308       
01309     default:
01310       TEUCHOS_TEST_FOR_EXCEPTION( !( (refPoints.rank() == 2) && (refPoints.rank() == 3) ), std::invalid_argument,
01311                              ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): rank 2 or 3 required for refPoints array. ");
01312   }
01313 }
01314 
01315 
01316 
01317 template<class Scalar>
01318 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
01319 void CellTools<Scalar>::mapToReferenceFrame(ArrayRefPoint        &        refPoints,
01320                                             const ArrayPhysPoint &        physPoints,
01321                                             const ArrayCell      &        cellWorkset,
01322                                             const shards::CellTopology &  cellTopo,
01323                                             const int &                   whichCell)
01324 {
01325   INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell) );
01326   
01327   int spaceDim  = (int)cellTopo.getDimension();
01328   int numPoints;
01329   int numCells;
01330   
01331   // Define initial guesses to be  the Cell centers of the reference cell topology
01332   FieldContainer<Scalar> cellCenter(spaceDim);
01333   switch( cellTopo.getKey() ){
01334     // Standard Base topologies (number of cellWorkset = number of vertices)
01335     case shards::Line<2>::key:
01336       cellCenter(0) = 0.0;    break;
01337       
01338     case shards::Triangle<3>::key:
01339     case shards::Triangle<6>::key:    
01340       cellCenter(0) = 1./3.;    cellCenter(1) = 1./3.;  break;
01341       
01342     case shards::Quadrilateral<4>::key:
01343     case shards::Quadrilateral<9>::key:
01344       cellCenter(0) = 0.0;      cellCenter(1) = 0.0;    break;
01345       
01346     case shards::Tetrahedron<4>::key:
01347     case shards::Tetrahedron<10>::key:
01348     case shards::Tetrahedron<11>::key:
01349       cellCenter(0) = 1./6.;    cellCenter(1) =  1./6.;    cellCenter(2) =  1./6.;  break;
01350       
01351     case shards::Hexahedron<8>::key:
01352     case shards::Hexahedron<27>::key:
01353       cellCenter(0) = 0.0;      cellCenter(1) =  0.0;       cellCenter(2) =  0.0;   break;
01354       
01355     case shards::Wedge<6>::key:
01356     case shards::Wedge<18>::key:
01357       cellCenter(0) = 1./3.;    cellCenter(1) =  1./3.;     cellCenter(2) = 0.0;    break;
01358 
01359     case shards::Pyramid<5>::key:
01360       cellCenter(0) = 0.;       cellCenter(1) = 0.;         cellCenter(2) = 0.25;    break;
01361       
01362       // These extended topologies are not used for mapping purposes
01363     case shards::Quadrilateral<8>::key:
01364     case shards::Hexahedron<20>::key:
01365     case shards::Wedge<15>::key:
01366       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01367                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
01368       break;
01369       
01370       // Base and Extended Line, Beam and Shell topologies  
01371     case shards::Line<3>::key:
01372     case shards::Beam<2>::key:
01373     case shards::Beam<3>::key:
01374     case shards::ShellLine<2>::key:
01375     case shards::ShellLine<3>::key:
01376     case shards::ShellTriangle<3>::key:
01377     case shards::ShellTriangle<6>::key:
01378     case shards::ShellQuadrilateral<4>::key:
01379     case shards::ShellQuadrilateral<8>::key:
01380     case shards::ShellQuadrilateral<9>::key:
01381       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01382                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
01383       break;
01384     default:
01385       TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, 
01386                           ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported.");        
01387   }// switch key 
01388   
01389   // Resize initial guess depending on the rank of the physical points array
01390   FieldContainer<Scalar> initGuess;
01391   
01392   // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) initial guess.
01393   if(whichCell == -1){
01394     numPoints = physPoints.dimension(1);
01395     numCells = cellWorkset.dimension(0);
01396     initGuess.resize(numCells, numPoints, spaceDim);
01397     // Set initial guess:
01398     for(int c = 0; c < numCells; c++){
01399       for(int p = 0; p < numPoints; p++){
01400         for(int d = 0; d < spaceDim; d++){
01401           initGuess(c, p, d) = cellCenter(d);
01402         }// d
01403       }// p
01404     }// c
01405   }
01406   // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) initial guess.
01407   else {
01408     numPoints = physPoints.dimension(0);
01409     initGuess.resize(numPoints, spaceDim);
01410     // Set initial guess:
01411     for(int p = 0; p < numPoints; p++){
01412       for(int d = 0; d < spaceDim; d++){
01413         initGuess(p, d) = cellCenter(d);
01414       }// d
01415     }// p
01416   }
01417   
01418   // Call method with initial guess
01419   mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell);  
01420 }
01421   
01422   
01423 
01424 template<class Scalar>
01425 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
01426 void CellTools<Scalar>::mapToReferenceFrameInitGuess(ArrayRefPoint        &        refPoints,
01427                                                      const ArrayInitGuess &        initGuess,
01428                                                      const ArrayPhysPoint &        physPoints,
01429                                                      const ArrayCell      &        cellWorkset,
01430                                                      const shards::CellTopology &  cellTopo,
01431                                                      const int &                   whichCell)
01432 {
01433   INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell) );
01434 
01435   int spaceDim  = (int)cellTopo.getDimension();
01436   int numPoints;
01437   int numCells=0;
01438   
01439   // Temp arrays for Newton iterates and Jacobians. Resize according to rank of ref. point array
01440   FieldContainer<Scalar> xOld;
01441   FieldContainer<Scalar> xTem;  
01442   FieldContainer<Scalar> jacobian;
01443   FieldContainer<Scalar> jacobInv;
01444   FieldContainer<Scalar> error; 
01445   FieldContainer<Scalar> cellCenter(spaceDim);
01446   
01447   // 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.
01448   if(whichCell == -1){
01449     numPoints = physPoints.dimension(1);
01450     numCells = cellWorkset.dimension(0);
01451     xOld.resize(numCells, numPoints, spaceDim);
01452     xTem.resize(numCells, numPoints, spaceDim);  
01453     jacobian.resize(numCells,numPoints, spaceDim, spaceDim);
01454     jacobInv.resize(numCells,numPoints, spaceDim, spaceDim);
01455     error.resize(numCells,numPoints); 
01456     // Set initial guess to xOld
01457     for(int c = 0; c < numCells; c++){
01458       for(int p = 0; p < numPoints; p++){
01459         for(int d = 0; d < spaceDim; d++){
01460           xOld(c, p, d) = initGuess(c, p, d);
01461         }// d
01462       }// p
01463     }// c
01464   }
01465   // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) temp arrays and (P,D,D) Jacobians.
01466   else {
01467     numPoints = physPoints.dimension(0);
01468     xOld.resize(numPoints, spaceDim);
01469     xTem.resize(numPoints, spaceDim);  
01470     jacobian.resize(numPoints, spaceDim, spaceDim);
01471     jacobInv.resize(numPoints, spaceDim, spaceDim);
01472     error.resize(numPoints); 
01473     // Set initial guess to xOld
01474     for(int p = 0; p < numPoints; p++){
01475       for(int d = 0; d < spaceDim; d++){
01476         xOld(p, d) = initGuess(p, d);
01477       }// d
01478     }// p
01479   }
01480   
01481   
01482   // Newton method to solve the equation F(refPoints) - physPoints = 0:
01483   // refPoints = xOld - DF^{-1}(xOld)*(F(xOld) - physPoints) = xOld + DF^{-1}(xOld)*(physPoints - F(xOld))
01484   for(int iter = 0; iter < INTREPID_MAX_NEWTON; ++iter) {
01485     
01486     // Jacobians at the old iterates and their inverses. 
01487     setJacobian(jacobian, xOld, cellWorkset, cellTopo, whichCell);
01488     setJacobianInv(jacobInv, jacobian);
01489         
01490     // The Newton step.
01491     mapToPhysicalFrame( xTem, xOld, cellWorkset, cellTopo, whichCell );      // xTem <- F(xOld)
01492     RealSpaceTools<Scalar>::subtract( xTem, physPoints, xTem );        // xTem <- physPoints - F(xOld)
01493     RealSpaceTools<Scalar>::matvec( refPoints, jacobInv, xTem);        // refPoints <- DF^{-1}( physPoints - F(xOld) )
01494     RealSpaceTools<Scalar>::add( refPoints, xOld );                    // refPoints <- DF^{-1}( physPoints - F(xOld) ) + xOld
01495     
01496     // l2 error (Euclidean distance) between old and new iterates: |xOld - xNew|
01497     RealSpaceTools<Scalar>::subtract( xTem, xOld, refPoints );
01498     RealSpaceTools<Scalar>::vectorNorm( error, xTem, NORM_TWO );
01499     
01500     // Average L2 error for a multiple sets of physical points: error is rank-2 (C,P) array 
01501     double totalError;
01502     if(whichCell == -1) {
01503       FieldContainer<Scalar> cellWiseError(numCells);
01504       // error(C,P) -> cellWiseError(P)
01505       RealSpaceTools<Scalar>::vectorNorm( cellWiseError, error, NORM_ONE );
01506       totalError = RealSpaceTools<Scalar>::vectorNorm( cellWiseError, NORM_ONE );
01507     }
01508     //Average L2 error for a single set of physical points: error is rank-1 (P) array
01509     else{
01510       totalError = RealSpaceTools<Scalar>::vectorNorm( error, NORM_ONE );   
01511       totalError = totalError;
01512     }
01513     
01514     // Stopping criterion:
01515     if (totalError < INTREPID_TOL) {
01516       break;
01517     } 
01518     else if ( iter > INTREPID_MAX_NEWTON) {
01519       INTREPID_VALIDATE(std::cout << " Intrepid::CellTools::mapToReferenceFrameInitGuess failed to converge to desired tolerance within " 
01520                       << INTREPID_MAX_NEWTON  << " iterations\n" );
01521       break;
01522     }
01523     
01524     // initialize next Newton step
01525     xOld = refPoints;
01526   } // for(iter)
01527 }
01528 
01529 
01530 
01531 template<class Scalar>
01532 template<class ArraySubcellPoint, class ArrayParamPoint>
01533 void CellTools<Scalar>::mapToReferenceSubcell(ArraySubcellPoint     &       refSubcellPoints,
01534                                               const ArrayParamPoint &       paramPoints,
01535                                               const int                     subcellDim,
01536                                               const int                     subcellOrd,
01537                                               const shards::CellTopology &  parentCell){
01538   
01539   int cellDim = parentCell.getDimension();
01540   int numPts  = paramPoints.dimension(0);
01541 
01542 #ifdef HAVE_INTREPID_DEBUG
01543   TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument, 
01544                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
01545   
01546   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
01547                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-dimensional subcells.");  
01548   
01549   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
01550                       ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
01551   
01552   // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
01553   std::string errmsg = ">>> ERROR (Intrepid::mapToReferenceSubcell):";
01554   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, refSubcellPoints, 2,2), std::invalid_argument, errmsg);
01555   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, refSubcellPoints, 1, cellDim, cellDim), std::invalid_argument, errmsg);
01556                     
01557   // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
01558   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, paramPoints, 2,2), std::invalid_argument, errmsg);
01559   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, paramPoints, 1, subcellDim, subcellDim), std::invalid_argument, errmsg);    
01560   
01561   // cross check: refSubcellPoints and paramPoints: dimension 0 must match
01562   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refSubcellPoints, 0,  paramPoints, 0), std::invalid_argument, errmsg);      
01563 #endif
01564   
01565   
01566   // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
01567   const FieldContainer<double>& subcellMap = getSubcellParametrization(subcellDim, parentCell);
01568 
01569   // Apply the parametrization map to every point in parameter domain
01570   if(subcellDim == 2) {
01571     for(int pt = 0; pt < numPts; pt++){
01572       double u = paramPoints(pt,0);
01573       double v = paramPoints(pt,1);
01574       
01575       // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
01576       for(int  dim = 0; dim < cellDim; dim++){
01577         refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + \
01578                                     subcellMap(subcellOrd, dim, 1)*u + \
01579                                     subcellMap(subcellOrd, dim, 2)*v;
01580       }
01581     }
01582   }
01583   else if(subcellDim == 1) {    
01584     for(int pt = 0; pt < numPts; pt++){
01585       for(int dim = 0; dim < cellDim; dim++) {
01586         refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + subcellMap(subcellOrd, dim, 1)*paramPoints(pt,0);
01587       }
01588     }
01589   }
01590   else{
01591     TEUCHOS_TEST_FOR_EXCEPTION( !( (subcellDim == 1) || (subcellDim == 2) ), std::invalid_argument, 
01592                         ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
01593   }
01594 }
01595 
01596 
01597 
01598 template<class Scalar>
01599 template<class ArrayEdgeTangent>
01600 void CellTools<Scalar>::getReferenceEdgeTangent(ArrayEdgeTangent &            refEdgeTangent,
01601                                                 const int &                   edgeOrd,
01602                                                 const shards::CellTopology &  parentCell){
01603   
01604   int spaceDim  = parentCell.getDimension();
01605   
01606 #ifdef HAVE_INTREPID_DEBUG
01607   
01608   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01609                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): two or three-dimensional parent cell required");
01610   
01611   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= edgeOrd) && (edgeOrd < (int)parentCell.getSubcellCount(1) ) ), std::invalid_argument,
01612                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): edge ordinal out of bounds");  
01613   
01614   TEUCHOS_TEST_FOR_EXCEPTION( !( refEdgeTangent.size() == spaceDim ), std::invalid_argument,
01615                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): output array size is required to match space dimension");  
01616 #endif
01617   
01618   // Edge parametrizations are computed in setSubcellParametrization and stored in rank-3 array 
01619   // (subcOrd, coordinate, coefficient)
01620   const FieldContainer<double>& edgeMap = getSubcellParametrization(1, parentCell);
01621   
01622   // All ref. edge maps have affine coordinate functions: f_dim(u) = C_0(dim) + C_1(dim)*u, 
01623   //                                     => edge Tangent: -> C_1(*)
01624   refEdgeTangent(0) = edgeMap(edgeOrd, 0, 1);
01625   refEdgeTangent(1) = edgeMap(edgeOrd, 1, 1);
01626   
01627   // Skip last coordinate for 2D parent cells
01628   if(spaceDim == 3) {
01629     refEdgeTangent(2) = edgeMap(edgeOrd, 2, 1);  
01630   }
01631 }
01632 
01633 
01634 
01635 template<class Scalar>
01636 template<class ArrayFaceTangentU, class ArrayFaceTangentV>
01637 void CellTools<Scalar>::getReferenceFaceTangents(ArrayFaceTangentU &           uTan,
01638                                                  ArrayFaceTangentV &           vTan,
01639                                                  const int &                   faceOrd,
01640                                                  const shards::CellTopology &  parentCell){
01641   
01642 #ifdef HAVE_INTREPID_DEBUG
01643   int spaceDim  = parentCell.getDimension();
01644   TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument, 
01645                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): three-dimensional parent cell required");  
01646   
01647   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
01648                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): face ordinal out of bounds");  
01649 
01650   TEUCHOS_TEST_FOR_EXCEPTION( !( (uTan.rank() == 1)  && (vTan.rank() == 1) ), std::invalid_argument,  
01651                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): rank = 1 required for output arrays"); 
01652   
01653   TEUCHOS_TEST_FOR_EXCEPTION( !( uTan.dimension(0) == spaceDim ), std::invalid_argument,
01654                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");  
01655 
01656   TEUCHOS_TEST_FOR_EXCEPTION( !( vTan.dimension(0) == spaceDim ), std::invalid_argument,
01657                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");  
01658 #endif
01659   
01660   // Face parametrizations are computed in setSubcellParametrization and stored in rank-3 array 
01661   // (subcOrd, coordinate, coefficient): retrieve this array
01662   const FieldContainer<double>& faceMap = getSubcellParametrization(2, parentCell);
01663   
01664   /*  All ref. face maps have affine coordinate functions:  f_dim(u,v) = C_0(dim) + C_1(dim)*u + C_2(dim)*v
01665    *                           `   => Tangent vectors are:  uTan -> C_1(*);    vTan -> C_2(*)
01666    */
01667     // set uTan -> C_1(*)
01668     uTan(0) = faceMap(faceOrd, 0, 1);
01669     uTan(1) = faceMap(faceOrd, 1, 1);
01670     uTan(2) = faceMap(faceOrd, 2, 1);
01671     
01672      // set vTan -> C_2(*)
01673     vTan(0) = faceMap(faceOrd, 0, 2);
01674     vTan(1) = faceMap(faceOrd, 1, 2);
01675     vTan(2) = faceMap(faceOrd, 2, 2);
01676 }
01677 
01678 
01679 
01680 template<class Scalar>
01681 template<class ArraySideNormal>
01682 void CellTools<Scalar>::getReferenceSideNormal(ArraySideNormal &             refSideNormal,
01683                                                const int &                   sideOrd,
01684                                                const shards::CellTopology &  parentCell){
01685   int spaceDim  = parentCell.getDimension();
01686 #ifdef HAVE_INTREPID_DEBUG
01687   
01688   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01689                       ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): two or three-dimensional parent cell required");
01690   
01691   // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
01692   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= sideOrd) && (sideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
01693                       ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): side ordinal out of bounds");    
01694 #endif  
01695   
01696   if(spaceDim == 2){
01697     
01698     // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
01699     getReferenceEdgeTangent(refSideNormal, sideOrd, parentCell);
01700     
01701     // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
01702     Scalar temp = refSideNormal(0);
01703     refSideNormal(0) = refSideNormal(1);
01704     refSideNormal(1) = -temp;
01705   }
01706   else{
01707     // 3D parent cell: side = 2D subcell (face), call the face normal method.
01708     getReferenceFaceNormal(refSideNormal, sideOrd, parentCell);
01709   }
01710 }
01711   
01712 
01713 
01714 template<class Scalar>
01715 template<class ArrayFaceNormal>
01716 void CellTools<Scalar>::getReferenceFaceNormal(ArrayFaceNormal &             refFaceNormal,
01717                                                const int &                   faceOrd,
01718                                                const shards::CellTopology &  parentCell){
01719   int spaceDim  = parentCell.getDimension();
01720   
01721 #ifdef HAVE_INTREPID_DEBUG
01722   
01723   TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument, 
01724                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): three-dimensional parent cell required");  
01725   
01726   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
01727                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): face ordinal out of bounds");  
01728   
01729   TEUCHOS_TEST_FOR_EXCEPTION( !( refFaceNormal.rank() == 1 ), std::invalid_argument,  
01730                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): rank = 1 required for output array"); 
01731     
01732   TEUCHOS_TEST_FOR_EXCEPTION( !( refFaceNormal.dimension(0) == spaceDim ), std::invalid_argument,
01733                       ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): dim0 (spatial dim) must match that of parent cell");  
01734 #endif
01735 
01736   // Reference face normal = vector product of reference face tangents. Allocate temp FC storage:
01737   FieldContainer<Scalar> uTan(spaceDim);
01738   FieldContainer<Scalar> vTan(spaceDim);
01739   getReferenceFaceTangents(uTan, vTan, faceOrd, parentCell);
01740   
01741   // Compute the vector product of the reference face tangents:
01742   RealSpaceTools<Scalar>::vecprod(refFaceNormal, uTan, vTan);
01743 }
01744 
01745 
01746 
01747 template<class Scalar>
01748 template<class ArrayEdgeTangent, class ArrayJac>
01749 void CellTools<Scalar>::getPhysicalEdgeTangents(ArrayEdgeTangent &            edgeTangents,
01750                                                 const ArrayJac &              worksetJacobians,
01751                                                 const int &                   worksetEdgeOrd,
01752                                                 const shards::CellTopology &  parentCell){
01753   int worksetSize = worksetJacobians.dimension(0);
01754   int edgePtCount = worksetJacobians.dimension(1); 
01755   int pCellDim    = parentCell.getDimension();
01756   
01757 #ifdef HAVE_INTREPID_DEBUG
01758   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents):";
01759   
01760   TEUCHOS_TEST_FOR_EXCEPTION( !( (pCellDim == 3) || (pCellDim == 2) ), std::invalid_argument, 
01761                       ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required");  
01762   
01763   // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
01764   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, edgeTangents, 3,3), std::invalid_argument, errmsg);
01765   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, edgeTangents, 2, 2,3), std::invalid_argument, errmsg);
01766  
01767   // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
01768   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01769   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 2,3), std::invalid_argument, errmsg);
01770   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 2,3), std::invalid_argument, errmsg);
01771   
01772   // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
01773   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, edgeTangents, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);      
01774   
01775 #endif
01776   
01777   // Temp storage for constant reference edge tangent: rank-1 (D) arrays
01778   FieldContainer<double> refEdgeTan(pCellDim);
01779   getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
01780   
01781   // Loop over workset faces and edge points
01782   for(int pCell = 0; pCell < worksetSize; pCell++){
01783     for(int pt = 0; pt < edgePtCount; pt++){
01784       
01785       // Apply parent cell Jacobian to ref. edge tangent
01786       for(int i = 0; i < pCellDim; i++){
01787         edgeTangents(pCell, pt, i) = 0.0;
01788         for(int j = 0; j < pCellDim; j++){
01789           edgeTangents(pCell, pt, i) +=  worksetJacobians(pCell, pt, i, j)*refEdgeTan(j);
01790         }// for j
01791       }// for i
01792     }// for pt
01793   }// for pCell
01794 }
01795 
01796 
01797 
01798 template<class Scalar>
01799 template<class ArrayFaceTangentU, class ArrayFaceTangentV, class ArrayJac>
01800 void CellTools<Scalar>::getPhysicalFaceTangents(ArrayFaceTangentU &           faceTanU,
01801                                                 ArrayFaceTangentV &           faceTanV,
01802                                                 const ArrayJac &              worksetJacobians,
01803                                                 const int &                   worksetFaceOrd,
01804                                                 const shards::CellTopology &  parentCell){
01805   int worksetSize = worksetJacobians.dimension(0);
01806   int facePtCount = worksetJacobians.dimension(1); 
01807   int pCellDim    = parentCell.getDimension();
01808   
01809 #ifdef HAVE_INTREPID_DEBUG
01810   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents):";
01811 
01812   TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument, 
01813                       ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");  
01814   
01815   // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
01816   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanU, 3,3), std::invalid_argument, errmsg);
01817   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanU, 2, 3,3), std::invalid_argument, errmsg);
01818 
01819   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanV, 3,3), std::invalid_argument, errmsg);
01820   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanV, 2, 3,3), std::invalid_argument, errmsg);
01821 
01822   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU,  faceTanV), std::invalid_argument, errmsg);      
01823 
01824   // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
01825   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01826   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
01827   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
01828 
01829   // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
01830   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);      
01831 
01832 #endif
01833     
01834   // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
01835   FieldContainer<double> refFaceTanU(pCellDim);
01836   FieldContainer<double> refFaceTanV(pCellDim);
01837   getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
01838 
01839   // Loop over workset faces and face points
01840   for(int pCell = 0; pCell < worksetSize; pCell++){
01841     for(int pt = 0; pt < facePtCount; pt++){
01842       
01843       // Apply parent cell Jacobian to ref. face tangents
01844       for(int dim = 0; dim < pCellDim; dim++){
01845         faceTanU(pCell, pt, dim) = 0.0;
01846         faceTanV(pCell, pt, dim) = 0.0;
01847         
01848         // Unroll loops: parent cell dimension can only be 3
01849         faceTanU(pCell, pt, dim) = \
01850           worksetJacobians(pCell, pt, dim, 0)*refFaceTanU(0) + \
01851           worksetJacobians(pCell, pt, dim, 1)*refFaceTanU(1) + \
01852           worksetJacobians(pCell, pt, dim, 2)*refFaceTanU(2);
01853         faceTanV(pCell, pt, dim) = \
01854           worksetJacobians(pCell, pt, dim, 0)*refFaceTanV(0) + \
01855           worksetJacobians(pCell, pt, dim, 1)*refFaceTanV(1) + \
01856           worksetJacobians(pCell, pt, dim, 2)*refFaceTanV(2);
01857       }// for dim
01858     }// for pt
01859   }// for pCell
01860 }
01861 
01862 
01863 template<class Scalar>
01864 template<class ArraySideNormal, class ArrayJac>
01865 void CellTools<Scalar>::getPhysicalSideNormals(ArraySideNormal &             sideNormals,
01866                                                const ArrayJac &              worksetJacobians,
01867                                                const int &                   worksetSideOrd,
01868                                                const shards::CellTopology &  parentCell){
01869   int worksetSize = worksetJacobians.dimension(0);
01870   int sidePtCount = worksetJacobians.dimension(1);   
01871   int spaceDim  = parentCell.getDimension();
01872   
01873 #ifdef HAVE_INTREPID_DEBUG
01874   TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument, 
01875                       ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
01876   
01877   // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
01878   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= worksetSideOrd) && (worksetSideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
01879                       ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): side ordinal out of bounds");  
01880 #endif  
01881   
01882   if(spaceDim == 2){
01883 
01884     // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
01885     getPhysicalEdgeTangents(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
01886     
01887     // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
01888     for(int cell = 0; cell < worksetSize; cell++){
01889       for(int pt = 0; pt < sidePtCount; pt++){
01890         Scalar temp = sideNormals(cell, pt, 0);
01891         sideNormals(cell, pt, 0) = sideNormals(cell, pt, 1);
01892         sideNormals(cell, pt, 1) = -temp;
01893       }// for pt
01894     }// for cell
01895   }
01896   else{
01897     // 3D parent cell: side = 2D subcell (face), call the face normal method.
01898     getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
01899   }
01900 }
01901   
01902   
01903   
01904 template<class Scalar>
01905 template<class ArrayFaceNormal, class ArrayJac>
01906 void CellTools<Scalar>::getPhysicalFaceNormals(ArrayFaceNormal &             faceNormals,
01907                                                const ArrayJac &              worksetJacobians,
01908                                                const int &                   worksetFaceOrd,
01909                                                const shards::CellTopology &  parentCell){
01910   int worksetSize = worksetJacobians.dimension(0);
01911   int facePtCount = worksetJacobians.dimension(1); 
01912   int pCellDim    = parentCell.getDimension();
01913   
01914 #ifdef HAVE_INTREPID_DEBUG
01915   std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals):";
01916   
01917   TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument, 
01918                       ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required");  
01919   
01920   // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
01921   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceNormals, 3,3), std::invalid_argument, errmsg);
01922   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceNormals, 2, 3,3), std::invalid_argument, errmsg);
01923   
01924   // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
01925   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
01926   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
01927   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
01928   
01929   // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
01930   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceNormals, 0,1,2,2,  worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);        
01931 #endif
01932   
01933   // Temp storage for physical face tangents: rank-3 (C,P,D) arrays
01934   FieldContainer<Scalar> faceTanU(worksetSize, facePtCount, pCellDim);
01935   FieldContainer<Scalar> faceTanV(worksetSize, facePtCount, pCellDim);
01936   getPhysicalFaceTangents(faceTanU, faceTanV, worksetJacobians, worksetFaceOrd, parentCell);
01937   
01938   // Compute the vector product of the physical face tangents:
01939   RealSpaceTools<Scalar>::vecprod(faceNormals, faceTanU, faceTanV);
01940   
01941   
01942 }
01943 
01944 //============================================================================================//
01945 //                                                                                            //
01946 //                                        Inclusion tests                                     //
01947 //                                                                                            //
01948 //============================================================================================//
01949 
01950 
01951 template<class Scalar>
01952 int CellTools<Scalar>::checkPointInclusion(const Scalar*                 point,
01953                                            const int                     pointDim,
01954                                            const shards::CellTopology &  cellTopo,
01955                                            const double &                threshold) {
01956 #ifdef HAVE_INTREPID_DEBUG
01957   TEUCHOS_TEST_FOR_EXCEPTION( !(pointDim == (int)cellTopo.getDimension() ), std::invalid_argument,
01958                       ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Point and cell dimensions do not match. ");
01959 #endif
01960   int testResult = 1;
01961   
01962   // Using these values in the tests effectievly inflates the reference element to a larger one
01963   double minus_one = -1.0 - threshold;
01964   double plus_one  =  1.0 + threshold;
01965   double minus_zero = - threshold;
01966   
01967   // A cell with extended topology has the same reference cell as a cell with base topology. 
01968   // => testing for inclusion in a reference Triangle<> and a reference Triangle<6> relies on 
01969   // on the same set of inequalities. To eliminate unnecessary cases we switch on the base topology
01970   unsigned key = cellTopo.getBaseCellTopologyData() -> key ;
01971   switch( key ) {
01972     
01973     case shards::Line<>::key :
01974       if( !(minus_one <= point[0] && point[0] <= plus_one))  testResult = 0;
01975       break;
01976       
01977     case shards::Triangle<>::key : {
01978       Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1.0 );
01979       if( distance > threshold ) testResult = 0;
01980       break;
01981     }
01982       
01983     case shards::Quadrilateral<>::key :
01984       if(!( (minus_one <= point[0] && point[0] <= plus_one) && \
01985             (minus_one <= point[1] && point[1] <= plus_one) ) ) testResult = 0;   
01986       break;
01987       
01988     case shards::Tetrahedron<>::key : {
01989       Scalar distance = std::max(  std::max(-point[0],-point[1]), \
01990                                    std::max(-point[2], point[0] + point[1] + point[2] - 1)  );
01991       if( distance > threshold ) testResult = 0;
01992       break;
01993     }
01994       
01995     case shards::Hexahedron<>::key :
01996       if(!((minus_one <= point[0] && point[0] <= plus_one) && \
01997            (minus_one <= point[1] && point[1] <= plus_one) && \
01998            (minus_one <= point[2] && point[2] <= plus_one)))  \
01999              testResult = 0;
02000       break;
02001       
02002     // The base of the reference prism is the same as the reference triangle => apply triangle test
02003     // to X and Y coordinates and test whether Z is in [-1,1]
02004     case shards::Wedge<>::key : {
02005       Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1 );
02006       if( distance > threshold  || \
02007           point[2] < minus_one || point[2] > plus_one) \
02008             testResult = 0;
02009       break;
02010     }
02011 
02012     // The base of the reference pyramid is the same as the reference quad cell => a horizontal plane
02013     // through a point P(x,y,z) intersects the pyramid at a quadrilateral that equals the base quad 
02014     // scaled by (1-z) => P(x,y,z) is inside the pyramid <=> (x,y) is in [-1+z,1-z]^2 && 0 <= Z <= 1 
02015     case shards::Pyramid<>::key : {
02016       Scalar left  = minus_one + point[2];
02017       Scalar right = plus_one  - point[2];
02018       if(!( (left       <= point[0] && point[0] <= right) && \
02019             (left       <= point[1] && point[1] <= right) && 
02020             (minus_zero <= point[2] && point[2] <= plus_one) ) )  \
02021              testResult = 0;  
02022       break;
02023     }
02024       
02025     default:
02026       TEUCHOS_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) ||
02027                              (key == shards::Triangle<>::key)  ||
02028                              (key == shards::Quadrilateral<>::key) ||
02029                              (key == shards::Tetrahedron<>::key)  ||
02030                              (key == shards::Hexahedron<>::key)  ||
02031                              (key == shards::Wedge<>::key)  ||
02032                              (key == shards::Pyramid<>::key) ),
02033                           std::invalid_argument,
02034                           ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Invalid cell topology. ");
02035   }
02036   return testResult;
02037 }
02038 
02039 
02040 
02041 template<class Scalar>
02042 template<class ArrayPoint>
02043 int CellTools<Scalar>::checkPointsetInclusion(const ArrayPoint&             points,
02044                                               const shards::CellTopology &  cellTopo, 
02045                                               const double &                threshold) {
02046   
02047   int rank = points.rank();  
02048   
02049 #ifdef HAVE_INTREPID_DEBUG
02050   TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= points.rank() ) && (points.rank() <= 3) ), std::invalid_argument,
02051                       ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): rank-1, 2 or 3 required for input points array. ");
02052 
02053   // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
02054   TEUCHOS_TEST_FOR_EXCEPTION( !( points.dimension(rank - 1) == (int)cellTopo.getDimension() ), std::invalid_argument,
02055                       ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): Point and cell dimensions do not match. ");
02056 #endif
02057   
02058   // create temp output array depending on the rank of the input array 
02059   FieldContainer<int> inRefCell;
02060   switch(rank) {
02061     case 1: inRefCell.resize(1); break;
02062     case 2: inRefCell.resize( points.dimension(0) ); break;
02063     case 3: inRefCell.resize( points.dimension(0), points.dimension(1) ); break;
02064   }
02065 
02066   // Call the inclusion method which returns inclusion results for all points
02067   checkPointwiseInclusion(inRefCell, points, cellTopo, threshold);
02068   
02069   // Check if any points were outside, break when finding the first one
02070   int allInside = 1;
02071   for(int i = 0; i < inRefCell.size(); i++ ){
02072     if (inRefCell[i] == 0) {
02073       allInside = 0;
02074       break;
02075     }
02076   }
02077    return allInside;
02078 }
02079 
02080 
02081 
02082 template<class Scalar>
02083 template<class ArrayIncl, class ArrayPoint>
02084 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl &                   inRefCell,
02085                                                 const ArrayPoint &            points,
02086                                                 const shards::CellTopology &  cellTopo, 
02087                                                 const double &                threshold) {
02088   int apRank   = points.rank();
02089   
02090 #ifdef HAVE_INTREPID_DEBUG
02091   
02092   // Verify that points and inRefCell have correct ranks and dimensions
02093   std::string errmsg = ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion):";
02094   if(points.rank() == 1) {
02095     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 1 ), std::invalid_argument, 
02096                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires rank-1 output array.");  
02097     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.dimension(0) == 1), std::invalid_argument,
02098                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires dim0 = 1 for output array.");  
02099   }
02100   else if(points.rank() == 2){
02101     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 1 ), std::invalid_argument, 
02102                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-2 input array requires rank-1 output array.");  
02103     // dimension 0 of the arrays must match
02104     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0,  points, 0), std::invalid_argument, errmsg);
02105   }
02106   else if (points.rank() == 3) {
02107     TEUCHOS_TEST_FOR_EXCEPTION( !(inRefCell.rank() == 2 ), std::invalid_argument, 
02108                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-3 input array requires rank-2 output array.");  
02109     // dimensions 0 and 1 of the arrays must match
02110     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0,1,  points, 0,1), std::invalid_argument, errmsg);
02111   }
02112   else{
02113     TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 1) || (points.rank() == 2) || (points.rank() == 3) ), std::invalid_argument,
02114                         ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");      
02115   }    
02116   
02117   // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
02118   TEUCHOS_TEST_FOR_EXCEPTION( !( points.dimension(apRank - 1) == (int)cellTopo.getDimension() ), std::invalid_argument,
02119                       ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Point and cell dimensions do not match. ");
02120   
02121 #endif
02122   
02123   // Initializations
02124   int dim0     = 1;
02125   int dim1     = 1;
02126   int pointDim = 0;
02127   switch(apRank) {
02128     case 1:
02129       pointDim = points.dimension(0);
02130       break;
02131     case 2:
02132       dim1     = points.dimension(0);
02133       pointDim = points.dimension(1);
02134       break;
02135     case 3:
02136       dim0     = points.dimension(0);
02137       dim1     = points.dimension(1);
02138       pointDim = points.dimension(2);
02139       break;
02140     default:
02141       TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= points.rank() ) && (points.rank() <= 3) ), std::invalid_argument,
02142                           ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");      
02143   }// switch
02144   
02145   
02146   // This method can handle up to rank-3 input arrays. The spatial dim must be the last dimension. 
02147   // The method uses [] accessor because array rank is determined at runtime and the appropriate
02148   // (i,j,..,k) accessor is not known. Use of [] requires the following offsets:
02149   //    for input array  = i0*dim1*pointDim + i1*dim1  (computed in 2 pieces: inPtr0 and inPtr1, resp)
02150   //    for output array = i0*dim1                     (computed in one piece: outPtr0)
02151   int inPtr0  = 0;
02152   int inPtr1  = 0;
02153   int outPtr0 = 0;
02154   Scalar point[3] = {0.0, 0.0, 0.0};
02155   
02156   for(int i0 = 0; i0 < dim0; i0++){
02157     outPtr0 = i0*dim1;
02158     inPtr0  = outPtr0*pointDim;
02159     
02160     for(int i1 = 0; i1 < dim1; i1++) {
02161       inPtr1 = inPtr0 + i1*pointDim;      
02162       point[0] = points[inPtr1];
02163       if(pointDim > 1) {
02164         point[1] = points[inPtr1 + 1];
02165         if(pointDim > 2) {
02166           point[2] = points[inPtr1 + 2];
02167           if(pointDim > 3) {
02168             TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= pointDim) && (pointDim <= 3)), std::invalid_argument, 
02169                                 ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Input array specifies invalid point dimension ");      
02170           }
02171         }
02172       } //if(pointDim > 1)
02173       inRefCell[outPtr0 + i1] = checkPointInclusion(point, pointDim, cellTopo, threshold);
02174     } // for (i1)
02175   } // for(i2)
02176 
02177 }  
02178 
02179 
02180 template<class Scalar>
02181 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
02182 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl &                   inCell,
02183                                                 const ArrayPoint &            points,
02184                                                 const ArrayCell &             cellWorkset,
02185                                                 const shards::CellTopology &  cell,
02186                                                 const int &                   whichCell, 
02187                                                 const double &                threshold)
02188 {
02189   INTREPID_VALIDATE( validateArguments_checkPointwiseInclusion(inCell, points, cellWorkset, whichCell, cell) );
02190   
02191   // For cell topologies with reference cells this test maps the points back to the reference cell
02192   // and uses the method for reference cells
02193   unsigned baseKey = cell.getBaseCellTopologyData() -> key;
02194   
02195   switch(baseKey){
02196     
02197     case shards::Line<>::key :
02198     case shards::Triangle<>::key:
02199     case shards::Quadrilateral<>::key :
02200     case shards::Tetrahedron<>::key :
02201     case shards::Hexahedron<>::key :
02202     case shards::Wedge<>::key :
02203     case shards::Pyramid<>::key :
02204       {
02205         FieldContainer<Scalar> refPoints;
02206         
02207         if(points.rank() == 2){
02208           refPoints.resize(points.dimension(0), points.dimension(1) );
02209           mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
02210           checkPointwiseInclusion(inCell, refPoints, cell, threshold );
02211         }
02212         else if(points.rank() == 3){
02213           refPoints.resize(points.dimension(0), points.dimension(1), points.dimension(2) );
02214           mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
02215           checkPointwiseInclusion(inCell, refPoints, cell, threshold );          
02216         }
02217         break;
02218       }
02219     default: 
02220       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
02221                           ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): cell topology not supported");
02222   }// switch
02223   
02224 }
02225 
02226 
02227 //============================================================================================//
02228 //                                                                                            //
02229 //                  Validation of input/output arguments for CellTools methods                //
02230 //                                                                                            //
02231 //============================================================================================//
02232 
02233 template<class Scalar>
02234 template<class ArrayJac, class ArrayPoint, class ArrayCell>
02235 void CellTools<Scalar>::validateArguments_setJacobian(const ArrayJac    &          jacobian,
02236                                                       const ArrayPoint  &          points,
02237                                                       const ArrayCell   &          cellWorkset,
02238                                                       const int &                  whichCell,
02239                                                       const shards::CellTopology & cellTopo){
02240   
02241   // Validate cellWorkset array
02242   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02243                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for cellWorkset array");
02244   
02245   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02246                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) >= 1 required for cellWorkset array");
02247   
02248   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02249                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02250   
02251   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02252                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02253     
02254   // validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02255   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02256                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): whichCell = -1 or a valid cell ordinal is required.");
02257   
02258   
02259   // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
02260   // If rank-2: admissible jacobians: rank-3 (P,D,D) or rank-4 (C,P,D,D); admissible whichCell: -1 (default) or cell ordinal.
02261   if(points.rank() == 2) {
02262     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(0) <= 0), std::invalid_argument,
02263                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) >= 1 required for points array ");
02264     
02265     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(1) != (int)cellTopo.getDimension() ), std::invalid_argument,
02266                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of points array does not match cell dimension");
02267     
02268     // Validate the output array for the Jacobian: if whichCell == -1 all Jacobians are computed, rank-4 (C,P,D,D) required
02269     if(whichCell == -1) {
02270       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.rank() != 4), std::invalid_argument, 
02271                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 4 required for jacobian array");
02272       
02273       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != cellWorkset.dimension(0)), std::invalid_argument,
02274                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of cellWorkset array");
02275       
02276       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(0)), std::invalid_argument,
02277                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 0 of points array");
02278 
02279       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(2) != points.dimension(1)), std::invalid_argument,
02280                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 1 of points array");
02281       
02282       TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(2) == jacobian.dimension(3) ), std::invalid_argument,
02283                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
02284       
02285       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(3) ) && (jacobian.dimension(3) < 4) ), std::invalid_argument,
02286                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
02287     }     
02288     // A single cell Jacobian is computed when whichCell != -1 (whichCell has been already validated), rank-3 (P,D,D) required
02289     else {
02290       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.rank() != 3), std::invalid_argument, 
02291                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for jacobian array");
02292       
02293       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != points.dimension(0)), std::invalid_argument,
02294                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) of jacobian array must equal dim 0 of points array");
02295 
02296       TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(1)), std::invalid_argument,
02297                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of jacobian array must equal dim 1 of points array");
02298       
02299       TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(1) == jacobian.dimension(2) ), std::invalid_argument,
02300                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 = dim 2 (same spatial dimensions) required for jacobian array. ");
02301       
02302       TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(1) ) && (jacobian.dimension(1) < 4) ), std::invalid_argument,
02303                           ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 and dim 2 (spatial dimensions) must be between 1 and 3. ");
02304     }
02305   }
02306   // Point array is rank-3 (C,P,D): requires whichCell = -1 and rank-4 (C,P,D,D) jacobians
02307   else if(points.rank() ==3){
02308     std::string errmsg  = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian):";
02309     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02310                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of points array must equal dim 0 of cellWorkset array");
02311 
02312     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(1) <= 0), std::invalid_argument,
02313                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) >= 1 required for points array ");
02314     
02315     TEUCHOS_TEST_FOR_EXCEPTION( (points.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02316                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of points array does not match cell dimension");
02317     
02318     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument,
02319                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): default value whichCell=-1 required for rank-3 input points");
02320     
02321     // rank-4 (C,P,D,D) jacobian required for rank-3 (C,P,D) input points
02322     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, jacobian,  4, 4), std::invalid_argument,errmsg);
02323     
02324     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(0) != points.dimension(0)), std::invalid_argument,
02325                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of points array");
02326     
02327     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(1) != points.dimension(1)), std::invalid_argument,
02328                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 1 of points array");
02329   
02330     TEUCHOS_TEST_FOR_EXCEPTION( (jacobian.dimension(2) != points.dimension(2)), std::invalid_argument,
02331                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 2 of points array");
02332     
02333     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(2) == jacobian.dimension(3) ), std::invalid_argument,
02334                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
02335     
02336     TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(3) ) && (jacobian.dimension(3) < 4) ), std::invalid_argument,
02337                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
02338   }
02339   else {
02340     TEUCHOS_TEST_FOR_EXCEPTION( !( (points.rank() == 2) && (points.rank() ==3) ), std::invalid_argument,
02341                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 2 or 3 required for points array");
02342   }  
02343 }
02344 
02345 
02346 
02347 template<class Scalar>
02348 template<class ArrayJacInv, class ArrayJac>
02349 void CellTools<Scalar>::validateArguments_setJacobianInv(const ArrayJacInv & jacobianInv,
02350                                                          const ArrayJac &    jacobian)
02351 {
02352   // Validate input jacobian array: admissible ranks & dimensions are: 
02353   // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
02354   int jacobRank = jacobian.rank();
02355   TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
02356                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
02357   
02358   // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
02359   TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
02360                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
02361   
02362   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
02363                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
02364   
02365   // Validate output jacobianInv array: must have the same rank and dimensions as the input array.
02366   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv):";
02367 
02368   TEUCHOS_TEST_FOR_EXCEPTION( !(requireRankMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
02369   
02370   TEUCHOS_TEST_FOR_EXCEPTION( !(requireDimensionMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
02371 }
02372 
02373 
02374 
02375 template<class Scalar>
02376 template<class ArrayJacDet, class ArrayJac>
02377 void CellTools<Scalar>::validateArguments_setJacobianDetArgs(const ArrayJacDet &  jacobianDet,
02378                                                              const ArrayJac    &  jacobian)
02379 {
02380   // Validate input jacobian array: admissible ranks & dimensions are: 
02381   // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
02382   int jacobRank = jacobian.rank();
02383   TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
02384                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
02385   
02386   // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
02387   TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
02388                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
02389   
02390   TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
02391                       ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
02392 
02393   
02394   // Validate output jacobianDet array: must be rank-2 with dimensions (C,P) if jacobian was rank-4:
02395   if(jacobRank == 4){
02396     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.rank() == 2), std::invalid_argument,
02397                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 2 required for jacobianDet if jacobian is rank-4. ");
02398     
02399     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(0) == jacobian.dimension(0) ), std::invalid_argument,
02400                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of cells) of jacobianDet array must equal dim 0 of jacobian array. ");
02401     
02402     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(1) == jacobian.dimension(1) ), std::invalid_argument,
02403                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 1 (number of points) of jacobianDet array must equal dim 1 of jacobian array.");  
02404   }
02405   
02406   // must be rank-1 with dimension (P) if jacobian was rank-3
02407   else {
02408     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.rank() == 1), std::invalid_argument,
02409                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 1 required for jacobianDet if jacobian is rank-3. ");
02410     
02411     TEUCHOS_TEST_FOR_EXCEPTION( !(jacobianDet.dimension(0) == jacobian.dimension(0) ), std::invalid_argument,
02412                         ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of points) of jacobianDet array must equal dim 0 of jacobian array.");  
02413   }
02414 }
02415 
02416 
02417 
02418 template<class Scalar>
02419 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
02420 void CellTools<Scalar>::validateArguments_mapToPhysicalFrame(const ArrayPhysPoint &        physPoints,
02421                                                              const ArrayRefPoint  &        refPoints,
02422                                                              const ArrayCell      &        cellWorkset,
02423                                                              const shards::CellTopology &  cellTopo,
02424                                                              const int&                    whichCell)
02425 {
02426   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame):";
02427   
02428   // Validate cellWorkset array
02429   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02430                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for cellWorkset array");
02431   
02432   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02433                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
02434   
02435   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02436                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02437   
02438   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02439                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02440   
02441     
02442   // validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02443   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02444                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): whichCell = -1 or a valid cell ordinal is required.");
02445   
02446   // Validate refPoints array: can be rank-2 (P,D) or rank-3 (C,P,D) array
02447   // If rank-2: admissible output array is (P,D) or (C,P,D); admissible whichCell: -1 (default) or cell ordinal
02448   if(refPoints.rank() == 2) {
02449     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(0) <= 0), std::invalid_argument,
02450                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) >= 1 required for refPoints array ");
02451     
02452     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(1) != (int)cellTopo.getDimension() ), std::invalid_argument,
02453                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) of refPoints array does not match cell dimension");
02454 
02455     // Validate output array: whichCell = -1 requires rank-3 array with dimensions (C,P,D)  
02456     if(whichCell == -1) {
02457       TEUCHOS_TEST_FOR_EXCEPTION( ( (physPoints.rank() != 3) && (whichCell == -1) ), std::invalid_argument,
02458                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for physPoints array for the default whichCell value");
02459       
02460       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != cellWorkset.dimension(0)), std::invalid_argument,
02461                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of physPoints array must equal dim 0 of cellWorkset array");
02462       
02463       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != refPoints.dimension(0)), std::invalid_argument,
02464                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) of physPoints array must equal dim 0 of refPoints array"); 
02465       
02466       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(2) != (int)cellTopo.getDimension()), std::invalid_argument,
02467                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) does not match cell dimension ");  
02468     }
02469     // 0 <= whichCell < num cells requires rank-2 (P,D) arrays for both refPoints and physPoints
02470     else{
02471       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.rank() != 2), std::invalid_argument,
02472                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 required for physPoints array");
02473       
02474       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != refPoints.dimension(0)), std::invalid_argument,
02475                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) of physPoints array must equal dim 0 of refPoints array"); 
02476       
02477       TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != (int)cellTopo.getDimension()), std::invalid_argument,
02478                           ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) does not match cell dimension ");      
02479     }
02480   }
02481   // refPoints is (C,P,D): requires physPoints to be (C,P,D) and whichCell=-1  (because all cell mappings are applied)
02482   else if(refPoints.rank() == 3) {
02483     
02484     // 1. validate refPoints dimensions and rank
02485     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02486                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of refPoints and cellWorkset arraya are required to match ");
02487 
02488     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(1) <= 0), std::invalid_argument,
02489                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) >= 1 required for refPoints array ");
02490     
02491     TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02492                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of refPoints array does not match cell dimension");
02493     
02494     // 2. whichCell  must be -1
02495     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument, 
02496                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): default value is required for rank-3 refPoints array");
02497 
02498     // 3.  physPoints must match rank and dimensions of refPoints
02499     TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg );
02500     TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg);
02501   }
02502   // if rank is not 2 or 3 throw exception
02503   else {
02504     TEUCHOS_TEST_FOR_EXCEPTION( !( (refPoints.rank() == 2) || (refPoints.rank() == 3) ), std::invalid_argument,
02505                         ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 or 3 required for refPoints array");
02506   }
02507 }
02508 
02509 
02510 
02511 template<class Scalar>
02512 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
02513 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint  &        refPoints,
02514                                                               const ArrayPhysPoint &        physPoints,
02515                                                               const ArrayCell      &        cellWorkset,
02516                                                               const shards::CellTopology &  cellTopo,
02517                                                               const int&                    whichCell)
02518 {
02519   std::string errmsg  = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02520   std::string errmsg1 = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02521   
02522   // Validate cellWorkset array
02523   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02524                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): rank = 3 required for cellWorkset array");
02525   
02526   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02527                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
02528   
02529   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cellTopo.getSubcellCount(0) ), std::invalid_argument,
02530                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02531   
02532   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cellTopo.getDimension() ), std::invalid_argument,
02533                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02534     
02535   // Validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
02536   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02537                       ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): whichCell = -1 or a valid cell ordinal is required.");
02538   
02539   // Admissible ranks and dimensions of refPoints and physPoints depend on whichCell value:
02540   // default is to map multiple sets of points to multiple sets of points. (C,P,D) arrays required
02541   int validRank;
02542   if(whichCell == -1) {
02543     validRank = 3;
02544     errmsg1 += " default value of whichCell requires rank-3 arrays:";
02545   }
02546   // whichCell is valid cell ordinal => we map single set of pts to a single set of pts. (P,D) arrays required
02547   else{
02548     errmsg1 += " rank-2 arrays required when whichCell is valid cell ordinal";
02549     validRank = 2;
02550   }
02551   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg1, refPoints,  validRank,validRank), std::invalid_argument, errmsg1);
02552   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg1, physPoints, refPoints),           std::invalid_argument, errmsg1);
02553   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg1, refPoints, physPoints),      std::invalid_argument, errmsg1);
02554 }
02555 
02556 
02557 
02558 template<class Scalar>
02559 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
02560 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint  &        refPoints,
02561                                                               const ArrayInitGuess &        initGuess,
02562                                                               const ArrayPhysPoint &        physPoints,
02563                                                               const ArrayCell      &        cellWorkset,
02564                                                               const shards::CellTopology &  cellTopo,
02565                                                               const int&                    whichCell)
02566 {
02567   // Call the method that validates arguments with the default initial guess selection
02568   validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell);
02569   
02570   // Then check initGuess: its rank and dimensions must match those of physPoints.
02571   std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
02572   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, initGuess, physPoints), std::invalid_argument, errmsg);  
02573 }
02574 
02575 
02576 template<class Scalar>
02577 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
02578 void CellTools<Scalar>::validateArguments_checkPointwiseInclusion(ArrayIncl &                   inCell,
02579                                                                   const ArrayPoint &            physPoints,
02580                                                                   const ArrayCell &             cellWorkset,
02581                                                                   const int &                   whichCell,
02582                                                                   const shards::CellTopology &  cell)
02583 {
02584   // Validate cellWorkset array
02585   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.rank() != 3), std::invalid_argument,
02586                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 3 required for cellWorkset array");
02587   
02588   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(0) <= 0), std::invalid_argument,
02589                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) >= 1 required for cellWorkset array");
02590   
02591   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(1) != (int)cell.getSubcellCount(0) ), std::invalid_argument,
02592                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
02593   
02594   TEUCHOS_TEST_FOR_EXCEPTION( (cellWorkset.dimension(2) != (int)cell.getDimension() ), std::invalid_argument,
02595                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of cellWorkset array  does not match cell dimension");
02596   
02597   
02598   // Validate whichCell It can be either -1 (default value) or a valid cell ordinal.
02599   TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
02600                       ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 or a valid cell ordinal is required.");
02601   
02602   
02603   // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
02604   // If rank-2: admissible inCell is rank-1 (P); admissible whichCell is valid cell ordinal but not -1.
02605   if(physPoints.rank() == 2) {
02606     
02607     TEUCHOS_TEST_FOR_EXCEPTION( (whichCell == -1), std::invalid_argument,
02608                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = a valid cell ordinal is required with rank-2 input array.");
02609 
02610     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) <= 0), std::invalid_argument,
02611                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) >= 1 required for physPoints array ");
02612     
02613     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) != (int)cell.getDimension() ), std::invalid_argument,
02614                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (spatial dimension) of physPoints array does not match cell dimension");
02615     
02616     // Validate inCell
02617     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.rank() != 1), std::invalid_argument, 
02618                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 1 required for inCell array");
02619     
02620     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(0) != physPoints.dimension(0)), std::invalid_argument,
02621                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) of inCell array must equal dim 0 of physPoints array");
02622   }
02623   // If rank-3: admissible inCell is rank-2 (C,P); admissible whichCell = -1.
02624   else if (physPoints.rank() == 3){
02625     
02626     TEUCHOS_TEST_FOR_EXCEPTION( !(whichCell == -1), std::invalid_argument,
02627                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 is required with rank-3 input array.");
02628     
02629     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(0) != cellWorkset.dimension(0) ), std::invalid_argument,
02630                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells)  of physPoints array must equal dim 0 of cellWorkset array ");
02631 
02632     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(1) <= 0), std::invalid_argument,
02633                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) >= 1 required for physPoints array ");
02634     
02635     TEUCHOS_TEST_FOR_EXCEPTION( (physPoints.dimension(2) != (int)cell.getDimension() ), std::invalid_argument,
02636                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of physPoints array does not match cell dimension");
02637     
02638     // Validate inCell
02639     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.rank() != 2), std::invalid_argument, 
02640                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 required for inCell array");
02641     
02642     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(0) != physPoints.dimension(0)), std::invalid_argument,
02643                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) of inCell array must equal dim 0 of physPoints array");    
02644 
02645     TEUCHOS_TEST_FOR_EXCEPTION( (inCell.dimension(1) != physPoints.dimension(1)), std::invalid_argument,
02646                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) of inCell array must equal dim 1 of physPoints array");    
02647   }
02648   else {
02649     TEUCHOS_TEST_FOR_EXCEPTION( !( (physPoints.rank() == 2) && (physPoints.rank() ==3) ), std::invalid_argument,
02650                         ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 or 3 required for points array");
02651   }
02652 }
02653 
02654 
02655 
02656 //============================================================================================//
02657 //                                                                                            //
02658 //                                           Debug                                            //
02659 //                                                                                            //
02660 //============================================================================================//
02661 
02662 
02663 template<class Scalar>
02664 void CellTools<Scalar>::printSubcellVertices(const int subcellDim,
02665                                              const int subcellOrd,
02666                                              const shards::CellTopology & parentCell){
02667   
02668   // Get number of vertices for the specified subcell and parent cell dimension
02669   int subcVertexCount = parentCell.getVertexCount(subcellDim, subcellOrd);
02670   int cellDim         = parentCell.getDimension();
02671   
02672   // Allocate space for the subcell vertex coordinates
02673   FieldContainer<double> subcellVertices(subcVertexCount, cellDim);
02674   
02675   // Retrieve the vertex coordinates
02676   getReferenceSubcellVertices(subcellVertices,
02677                               subcellDim,
02678                               subcellOrd,
02679                               parentCell);
02680   
02681   // Print the vertices
02682   std::cout 
02683     << " Subcell " << std::setw(2) << subcellOrd 
02684     <<  " is " << parentCell.getName(subcellDim, subcellOrd) << " with vertices = {";
02685   
02686   // Loop over subcell vertices
02687   for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
02688     std::cout<< "(";
02689     
02690     // Loop over vertex Cartesian coordinates
02691     for(int dim = 0; dim < (int)parentCell.getDimension(); dim++){
02692       std::cout << subcellVertices(subcVertOrd, dim);
02693       if(dim < (int)parentCell.getDimension()-1 ) { std::cout << ","; }
02694     }
02695     std::cout<< ")";
02696     if(subcVertOrd < subcVertexCount - 1) { std::cout << ", "; }
02697   }
02698   std::cout << "}\n";
02699 }
02700   
02701 
02702 template<class Scalar>
02703 template<class ArrayCell>
02704 void CellTools<Scalar>::printWorksetSubcell(const ArrayCell &             cellWorkset,
02705                                             const shards::CellTopology &  parentCell,
02706                                             const int&                    pCellOrd,
02707                                             const int&                    subcellDim,
02708                                             const int&                    subcellOrd,
02709                                             const int&                    fieldWidth){
02710   
02711   // Get the ordinals, relative to reference cell, of subcell cellWorkset
02712   int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
02713   int pCellDim      = parentCell.getDimension();
02714   std::vector<int> subcNodeOrdinals(subcNodeCount);
02715   
02716   for(int i = 0; i < subcNodeCount; i++){
02717     subcNodeOrdinals[i] = parentCell.getNodeMap(subcellDim, subcellOrd, i);
02718   }
02719   
02720   // Loop over parent cells and print subcell cellWorkset
02721   
02722   std::cout 
02723     << " Subcell " << subcellOrd << " on parent cell " << pCellOrd << " is " 
02724     << parentCell.getName(subcellDim, subcellOrd) << " with node(s) \n ({";
02725   
02726   for(int i = 0; i < subcNodeCount; i++){
02727     
02728     // print Cartesian coordinates of the node
02729     for(int dim = 0; dim < pCellDim; dim++){
02730       std::cout
02731       << std::setw(fieldWidth) << std::right << cellWorkset(pCellOrd, subcNodeOrdinals[i], dim); 
02732       if(dim < pCellDim - 1){ std::cout << ","; }
02733     }
02734     std::cout << "}";
02735     if(i < subcNodeCount - 1){ std::cout <<", {"; }
02736   }
02737   std::cout << ")\n\n";
02738 }
02739 
02740 
02741 
02742 } // namespace Intrepid
02743 #endif