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