Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TRI_C2_FEMDef.hpp
00001 #ifndef INTREPID_HGRAD_TRI_C2_FEMDEF_HPP
00002 #define INTREPID_HGRAD_TRI_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   
00054 template<class Scalar, class ArrayScalar>
00055 Basis_HGRAD_TRI_C2_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TRI_C2_FEM()
00056   {
00057     this -> basisCardinality_  = 6;
00058     this -> basisDegree_       = 2;    
00059     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
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_TRI_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
00069   
00070   // Basis-dependent initializations
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 the 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                   1, 0, 0, 1,
00081                   1, 1, 0, 1,
00082                   1, 2, 0, 1};
00083   
00084   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00085   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00086                               this -> ordinalToTag_,
00087                               tags,
00088                               this -> basisCardinality_,
00089                               tagSize,
00090                               posScDim,
00091                               posScOrd,
00092                               posDfOrd);
00093 }  
00094 
00095 
00096 
00097 template<class Scalar, class ArrayScalar> 
00098 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00099                                                             const ArrayScalar &  inputPoints,
00100                                                             const EOperator      operatorType) const {
00101   
00102   // Verify arguments
00103 #ifdef HAVE_INTREPID_DEBUG
00104   Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00105                                                       inputPoints,
00106                                                       operatorType,
00107                                                       this -> getBaseCellTopology(),
00108                                                       this -> getCardinality() );
00109 #endif
00110   
00111   // Number of evaluation points = dim 0 of inputPoints
00112   int dim0 = inputPoints.dimension(0);  
00113   
00114   // Temporaries: (x,y) coordinates of the evaluation point
00115   Scalar x = 0.0;                                    
00116   Scalar y = 0.0;   
00117   
00118   switch (operatorType) {
00119     
00120     case OPERATOR_VALUE:
00121       for (int i0 = 0; i0 < dim0; i0++) {
00122         x = inputPoints(i0, 0);
00123         y = inputPoints(i0, 1);
00124           
00125         // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
00126         outputValues(0, i0) = (x + y - 1.0)*(2.0*x + 2.0*y - 1.0);
00127         outputValues(1, i0) = x*(2.0*x - 1.0);
00128         outputValues(2, i0) = y*(2.0*y - 1.0);
00129         outputValues(3, i0) = -4.0*x*(x + y - 1.0);
00130         outputValues(4, i0) =  4.0*x*y;
00131         outputValues(5, i0) = -4.0*y*(x + y - 1.0);
00132         
00133       }
00134       break;
00135 
00136     case OPERATOR_GRAD:
00137     case OPERATOR_D1:
00138       for (int i0 = 0; i0 < dim0; i0++) {
00139         x = inputPoints(i0, 0);
00140         y = inputPoints(i0, 1);
00141 
00142         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00143         outputValues(0, i0, 0) =  4.0*x + 4.0*y - 3.0;
00144         outputValues(0, i0, 1) =  4.0*x + 4.0*y - 3.0;
00145 
00146         outputValues(1, i0, 0) =  4.0*x - 1.0;
00147         outputValues(1, i0, 1) =  0.0;
00148 
00149         outputValues(2, i0, 0) =  0.0;
00150         outputValues(2, i0, 1) =  4.0*y - 1.0;
00151         
00152         outputValues(3, i0, 0) = -4.0*(2.0*x + y - 1.0);
00153         outputValues(3, i0, 1) = -4.0*x;
00154 
00155         outputValues(4, i0, 0) =  4.0*y;
00156         outputValues(4, i0, 1) =  4.0*x;
00157 
00158         outputValues(5, i0, 0) = -4.0*y;
00159         outputValues(5, i0, 1) = -4.0*(x + 2.0*y - 1.0);
00160       }
00161       break;
00162 
00163     case OPERATOR_CURL:
00164       for (int i0 = 0; i0 < dim0; i0++) {
00165         x = inputPoints(i0, 0);
00166         y = inputPoints(i0, 1);
00167         
00168         // CURL(u) = (u_y, -u_x), is rotated GRAD
00169         outputValues(0, i0, 1) =-(4.0*x + 4.0*y - 3.0);
00170         outputValues(0, i0, 0) =  4.0*x + 4.0*y - 3.0;
00171         
00172         outputValues(1, i0, 1) =-(4.0*x - 1.0);
00173         outputValues(1, i0, 0) =  0.0;
00174         
00175         outputValues(2, i0, 1) =  0.0;
00176         outputValues(2, i0, 0) =  4.0*y - 1.0;
00177         
00178         outputValues(3, i0, 1) =  4.0*(2.0*x + y - 1.0);
00179         outputValues(3, i0, 0) = -4.0*x;
00180         
00181         outputValues(4, i0, 1) = -4.0*y;
00182         outputValues(4, i0, 0) =  4.0*x;
00183         
00184         outputValues(5, i0, 1) =  4.0*y;
00185         outputValues(5, i0, 0) = -4.0*(x + 2.0*y - 1.0);
00186       }
00187       break;
00188 
00189     case OPERATOR_DIV:
00190       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00191                           ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D.");
00192       break;
00193 
00194     case OPERATOR_D2:
00195       for (int i0 = 0; i0 < dim0; i0++) {
00196         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
00197         // D2 -> (2,0) -> dx^2. 
00198         outputValues(0, i0, 0) = 4.0;
00199         outputValues(1, i0, 0) = 4.0;
00200         outputValues(2, i0, 0) = 0.0;
00201         outputValues(3, i0, 0) =-8.0;
00202         outputValues(4, i0, 0) = 0.0;
00203         outputValues(5, i0, 0) = 0.0;
00204         
00205         // D2 -> (1,1) -> dx dy
00206         outputValues(0, i0, 1) = 4.0;
00207         outputValues(1, i0, 1) = 0.0;
00208         outputValues(2, i0, 1) = 0.0;
00209         outputValues(3, i0, 1) =-4.0;
00210         outputValues(4, i0, 1) = 4.0;
00211         outputValues(5, i0, 1) =-4.0;
00212         
00213         // D2 -> (0,2) -> dy^2
00214         outputValues(0, i0, 2) = 4.0;
00215         outputValues(1, i0, 2) = 0.0;
00216         outputValues(2, i0, 2) = 4.0;
00217         outputValues(3, i0, 2) = 0.0;
00218         outputValues(4, i0, 2) = 0.0;
00219         outputValues(5, i0, 2) =-8.0;
00220       }// for i0
00221       break;
00222       
00223     case OPERATOR_D3:
00224     case OPERATOR_D4:
00225     case OPERATOR_D5:
00226     case OPERATOR_D6:
00227     case OPERATOR_D7:
00228     case OPERATOR_D8:
00229     case OPERATOR_D9:
00230     case OPERATOR_D10:
00231       {
00232         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
00233         int DkCardinality = Intrepid::getDkCardinality(operatorType, 
00234                                                        this -> basisCellTopology_.getDimension() );
00235         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00236           for (int i0 = 0; i0 < dim0; i0++) {
00237             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00238               outputValues(dofOrd, i0, dkOrd) = 0.0;
00239             }
00240           }
00241         }
00242       }
00243       break;
00244 
00245     default:
00246       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00247                           ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type");
00248   }
00249 }
00250   
00251 
00252   
00253 template<class Scalar, class ArrayScalar>
00254 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00255                                                             const ArrayScalar &    inputPoints,
00256                                                             const ArrayScalar &    cellVertices,
00257                                                             const EOperator        operatorType) const {
00258   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00259                       ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): FEM Basis calling an FVD member function");
00260 }
00261 
00262 
00263 }// namespace Intrepid
00264 #endif