Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_TET_In_FEMDef.hpp
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 
00050 namespace Intrepid {
00051 
00052   template<class Scalar, class ArrayScalar>
00053   Basis_HCURL_TET_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_TET_In_FEM( const int n ,
00054                                                                       const EPointType pointType ):
00055     Phis_( n ),
00056     coeffs_( (n+1)*(n+2)*(n+3)/2 , n*(n+2)*(n+3)/2  )
00057   {
00058      const int N = n*(n+2)*(n+3)/2;
00059      this -> basisCardinality_  = N;
00060      this -> basisDegree_       = n;
00061      this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00062      this -> basisType_         = BASIS_FEM_FIAT;
00063      this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00064      this -> basisTagsAreSet_   = false;
00065 
00066      const int littleN = n*(n+1)*(n+2)/2;    // dim of (P_{n-1})^3 -- smaller space
00067      const int bigN = (n+1)*(n+2)*(n+3)/2;   // dim of (P_{n})^3 -- larger space
00068      const int start_PkH = (n-1)*n*(n+1)/6;  // dim of P({n-2}), offset into 
00069      const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH;
00070      const int scalarLittleN = littleN/3;
00071      const int scalarBigN = bigN/3;
00072 
00073      // first, need to project the basis for Nedelec space onto the
00074      // orthogonal basis of degree n
00075      // get coefficients of PkHx
00076 
00077      Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, littleN + 3 * dim_PkH);
00078 
00079      // these two loops get the first three sets of basis functions
00080      for (int i=0;i<scalarLittleN;i++) {
00081        for (int k=0;k<3;k++) {
00082          V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0;
00083        }
00084      }
00085 
00086      // first 3*scalarLittleN columns are (P_{n-1})^3 space
00087 
00088 
00089      // now I need to integrate { (x,y,z) \times } against the big basis
00090      // first, get a cubature rule.
00091      CubatureDirectTetDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
00092      FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 );
00093      FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
00094      myCub.getCubature( cubPoints , cubWeights );
00095 
00096      // tabulate the scalar orthonormal basis at cubature points
00097      FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
00098      Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
00099 
00100 
00101 
00102      // first set of these functions will write into the first dimPkH columns of remainder
00103 
00104      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00105        // write into second spatial component, where
00106        // I integrate z phi_j phi_i
00107        for (int i=0;i<scalarBigN;i++) {
00108          V1(scalarBigN+i,littleN+j) = 0.0;
00109          for (int k=0;k<myCub.getNumPoints();k++) {
00110            V1(scalarBigN+i,littleN+j) -= cubWeights(k) * cubPoints(k,2) 
00111              * phisAtCubPoints(start_PkH+j,k)
00112              * phisAtCubPoints(i,k);
00113          }
00114        }  
00115        // write into third spatial component (-y phi_j, phi_i)
00116        for (int i=0;i<scalarBigN;i++) {
00117          V1(2*scalarBigN+i,littleN+j) = 0.0;
00118          for (int k=0;k<myCub.getNumPoints();k++) {
00119            V1(2*scalarBigN+i,littleN+j) += cubWeights(k) * cubPoints(k,1) 
00120              * phisAtCubPoints(start_PkH+j,k)
00121              * phisAtCubPoints(i,k);
00122          }
00123        }  
00124      }
00125 
00126      // second set of basis functions, write into second set of dimPkH columns
00127      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00128        // write into first spatial component, where
00129        // I integrate -z phi_j phi_i
00130        for (int i=0;i<scalarBigN;i++) {
00131          V1(i,littleN+dim_PkH+j) = 0.0;
00132          for (int k=0;k<myCub.getNumPoints();k++) {
00133            V1(i,littleN+dim_PkH+j) += cubWeights(k) * cubPoints(k,2) 
00134              * phisAtCubPoints(start_PkH+j,k)
00135              * phisAtCubPoints(i,k);
00136          }
00137        } 
00138      
00139        // third spatial component, x phi_j phi_i
00140        for (int i=0;i<scalarBigN;i++) {
00141          V1(2*scalarBigN+i,littleN+dim_PkH+j) = 0.0;
00142          for (int k=0;k<myCub.getNumPoints();k++) {
00143            V1(2*scalarBigN+i,littleN+dim_PkH+j) -= cubWeights(k) * cubPoints(k,0) 
00144              * phisAtCubPoints(start_PkH+j,k)
00145              * phisAtCubPoints(i,k);
00146          }
00147        }  
00148      }
00149     
00150      // third clump of dimPkH columns
00151      for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
00152        // write into first spatial component, where
00153        // I integrate y phi_j phi_i
00154        for (int i=0;i<scalarBigN;i++) {
00155          V1(i,littleN+2*dim_PkH+j) = 0.0;
00156          for (int k=0;k<myCub.getNumPoints();k++) {
00157            V1(i,littleN+2*dim_PkH+j) -= cubWeights(k) * cubPoints(k,1) 
00158              * phisAtCubPoints(start_PkH+j,k)
00159              * phisAtCubPoints(i,k);
00160          }
00161        }  
00162        // second spatial component, -x phi_j phi_i
00163        for (int i=0;i<scalarBigN;i++) {
00164          V1(scalarBigN+i,littleN+2*dim_PkH+j) = 0.0;
00165          for (int k=0;k<myCub.getNumPoints();k++) {
00166            V1(scalarBigN+i,littleN+2*dim_PkH+j) += cubWeights(k) * cubPoints(k,0) 
00167              * phisAtCubPoints(start_PkH+j,k)
00168              * phisAtCubPoints(i,k);
00169          }
00170        }  
00171      }
00172 
00173      // now I need to set up an SVD to get a basis for the space
00174      Teuchos::SerialDenseMatrix<int,Scalar> S(bigN,1);
00175      Teuchos::SerialDenseMatrix<int,Scalar> U(bigN, bigN);
00176      Teuchos::SerialDenseMatrix<int,Scalar> Vt(bigN,bigN);
00177      Teuchos::SerialDenseMatrix<int,Scalar> work(5*bigN,1);
00178      Teuchos::SerialDenseMatrix<int,Scalar> rWork(1,1);
00179      int info;
00180 
00181      Teuchos::LAPACK<int,Scalar> lala;
00182 
00183      lala.GESVD( 'A',  
00184                  'N',
00185                  V1.numRows() ,
00186                  V1.numCols() ,
00187                  V1.values() ,
00188                  V1.stride() ,
00189                  S.values() ,
00190                  U.values() ,
00191                  U.stride() ,
00192                  Vt.values() ,
00193                  Vt.stride() ,
00194                  work.values() ,
00195                  5*bigN ,
00196                  rWork.values() ,
00197                  &info );
00198                         
00199      int num_nonzero_sv = 0;
00200      for (int i=0;i<bigN;i++) {
00201        if (S(i,0) > INTREPID_TOL) {
00202          num_nonzero_sv++;
00203        }
00204      }
00205      
00206      Teuchos::SerialDenseMatrix<int,Scalar> Uslender(bigN, num_nonzero_sv);
00207      for (int j=0;j<num_nonzero_sv;j++) {
00208        for (int i=0;i<bigN;i++) {
00209          Uslender(i,j) = U(i,j);
00210        }
00211      }
00212 
00213      // apply nodes to big space
00214      Teuchos::SerialDenseMatrix<int,Scalar> V2(N, bigN);
00215 
00216      shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
00217      shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
00218 
00219 
00220      const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
00221                                                            n+1 ,
00222                                                            1 );
00223 
00224      const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
00225                                                            n+1 ,
00226                                                            1 );
00227 
00228      const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
00229                                                            n+1 ,
00230                                                            1 );
00231     
00232      // these hold the reference domain points that will be mapped to each edge or face
00233      FieldContainer<Scalar> oneDPts( numPtsPerEdge , 1 );
00234      FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 );
00235 
00236      if (pointType == POINTTYPE_WARPBLEND) {
00237        CubatureDirectLineGauss<Scalar> edgeRule( numPtsPerEdge );
00238        FieldContainer<Scalar> edgeCubWts( numPtsPerEdge );
00239        edgeRule.getCubature( oneDPts , edgeCubWts );
00240      }
00241      else if (pointType == POINTTYPE_EQUISPACED ) {
00242        PointTools::getLattice<Scalar,FieldContainer<Scalar> >( oneDPts , 
00243                                                                edgeTop ,
00244                                                                n+1 , 
00245                                                                1 ,
00246                                                                pointType );
00247      }
00248 
00249      PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
00250                                                              faceTop ,
00251                                                              n+1 ,
00252                                                              1 ,
00253                                                              pointType );
00254 
00255      FieldContainer<Scalar> edgePts( numPtsPerEdge , 3 );
00256      FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , numPtsPerEdge );
00257 
00258      FieldContainer<Scalar> facePts( numPtsPerFace , 3 );
00259      FieldContainer<Scalar> phisAtFacePoints( scalarBigN , 
00260                                              numPtsPerFace );   
00261 
00262      FieldContainer<Scalar> edgeTan( 3 );
00263     
00264      // loop over the edges
00265      for (int edge=0;edge<6;edge++) {
00266        CellTools<Scalar>::getReferenceEdgeTangent( edgeTan ,
00267                                                    edge ,
00268                                                    this->basisCellTopology_ );
00269        /* multiply by 2.0 to account for a scaling in Pavel's definition */
00270        for (int j=0;j<3;j++) {
00271          edgeTan(j) *= 2.0;
00272        }
00273 
00274        CellTools<Scalar>::mapToReferenceSubcell( edgePts ,
00275                                                  oneDPts ,
00276                                                  1 ,
00277                                                  edge ,
00278                                                  this->basisCellTopology_ );
00279       
00280        Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
00281    
00282        // loop over points (rows of V2)
00283        for (int j=0;j<numPtsPerEdge;j++) {
00284          // loop over orthonormal basis functions (columns of V2)
00285          for (int k=0;k<scalarBigN;k++) {
00286            for (int d=0;d<3;d++) {
00287              V2(edge*numPtsPerEdge+j,k+scalarBigN*d) = edgeTan(d) * phisAtEdgePoints(k,j);
00288            }
00289          }
00290        }
00291      }
00292 
00293    // handle the faces, if needed
00294     if (n > 1) {
00295       FieldContainer<Scalar> refFaceTanU(3);
00296       FieldContainer<Scalar> refFaceTanV(3);
00297       for (int face=0;face<4;face++) {
00298         CellTools<Scalar>::getReferenceFaceTangents( refFaceTanU ,
00299                                                     refFaceTanV ,
00300                                                     face ,
00301                                                     this->basisCellTopology_ );
00302         CellTools<Scalar>::mapToReferenceSubcell( facePts ,
00303                                                   twoDPts ,
00304                                                   2 ,
00305                                                   face ,
00306                                                   this->basisCellTopology_ );
00307         Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE );
00308         for (int j=0;j<numPtsPerFace;j++) {
00309           for (int k=0;k<scalarBigN;k++) {
00310             for (int d=0;d<3;d++) {
00311               V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j,k+scalarBigN*d) =
00312                 refFaceTanU(d) * phisAtFacePoints(k,j);
00313               V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j+1,k+scalarBigN*d) =
00314                 refFaceTanV(d) * phisAtFacePoints(k,j);
00315             }
00316           }
00317         }
00318       }
00319    }
00320 
00321      // internal dof, if needed
00322      if (n > 2) {
00323        FieldContainer<Scalar> cellPoints( numPtsPerCell , 3 );
00324        PointTools::getLattice<Scalar,FieldContainer<Scalar> >( cellPoints ,
00325                                                                this->getBaseCellTopology() , 
00326                                                                n + 1 ,
00327                                                                1 ,
00328                                                                pointType );
00329        FieldContainer<Scalar> phisAtCellPoints( scalarBigN , numPtsPerCell );
00330        Phis_.getValues( phisAtCellPoints , cellPoints , OPERATOR_VALUE );
00331        for (int i=0;i<numPtsPerCell;i++) {
00332          for (int j=0;j<scalarBigN;j++) {
00333            for (int k=0;k<3;k++) {
00334              V2(6*numPtsPerEdge+8*numPtsPerFace+k*numPtsPerCell+i,k*scalarBigN+j) = phisAtCellPoints(j,i);
00335            }
00336          }
00337        }
00338      }
00339 
00340     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
00341     
00342     // multiply V2 * U --> V
00343     Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , Uslender , 0.0 );
00344 
00345     Teuchos::SerialDenseSolver<int,Scalar> solver;
00346     solver.setMatrix( rcp( &Vsdm , false ) );
00347 
00348     solver.invert( );
00349 
00350 
00351     Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
00352     Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , Uslender , Vsdm , 0.0 );
00353 
00354     //std::cout << Csdm << "\n";
00355 
00356     for (int i=0;i<bigN;i++) {
00357       for (int j=0;j<N;j++) {
00358         coeffs_(i,j) = Csdm(i,j);
00359       }
00360     }
00361 
00362     //std::cout << coeffs_ << std::endl;
00363 
00364   }
00365     
00366   template<class Scalar, class ArrayScalar>
00367   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::initializeTags() {
00368     // Basis-dependent initializations
00369     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00370     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00371     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00372     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00373     
00374     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00375     
00376     int *tags = new int[ tagSize * this->getCardinality() ];
00377     int *tag_cur = tags;
00378     const int deg = this->getDegree();
00379 
00380     shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
00381     shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
00382 
00383 
00384     const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
00385                                                           deg+1 ,
00386                                                           1 );
00387 
00388     const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
00389                                                           deg+1 ,
00390                                                           1 );
00391 
00392     const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
00393                                                           deg+1 ,
00394                                                           1 );
00395 
00396     // edge dof first
00397     for (int e=0;e<6;e++) {
00398       for (int i=0;i<numPtsPerEdge;i++) {
00399         tag_cur[0] = 1;  tag_cur[1] = e;  tag_cur[2] = i;  tag_cur[3] = numPtsPerEdge;
00400         tag_cur += tagSize;
00401       }
00402     }
00403 
00404     // face dof, 2 * numPtsPerFace dof per face
00405     for (int f=0;f<4;f++) {
00406       for (int i=0;i<2*numPtsPerFace;i++) {
00407         tag_cur[0] = 2;  tag_cur[1] = f;  tag_cur[2] = i;  tag_cur[3] = 2*numPtsPerFace;
00408         tag_cur+= tagSize;
00409       }
00410     }
00411     
00412     // internal dof, 3 * numPtsPerCell
00413     for (int i=0;i<3*numPtsPerCell;i++) {
00414       tag_cur[0] = 3;  tag_cur[1] = 0;  tag_cur[2] = i;  tag_cur[3] = 3*numPtsPerCell;
00415       tag_cur += tagSize;
00416     }
00417     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00418                                 this -> ordinalToTag_,
00419                                 tags,
00420                                 this -> basisCardinality_,
00421                                 tagSize,
00422                                 posScDim,
00423                                 posScOrd,
00424                                 posDfOrd);
00425 
00426     delete []tags;   
00427     
00428   }  
00429 
00430 
00431 
00432   template<class Scalar, class ArrayScalar> 
00433   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00434                                                               const ArrayScalar &  inputPoints,
00435                                                               const EOperator      operatorType) const {
00436   
00437     // Verify arguments
00438 #ifdef HAVE_INTREPID_DEBUG
00439     Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
00440                                                         inputPoints,
00441                                                         operatorType,
00442                                                         this -> getBaseCellTopology(),
00443                                                         this -> getCardinality() );
00444 #endif
00445     const int numPts = inputPoints.dimension(0);
00446     const int deg = this -> getDegree();
00447     const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6;
00448 
00449     try {
00450       switch (operatorType) {
00451       case OPERATOR_VALUE:
00452         {
00453           FieldContainer<Scalar> phisCur( scalarBigN , numPts );
00454           Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
00455 
00456           for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
00457             for (int j=0;j<outputValues.dimension(1);j++) {  // point
00458               for (int d=0;d<3;d++) {
00459                 outputValues(i,j,d) = 0.0;
00460               }
00461               for (int k=0;k<scalarBigN;k++) { // Dubiner bf
00462                 for (int d=0;d<3;d++) {
00463                   outputValues(i,j,d) += coeffs_(k+d*scalarBigN,i) * phisCur(k,j);
00464                 }
00465               }
00466             }
00467           }
00468         }
00469         break;
00470       case OPERATOR_CURL:
00471         {
00472           FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 );
00473           Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
00474           for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
00475             for (int j=0;j<outputValues.dimension(1);j++) { // point loop
00476               outputValues(i,j,0) = 0.0;
00477               for (int k=0;k<scalarBigN;k++) {
00478                 outputValues(i,j,0) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,1);
00479                 outputValues(i,j,0) -= coeffs_(k+scalarBigN,i) * phisCur(k,j,2);
00480               }
00481               
00482               outputValues(i,j,1) = 0.0;
00483               for (int k=0;k<scalarBigN;k++) {
00484                 outputValues(i,j,1) += coeffs_(k,i) * phisCur(k,j,2);
00485                 outputValues(i,j,1) -= coeffs_(k+2*scalarBigN,i) * phisCur(k,j,0);
00486               }
00487 
00488               outputValues(i,j,2) = 0.0;
00489               for (int k=0;k<scalarBigN;k++) {
00490                 outputValues(i,j,2) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0);
00491                 outputValues(i,j,2) -= coeffs_(k,i) * phisCur(k,j,1);
00492               }
00493             }
00494           }
00495         }
00496         break;
00497       default:
00498         TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00499                             ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
00500         break;
00501       }
00502     }
00503     catch (std::invalid_argument &exception){
00504       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00505                           ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");    
00506     }
00507 
00508   }
00509   
00510 
00511   
00512   template<class Scalar, class ArrayScalar>
00513   void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00514                                                               const ArrayScalar &    inputPoints,
00515                                                               const ArrayScalar &    cellVertices,
00516                                                               const EOperator        operatorType) const {
00517     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00518                         ">>> ERROR (Basis_HCURL_TET_In_FEM): FEM Basis calling an FVD member function");
00519   }
00520 
00521 
00522 }// namespace Intrepid