Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_QUAD_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_QUAD_C2_FEM<Scalar, ArrayScalar>::Basis_HGRAD_QUAD_C2_FEM()
00039   {
00040     this -> basisCardinality_  = 9;
00041     this -> basisDegree_       = 2;    
00042     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<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_QUAD_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
00051   
00052   // Basis-dependent intializations
00053   int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the 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                   // edge midpoints
00064                   1, 0, 0, 1,
00065                   1, 1, 0, 1,
00066                   1, 2, 0, 1,
00067                   1, 3, 0, 1,
00068                   // quad center
00069                   2, 0, 0, 1};
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_QUAD_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) coordinates of the evaluation point
00102   Scalar x = 0.0;                                    
00103   Scalar y = 0.0;                                    
00104   
00105   switch (operatorType) {
00106     
00107     case OPERATOR_VALUE:
00108       for (int i0 = 0; i0 < dim0; i0++) {
00109         x = inputPoints(i0, 0);
00110         y = inputPoints(i0, 1);
00111         
00112         // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
00113         outputValues(0, i0) = x*(x - 1.0)*y*(y - 1.0)/4.0;
00114         outputValues(1, i0) = x*(x + 1.0)*y*(y - 1.0)/4.0;
00115         outputValues(2, i0) = x*(x + 1.0)*y*(y + 1.0)/4.0;
00116         outputValues(3, i0) = x*(x - 1.0)*y*(y + 1.0)/4.0;
00117         // edge midpoints basis functions
00118         outputValues(4, i0) = (1.0 - x)*(1.0 + x)*y*(y - 1.0)/2.0;
00119         outputValues(5, i0) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0;
00120         outputValues(6, i0) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0;
00121         outputValues(7, i0) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0;
00122         // quad bubble basis function
00123         outputValues(8, i0) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y); 
00124       }
00125       break;
00126       
00127     case OPERATOR_GRAD:
00128     case OPERATOR_D1:
00129       for (int i0 = 0; i0 < dim0; i0++) {
00130         x = inputPoints(i0,0);
00131         y = inputPoints(i0,1);
00132         
00133         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00134         outputValues(0, i0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y;
00135         outputValues(0, i0, 1) = (-1.0 + x)*x*(-0.25 + 0.5*y);
00136         
00137         outputValues(1, i0, 0) = (0.25 + 0.5*x)*(-1. + y)*y;
00138         outputValues(1, i0, 1) = x*(1. + x)*(-0.25 + 0.5*y);
00139         
00140         outputValues(2, i0, 0) = (0.25 + 0.5*x)*y*(1. + y);
00141         outputValues(2, i0, 1) = x*(1. + x)*(0.25 + 0.5*y);
00142  
00143         outputValues(3, i0, 0) = (-0.25 + 0.5*x)*y*(1. + y);
00144         outputValues(3, i0, 1) = (-1. + x)*x*(0.25 + 0.5*y);
00145 
00146         outputValues(4, i0, 0) = x*(1.0 - y)*y;
00147         outputValues(4, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
00148           
00149         outputValues(5, i0, 0) = 0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
00150         outputValues(5, i0, 1) =-x*(1.0 + x)*y;
00151           
00152         outputValues(6, i0, 0) =-y*(1.0 + y)*x;
00153         outputValues(6, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
00154           
00155         outputValues(7, i0, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x);
00156         outputValues(7, i0, 1) = (1.0 - x)*x*y;
00157  
00158         outputValues(8, i0, 0) =-2.0*(1.0 - y)*(1.0 + y)*x;
00159         outputValues(8, i0, 1) =-2.0*(1.0 - x)*(1.0 + x)*y;          
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         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00169         // CURL(u) = (u_y, -u_x), is rotated GRAD
00170         outputValues(0, i0, 1) =-(-0.25 + 0.5*x)*(-1. + y)*y;
00171         outputValues(0, i0, 0) = (-1.0 + x)*x*(-0.25 + 0.5*y);
00172         
00173         outputValues(1, i0, 1) =-(0.25 + 0.5*x)*(-1. + y)*y;
00174         outputValues(1, i0, 0) = x*(1. + x)*(-0.25 + 0.5*y);
00175         
00176         outputValues(2, i0, 1) =-(0.25 + 0.5*x)*y*(1. + y);
00177         outputValues(2, i0, 0) = x*(1. + x)*(0.25 + 0.5*y);
00178         
00179         outputValues(3, i0, 1) =-(-0.25 + 0.5*x)*y*(1. + y);
00180         outputValues(3, i0, 0) = (-1. + x)*x*(0.25 + 0.5*y);
00181         
00182         outputValues(4, i0, 1) =-x*(1.0 - y)*y;
00183         outputValues(4, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
00184         
00185         outputValues(5, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
00186         outputValues(5, i0, 0) =-x*(1.0 + x)*y;
00187         
00188         outputValues(6, i0, 1) = y*(1.0 + y)*x;
00189         outputValues(6, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
00190         
00191         outputValues(7, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x);
00192         outputValues(7, i0, 0) = (1.0 - x)*x*y;
00193         
00194         outputValues(8, i0, 1) = 2.0*(1.0 - y)*(1.0 + y)*x;
00195         outputValues(8, i0, 0) =-2.0*(1.0 - x)*(1.0 + x)*y;          
00196       }
00197       break;
00198       
00199     case OPERATOR_DIV:
00200       TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00201                           ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
00202       break;
00203       
00204     case OPERATOR_D2:
00205       for (int i0 = 0; i0 < dim0; i0++) {
00206         x = inputPoints(i0,0);
00207         y = inputPoints(i0,1);
00208        
00209         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality=3) 
00210         outputValues(0, i0, 0) = 0.5*(-1.0 + y)*y;
00211         outputValues(0, i0, 1) = 0.25 - 0.5*y + x*(-0.5 + 1.*y);
00212         outputValues(0, i0, 2) = 0.5*(-1.0 + x)*x;
00213 
00214         outputValues(1, i0, 0) = 0.5*(-1.0 + y)*y;
00215         outputValues(1, i0, 1) =-0.25 + 0.5*y + x*(-0.5 + 1.*y);
00216         outputValues(1, i0, 2) = 0.5*x*(1.0 + x);
00217         
00218         outputValues(2, i0, 0) = 0.5*y*(1.0 + y);
00219         outputValues(2, i0, 1) = 0.25 + 0.5*y + x*(0.5 + 1.*y);
00220         outputValues(2, i0, 2) = 0.5*x*(1.0 + x);
00221         
00222         outputValues(3, i0, 0) = 0.5*y*(1.0 + y);
00223         outputValues(3, i0, 1) =-0.25 - 0.5*y + x*(0.5 + 1.*y);
00224         outputValues(3, i0, 2) = 0.5*(-1.0 + x)*x;
00225         
00226         outputValues(4, i0, 0) = (1.0 - y)*y;
00227         outputValues(4, i0, 1) = x*(1. - 2.*y);
00228         outputValues(4, i0, 2) = (1.0 - x)*(1.0 + x);
00229 
00230         outputValues(5, i0, 0) = (1.0 - y)*(1.0 + y);
00231         outputValues(5, i0, 1) = x*(0. - 2.*y) - 1.*y;
00232         outputValues(5, i0, 2) =-x*(1.0 + x);
00233 
00234         outputValues(6, i0, 0) =-y*(1.0 + y);
00235         outputValues(6, i0, 1) = x*(-1. - 2.*y);
00236         outputValues(6, i0, 2) = (1.0 - x)*(1.0 + x);
00237 
00238         outputValues(7, i0, 0) = (1.0 - y)*(1.0 + y);
00239         outputValues(7, i0, 1) = x*(0. - 2.*y) + 1.*y;
00240         outputValues(7, i0, 2) = (1.0 - x)*x;
00241 
00242         outputValues(8, i0, 0) =-2.0 + 2.0*y*y;
00243         outputValues(8, i0, 1) = 4*x*y;
00244         outputValues(8, i0, 2) =-2.0 + 2.0*x*x;
00245         
00246       }
00247       break;
00248       
00249     case OPERATOR_D3:
00250       // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D3Cardinality=4) 
00251       for (int i0 = 0; i0 < dim0; i0++) {
00252         x = inputPoints(i0,0);
00253         y = inputPoints(i0,1);
00254 
00255         outputValues(0, i0, 0) = 0.0;
00256         outputValues(0, i0, 1) =-0.5 + y;
00257         outputValues(0, i0, 2) =-0.5 + x;
00258         outputValues(0, i0, 3) = 0.0;
00259 
00260         outputValues(1, i0, 0) = 0.0;
00261         outputValues(1, i0, 1) =-0.5 + y;
00262         outputValues(1, i0, 2) = 0.5 + x;
00263         outputValues(1, i0, 3) = 0.0;
00264 
00265         outputValues(2, i0, 0) = 0.0;
00266         outputValues(2, i0, 1) = 0.5 + y;
00267         outputValues(2, i0, 2) = 0.5 + x;
00268         outputValues(2, i0, 3) = 0.0;
00269 
00270         outputValues(3, i0, 0) = 0.0;
00271         outputValues(3, i0, 1) = 0.5 + y;
00272         outputValues(3, i0, 2) =-0.5 + x;
00273         outputValues(3, i0, 3) = 0.0;
00274 
00275         outputValues(4, i0, 0) = 0.0;
00276         outputValues(4, i0, 1) = 1.0 - 2.0*y;
00277         outputValues(4, i0, 2) =-2.0*x;
00278         outputValues(4, i0, 3) = 0.0;
00279 
00280         outputValues(5, i0, 0) = 0.0;
00281         outputValues(5, i0, 1) =-2.0*y;
00282         outputValues(5, i0, 2) =-1.0 - 2.0*x;
00283         outputValues(5, i0, 3) = 0.0;
00284 
00285         outputValues(6, i0, 0) = 0.0;
00286         outputValues(6, i0, 1) =-1.0 - 2.0*y;
00287         outputValues(6, i0, 2) =-2.0*x;
00288         outputValues(6, i0, 3) = 0.0;
00289 
00290         outputValues(7, i0, 0) = 0.0;
00291         outputValues(7, i0, 1) =-2.0*y;
00292         outputValues(7, i0, 2) = 1.0 - 2.0*x;
00293         outputValues(7, i0, 3) = 0.0;        
00294         
00295         outputValues(8, i0, 0) = 0.0;
00296         outputValues(8, i0, 1) = 4.0*y;
00297         outputValues(8, i0, 2) = 4.0*x;
00298         outputValues(8, i0, 3) = 0.0;                
00299       }
00300       break;
00301     
00302     case OPERATOR_D4:
00303       // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D4Cardinality=5) 
00304       for (int i0 = 0; i0 < dim0; i0++) {
00305         
00306         outputValues(0, i0, 0) = 0.0;
00307         outputValues(0, i0, 1) = 0.0;
00308         outputValues(0, i0, 2) = 1.0;
00309         outputValues(0, i0, 3) = 0.0;                
00310         outputValues(0, i0, 4) = 0.0;                
00311 
00312         outputValues(1, i0, 0) = 0.0;
00313         outputValues(1, i0, 1) = 0.0;
00314         outputValues(1, i0, 2) = 1.0;
00315         outputValues(1, i0, 3) = 0.0;                
00316         outputValues(1, i0, 4) = 0.0;                
00317         
00318         outputValues(2, i0, 0) = 0.0;
00319         outputValues(2, i0, 1) = 0.0;
00320         outputValues(2, i0, 2) = 1.0;
00321         outputValues(2, i0, 3) = 0.0;                
00322         outputValues(2, i0, 4) = 0.0;                
00323         
00324         outputValues(3, i0, 0) = 0.0;
00325         outputValues(3, i0, 1) = 0.0;
00326         outputValues(3, i0, 2) = 1.0;
00327         outputValues(3, i0, 3) = 0.0;                
00328         outputValues(3, i0, 4) = 0.0;                
00329         
00330         outputValues(4, i0, 0) = 0.0;
00331         outputValues(4, i0, 1) = 0.0;
00332         outputValues(4, i0, 2) =-2.0;
00333         outputValues(4, i0, 3) = 0.0;                
00334         outputValues(4, i0, 4) = 0.0;                
00335 
00336         outputValues(5, i0, 0) = 0.0;
00337         outputValues(5, i0, 1) = 0.0;
00338         outputValues(5, i0, 2) =-2.0;
00339         outputValues(5, i0, 3) = 0.0;                
00340         outputValues(5, i0, 4) = 0.0;                
00341         
00342         outputValues(6, i0, 0) = 0.0;
00343         outputValues(6, i0, 1) = 0.0;
00344         outputValues(6, i0, 2) =-2.0;
00345         outputValues(6, i0, 3) = 0.0;                
00346         outputValues(6, i0, 4) = 0.0;                
00347         
00348         outputValues(7, i0, 0) = 0.0;
00349         outputValues(7, i0, 1) = 0.0;
00350         outputValues(7, i0, 2) =-2.0;
00351         outputValues(7, i0, 3) = 0.0;                
00352         outputValues(7, i0, 4) = 0.0;                
00353         
00354         outputValues(8, i0, 0) = 0.0;
00355         outputValues(8, i0, 1) = 0.0;
00356         outputValues(8, i0, 2) = 4.0;
00357         outputValues(8, i0, 3) = 0.0;                
00358         outputValues(8, i0, 4) = 0.0;                
00359       }
00360       break;
00361       
00362     case OPERATOR_D5:
00363     case OPERATOR_D6:
00364     case OPERATOR_D7:
00365     case OPERATOR_D8:
00366     case OPERATOR_D9:
00367     case OPERATOR_D10:
00368       {
00369         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
00370         int DkCardinality = Intrepid::getDkCardinality(operatorType, 
00371                                                        this -> basisCellTopology_.getDimension() );
00372         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00373           for (int i0 = 0; i0 < dim0; i0++) {
00374             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00375               outputValues(dofOrd, i0, dkOrd) = 0.0;
00376             }
00377           }
00378         }
00379       }
00380       break;
00381       
00382     default:
00383       TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00384                           ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
00385   }
00386 }
00387 
00388 
00389   
00390 template<class Scalar, class ArrayScalar>
00391 void Basis_HGRAD_QUAD_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00392                                                              const ArrayScalar &    inputPoints,
00393                                                              const ArrayScalar &    cellVertices,
00394                                                              const EOperator        operatorType) const {
00395   TEST_FOR_EXCEPTION( (true), std::logic_error,
00396                       ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): FEM Basis calling an FVD member function");
00397 }
00398 
00399 }// namespace Intrepid
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_QUAD_C2_FEMDef.hpp
