Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_TET_C2

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