Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/src/Discretization/Basis/Intrepid_HDIV_HEX_In_FEMDef.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 
00049 namespace Intrepid {
00050 
00051   template<class Scalar, class ArrayScalar>
00052   Basis_HDIV_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_HEX_In_FEM( int order ,
00053                                                                       const ArrayScalar & ptsClosed ,
00054                                                                       const ArrayScalar & ptsOpen):
00055     closedBasis_( order , ptsClosed ),
00056     openBasis_( order-1 , ptsOpen ),
00057     closedPts_( ptsClosed ),
00058     openPts_( ptsOpen )
00059   {
00060     this -> basisDegree_       = order;
00061     this -> basisCardinality_  = 3 * closedBasis_.getCardinality() * openBasis_.getCardinality() * openBasis_.getCardinality();
00062     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00063     this -> basisType_         = BASIS_FEM_FIAT;
00064     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00065     this -> basisTagsAreSet_   = false;
00066 
00067     Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3);
00068     bases[0].resize(3); bases[1].resize(3); bases[2].resize(3);
00069     bases[0][0] = rcp( &closedBasis_ , false );
00070     bases[0][1] = rcp( &openBasis_ , false );
00071     bases[0][2] = rcp( &openBasis_ , false );
00072     bases[1][0] = rcp( &openBasis_ , false );
00073     bases[1][1] = rcp( &closedBasis_ , false );
00074     bases[1][2] = rcp( &openBasis_ , false );
00075     bases[2][0] = rcp( &openBasis_ , false );
00076     bases[2][1] = rcp( &openBasis_ , false );
00077     bases[2][2] = rcp( &closedBasis_ , false );
00078     this->setBases( bases );
00079 
00080   }
00081 
00082   template<class Scalar, class ArrayScalar>
00083   Basis_HDIV_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_HEX_In_FEM( int order , const EPointType &pointType ):
00084     closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
00085     openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
00086     closedPts_( order+1 , 1 ),
00087     openPts_( order , 1 )
00088   {
00089     this -> basisDegree_       = order;
00090     this -> basisCardinality_  = 3 * closedBasis_.getCardinality() * openBasis_.getCardinality() * openBasis_.getCardinality();
00091     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00092     this -> basisType_         = BASIS_FEM_FIAT;
00093     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00094     this -> basisTagsAreSet_   = false;
00095 
00096     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
00097                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00098                                                             order ,
00099                                                             0 ,
00100                                                             pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
00101 
00102     if (pointType == POINTTYPE_SPECTRAL)
00103       {
00104         PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
00105                                                                     order - 1 );
00106       }
00107     else
00108       {
00109         PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
00110                                                                 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00111                                                                 order - 1,
00112                                                                 0 ,
00113                                                                 POINTTYPE_EQUISPACED );
00114 
00115       }
00116 
00117 
00118     Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3);
00119     bases[0].resize(3); bases[1].resize(3); bases[2].resize(3);
00120     bases[0][0] = rcp( &closedBasis_ , false );
00121     bases[0][1] = rcp( &openBasis_ , false );
00122     bases[0][2] = rcp( &openBasis_ , false );
00123     bases[1][0] = rcp( &openBasis_ , false );
00124     bases[1][1] = rcp( &closedBasis_ , false );
00125     bases[1][2] = rcp( &openBasis_ , false );
00126     bases[2][0] = rcp( &openBasis_ , false );
00127     bases[2][1] = rcp( &openBasis_ , false );
00128     bases[2][2] = rcp( &closedBasis_ , false );
00129     this->setBases( bases );
00130 
00131   }
00132   
00133 
00134   
00135   template<class Scalar, class ArrayScalar>
00136   void Basis_HDIV_HEX_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00137     
00138     // Basis-dependent intializations
00139     int tagSize  = 4;        // size of DoF tag
00140     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00141     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00142     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00143     
00144     std::vector<int> tags( tagSize * this->getCardinality() );
00145     
00146     const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
00147     const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
00148 
00149     std::map<int,std::map<int,int> > total_dof_per_entity;
00150     std::map<int,std::map<int,int> > current_dof_per_entity;
00151 
00152     // four vertices
00153     for (int i=0;i<4;i++) {
00154       total_dof_per_entity[0][i] = 0;
00155       current_dof_per_entity[0][i] = 0;
00156     }
00157     // twelve edges
00158     for (int i=0;i<12;i++) {
00159       total_dof_per_entity[1][i] = 0;
00160       current_dof_per_entity[1][i] = 0;
00161     }
00162     // six faces
00163     for (int i=0;i<6;i++) {
00164       total_dof_per_entity[2][i] = 0;
00165       current_dof_per_entity[2][i] = 0;
00166     }
00167     total_dof_per_entity[3][0] = 0;
00168     current_dof_per_entity[3][0] = 0;
00169 
00170     // tally dof on each facet.  none on vertex or edge
00171     for (int i=0;i<6;i++) {
00172       total_dof_per_entity[2][i] = openBasis_.getCardinality() * openBasis_.getCardinality();
00173     }
00174 
00175     total_dof_per_entity[2][0] = this->getCardinality() - 6 * openBasis_.getCardinality()* openBasis_.getCardinality();
00176 
00177     int tagcur = 0;
00178     // loop over the x-component basis functions, which are (psi(x)phi(y)phi(z),0,0)
00179     // for psi in the closed basis and phi in the open
00180     for (int k=0;k<openBasis_.getCardinality();k++) {
00181       const int odimk = openDofTags[k][0];
00182       const int oentk = openDofTags[k][1];
00183       for (int j=0;j<openBasis_.getCardinality();j++) {
00184         const int odimj = openDofTags[j][0];
00185         const int oentj = openDofTags[j][1];
00186         for (int i=0;i<closedBasis_.getCardinality();i++) {
00187           const int cdim = closedDofTags[i][0];
00188           const int cent = closedDofTags[i][1];
00189           int dofdim;
00190           int dofent;
00191           ProductTopology::lineProduct3d(cdim,cent,odimj,oentj,odimk,oentk,dofdim,dofent);
00192           tags[4*tagcur] = dofdim;
00193           tags[4*tagcur+1] = dofent;
00194           tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00195           current_dof_per_entity[dofdim][dofent]++;
00196           tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00197           tagcur++;
00198         }
00199       }
00200     }
00201 
00202     // now we have to do it for the y-component basis functions, which are
00203     // (0,phi(x)psi(y)phi(z),0) for psi in the closed basis and phi in the open
00204     for (int k=0;k<openBasis_.getCardinality();k++) {
00205       const int odimk = openDofTags[k][0];
00206       const int oentk = openDofTags[k][1];
00207       for (int j=0;j<closedBasis_.getCardinality();j++) {
00208         const int cdim = closedDofTags[j][0];
00209         const int cent = closedDofTags[j][1];
00210         for (int i=0;i<openBasis_.getCardinality();i++) {
00211           const int odimi = openDofTags[i][0];
00212           const int oenti = openDofTags[i][1];
00213           int dofdim;
00214           int dofent;
00215           ProductTopology::lineProduct3d(odimi,oenti,cdim,cent,odimk,oentk,dofdim,dofent);
00216           tags[4*tagcur] = dofdim;
00217           tags[4*tagcur+1] = dofent;
00218           tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00219           current_dof_per_entity[dofdim][dofent]++;
00220           tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00221           tagcur++;
00222         }
00223       }
00224     }
00225 
00226     // now we have to do it for the z-component basis functions, which are
00227     // (0,0,phi(x)phi(y)psi(z)) for psi in the closed basis and phi in the open
00228     for (int k=0;k<closedBasis_.getCardinality();k++) {
00229       const int cdim = closedDofTags[k][0];
00230       const int cent = closedDofTags[k][1];
00231       for (int j=0;j<openBasis_.getCardinality();j++) {
00232         const int odimj = openDofTags[j][0];
00233         const int oentj = openDofTags[j][1];
00234         for (int i=0;i<openBasis_.getCardinality();i++) {
00235           const int odimi = openDofTags[i][0];
00236           const int oenti = openDofTags[i][1];
00237           int dofdim;
00238           int dofent;
00239           ProductTopology::lineProduct3d(odimi,oenti,odimj,oentj,cdim,cent,dofdim,dofent);
00240           tags[4*tagcur] = dofdim;
00241           tags[4*tagcur+1] = dofent;
00242           tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00243           current_dof_per_entity[dofdim][dofent]++;
00244           tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00245           tagcur++;
00246         }
00247       }
00248     }
00249   
00250 //     for (int i=0;i<this->getCardinality();i++) {
00251 //       for (int j=0;j<4;j++) {
00252 //      std::cout << tags[4*i+j] << " ";
00253 //       }
00254 //       std::cout << std::endl;
00255 //     }
00256   
00257     // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00258     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00259                                 this -> ordinalToTag_,
00260                                 &(tags[0]),
00261                                 this -> basisCardinality_,
00262                                 tagSize,
00263                                 posScDim,
00264                                 posScOrd,
00265                                 posDfOrd);
00266   }
00267 
00268 
00269   template<class Scalar, class ArrayScalar>
00270   void Basis_HDIV_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00271                                                               const ArrayScalar &  inputPoints,
00272                                                               const EOperator      operatorType) const {
00273     
00274     // Verify arguments
00275 #ifdef HAVE_INTREPID_DEBUG
00276     Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
00277                                                        inputPoints,
00278                                                        operatorType,
00279                                                        this -> getBaseCellTopology(),
00280                                                        this -> getCardinality() );
00281 #endif
00282     
00283     // Number of evaluation points = dim 0 of inputPoints
00284     int dim0 = inputPoints.dimension(0);
00285     
00286     // separate out points
00287     FieldContainer<Scalar> xPoints(dim0,1);
00288     FieldContainer<Scalar> yPoints(dim0,1);
00289     FieldContainer<Scalar> zPoints(dim0,1);
00290 
00291     
00292     for (int i=0;i<dim0;i++) {
00293       xPoints(i,0) = inputPoints(i,0);
00294       yPoints(i,0) = inputPoints(i,1);
00295       zPoints(i,0) = inputPoints(i,2);
00296     }
00297     
00298     switch (operatorType) {
00299     case OPERATOR_VALUE:
00300       {
00301         FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
00302         FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
00303         FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 );
00304         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00305         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00306         FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 );
00307         
00308         closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
00309         closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
00310         closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE );
00311         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00312         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00313         openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
00314         
00315         int bfcur = 0;
00316         // x component bfs are (closed(x) open(y) open(z),0,0)
00317         for (int k=0;k<openBasis_.getCardinality();k++) {
00318           for (int j=0;j<openBasis_.getCardinality();j++) {
00319             for (int i=0;i<closedBasis_.getCardinality();i++) {
00320               for (int l=0;l<dim0;l++) {
00321                 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l);
00322                 outputValues(bfcur,l,1) = 0.0;
00323                 outputValues(bfcur,l,2) = 0.0;
00324               }
00325               bfcur++;
00326             }
00327           }
00328         }
00329         
00330         // y component bfs are (0,open(x) closed(y) open(z),0)
00331         for (int k=0;k<openBasis_.getCardinality();k++) {
00332           for (int j=0;j<closedBasis_.getCardinality();j++) {
00333             for (int i=0;i<openBasis_.getCardinality();i++) {
00334               for (int l=0;l<dim0;l++) {
00335                 outputValues(bfcur,l,0) = 0.0;
00336                 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l);
00337                 outputValues(bfcur,l,2) = 0.0;
00338               }
00339               bfcur++;
00340             }
00341           }
00342         }
00343 
00344         // z component bfs are (0,0,open(x) open(y) closed(z))
00345         for (int k=0;k<closedBasis_.getCardinality();k++) {
00346           for (int j=0;j<openBasis_.getCardinality();j++) {
00347             for (int i=0;i<openBasis_.getCardinality();i++) {
00348               for (int l=0;l<dim0;l++) {
00349                 outputValues(bfcur,l,0) = 0.0;
00350                 outputValues(bfcur,l,1) = 0.0;
00351                 outputValues(bfcur,l,2) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l);
00352               }
00353               bfcur++;
00354             }
00355           }
00356         }
00357 
00358 
00359       }
00360       break;
00361     case OPERATOR_DIV:
00362       {
00363         FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
00364         FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
00365         FieldContainer<Scalar> closedBasisDerivsZPts( closedBasis_.getCardinality() , dim0 , 1 );
00366         FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
00367         FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
00368         FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 );
00369         
00370         closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
00371         closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
00372         closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 );
00373         openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
00374         openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
00375         openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
00376         
00377         int bfcur = 0;
00378         
00379         // x component basis functions first
00380         for (int k=0;k<openBasis_.getCardinality();k++) {
00381           for (int j=0;j<openBasis_.getCardinality();j++) {
00382             for (int i=0;i<closedBasis_.getCardinality();i++) {
00383               for (int l=0;l<dim0;l++) {
00384                 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l);
00385               }
00386               bfcur++;
00387             }
00388           }
00389         }
00390         
00391         // now y component basis functions
00392         for (int k=0;k<openBasis_.getCardinality();k++) {
00393           for (int j=0;j<closedBasis_.getCardinality();j++) {
00394             for (int i=0;i<openBasis_.getCardinality();i++) {
00395               for (int l=0;l<dim0;l++) {
00396                 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l);
00397               }
00398               bfcur++;
00399             }
00400           }
00401         }
00402 
00403         // now z component basis functions
00404         for (int k=0;k<closedBasis_.getCardinality();k++) {
00405           for (int j=0;j<openBasis_.getCardinality();j++) {
00406             for (int i=0;i<openBasis_.getCardinality();i++) {
00407               for (int l=0;l<dim0;l++) {
00408                 outputValues(bfcur,l) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0);
00409               }
00410               bfcur++;
00411             }
00412           }
00413         }
00414       }
00415       break;
00416     case OPERATOR_CURL:
00417       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00418                           ">>> ERROR (Basis_HDIV_HEX_In_FEM): CURL is invalid operator for HDIV Basis Functions");
00419       break;
00420       
00421     case OPERATOR_GRAD:
00422       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
00423                           ">>> ERROR (Basis_HDIV_HEX_In_FEM): GRAD is invalid operator for HDIV Basis Functions");
00424       break;
00425       
00426     case OPERATOR_D1:
00427     case OPERATOR_D2:
00428     case OPERATOR_D3:
00429     case OPERATOR_D4:
00430     case OPERATOR_D5:
00431     case OPERATOR_D6:
00432     case OPERATOR_D7:
00433     case OPERATOR_D8:
00434     case OPERATOR_D9:
00435     case OPERATOR_D10:
00436       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1)    ||
00437                             (operatorType == OPERATOR_D2)    ||
00438                             (operatorType == OPERATOR_D3)    ||
00439                             (operatorType == OPERATOR_D4)    ||
00440                             (operatorType == OPERATOR_D5)    ||
00441                             (operatorType == OPERATOR_D6)    ||
00442                             (operatorType == OPERATOR_D7)    ||
00443                             (operatorType == OPERATOR_D8)    ||
00444                             (operatorType == OPERATOR_D9)    ||
00445                             (operatorType == OPERATOR_D10) ),
00446                           std::invalid_argument,
00447                           ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type");
00448       break;
00449       
00450     default:
00451       TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
00452                             (operatorType != OPERATOR_GRAD)  &&
00453                             (operatorType != OPERATOR_CURL)  &&
00454                             (operatorType != OPERATOR_DIV)   &&
00455                             (operatorType != OPERATOR_D1)    &&
00456                             (operatorType != OPERATOR_D2)    &&
00457                             (operatorType != OPERATOR_D3)    &&
00458                             (operatorType != OPERATOR_D4)    &&
00459                             (operatorType != OPERATOR_D5)    &&
00460                             (operatorType != OPERATOR_D6)    &&
00461                             (operatorType != OPERATOR_D7)    &&
00462                             (operatorType != OPERATOR_D8)    &&
00463                             (operatorType != OPERATOR_D9)    &&
00464                             (operatorType != OPERATOR_D10) ),
00465                           std::invalid_argument,
00466                           ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type");
00467     }
00468   }
00469   
00470   
00471   
00472   template<class Scalar, class ArrayScalar>
00473   void Basis_HDIV_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00474                                                               const ArrayScalar &    inputPoints,
00475                                                               const ArrayScalar &    cellVertices,
00476                                                               const EOperator        operatorType) const {
00477     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00478                         ">>> ERROR (Basis_HDIV_HEX_In_FEM): FEM Basis calling an FVD member function");
00479   }
00480 
00481   template<class Scalar, class ArrayScalar>
00482   void Basis_HDIV_HEX_In_FEM<Scalar,ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const
00483   {
00484     // x-component basis functions
00485     int cur = 0;
00486 
00487     for (int k=0;k<openPts_.dimension(0);k++)
00488       {
00489         for (int j=0;j<openPts_.dimension(0);j++)
00490           {
00491             for (int i=0;i<closedPts_.dimension(0);i++)
00492               {
00493                 DofCoords(cur,0) = closedPts_(i,0);
00494                 DofCoords(cur,1) = openPts_(j,0);
00495                 DofCoords(cur,2) = openPts_(k,0);
00496                 cur++;
00497               }
00498           }
00499       }
00500 
00501     // y-component basis functions
00502     for (int k=0;k<openPts_.dimension(0);k++)
00503       {
00504         for (int j=0;j<closedPts_.dimension(0);j++)
00505           {
00506             for (int i=0;i<openPts_.dimension(0);i++)
00507               {
00508                 DofCoords(cur,0) = openPts_(i,0);
00509                 DofCoords(cur,1) = closedPts_(j,0);
00510                 DofCoords(cur,2) = openPts_(k,0);
00511                 cur++;
00512               }
00513           }
00514       }
00515 
00516     // z-component basis functions
00517     for (int k=0;k<closedPts_.dimension(0);k++)
00518       {
00519         for (int j=0;j<openPts_.dimension(0);j++)
00520           {
00521             for (int i=0;i<openPts_.dimension(0);i++)
00522               {
00523                 DofCoords(cur,0) = openPts_(i,0);
00524                 DofCoords(cur,1) = openPts_(j,0);
00525                 DofCoords(cur,2) = closedPts_(k,0);
00526                 cur++;
00527               }
00528           }
00529       }
00530 
00531 
00532 
00533     return;
00534   }
00535   
00536 }// namespace Intrepid