Intrepid
http://trilinos.sandia.gov/packages/docs/r10.8/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TRI_Cn_FEMDef.hpp
00001 #ifndef INTREPID_HGRAD_TRI_CN_FEMDEF_HPP
00002 #define INTREPID_HGRAD_TRI_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_TRI_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TRI_Cn_FEM( const int n ,
00055                                                                       const EPointType pointType ):
00056     Phis( n ),
00057     V((n+1)*(n+2)/2,(n+1)*(n+2)/2),
00058     Vinv((n+1)*(n+2)/2,(n+1)*(n+2)/2),
00059     latticePts( (n+1)*(n+2)/2 , 2 )
00060   {
00061     const int N = (n+1)*(n+2)/2;
00062     this -> basisCardinality_  = N;
00063     this -> basisDegree_       = n;
00064     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
00065     this -> basisType_         = BASIS_FEM_FIAT;
00066     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00067     this -> basisTagsAreSet_   = false;
00068 
00069     // construct lattice
00070 
00071     shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );  
00072     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
00073                                                             myTri_3 ,
00074                                                             n ,
00075                                                             0 ,
00076                                                             pointType );
00077 
00078     
00079     // form Vandermonde matrix.  Actually, this is the transpose of the VDM,
00080     // so we transpose on copy below.
00081   
00082     Phis.getValues( V , latticePts , OPERATOR_VALUE );
00083 
00084     // now I need to copy V into a Teuchos array to do the inversion
00085     Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
00086     for (int i=0;i<N;i++) {
00087       for (int j=0;j<N;j++) {
00088         Vsdm(i,j) = V(i,j);
00089       }
00090     }
00091 
00092     // invert the matrix
00093     Teuchos::SerialDenseSolver<int,Scalar> solver;
00094     solver.setMatrix( rcp( &Vsdm , false ) );
00095     solver.invert( );
00096 
00097     // now I need to copy the inverse into Vinv
00098     for (int i=0;i<N;i++) {
00099       for (int j=0;j<N;j++) {
00100         Vinv(i,j) = Vsdm(j,i);
00101       }
00102     }
00103 
00104   }  
00105   
00106   
00107   template<class Scalar, class ArrayScalar>
00108   void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
00109   
00110     // Basis-dependent initializations
00111     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00112     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00113     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00114     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00115   
00116     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00117 
00118     int *tags = new int[ tagSize * this->getCardinality() ];
00119     int *tag_cur = tags;
00120     const int degree = this->getDegree();
00121 
00122     // BEGIN DOF ALONG BOTTOM EDGE
00123 
00124     // the first dof is on vertex 0
00125     tag_cur[0] = 0;  tag_cur[1] = 0;  tag_cur[2] = 0;  tag_cur[3] = 1;
00126     tag_cur += tagSize;
00127 
00128     // next degree-1 dof are on edge 0
00129     for (int i=1;i<degree;i++) {
00130       tag_cur[0] = 1;  tag_cur[1] = 0; tag_cur[2] = i-1;  tag_cur[3] = degree-1;
00131       tag_cur += tagSize;
00132     }
00133 
00134     // last dof is on vertex 1
00135     tag_cur[0] = 0;  tag_cur[1] = 1;  tag_cur[2] = 0;  tag_cur[3] = 1;
00136     tag_cur += tagSize;
00137 
00138     // END DOF ALONG BOTTOM EDGE
00139 
00140     int num_internal_dof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
00141                                                       this->getDegree() ,
00142                                                       1 );
00143 
00144     int internal_dof_cur = 0;
00145   
00146     // BEGIN DOF ALONG INTERNAL HORIZONTAL LINES
00147     for (int i=1;i<degree;i++) {      
00148       // first dof along line is on edge #2
00149       tag_cur[0] = 1;  tag_cur[1] = 2;  tag_cur[2] = i-1;  tag_cur[3] = degree-1;
00150       tag_cur += tagSize;
00151 
00152       // next dof are internal
00153       for (int j=1;j<degree-i;j++) {
00154         tag_cur[0] = 2;  tag_cur[1] = 0;  tag_cur[2] = internal_dof_cur;  tag_cur[3] = num_internal_dof;
00155         internal_dof_cur++;
00156         tag_cur += tagSize;
00157       }
00158 
00159       // last dof along line is on edge 1
00160       tag_cur[0] = 1;  tag_cur[1] = 1;  tag_cur[2] = i-1;  tag_cur[3] = degree-1;
00161       tag_cur += tagSize;
00162 
00163     }
00164     // END DOF ALONG INTERNAL HORIZONTAL LINES
00165   
00166     // LAST DOF IS ON VERTEX 2
00167     tag_cur[0] = 0;  tag_cur[1] = 2;  tag_cur[2] = 0;  tag_cur[3] = 1;
00168     // END LAST DOF
00169 
00170     
00171     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00172                                 this -> ordinalToTag_,
00173                                 tags,
00174                                 this -> basisCardinality_,
00175                                 tagSize,
00176                                 posScDim,
00177                                 posScOrd,
00178                                 posDfOrd);
00179 
00180     delete []tags;
00181   
00182   }  
00183 
00184 
00185 
00186   template<class Scalar, class ArrayScalar> 
00187   void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00188                                                               const ArrayScalar &  inputPoints,
00189                                                               const EOperator      operatorType) const {
00190   
00191     // Verify arguments
00192 #ifdef HAVE_INTREPID_DEBUG
00193     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00194                                                         inputPoints,
00195                                                         operatorType,
00196                                                         this -> getBaseCellTopology(),
00197                                                         this -> getCardinality() );
00198 #endif
00199     const int numPts = inputPoints.dimension(0);
00200     const int numBf = this->getCardinality();
00201 
00202     try {
00203       switch (operatorType) {
00204       case OPERATOR_VALUE:
00205         {
00206           FieldContainer<Scalar> phisCur( numBf , numPts );
00207           Phis.getValues( phisCur , inputPoints , operatorType );
00208           for (int i=0;i<outputValues.dimension(0);i++) {
00209             for (int j=0;j<outputValues.dimension(1);j++) {
00210               outputValues(i,j) = 0.0;
00211               for (int k=0;k<this->getCardinality();k++) {
00212                 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
00213               }
00214             }
00215           }
00216         }
00217         break;
00218       case OPERATOR_GRAD:
00219       case OPERATOR_D1:
00220       case OPERATOR_D2:
00221       case OPERATOR_D3:
00222       case OPERATOR_D4:
00223       case OPERATOR_D5:
00224       case OPERATOR_D6:
00225       case OPERATOR_D7:
00226       case OPERATOR_D8:
00227       case OPERATOR_D9:
00228       case OPERATOR_D10:
00229         {
00230           const int dkcard = 
00231             (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,2): getDkCardinality(operatorType,2);
00232           
00233           FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
00234           Phis.getValues( phisCur , inputPoints , operatorType );
00235 
00236           for (int i=0;i<outputValues.dimension(0);i++) {
00237             for (int j=0;j<outputValues.dimension(1);j++) {
00238               for (int k=0;k<outputValues.dimension(2);k++) {
00239                 outputValues(i,j,k) = 0.0;
00240                 for (int l=0;l<this->getCardinality();l++) {
00241                   outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
00242                 }
00243               }
00244             }
00245           }
00246         }
00247         break;
00248       case OPERATOR_CURL:  // only works in 2d. first component is -d/dy, second is d/dx
00249         {
00250           FieldContainer<Scalar> phisCur( numBf , numPts , getDkCardinality( OPERATOR_D1 , 2 ) );
00251           Phis.getValues( phisCur , inputPoints , OPERATOR_D1 );
00252 
00253           for (int i=0;i<outputValues.dimension(0);i++) {
00254             for (int j=0;j<outputValues.dimension(1);j++) {
00255               outputValues(i,j,0) = 0.0;
00256               outputValues(i,j,1) = 0.0;
00257               for (int k=0;k<this->getCardinality();k++) {
00258                 outputValues(i,j,0) = this->Vinv(k,i) * phisCur(k,j,1);
00259               }
00260               for (int k=0;k<this->getCardinality();k++) {
00261                 outputValues(i,j,1) -= this->Vinv(k,i) * phisCur(k,j,0);
00262               }
00263             }
00264           }
00265         }
00266         break;
00267       default:
00268         TEST_FOR_EXCEPTION( true , std::invalid_argument,
00269                             ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");
00270         break;
00271       }
00272     }
00273     catch (std::invalid_argument &exception){
00274       TEST_FOR_EXCEPTION( true , std::invalid_argument,
00275                           ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");    
00276     }
00277 
00278   }
00279   
00280 
00281   
00282   template<class Scalar, class ArrayScalar>
00283   void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00284                                                               const ArrayScalar &    inputPoints,
00285                                                               const ArrayScalar &    cellVertices,
00286                                                               const EOperator        operatorType) const {
00287     TEST_FOR_EXCEPTION( (true), std::logic_error,
00288                         ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function");
00289   }
00290 
00291 
00292 }// namespace Intrepid
00293 #endif