Intrepid
http://trilinos.sandia.gov/packages/docs/r10.6/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_Cn_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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or
00025 //                    Denis Ridzal (dridzal@sandia.gov).
00026 //
00027 // ************************************************************************
00028 // @HEADER
00029 
00035 namespace Intrepid {
00036 
00037   template<class Scalar, class ArrayScalar>
00038   Basis_HGRAD_TET_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TET_Cn_FEM( const int n ,
00039                                                                       const EPointType pointType ):
00040     Phis( n ),
00041     V((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
00042     Vinv((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
00043     latticePts( (n+1)*(n+2)*(n+3)/6 , 3 )
00044   {
00045     const int N = (n+1)*(n+2)*(n+3)/6;
00046     this -> basisCardinality_  = N;
00047     this -> basisDegree_       = n;
00048     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00049     this -> basisType_         = BASIS_FEM_FIAT;
00050     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00051     this -> basisTagsAreSet_   = false;
00052 
00053     // construct lattice
00054 
00055     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
00056                                                             this->getBaseCellTopology() ,
00057                                                             n ,
00058                                                             0 ,
00059                                                             pointType );
00060 
00061     
00062     // form Vandermonde matrix.  Actually, this is the transpose of the VDM,
00063     // so we transpose on copy below.
00064   
00065     Phis.getValues( V , latticePts , OPERATOR_VALUE );
00066 
00067     // now I need to copy V into a Teuchos array to do the inversion
00068     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
00069     for (int i=0;i<N;i++) {
00070       for (int j=0;j<N;j++) {
00071         Vsdm(i,j) = V(i,j);
00072       }
00073     }
00074 
00075     // invert the matrix
00076     Teuchos::SerialDenseSolver<int,Scalar> solver;
00077     solver.setMatrix( rcp( &Vsdm , false ) );
00078     solver.invert( );
00079 
00080     // now I need to copy the inverse into Vinv
00081     for (int i=0;i<N;i++) {
00082       for (int j=0;j<N;j++) {
00083         Vinv(i,j) = Vsdm(j,i);
00084       }
00085     }
00086 
00087   }  
00088   
00089   
00090   template<class Scalar, class ArrayScalar>
00091   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
00092   
00093     // Basis-dependent initializations
00094     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00095     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00096     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00097     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00098   
00099     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00100 
00101     int *tags = new int[ tagSize * this->getCardinality() ];
00102     int *tag_cur = tags;
00103     const int degree = this->getDegree();
00104     const int numEdgeDof = degree - 1;
00105     const int numFaceDof = PointTools::getLatticeSize( shards::CellTopology( shards::getCellTopologyData<shards::Triangle<3> >() ) ,
00106                                                       degree , 
00107                                                       1);
00108     const int numCellDof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
00109                                                       degree ,
00110                                                       1 );
00111     int edge_dof_cur[] = {0,0,0,0,0,0};
00112     int face_dof_cur[] = {0,0,0,0};
00113     int cell_dof_cur = 0;
00114 
00115     // this is the really big mess :(
00116     // BEGIN DOF ON BOTTOM FACE
00117     // first vertex: 0
00118     tag_cur[0] = 0;  tag_cur[1] = 0;  tag_cur[2] = 0;  tag_cur[3] = 1;
00119     tag_cur += tagSize;
00120     // end first vertex
00121 
00122     // internal points on line from vertex 0 to vertex 1.  This is edge 0
00123     for (int i=1;i<degree;i++) {
00124       tag_cur[0] = 1;  tag_cur[1] = 0;  tag_cur[2] = edge_dof_cur[0];  tag_cur[3] = numEdgeDof;
00125       edge_dof_cur[0]++;
00126       tag_cur += tagSize;
00127     }
00128     // end line from vertex 0 to vertex 1
00129 
00130     // begin vertex 1
00131     tag_cur[0] = 0;  tag_cur[1] = 1;  tag_cur[2] = 0;  tag_cur[3] = 1;
00132     tag_cur += tagSize;
00133     // end vertex 1
00134 
00135     // internal lines on bottom face
00136     for (int i=1;i<degree;i++) {
00137       // first dof is on edge 2
00138       tag_cur[0] = 1;  tag_cur[1] = 2;  tag_cur[2] = edge_dof_cur[2];  tag_cur[3] = numEdgeDof;
00139       edge_dof_cur[2]++;
00140       tag_cur += tagSize;
00141       // end dof on edge 2
00142 
00143       // internal points are on bottom face, which is face 3
00144       for (int j=1;j<degree-i;j++) {
00145         tag_cur[0] = 2;  tag_cur[1] = 3;  tag_cur[2] = face_dof_cur[3];  tag_cur[3] = numFaceDof;
00146         face_dof_cur[3]++;
00147         tag_cur += tagSize;
00148       }
00149       // end internal points on face 
00150 
00151       // last dof is on edge 1
00152       tag_cur[0] = 1;  tag_cur[1] = 1;  tag_cur[2] = edge_dof_cur[1];  tag_cur[3] = numEdgeDof;
00153       edge_dof_cur[1]++;
00154       tag_cur += tagSize;
00155       // end dof on edge 1
00156     }
00157     // end internal lines on bottom face
00158 
00159     // vertex 2 on bottom face
00160     tag_cur[0] = 0;  tag_cur[1] = 2;  tag_cur[2] = 0;  tag_cur[3] = 1;
00161     tag_cur += tagSize;
00162     // end vertex 2 on bottom face
00163 
00164     // END DOF ON BOTTOM FACE
00165     
00166     // BEGIN DOF ON INTERNAL FACE SLICES (ascending z)
00167     for (int i=1;i<degree;i++) {
00168       //   bottom line of internal face
00169       //     first point is on edge 3 (from vertex 0 to 3)
00170       tag_cur[0] = 1;  tag_cur[1] = 3;  tag_cur[2] = edge_dof_cur[3];  tag_cur[3] = numEdgeDof;
00171       edge_dof_cur[3]++;
00172       tag_cur += tagSize;
00173       //     end first point
00174       //     points internal to face of vertices (0,1,3), which is face 0
00175       for (int j=1;j<degree-i;j++) {
00176         tag_cur[0] = 2;  tag_cur[1] = 0;  tag_cur[2] = face_dof_cur[0];  tag_cur[3] = numFaceDof;
00177         face_dof_cur[0]++;
00178         tag_cur += tagSize;
00179       }
00180       //     end points internal to face 0
00181       //     last point on bottom line is on edge 4
00182       tag_cur[0] = 1;  tag_cur[1] = 4;  tag_cur[2] = edge_dof_cur[4];  tag_cur[3] = numEdgeDof;
00183       edge_dof_cur[4]++;
00184       tag_cur += tagSize;
00185       //     end last point on bottom edge
00186       //  end bottom line of internal face
00187 
00188       //  begin internal lines of internal face
00189       for (int j=1;j<degree-i;j++) {
00190         //    first point on line is on face of vertices (0,3,2), which is face 2
00191         tag_cur[0] = 2;  tag_cur[1] = 2;  tag_cur[2] = edge_dof_cur[2];  tag_cur[3] = numFaceDof;
00192         edge_dof_cur[2]++;
00193         tag_cur += tagSize;
00194         //    end first point of line
00195         //    begin internal points on the cell
00196         for (int k=1;k<degree-i-j;k++) {
00197           tag_cur[0] = 3;  tag_cur[1] = 0;  tag_cur[2] = cell_dof_cur;  tag_cur[3] = numCellDof;
00198           cell_dof_cur++;
00199           tag_cur += tagSize;
00200         }
00201         //    end internal points on the cell
00202         //    last point on the line is on face with vertices (1,2,3) , which is face 1
00203         tag_cur[0] = 2;  tag_cur[1] = 1;  tag_cur[2] = face_dof_cur[1];  tag_cur[3] = numFaceDof;
00204         face_dof_cur[1]++;
00205         tag_cur += tagSize;
00206         //    end last point of line
00207       }
00208       //  end internal lines of internal face
00209       // begin top point on current face slice:  on edge 5
00210       tag_cur[0] = 1;  tag_cur[1] = 5;  tag_cur[2] = edge_dof_cur[5];  tag_cur[3] = numEdgeDof;
00211       edge_dof_cur[5]++;
00212       tag_cur += 4;
00213       // end top point on current face slice
00214     }
00215     // END DOF ON INTERNAL FACE SLICES
00216 
00217     // TOP VERTEX: 3
00218     tag_cur[0] = 0;  tag_cur[1] = 3;  tag_cur[2] = 0;  tag_cur[3] = 1;
00219     // END TOP VERTEX:3
00220 
00221     // end of really big mess :)
00222     
00223 
00224     
00225     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00226                                 this -> ordinalToTag_,
00227                                 tags,
00228                                 this -> basisCardinality_,
00229                                 tagSize,
00230                                 posScDim,
00231                                 posScOrd,
00232                                 posDfOrd);
00233 
00234     delete []tags;
00235   
00236   }  
00237 
00238 
00239 
00240   template<class Scalar, class ArrayScalar> 
00241   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00242                                                               const ArrayScalar &  inputPoints,
00243                                                               const EOperator      operatorType) const {
00244   
00245     // Verify arguments
00246 #ifdef HAVE_INTREPID_DEBUG
00247     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00248                                                         inputPoints,
00249                                                         operatorType,
00250                                                         this -> getBaseCellTopology(),
00251                                                         this -> getCardinality() );
00252 #endif
00253     const int numPts = inputPoints.dimension(0);
00254     const int numBf = this->getCardinality();
00255 
00256     try {
00257       switch (operatorType) {
00258       case OPERATOR_VALUE:
00259         {
00260           FieldContainer<Scalar> phisCur( numBf , numPts );
00261           Phis.getValues( phisCur , inputPoints , operatorType );
00262           for (int i=0;i<outputValues.dimension(0);i++) {
00263             for (int j=0;j<outputValues.dimension(1);j++) {
00264               outputValues(i,j) = 0.0;
00265               for (int k=0;k<this->getCardinality();k++) {
00266                 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
00267               }
00268             }
00269           }
00270         }
00271         break;
00272       case OPERATOR_GRAD:
00273       case OPERATOR_D1:
00274       case OPERATOR_D2:
00275       case OPERATOR_D3:
00276       case OPERATOR_D4:
00277       case OPERATOR_D5:
00278       case OPERATOR_D6:
00279       case OPERATOR_D7:
00280       case OPERATOR_D8:
00281       case OPERATOR_D9:
00282       case OPERATOR_D10:
00283         {
00284           const int dkcard = 
00285             (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,3): getDkCardinality(operatorType,3);
00286           
00287           FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
00288           Phis.getValues( phisCur , inputPoints , operatorType );
00289 
00290           for (int i=0;i<outputValues.dimension(0);i++) {
00291             for (int j=0;j<outputValues.dimension(1);j++) {
00292               for (int k=0;k<outputValues.dimension(2);k++) {
00293                 outputValues(i,j,k) = 0.0;
00294                 for (int l=0;l<this->getCardinality();l++) {
00295                   outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
00296                 }
00297               }
00298             }
00299           }
00300 
00301         }
00302         break;
00303       default:
00304         TEST_FOR_EXCEPTION( true , std::invalid_argument,
00305                             ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");
00306         break;
00307       }
00308     }
00309     catch (std::invalid_argument &exception){
00310       TEST_FOR_EXCEPTION( true , std::invalid_argument,
00311                           ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");    
00312     }
00313 
00314   }
00315   
00316 
00317   
00318   template<class Scalar, class ArrayScalar>
00319   void Basis_HGRAD_TET_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00320                                                               const ArrayScalar &    inputPoints,
00321                                                               const ArrayScalar &    cellVertices,
00322                                                               const EOperator        operatorType) const {
00323     TEST_FOR_EXCEPTION( (true), std::logic_error,
00324                         ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): FEM Basis calling an FVD member function");
00325   }
00326 
00327 
00328 }// namespace Intrepid