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