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