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