Intrepid
http://trilinos.sandia.gov/packages/docs/r11.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_Cn_FEMDef.hpp
00001 #ifndef INTREPID_HGRAD_TET_CN_FEMDEF_HPP
00002 #define INTREPID_HGRAD_TET_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 
00053   template<class Scalar, class ArrayScalar>
00054   Basis_HGRAD_TET_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TET_Cn_FEM( const int n ,
00055                                                                       const EPointType pointType ):
00056     Phis( n ),
00057     V((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
00058     Vinv((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
00059     latticePts( (n+1)*(n+2)*(n+3)/6 , 3 )
00060   {
00061     const int N = (n+1)*(n+2)*(n+3)/6;
00062     this -> basisCardinality_  = N;
00063     this -> basisDegree_       = n;
00064     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00065     this -> basisType_         = BASIS_FEM_FIAT;
00066     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00067     this -> basisTagsAreSet_   = false;
00068 
00069     // construct lattice
00070 
00071     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
00072                                                             this->getBaseCellTopology() ,
00073                                                             n ,
00074                                                             0 ,
00075                                                             pointType );
00076 
00077     
00078     // form Vandermonde matrix.  Actually, this is the transpose of the VDM,
00079     // so we transpose on copy below.
00080   
00081     Phis.getValues( V , latticePts , OPERATOR_VALUE );
00082 
00083     // now I need to copy V into a Teuchos array to do the inversion
00084     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
00085     for (int i=0;i<N;i++) {
00086       for (int j=0;j<N;j++) {
00087         Vsdm(i,j) = V(i,j);
00088       }
00089     }
00090 
00091     // invert the matrix
00092     Teuchos::SerialDenseSolver<int,Scalar> solver;
00093     solver.setMatrix( rcp( &Vsdm , false ) );
00094     solver.invert( );
00095 
00096     // now I need to copy the inverse into Vinv
00097     for (int i=0;i<N;i++) {
00098       for (int j=0;j<N;j++) {
00099         Vinv(i,j) = Vsdm(j,i);
00100       }
00101     }
00102 
00103   }  
00104   
00105   
00106   template<class Scalar, class ArrayScalar>
00107   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
00108   
00109     // Basis-dependent initializations
00110     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00111     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00112     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00113     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00114   
00115     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00116 
00117     int *tags = new int[ tagSize * this->getCardinality() ];
00118     int *tag_cur = tags;
00119     const int degree = this->getDegree();
00120     const int numEdgeDof = degree - 1;
00121     const int numFaceDof = PointTools::getLatticeSize( shards::CellTopology( shards::getCellTopologyData<shards::Triangle<3> >() ) ,
00122                                                       degree , 
00123                                                       1);
00124     const int numCellDof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
00125                                                       degree ,
00126                                                       1 );
00127     int edge_dof_cur[] = {0,0,0,0,0,0};
00128     int face_dof_cur[] = {0,0,0,0};
00129     int cell_dof_cur = 0;
00130 
00131     // this is the really big mess :(
00132     // BEGIN DOF ON BOTTOM FACE
00133     // first vertex: 0
00134     tag_cur[0] = 0;  tag_cur[1] = 0;  tag_cur[2] = 0;  tag_cur[3] = 1;
00135     tag_cur += tagSize;
00136     // end first vertex
00137 
00138     // internal points on line from vertex 0 to vertex 1.  This is edge 0
00139     for (int i=1;i<degree;i++) {
00140       tag_cur[0] = 1;  tag_cur[1] = 0;  tag_cur[2] = edge_dof_cur[0];  tag_cur[3] = numEdgeDof;
00141       edge_dof_cur[0]++;
00142       tag_cur += tagSize;
00143     }
00144     // end line from vertex 0 to vertex 1
00145 
00146     // begin vertex 1
00147     tag_cur[0] = 0;  tag_cur[1] = 1;  tag_cur[2] = 0;  tag_cur[3] = 1;
00148     tag_cur += tagSize;
00149     // end vertex 1
00150 
00151     // internal lines on bottom face
00152     for (int i=1;i<degree;i++) {
00153       // first dof is on edge 2
00154       tag_cur[0] = 1;  tag_cur[1] = 2;  tag_cur[2] = edge_dof_cur[2];  tag_cur[3] = numEdgeDof;
00155       edge_dof_cur[2]++;
00156       tag_cur += tagSize;
00157       // end dof on edge 2
00158 
00159       // internal points are on bottom face, which is face 3
00160       for (int j=1;j<degree-i;j++) {
00161         tag_cur[0] = 2;  tag_cur[1] = 3;  tag_cur[2] = face_dof_cur[3];  tag_cur[3] = numFaceDof;
00162         face_dof_cur[3]++;
00163         tag_cur += tagSize;
00164       }
00165       // end internal points on face 
00166 
00167       // last dof is on edge 1
00168       tag_cur[0] = 1;  tag_cur[1] = 1;  tag_cur[2] = edge_dof_cur[1];  tag_cur[3] = numEdgeDof;
00169       edge_dof_cur[1]++;
00170       tag_cur += tagSize;
00171       // end dof on edge 1
00172     }
00173     // end internal lines on bottom face
00174 
00175     // vertex 2 on bottom face
00176     tag_cur[0] = 0;  tag_cur[1] = 2;  tag_cur[2] = 0;  tag_cur[3] = 1;
00177     tag_cur += tagSize;
00178     // end vertex 2 on bottom face
00179 
00180     // END DOF ON BOTTOM FACE
00181     
00182     // BEGIN DOF ON INTERNAL FACE SLICES (ascending z)
00183     for (int i=1;i<degree;i++) {
00184       //   bottom line of internal face
00185       //     first point is on edge 3 (from vertex 0 to 3)
00186       tag_cur[0] = 1;  tag_cur[1] = 3;  tag_cur[2] = edge_dof_cur[3];  tag_cur[3] = numEdgeDof;
00187       edge_dof_cur[3]++;
00188       tag_cur += tagSize;
00189       //     end first point
00190       //     points internal to face of vertices (0,1,3), which is face 0
00191       for (int j=1;j<degree-i;j++) {
00192         tag_cur[0] = 2;  tag_cur[1] = 0;  tag_cur[2] = face_dof_cur[0];  tag_cur[3] = numFaceDof;
00193         face_dof_cur[0]++;
00194         tag_cur += tagSize;
00195       }
00196       //     end points internal to face 0
00197       //     last point on bottom line is on edge 4
00198       tag_cur[0] = 1;  tag_cur[1] = 4;  tag_cur[2] = edge_dof_cur[4];  tag_cur[3] = numEdgeDof;
00199       edge_dof_cur[4]++;
00200       tag_cur += tagSize;
00201       //     end last point on bottom edge
00202       //  end bottom line of internal face
00203 
00204       //  begin internal lines of internal face
00205       for (int j=1;j<degree-i;j++) {
00206         //    first point on line is on face of vertices (0,3,2), which is face 2
00207         tag_cur[0] = 2;  tag_cur[1] = 2;  tag_cur[2] = edge_dof_cur[2];  tag_cur[3] = numFaceDof;
00208         edge_dof_cur[2]++;
00209         tag_cur += tagSize;
00210         //    end first point of line
00211         //    begin internal points on the cell
00212         for (int k=1;k<degree-i-j;k++) {
00213           tag_cur[0] = 3;  tag_cur[1] = 0;  tag_cur[2] = cell_dof_cur;  tag_cur[3] = numCellDof;
00214           cell_dof_cur++;
00215           tag_cur += tagSize;
00216         }
00217         //    end internal points on the cell
00218         //    last point on the line is on face with vertices (1,2,3) , which is face 1
00219         tag_cur[0] = 2;  tag_cur[1] = 1;  tag_cur[2] = face_dof_cur[1];  tag_cur[3] = numFaceDof;
00220         face_dof_cur[1]++;
00221         tag_cur += tagSize;
00222         //    end last point of line
00223       }
00224       //  end internal lines of internal face
00225       // begin top point on current face slice:  on edge 5
00226       tag_cur[0] = 1;  tag_cur[1] = 5;  tag_cur[2] = edge_dof_cur[5];  tag_cur[3] = numEdgeDof;
00227       edge_dof_cur[5]++;
00228       tag_cur += 4;
00229       // end top point on current face slice
00230     }
00231     // END DOF ON INTERNAL FACE SLICES
00232 
00233     // TOP VERTEX: 3
00234     tag_cur[0] = 0;  tag_cur[1] = 3;  tag_cur[2] = 0;  tag_cur[3] = 1;
00235     // END TOP VERTEX:3
00236 
00237     // end of really big mess :)
00238     
00239 
00240     
00241     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00242                                 this -> ordinalToTag_,
00243                                 tags,
00244                                 this -> basisCardinality_,
00245                                 tagSize,
00246                                 posScDim,
00247                                 posScOrd,
00248                                 posDfOrd);
00249 
00250     delete []tags;
00251   
00252   }  
00253 
00254 
00255 
00256   template<class Scalar, class ArrayScalar> 
00257   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00258                                                               const ArrayScalar &  inputPoints,
00259                                                               const EOperator      operatorType) const {
00260   
00261     // Verify arguments
00262 #ifdef HAVE_INTREPID_DEBUG
00263     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00264                                                         inputPoints,
00265                                                         operatorType,
00266                                                         this -> getBaseCellTopology(),
00267                                                         this -> getCardinality() );
00268 #endif
00269     const int numPts = inputPoints.dimension(0);
00270     const int numBf = this->getCardinality();
00271 
00272     try {
00273       switch (operatorType) {
00274       case OPERATOR_VALUE:
00275         {
00276           FieldContainer<Scalar> phisCur( numBf , numPts );
00277           Phis.getValues( phisCur , inputPoints , operatorType );
00278           for (int i=0;i<outputValues.dimension(0);i++) {
00279             for (int j=0;j<outputValues.dimension(1);j++) {
00280               outputValues(i,j) = 0.0;
00281               for (int k=0;k<this->getCardinality();k++) {
00282                 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
00283               }
00284             }
00285           }
00286         }
00287         break;
00288       case OPERATOR_GRAD:
00289       case OPERATOR_D1:
00290       case OPERATOR_D2:
00291       case OPERATOR_D3:
00292       case OPERATOR_D4:
00293       case OPERATOR_D5:
00294       case OPERATOR_D6:
00295       case OPERATOR_D7:
00296       case OPERATOR_D8:
00297       case OPERATOR_D9:
00298       case OPERATOR_D10:
00299         {
00300           const int dkcard = 
00301             (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,3): getDkCardinality(operatorType,3);
00302           
00303           FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
00304           Phis.getValues( phisCur , inputPoints , operatorType );
00305 
00306           for (int i=0;i<outputValues.dimension(0);i++) {
00307             for (int j=0;j<outputValues.dimension(1);j++) {
00308               for (int k=0;k<outputValues.dimension(2);k++) {
00309                 outputValues(i,j,k) = 0.0;
00310                 for (int l=0;l<this->getCardinality();l++) {
00311                   outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
00312                 }
00313               }
00314             }
00315           }
00316 
00317         }
00318         break;
00319       default:
00320         TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00321                             ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");
00322         break;
00323       }
00324     }
00325     catch (std::invalid_argument &exception){
00326       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00327                           ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");    
00328     }
00329 
00330   }
00331   
00332 
00333   
00334   template<class Scalar, class ArrayScalar>
00335   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00336                                                               const ArrayScalar &    inputPoints,
00337                                                               const ArrayScalar &    cellVertices,
00338                                                               const EOperator        operatorType) const {
00339     TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00340                         ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): FEM Basis calling an FVD member function");
00341   }
00342 
00343 
00344 }// namespace Intrepid
00345 #endif