Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_HEX_Cn_FEMDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_HEX_CN_FEMDEF_HPP
00002 #define INTREPID_HGRAD_HEX_CN_FEMDEF_HPP
00003 // @HEADER
00004 // ************************************************************************
00005 //
00006 //                           Intrepid Package
00007 //                 Copyright (2007) Sandia Corporation
00008 //
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 //
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00040 //                    Denis Ridzal  (dridzal@sandia.gov), or
00041 //                    Kara Peterson (kjpeter@sandia.gov)
00042 //
00043 // ************************************************************************
00044 // @HEADER
00045 
00051 namespace Intrepid {
00052   template<class Scalar, class ArrayScalar>
00053   Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int orderx , 
00054                                                                       const int ordery ,
00055                                                                       const int orderz ,
00056                                                                       const ArrayScalar &pts_x ,
00057                                                                       const ArrayScalar &pts_y ,
00058                                                                       const ArrayScalar &pts_z ):
00059     ptsx_( pts_x.dimension(0) , 1 ),
00060     ptsy_( pts_y.dimension(0) , 1 ),
00061     ptsz_( pts_z.dimension(0) , 1 )
00062   {
00063     for (int i=0;i<pts_x.dimension(0);i++)
00064       {
00065         ptsx_(i,0) = pts_x(i,0);
00066       }
00067     for (int i=0;i<pts_y.dimension(0);i++)
00068       {
00069         ptsy_(i,0) = pts_y(i,0);
00070       }
00071     for (int i=0;i<pts_z.dimension(0);i++)
00072       {
00073         ptsz_(i,0) = pts_z(i,0);
00074       }
00075 
00076     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
00077     bases[0].resize(3);
00078 
00079     bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderx , pts_x ) );
00080     bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( ordery , pts_y ) );
00081     bases[0][2] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderz , pts_z ) );
00082 
00083     this->setBases( bases );
00084 
00085     this->basisCardinality_ = (orderx+1)*(ordery+1)*(orderz+1);
00086     if (orderx >= ordery && orderx >= orderz ) {
00087       this->basisDegree_ = orderx;
00088     }
00089     else if (ordery >= orderx && ordery >= orderz) {
00090       this->basisDegree_ = ordery;
00091     }
00092     else {
00093       this->basisDegree_ = orderz;
00094     }
00095     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00096     this -> basisType_         = BASIS_FEM_FIAT;
00097     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00098     this -> basisTagsAreSet_   = false;
00099   }
00100 
00101   template<class Scalar, class ArrayScalar>
00102   Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int order , 
00103                                                                       const EPointType & pointType ):
00104     ptsx_( order+1 , 1 ),
00105     ptsy_( order+1 , 1 ),
00106     ptsz_( order+1 , 1 )
00107   {
00108     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
00109     bases[0].resize(3);
00110 
00111     bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( order , pointType ) );
00112     // basis is same in each direction, so I only need to instantiate it once!
00113     bases[0][1] = bases[0][0];
00114     bases[0][2] = bases[0][0];
00115 
00116     this->setBases( bases );
00117 
00118     this->basisCardinality_ = (order+1)*(order+1)*(order+1);
00119     this->basisDegree_ = order;
00120     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00121     this -> basisType_         = BASIS_FEM_FIAT;
00122     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00123     this -> basisTagsAreSet_   = false;
00124 
00125     // get points
00126     EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
00127     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
00128                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
00129                                                             order ,
00130                                                             0 ,
00131                                                             pt );
00132     for (int i=0;i<order+1;i++)
00133       {
00134         ptsy_(i,0) = ptsx_(i,0);
00135         ptsz_(i,0) = ptsx_(i,0);
00136       }
00137 
00138   }
00139 
00140   template<class Scalar, class ArrayScalar>
00141   void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::initializeTags()
00142   {
00143     // Basis-dependent initializations
00144     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00145     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00146     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00147     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00148   
00149     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00150 
00151     std::vector<int> tags( tagSize * this->getCardinality() );
00152 
00153     // temporarily just put everything on the cell itself
00154     for (int i=0;i<this->getCardinality();i++) {
00155        tags[tagSize*i] = 3;
00156        tags[tagSize*i+1] = 0;
00157        tags[tagSize*i+2] = i;
00158        tags[tagSize*i+3] = this->getCardinality();
00159      }
00160 
00161     Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
00162     Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
00163     Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2];
00164 
00165 
00166     // now let's try to do it "right"
00167     // let's get the x, y, z bases and their dof
00168     const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags();
00169     const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags();
00170     const std::vector<std::vector<int> >& zdoftags = zBasis_.getAllDofTags();
00171 
00172     std::map<int,std::map<int,int> > total_dof_per_entity;
00173     std::map<int,std::map<int,int> > current_dof_per_entity;
00174 
00175     // vertex dof
00176     for (int i=0;i<8;i++) {
00177       total_dof_per_entity[0][i] = 0;
00178       current_dof_per_entity[0][i] = 0;
00179     }
00180     // edge dof (12 edges)
00181     for (int i=0;i<12;i++) {
00182       total_dof_per_entity[1][i] = 0;
00183       current_dof_per_entity[1][i] = 0;
00184     }
00185     // face dof (6 faces)
00186     for (int i=0;i<6;i++) {
00187       total_dof_per_entity[2][i] = 0;
00188       current_dof_per_entity[2][i] = 0;
00189     }
00190     // internal dof
00191     total_dof_per_entity[3][0] = 0;
00192     //total_dof_per_entity[3][1] = 0;
00193 
00194     // let's tally the total degrees of freedom on each entity
00195     for (int k=0;k<zBasis_.getCardinality();k++) {
00196       const int zdim = zdoftags[k][0];
00197       const int zent = zdoftags[k][1];
00198       for (int j=0;j<yBasis_.getCardinality();j++) {
00199         const int ydim = ydoftags[j][0];
00200         const int yent = ydoftags[j][1];
00201         for (int i=0;i<xBasis_.getCardinality();i++) {
00202           const int xdim = xdoftags[i][0];
00203           const int xent = xdoftags[i][1];
00204           int dofdim;
00205           int dofent;
00206           ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent );
00207           std::cout << i << " " << j << " " << k << "\t" << dofdim << " " << dofent << std::endl;
00208           total_dof_per_entity[dofdim][dofent] += 1;
00209 
00210         }
00211       }
00212     }
00213 
00214     int tagcur = 0;
00215     for (int k=0;k<zBasis_.getCardinality();k++) {
00216       const int zdim = zdoftags[k][0];
00217       const int zent = zdoftags[k][1];
00218       for (int j=0;j<yBasis_.getCardinality();j++) {
00219         const int ydim = ydoftags[j][0];
00220         const int yent = ydoftags[j][1];
00221         for (int i=0;i<xBasis_.getCardinality();i++) {
00222           const int xdim = xdoftags[i][0];
00223           const int xent = xdoftags[i][1];
00224           int dofdim;
00225           int dofent;
00226           ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent );
00227           tags[4*tagcur] = dofdim;
00228           tags[4*tagcur+1] = dofent;
00229           tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
00230           current_dof_per_entity[dofdim][dofent]++;
00231           tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
00232           tagcur++;
00233         }
00234       }
00235     }
00236 
00237     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00238                                 this -> ordinalToTag_,
00239                                 &(tags[0]),
00240                                 this -> basisCardinality_,
00241                                 tagSize,
00242                                 posScDim,
00243                                 posScOrd,
00244                                 posDfOrd);
00245 
00246   }
00247 
00248   template<class Scalar, class ArrayScalar>
00249   void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues ,
00250                                                                const ArrayScalar &inputPoints ,
00251                                                                const EOperator operatorType ) const 
00252   {
00253 #ifdef HAVE_INTREPID_DEBUG
00254     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00255                                                         inputPoints,
00256                                                         operatorType,
00257                                                         this -> getBaseCellTopology(),
00258                                                         this -> getCardinality() );
00259 #endif
00260     
00261     Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
00262     Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
00263     Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2];
00264 
00265 
00266     FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1);
00267     FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1);
00268     FieldContainer<Scalar> zInputPoints(inputPoints.dimension(0),1);
00269 
00270     for (int i=0;i<inputPoints.dimension(0);i++) {
00271       xInputPoints(i,0) = inputPoints(i,0);
00272       yInputPoints(i,0) = inputPoints(i,1);
00273       zInputPoints(i,0) = inputPoints(i,2);
00274     }
00275 
00276     switch (operatorType) {
00277     case OPERATOR_VALUE:
00278       {
00279         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00280         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00281         FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0));
00282 
00283         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
00284         yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
00285         zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
00286 
00287         int bfcur = 0;
00288         for (int k=0;k<zBasis_.getCardinality();k++) {
00289           for (int j=0;j<yBasis_.getCardinality();j++) {
00290             for (int i=0;i<xBasis_.getCardinality();i++) {
00291               for (int l=0;l<inputPoints.dimension(0);l++) {
00292                 outputValues(bfcur,l) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
00293               }
00294               bfcur++;
00295             }
00296           }
00297         }
00298       }
00299       break;
00300     case OPERATOR_D1:
00301     case OPERATOR_GRAD:
00302       {
00303         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00304         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00305         FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0));
00306         FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
00307         FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
00308         FieldContainer<Scalar> zBasisDerivs(zBasis_.getCardinality(),zInputPoints.dimension(0),1);
00309 
00310         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
00311         yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
00312         zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
00313         xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
00314         yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);       
00315         zBasis_.getValues(zBasisDerivs,zInputPoints,OPERATOR_D1);       
00316 
00317         int bfcur = 0;
00318         for (int k=0;k<zBasis_.getCardinality();k++) {
00319           for (int j=0;j<yBasis_.getCardinality();j++) {
00320             for (int i=0;i<xBasis_.getCardinality();i++) {
00321               for (int l=0;l<inputPoints.dimension(0);l++) {
00322                 outputValues(bfcur,l,0) = xBasisDerivs(i,l,0) * yBasisValues(j,l) * zBasisValues(k,l);
00323                 outputValues(bfcur,l,1) = xBasisValues(i,l) * yBasisDerivs(j,l,0) * zBasisValues(k,l);
00324                 outputValues(bfcur,l,2) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisDerivs(k,l,0);
00325               }
00326               bfcur++;
00327             }
00328           }
00329         }
00330       }
00331       break;
00332     case OPERATOR_D2:
00333     case OPERATOR_D3:
00334     case OPERATOR_D4:
00335     case OPERATOR_D5: 
00336     case OPERATOR_D6:
00337     case OPERATOR_D7:
00338     case OPERATOR_D8:
00339     case OPERATOR_D9:
00340     case OPERATOR_D10:
00341       {
00342         FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
00343         FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
00344         FieldContainer<Scalar> zBasisValues(yBasis_.getCardinality(),zInputPoints.dimension(0));
00345 
00346         Teuchos::Array<int> partialMult;
00347         
00348         for (int d=0;d<getDkCardinality(operatorType,3);d++) {
00349           getDkMultiplicities( partialMult , d , operatorType , 3 );
00350           if (partialMult[0] == 0) {
00351             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
00352             xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
00353           }
00354           else {
00355             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
00356             EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 );
00357             xBasis_.getValues( xBasisValues , xInputPoints, xop );
00358             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
00359           }
00360           if (partialMult[1] == 0) {
00361             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
00362             yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
00363           }
00364           else {
00365             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
00366             EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 );
00367             yBasis_.getValues( yBasisValues , yInputPoints, yop );
00368             yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
00369           }
00370           if (partialMult[2] == 0) {
00371             zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0));
00372             zBasis_.getValues( zBasisValues , zInputPoints, OPERATOR_VALUE );
00373           }
00374           else {
00375             zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0),1);
00376             EOperator zop = (EOperator) ( (int) OPERATOR_D1 + partialMult[2] - 1 );
00377             zBasis_.getValues( zBasisValues , zInputPoints, zop );
00378             zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0));
00379           }
00380 
00381 
00382           int bfcur = 0;
00383           for (int k=0;k<zBasis_.getCardinality();k++) {
00384             for (int j=0;j<yBasis_.getCardinality();j++) {
00385               for (int i=0;i<xBasis_.getCardinality();i++) {
00386                 for (int l=0;l<inputPoints.dimension(0);l++) {
00387                   outputValues(bfcur,l,d) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
00388                 }
00389                 bfcur++;
00390               }
00391             }
00392           }
00393         }
00394       }
00395       break;
00396     default:
00397         TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00398                             ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented");
00399         break;
00400     }
00401   }
00402 
00403   template<class Scalar,class ArrayScalar>
00404 void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00405                                                              const ArrayScalar &    inputPoints,
00406                                                              const ArrayScalar &    cellVertices,
00407                                                              const EOperator        operatorType) const {
00408   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00409                       ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): FEM Basis calling an FVD member function");
00410 }
00411 
00412   template<class Scalar,class ArrayScalar>
00413   void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
00414   {
00415     int cur = 0;
00416     for (int k=0;k<ptsz_.dimension(0);k++)
00417       {
00418         for (int j=0;j<ptsy_.dimension(0);j++)
00419           {
00420             for (int i=0;i<ptsx_.dimension(0);i++)
00421               {
00422                 dofCoords(cur,0) = ptsx_(i);
00423                 dofCoords(cur,1) = ptsy_(j);
00424                 dofCoords(cur,2) = ptsz_(k);
00425                 cur++;
00426               }
00427           }
00428       }
00429   }
00430   
00431 }// namespace Intrepid
00432 
00433 #endif