Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_C2_FEMDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_TET_C2_FEMDEF_HPP
00002 #define INTREPID_HGRAD_TET_C2_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_C2_FEM<Scalar, ArrayScalar>::Basis_HGRAD_TET_C2_FEM()
00055   {
00056     this -> basisCardinality_  = 10;
00057     this -> basisDegree_       = 2;    
00058     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
00059     this -> basisType_         = BASIS_FEM_DEFAULT;
00060     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00061     this -> basisTagsAreSet_   = false;
00062   }
00063   
00064   
00065 template<class Scalar, class ArrayScalar>
00066 void Basis_HGRAD_TET_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
00067   
00068   // Basis-dependent intializations
00069   int tagSize  = 4;        // size of DoF tag
00070   int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00071   int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00072   int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00073 
00074   // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 
00075   int tags[]  = { 0, 0, 0, 1,
00076                   0, 1, 0, 1,
00077                   0, 2, 0, 1,
00078                   0, 3, 0, 1,
00079                   1, 0, 0, 1,
00080                   1, 1, 0, 1,
00081                   1, 2, 0, 1,
00082                   1, 3, 0, 1,
00083                   1, 4, 0, 1,
00084                   1, 5, 0, 1,
00085   };
00086   
00087   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00088   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00089                               this -> ordinalToTag_,
00090                               tags,
00091                               this -> basisCardinality_,
00092                               tagSize,
00093                               posScDim,
00094                               posScOrd,
00095                               posDfOrd);
00096 }
00097 
00098 
00099 
00100 template<class Scalar, class ArrayScalar>
00101 void Basis_HGRAD_TET_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00102                                                              const ArrayScalar &  inputPoints,
00103                                                              const EOperator      operatorType) const {
00104   
00105   // Verify arguments
00106 #ifdef HAVE_INTREPID_DEBUG
00107   Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00108                                                       inputPoints,
00109                                                       operatorType,
00110                                                       this -> getBaseCellTopology(),
00111                                                       this -> getCardinality() );
00112 #endif
00113   
00114   // Number of evaluation points = dim 0 of inputPoints
00115   int dim0 = inputPoints.dimension(0);  
00116   
00117   // Temporaries: (x,y,z) coordinates of the evaluation point
00118   Scalar x = 0.0;                                    
00119   Scalar y = 0.0;   
00120   Scalar z = 0.0;
00121   
00122   switch (operatorType) {
00123     
00124     case OPERATOR_VALUE:
00125       for (int i0 = 0; i0 < dim0; i0++) {
00126         x = inputPoints(i0, 0);
00127         y = inputPoints(i0, 1);
00128         z = inputPoints(i0, 2);
00129         
00130         // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
00131         outputValues(0, i0) = (-1. + x + y + z)*(-1. + 2.*x + 2.*y + 2.*z);
00132         outputValues(1, i0) = x*(-1. + 2.*x);
00133         outputValues(2, i0) = y*(-1. + 2.*y);
00134         outputValues(3, i0) = z*(-1. + 2.*z);
00135 
00136         outputValues(4, i0) = -4.*x*(-1. + x + y + z);
00137         outputValues(5, i0) =  4.*x*y;
00138         outputValues(6, i0) = -4.*y*(-1. + x + y + z);
00139         outputValues(7, i0) = -4.*z*(-1. + x + y + z);
00140         outputValues(8, i0) =  4.*x*z;
00141         outputValues(9, i0) =  4.*y*z;
00142 
00143       }
00144       break;
00145       
00146     case OPERATOR_GRAD:
00147     case OPERATOR_D1:
00148       for (int i0 = 0; i0 < dim0; i0++) {
00149         x = inputPoints(i0, 0);
00150         y = inputPoints(i0, 1);
00151         z = inputPoints(i0, 2);
00152         
00153         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00154         outputValues(0, i0, 0) = -3.+ 4.*x + 4.*y + 4.*z;
00155         outputValues(0, i0, 1) = -3.+ 4.*x + 4.*y + 4.*z;
00156         outputValues(0, i0, 2) = -3.+ 4.*x + 4.*y + 4.*z; 
00157         
00158         outputValues(1, i0, 0) = -1.+ 4.*x; 
00159         outputValues(1, i0, 1) =  0.;
00160         outputValues(1, i0, 2) =  0.;
00161         
00162         outputValues(2, i0, 0) =  0.;        
00163         outputValues(2, i0, 1) = -1.+ 4.*y;
00164         outputValues(2, i0, 2) =  0.;
00165   
00166         outputValues(3, i0, 0) =  0.;         
00167         outputValues(3, i0, 1) =  0.;
00168         outputValues(3, i0, 2) = -1.+ 4.*z;
00169  
00170         
00171         outputValues(4, i0, 0) = -4.*(-1.+ 2*x + y + z);         
00172         outputValues(4, i0, 1) = -4.*x;
00173         outputValues(4, i0, 2) = -4.*x;
00174         
00175         outputValues(5, i0, 0) =  4.*y;       
00176         outputValues(5, i0, 1) =  4.*x; 
00177         outputValues(5, i0, 2) =  0.;
00178 
00179         outputValues(6, i0, 0) = -4.*y;          
00180         outputValues(6, i0, 1) = -4.*(-1.+ x + 2*y + z);
00181         outputValues(6, i0, 2) = -4.*y; 
00182 
00183         outputValues(7, i0, 0) = -4.*z;          
00184         outputValues(7, i0, 1) = -4.*z;
00185         outputValues(7, i0, 2) = -4.*(-1.+ x + y + 2*z);
00186 
00187         outputValues(8, i0, 0) =  4.*z;     
00188         outputValues(8, i0, 1) =  0.;
00189         outputValues(8, i0, 2) =  4.*x;
00190 
00191         outputValues(9, i0, 0) =  0.;         
00192         outputValues(9, i0, 1) =  4.*z;
00193         outputValues(9, i0, 2) =  4.*y;
00194         
00195       }
00196       break;
00197       
00198     case OPERATOR_CURL:
00199       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00200                           ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
00201       break;
00202       
00203     case OPERATOR_DIV:
00204       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00205                           ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
00206       break;
00207       
00208     case OPERATOR_D2:
00209       for (int i0 = 0; i0 < dim0; i0++) {
00210           
00211         outputValues(0, i0, 0) =  4.;
00212         outputValues(0, i0, 1) =  4.;
00213         outputValues(0, i0, 2) =  4.;
00214         outputValues(0, i0, 3) =  4.;
00215         outputValues(0, i0, 4) =  4.;
00216         outputValues(0, i0, 5) =  4.;
00217         
00218         outputValues(1, i0, 0) =  4.;
00219         outputValues(1, i0, 1) =  0.;
00220         outputValues(1, i0, 2) =  0.;
00221         outputValues(1, i0, 3) =  0.;
00222         outputValues(1, i0, 4) =  0.;
00223         outputValues(1, i0, 5) =  0.;
00224 
00225         outputValues(2, i0, 0) =  0.;
00226         outputValues(2, i0, 1) =  0.;
00227         outputValues(2, i0, 2) =  0.;
00228         outputValues(2, i0, 3) =  4.;
00229         outputValues(2, i0, 4) =  0.;
00230         outputValues(2, i0, 5) =  0.;
00231 
00232         outputValues(3, i0, 0) =  0.;
00233         outputValues(3, i0, 1) =  0.;
00234         outputValues(3, i0, 2) =  0.;
00235         outputValues(3, i0, 3) =  0.;
00236         outputValues(3, i0, 4) =  0.;
00237         outputValues(3, i0, 5) =  4.;
00238 
00239         outputValues(4, i0, 0) = -8.;
00240         outputValues(4, i0, 1) = -4.;
00241         outputValues(4, i0, 2) = -4.;
00242         outputValues(4, i0, 3) =  0.;
00243         outputValues(4, i0, 4) =  0.;
00244         outputValues(4, i0, 5) =  0.;
00245 
00246         outputValues(5, i0, 0) =  0.;
00247         outputValues(5, i0, 1) =  4.;
00248         outputValues(5, i0, 2) =  0.;
00249         outputValues(5, i0, 3) =  0.;
00250         outputValues(5, i0, 4) =  0.;
00251         outputValues(5, i0, 5) =  0.;
00252 
00253         outputValues(6, i0, 0) =  0.;
00254         outputValues(6, i0, 1) = -4.;
00255         outputValues(6, i0, 2) =  0.;
00256         outputValues(6, i0, 3) = -8.;
00257         outputValues(6, i0, 4) = -4.;
00258         outputValues(6, i0, 5) =  0;
00259 
00260         outputValues(7, i0, 0) =  0.;
00261         outputValues(7, i0, 1) =  0.;
00262         outputValues(7, i0, 2) = -4.;
00263         outputValues(7, i0, 3) =  0.;
00264         outputValues(7, i0, 4) = -4.;
00265         outputValues(7, i0, 5) = -8.;
00266 
00267         outputValues(8, i0, 0) =  0.;
00268         outputValues(8, i0, 1) =  0.;
00269         outputValues(8, i0, 2) =  4.;
00270         outputValues(8, i0, 3) =  0.;
00271         outputValues(8, i0, 4) =  0.;
00272         outputValues(8, i0, 5) =  0.;
00273 
00274         outputValues(9, i0, 0) =  0.;
00275         outputValues(9, i0, 1) =  0.;
00276         outputValues(9, i0, 2) =  0.;
00277         outputValues(9, i0, 3) =  0.;
00278         outputValues(9, i0, 4) =  4.;
00279         outputValues(9, i0, 5) =  0.;
00280       }
00281       break;
00282       
00283     case OPERATOR_D3:
00284     case OPERATOR_D4:
00285     case OPERATOR_D5:
00286     case OPERATOR_D6:
00287     case OPERATOR_D7:
00288     case OPERATOR_D8:
00289     case OPERATOR_D9:
00290     case OPERATOR_D10:
00291       {
00292         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
00293         int DkCardinality = Intrepid::getDkCardinality(operatorType, 
00294                                                        this -> basisCellTopology_.getDimension() );
00295         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00296           for (int i0 = 0; i0 < dim0; i0++) {
00297             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00298               outputValues(dofOrd, i0, dkOrd) = 0.0;
00299             }
00300           }
00301         }
00302       }
00303       break;
00304     default:
00305       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00306                           ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
00307   }
00308 }
00309 
00310 
00311   
00312 template<class Scalar, class ArrayScalar>
00313 void Basis_HGRAD_TET_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00314                                                              const ArrayScalar &    inputPoints,
00315                                                              const ArrayScalar &    cellVertices,
00316                                                              const EOperator        operatorType) const {
00317   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00318                       ">>> ERROR (Basis_HGRAD_TET_C2_FEM): FEM Basis calling an FVD member function");
00319 }
00320 }// namespace Intrepid
00321 #endif