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