Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_WEDGE_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_WEDGE_C2_FEM<Scalar, ArrayScalar>::Basis_HGRAD_WEDGE_C2_FEM()
00039   {
00040     this -> basisCardinality_  = 18;
00041     this -> basisDegree_       = 2;    
00042     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >() );
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_WEDGE_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                   0, 4, 0, 1,
00064                   0, 5, 0, 1,
00065                   1, 0, 0, 1,
00066                   1, 1, 0, 1,
00067                   1, 2, 0, 1,
00068                   1, 6, 0, 1,
00069                   1, 7, 0, 1,
00070                   1, 8, 0, 1,
00071                   1, 3, 0, 1,
00072                   1, 4, 0, 1,
00073                   1, 5, 0, 1,
00074                   2, 0, 0, 1,
00075                   2, 1, 0, 1,
00076                   2, 2, 0, 1
00077   };
00078   
00079   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00080   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00081                               this -> ordinalToTag_,
00082                               tags,
00083                               this -> basisCardinality_,
00084                               tagSize,
00085                               posScDim,
00086                               posScOrd,
00087                               posDfOrd);
00088 }
00089 
00090 
00091 
00092 template<class Scalar, class ArrayScalar>
00093 void Basis_HGRAD_WEDGE_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &    outputValues,
00094                                                              const ArrayScalar &  inputPoints,
00095                                                              const EOperator      operatorType) const {
00096   
00097   // Verify arguments
00098 #ifdef HAVE_INTREPID_DEBUG
00099   Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00100                                                       inputPoints,
00101                                                       operatorType,
00102                                                       this -> getBaseCellTopology(),
00103                                                       this -> getCardinality() );
00104 #endif
00105   
00106   // Number of evaluation points = dim 0 of inputPoints
00107   int dim0 = inputPoints.dimension(0);  
00108   
00109   // Temporaries: (x,y,z) coordinates of the evaluation point
00110   Scalar x = 0.0;                                    
00111   Scalar y = 0.0;   
00112   Scalar z = 0.0;
00113   
00114   switch (operatorType) {
00115     
00116     case OPERATOR_VALUE:
00117       for (int i0 = 0; i0 < dim0; i0++) {
00118         x = inputPoints(i0, 0);
00119         y = inputPoints(i0, 1);
00120         z = inputPoints(i0, 2);
00121         
00122         // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
00123         outputValues(0, i0) =  ((-1. + x + y)*(-1. + 2.*x + 2.*y)*(-1. + z)*z)/2.;
00124         outputValues(1, i0) =  (x*(-1. + 2.*x)*(-1. + z)*z)/2.;
00125         outputValues(2, i0) =  (y*(-1. + 2.*y)*(-1. + z)*z)/2.;
00126         outputValues(3, i0) =  ((-1. + x + y)*(-1. + 2.*x + 2.*y)*z*(1. + z))/2.;
00127         outputValues(4, i0) =  (x*(-1. + 2.*x)*z*(1. + z))/2.;
00128         outputValues(5, i0) =  (y*(-1. + 2.*y)*z*(1. + z))/2.;
00129         
00130         outputValues(6, i0) = -2.*x*(-1. + x + y)*(-1. + z)*z;
00131         outputValues(7, i0) =  2.*x*y*(-1. + z)*z;
00132         outputValues(8, i0) = -2.*y*(-1. + x + y)*(-1. + z)*z;
00133         outputValues(9, i0) = -((-1. + x + y)*(-1. + 2.*x + 2.*y)*(-1. + z)*(1. + z));
00134         outputValues(10,i0) = -(x*(-1. + 2.*x)*(-1. + z)*(1. + z));
00135         outputValues(11,i0) = -(y*(-1. + 2.*y)*(-1. + z)*(1. + z));
00136         outputValues(12,i0) = -2.*x*(-1. + x + y)*z*(1. + z);
00137         outputValues(13,i0) =  2.*x*y*z*(1. + z);
00138         outputValues(14,i0) = -2.*y*(-1. + x + y)*z*(1. + z);
00139         outputValues(15,i0) =  4.*x*(-1. + x + y)*(-1. + z)*(1. + z);
00140         outputValues(16,i0) = -4.*x*y*(-1. + z)*(1. + z);
00141         outputValues(17,i0) =  4.*y*(-1. + x + y)*(-1. + z)*(1. + z);
00142       }
00143       break;
00144       
00145     case OPERATOR_GRAD:
00146     case OPERATOR_D1:
00147       for (int i0 = 0; i0 < dim0; i0++) {
00148         x = inputPoints(i0,0);
00149         y = inputPoints(i0,1);
00150         z = inputPoints(i0,2);
00151         
00152         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00153         outputValues(0, i0, 0) = ((-3 + 4*x + 4*y)*(-1 + z)*z)/2.;
00154         outputValues(0, i0, 1) = ((-3 + 4*x + 4*y)*(-1 + z)*z)/2.;
00155         outputValues(0, i0, 2) = ((-1 + x + y)*(-1 + 2*x + 2*y)*(-1 + 2*z))/2.;
00156         
00157         outputValues(1, i0, 0) = ((-1 + 4*x)*(-1 + z)*z)/2.;
00158         outputValues(1, i0, 1) = 0.;
00159         outputValues(1, i0, 2) = (x*(-1 + 2*x)*(-1 + 2*z))/2.;
00160         
00161         outputValues(2, i0, 0) = 0.;
00162         outputValues(2, i0, 1) = ((-1 + 4*y)*(-1 + z)*z)/2.;
00163         outputValues(2, i0, 2) = (y*(-1 + 2*y)*(-1 + 2*z))/2.;
00164         
00165         outputValues(3, i0, 0) = ((-3 + 4*x + 4*y)*z*(1 + z))/2.; 
00166         outputValues(3, i0, 1) = ((-3 + 4*x + 4*y)*z*(1 + z))/2.;  
00167         outputValues(3, i0, 2) = ((-1 + x + y)*(-1 + 2*x + 2*y)*(1 + 2*z))/2.;
00168         
00169         outputValues(4, i0, 0) = ((-1 + 4*x)*z*(1 + z))/2.;
00170         outputValues(4, i0, 1) = 0.;
00171         outputValues(4, i0, 2) = (x*(-1 + 2*x)*(1 + 2*z))/2.;
00172         
00173         outputValues(5, i0, 0) = 0.;
00174         outputValues(5, i0, 1) = ((-1 + 4*y)*z*(1 + z))/2.; 
00175         outputValues(5, i0, 2) = (y*(-1 + 2*y)*(1 + 2*z))/2.;
00176         
00177         outputValues(6, i0, 0) = -2*(-1 + 2*x + y)*(-1 + z)*z;
00178         outputValues(6, i0, 1) = -2*x*(-1 + z)*z;
00179         outputValues(6, i0, 2) = 2*x*(-1 + x + y)*(1 - 2*z);   
00180         
00181         outputValues(7, i0, 0) = 2*y*(-1 + z)*z;
00182         outputValues(7, i0, 1) = 2*x*(-1 + z)*z;
00183         outputValues(7, i0, 2) = 2*x*y*(-1 + 2*z);   
00184         
00185         outputValues(8, i0, 0) = -2*y*(-1 + z)*z;
00186         outputValues(8, i0, 1) = -2*(-1 + x + 2*y)*(-1 + z)*z;
00187         outputValues(8, i0, 2) = 2*y*(-1 + x + y)*(1 - 2*z);   
00188         
00189         outputValues(9, i0, 0) = -(-3 + 4*x + 4*y)*(-1 + z*z);
00190         outputValues(9, i0, 1) = -(-3 + 4*x + 4*y)*(-1 + z*z);   
00191         outputValues(9, i0, 2) = -2*(1 + 2*x*x - 3*y + 2*y*y + x*(-3 + 4*y))*z;   
00192         
00193         outputValues(10,i0, 0) = -(-1 + 4*x)*(-1 + z*z);
00194         outputValues(10,i0, 1) =  0;
00195         outputValues(10,i0, 2) =  2*(1 - 2*x)*x*z;   
00196         
00197         outputValues(11,i0, 0) =  0;
00198         outputValues(11,i0, 1) =  -(-1 + 4*y)*(-1 + z*z);
00199         outputValues(11,i0, 2) =  2*(1 - 2*y)*y*z;   
00200         
00201         outputValues(12,i0, 0) = -2*(-1 + 2*x + y)*z*(1 + z);
00202         outputValues(12,i0, 1) = -2*x*z*(1 + z);
00203         outputValues(12,i0, 2) = -2*x*(-1 + x + y)*(1 + 2*z);   
00204         
00205         outputValues(13,i0, 0) =  2*y*z*(1 + z);
00206         outputValues(13,i0, 1) =  2*x*z*(1 + z);
00207         outputValues(13,i0, 2) =  2*x*y*(1 + 2*z);   
00208         
00209         outputValues(14,i0, 0) = -2*y*z*(1 + z);
00210         outputValues(14,i0, 1) = -2*(-1 + x + 2*y)*z*(1 + z);
00211         outputValues(14,i0, 2) = -2*y*(-1 + x + y)*(1 + 2*z);   
00212         
00213         outputValues(15,i0, 0) =  4*(-1 + 2*x + y)*(-1 + z*z);
00214         outputValues(15,i0, 1) =  4*x*(-1 + z)*(1 + z);
00215         outputValues(15,i0, 2) =  8*x*(-1 + x + y)*z;   
00216         
00217         outputValues(16,i0, 0) = -4*y*(-1 + z)*(1 + z);
00218         outputValues(16,i0, 1) = -4*x*(-1 + z)*(1 + z);
00219         outputValues(16,i0, 2) = -8*x*y*z;   
00220         
00221         outputValues(17,i0, 0) =  4*y*(-1 + z)*(1 + z);
00222         outputValues(17,i0, 1) =  4*(-1 + x + 2*y)*(-1 + z*z);
00223         outputValues(17,i0, 2) =  8*y*(-1 + x + y)*z;
00224         
00225       }
00226       break;
00227       
00228     case OPERATOR_CURL:
00229       TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00230                           ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
00231       break;
00232       
00233     case OPERATOR_DIV:
00234       TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00235                           ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
00236       break;
00237       
00238     case OPERATOR_D2:
00239       for (int i0 = 0; i0 < dim0; i0++) {
00240         x = inputPoints(i0,0);
00241         y = inputPoints(i0,1);
00242         z = inputPoints(i0,2);
00243         
00244         outputValues(0, i0, 0) =  2.*(-1. + z)*z;     
00245         outputValues(0, i0, 1) =  2.*(-1. + z)*z;     
00246         outputValues(0, i0, 2) =  ((-3. + 4.*x + 4.*y)*(-1. + 2.*z))/2.;       
00247         outputValues(0, i0, 3) =  2.*(-1. + z)*z;     
00248         outputValues(0, i0, 4) =  ((-3. + 4.*x + 4.*y)*(-1. + 2.*z))/2.;      
00249         outputValues(0, i0, 5) =  (-1. + x + y)*(-1. + 2.*x + 2.*y);     
00250         
00251         outputValues(1, i0, 0) =  2.*(-1. + z)*z;      
00252         outputValues(1, i0, 1) =  0.;     
00253         outputValues(1, i0, 2) =  ((-1. + 4.*x)*(-1. + 2.*z))/2.;     
00254         outputValues(1, i0, 3) =  0.;     
00255         outputValues(1, i0, 4) =  0.;     
00256         outputValues(1, i0, 5) =  x*(-1. + 2.*x);     
00257         
00258         outputValues(2, i0, 0) =  0.;     
00259         outputValues(2, i0, 1) =  0.;     
00260         outputValues(2, i0, 2) =  0.;     
00261         outputValues(2, i0, 3) =  2.*(-1. + z)*z;      
00262         outputValues(2, i0, 4) =  ((-1. + 4.*y)*(-1. + 2.*z))/2.;     
00263         outputValues(2, i0, 5) =  y*(-1. + 2.*y);     
00264         
00265         outputValues(3, i0, 0) =  2.*z*(1. + z); 
00266         outputValues(3, i0, 1) =  2.*z*(1. + z);
00267         outputValues(3, i0, 2) =  ((-3. + 4.*x + 4.*y)*(1. + 2.*z))/2.;  
00268         outputValues(3, i0, 3) =  2.*z*(1. + z);
00269         outputValues(3, i0, 4) =  ((-3. + 4.*x + 4.*y)*(1. + 2.*z))/2.;
00270         outputValues(3, i0, 5) =  (-1. + x + y)*(-1. + 2.*x + 2.*y);
00271         
00272         outputValues(4, i0, 0) =  2.*z*(1. + z);
00273         outputValues(4, i0, 1) =  0.; 
00274         outputValues(4, i0, 2) =  ((-1. + 4.*x)*(1. + 2.*z))/2.;;
00275         outputValues(4, i0, 3) =  0.;
00276         outputValues(4, i0, 4) =  0.;
00277         outputValues(4, i0, 5) =  x*(-1. + 2.*x);
00278         
00279         outputValues(5, i0, 0) =  0.;
00280         outputValues(5, i0, 1) =  0.;
00281         outputValues(5, i0, 2) =  0.;
00282         outputValues(5, i0, 3) =  2.*z*(1. + z);
00283         outputValues(5, i0, 4) =  ((-1. + 4.*y)*(1. + 2.*z))/2.;
00284         outputValues(5, i0, 5) =  y*(-1. + 2.*y);         
00285         
00286         outputValues(6, i0, 0) = -4.*(-1. + z)*z;
00287         outputValues(6, i0, 1) = -2.*(-1. + z)*z;
00288         outputValues(6, i0, 2) = -2.*(-1. + 2.*x + y)*(-1. + 2.*z);
00289         outputValues(6, i0, 3) =  0.;
00290         outputValues(6, i0, 4) =  x*(2. - 4.*z);     
00291         outputValues(6, i0, 5) = -4.*x*(-1. + x + y);     
00292         
00293         outputValues(7, i0, 0) =  0.;
00294         outputValues(7, i0, 1) =  2.*(-1. + z)*z;
00295         outputValues(7, i0, 2) =  2.*y*(-1. + 2.*z);
00296         outputValues(7, i0, 3) =  0.;
00297         outputValues(7, i0, 4) =  2.*x*(-1. + 2.*z);
00298         outputValues(7, i0, 5) =  4.*x*y;     
00299         
00300         outputValues(8, i0, 0) =  0.;
00301         outputValues(8, i0, 1) = -2.*(-1. + z)*z;
00302         outputValues(8, i0, 2) =  y*(2. - 4.*z);
00303         outputValues(8, i0, 3) = -4.*(-1. + z)*z;
00304         outputValues(8, i0, 4) = -2.*(-1. + x + 2.*y)*(-1. + 2.*z);     
00305         outputValues(8, i0, 5) = -4.*y*(-1. + x + y);
00306 
00307         outputValues(9, i0, 0) =  4. - 4.*z*z;
00308         outputValues(9, i0, 1) =  4. - 4.*z*z;
00309         outputValues(9, i0, 2) = -2.*(-3. + 4.*x + 4.*y)*z;     
00310         outputValues(9, i0, 3) =  4. - 4.*z*z;
00311         outputValues(9, i0, 4) = -2.*(-3. + 4.*x + 4.*y)*z;
00312         outputValues(9, i0, 5) = -2.*(-1. + x + y)*(-1. + 2.*x + 2.*y);     
00313          
00314         outputValues(10,i0, 0) =  4. - 4.*z*z;
00315         outputValues(10,i0, 1) =  0.;
00316         outputValues(10,i0, 2) =  (2. - 8.*x)*z;
00317         outputValues(10,i0, 3) =  0.;
00318         outputValues(10,i0, 4) =  0.;
00319         outputValues(10,i0, 5) = -2.*x*(-1. + 2.*x);     
00320         
00321         outputValues(11,i0, 0) =  0.;
00322         outputValues(11,i0, 1) =  0.;
00323         outputValues(11,i0, 2) =  0.;
00324         outputValues(11,i0, 3) =  4. - 4.*z*z;
00325         outputValues(11,i0, 4) =  (2. - 8.*y)*z;
00326         outputValues(11,i0, 5) = -2.*y*(-1. + 2.*y);     
00327         
00328         outputValues(12,i0, 0) = -4.*z*(1. + z);
00329         outputValues(12,i0, 1) = -2.*z*(1. + z);
00330         outputValues(12,i0, 2) = -2.*(-1. + 2.*x + y)*(1. + 2.*z);
00331         outputValues(12,i0, 3) =  0.;
00332         outputValues(12,i0, 4) = -2.*(x + 2.*x*z);     
00333         outputValues(12,i0, 5) = -4.*x*(-1. + x + y);     
00334         
00335         outputValues(13,i0, 0) =  0.;
00336         outputValues(13,i0, 1) =  2.*z*(1. + z);
00337         outputValues(13,i0, 2) =  2.*(y + 2.*y*z);
00338         outputValues(13,i0, 3) =  0.;
00339         outputValues(13,i0, 4) =  2.*(x + 2.*x*z);
00340         outputValues(13,i0, 5) =  4.*x*y;     
00341         
00342         outputValues(14,i0, 0) =  0.;
00343         outputValues(14,i0, 1) = -2.*z*(1. + z);
00344         outputValues(14,i0, 2) = -2.*(y + 2.*y*z);
00345         outputValues(14,i0, 3) = -4.*z*(1. + z);
00346         outputValues(14,i0, 4) = -2.*(-1. + x + 2.*y)*(1. + 2.*z);     
00347         outputValues(14,i0, 5) = -4.*y*(-1. + x + y);     
00348         
00349         outputValues(15,i0, 0) =  8.*(-1. + z*z);
00350         outputValues(15,i0, 1) =  4.*(-1. + z*z);
00351         outputValues(15,i0, 2) =  8.*(-1. + 2.*x + y)*z;
00352         outputValues(15,i0, 3) =  0.;     
00353         outputValues(15,i0, 4) =  8.*x*z;
00354         outputValues(15,i0, 5) =  8.*x*(-1. + x + y);     
00355         
00356         outputValues(16,i0, 0) =  0.;
00357         outputValues(16,i0, 1) =  4. - 4.*z*z;
00358         outputValues(16,i0, 2) = -8.*y*z;
00359         outputValues(16,i0, 3) =  0.;
00360         outputValues(16,i0, 4) = -8.*x*z;
00361         outputValues(16,i0, 5) = -8.*x*y;     
00362          
00363         
00364         outputValues(17,i0, 0) =  0.;
00365         outputValues(17,i0, 1) =  4.*(-1. + z*z);
00366         outputValues(17,i0, 2) =  8.*y*z;
00367         outputValues(17,i0, 3) =  8.*(-1. + z*z);
00368         outputValues(17,i0, 4) =  8.*(-1. + x + 2.*y)*z;     
00369         outputValues(17,i0, 5) =  8.*y*(-1. + x + y);
00370       }
00371       break;
00372       
00373     case OPERATOR_D3:
00374       for (int i0 = 0; i0 < dim0; i0++) {
00375         x = inputPoints(i0,0);
00376         y = inputPoints(i0,1);
00377         z = inputPoints(i0,2);
00378 
00379         outputValues(0, i0, 0) =  0.;
00380         outputValues(0, i0, 1) =  0.;
00381         outputValues(0, i0, 2) = -2. + 4.*z;
00382         outputValues(0, i0, 3) =  0.;
00383         outputValues(0, i0, 4) = -2. + 4.*z;
00384         outputValues(0, i0, 5) = -3. + 4.*x + 4.*y;
00385         outputValues(0, i0, 6) =  0.;
00386         outputValues(0, i0, 7) = -2. + 4.*z;
00387         outputValues(0, i0, 8) = -3. + 4.*x + 4.*y;
00388         outputValues(0, i0, 9) =  0.;  
00389           
00390         outputValues(1, i0, 0) =  0.;
00391         outputValues(1, i0, 1) =  0.;
00392         outputValues(1, i0, 2) = -2. + 4.*z;
00393         outputValues(1, i0, 3) =  0.;
00394         outputValues(1, i0, 4) =  0.;
00395         outputValues(1, i0, 5) = -1 + 4.*x;
00396         outputValues(1, i0, 6) =  0.;
00397         outputValues(1, i0, 7) =  0.;
00398         outputValues(1, i0, 8) =  0.;
00399         outputValues(1, i0, 9) =  0.;  
00400         
00401         outputValues(2, i0, 0) =  0.;
00402         outputValues(2, i0, 1) =  0.;
00403         outputValues(2, i0, 2) =  0.;
00404         outputValues(2, i0, 3) =  0.;
00405         outputValues(2, i0, 4) =  0.;
00406         outputValues(2, i0, 5) =  0.;
00407         outputValues(2, i0, 6) =  0.;
00408         outputValues(2, i0, 7) = -2. + 4.*z;
00409         outputValues(2, i0, 8) = -1 + 4.*y;
00410         outputValues(2, i0, 9) =  0.;  
00411         
00412         outputValues(3, i0, 0) =  0.;
00413         outputValues(3, i0, 1) =  0.;
00414         outputValues(3, i0, 2) =  2. + 4.*z;
00415         outputValues(3, i0, 3) =  0.;
00416         outputValues(3, i0, 4) =  2. + 4.*z;
00417         outputValues(3, i0, 5) = -3. + 4.*x + 4.*y;
00418         outputValues(3, i0, 6) =  0.;
00419         outputValues(3, i0, 7) =  2. + 4.*z;
00420         outputValues(3, i0, 8) = -3. + 4.*x + 4.*y;
00421         outputValues(3, i0, 9) =  0.;  
00422         
00423         outputValues(4, i0, 0) =  0.;
00424         outputValues(4, i0, 1) =  0.;
00425         outputValues(4, i0, 2) =  2. + 4.*z;
00426         outputValues(4, i0, 3) =  0.;
00427         outputValues(4, i0, 4) =  0.;
00428         outputValues(4, i0, 5) = -1 + 4.*x;
00429         outputValues(4, i0, 6) =  0.;
00430         outputValues(4, i0, 7) =  0.;
00431         outputValues(4, i0, 8) =  0.;
00432         outputValues(4, i0, 9) =  0.;  
00433         
00434         outputValues(5, i0, 0) =  0.;
00435         outputValues(5, i0, 1) =  0.;
00436         outputValues(5, i0, 2) =  0.;
00437         outputValues(5, i0, 3) =  0.;
00438         outputValues(5, i0, 4) =  0.;
00439         outputValues(5, i0, 5) =  0.;
00440         outputValues(5, i0, 6) =  0.;
00441         outputValues(5, i0, 7) =  2. + 4.*z;
00442         outputValues(5, i0, 8) = -1 + 4.*y;
00443         outputValues(5, i0, 9) =  0.;  
00444         
00445         outputValues(6, i0, 0) =  0.;
00446         outputValues(6, i0, 1) =  0.;
00447         outputValues(6, i0, 2) =  4. - 8.*z;
00448         outputValues(6, i0, 3) =  0.;
00449         outputValues(6, i0, 4) =  2. - 4.*z;
00450         outputValues(6, i0, 5) = -4.*(-1 + 2*x + y);
00451         outputValues(6, i0, 6) =  0.;
00452         outputValues(6, i0, 7) =  0.;
00453         outputValues(6, i0, 8) = -4.*x;
00454         outputValues(6, i0, 9) =  0.;  
00455         
00456         outputValues(7, i0, 0) =  0.;
00457         outputValues(7, i0, 1) =  0.;
00458         outputValues(7, i0, 2) =  0.;
00459         outputValues(7, i0, 3) =  0.;
00460         outputValues(7, i0, 4) = -2. + 4.*z;
00461         outputValues(7, i0, 5) =  4.*y;
00462         outputValues(7, i0, 6) =  0.;
00463         outputValues(7, i0, 7) =  0.;
00464         outputValues(7, i0, 8) =  4.*x;
00465         outputValues(7, i0, 9) =  0.;  
00466         
00467         outputValues(8, i0, 0) =  0.;
00468         outputValues(8, i0, 1) =  0.;
00469         outputValues(8, i0, 2) =  0.;
00470         outputValues(8, i0, 3) =  0.;
00471         outputValues(8, i0, 4) =  2. - 4.*z;
00472         outputValues(8, i0, 5) = -4.*y;
00473         outputValues(8, i0, 6) =  0.;
00474         outputValues(8, i0, 7) =  4. - 8.*z;
00475         outputValues(8, i0, 8) = -4.*(-1 + x + 2*y);
00476         outputValues(8, i0, 9) =  0.;  
00477         
00478         outputValues(9, i0, 0) =  0.;
00479         outputValues(9, i0, 1) =  0.;
00480         outputValues(9, i0, 2) = -8.*z;
00481         outputValues(9, i0, 3) =  0.;
00482         outputValues(9, i0, 4) = -8.*z;
00483         outputValues(9, i0, 5) =  6. - 8.*x - 8.*y;
00484         outputValues(9, i0, 6) =  0.;
00485         outputValues(9, i0, 7) = -8.*z;
00486         outputValues(9, i0, 8) =  6. - 8.*x - 8.*y;
00487         outputValues(9, i0, 9) =  0.;  
00488         
00489         outputValues(10,i0, 0) =  0.;
00490         outputValues(10,i0, 1) =  0.;
00491         outputValues(10,i0, 2) = -8.*z;
00492         outputValues(10,i0, 3) =  0.;
00493         outputValues(10,i0, 4) =  0.;
00494         outputValues(10,i0, 5) =  2. - 8.*x;
00495         outputValues(10,i0, 6) =  0.;
00496         outputValues(10,i0, 7) =  0.;
00497         outputValues(10,i0, 8) =  0.;
00498         outputValues(10,i0, 9) =  0.;  
00499         
00500         outputValues(11,i0, 0) =  0.;
00501         outputValues(11,i0, 1) =  0.;
00502         outputValues(11,i0, 2) =  0.;
00503         outputValues(11,i0, 3) =  0.;
00504         outputValues(11,i0, 4) =  0.;
00505         outputValues(11,i0, 5) =  0.;
00506         outputValues(11,i0, 6) =  0.;
00507         outputValues(11,i0, 7) = -8.*z;
00508         outputValues(11,i0, 8) =  2. - 8.*y;
00509         outputValues(11,i0, 9) =  0.;  
00510         
00511         outputValues(12,i0, 0) =  0.;
00512         outputValues(12,i0, 1) =  0.;
00513         outputValues(12,i0, 2) = -4. - 8.*z;
00514         outputValues(12,i0, 3) =  0.;
00515         outputValues(12,i0, 4) = -2. - 4.*z;
00516         outputValues(12,i0, 5) = -4.*(-1 + 2*x + y);
00517         outputValues(12,i0, 6) =  0.;
00518         outputValues(12,i0, 7) =  0.;
00519         outputValues(12,i0, 8) = -4.*x;
00520         outputValues(12,i0, 9) =  0.;  
00521         
00522         outputValues(13,i0, 0) =  0.;
00523         outputValues(13,i0, 1) =  0.;
00524         outputValues(13,i0, 2) =  0.;
00525         outputValues(13,i0, 3) =  0.;
00526         outputValues(13,i0, 4) =  2. + 4.*z;
00527         outputValues(13,i0, 5) =  4.*y;
00528         outputValues(13,i0, 6) =  0.;
00529         outputValues(13,i0, 7) =  0.;
00530         outputValues(13,i0, 8) =  4.*x;
00531         outputValues(13,i0, 9) =  0.;  
00532         
00533         outputValues(14,i0, 0) =  0.;
00534         outputValues(14,i0, 1) =  0.;
00535         outputValues(14,i0, 2) =  0.;
00536         outputValues(14,i0, 3) =  0.;
00537         outputValues(14,i0, 4) = -2. - 4.*z;
00538         outputValues(14,i0, 5) = -4.*y;
00539         outputValues(14,i0, 6) =  0.;
00540         outputValues(14,i0, 7) = -4. - 8.*z;
00541         outputValues(14,i0, 8) = -4.*(-1 + x + 2*y);
00542         outputValues(14,i0, 9) =  0.;  
00543         
00544         outputValues(15,i0, 0) =  0.;
00545         outputValues(15,i0, 1) =  0.;
00546         outputValues(15,i0, 2) =  16.*z;
00547         outputValues(15,i0, 3) =  0.;
00548         outputValues(15,i0, 4) =  8.*z;
00549         outputValues(15,i0, 5) =  8.*(-1 + 2*x + y);
00550         outputValues(15,i0, 6) =  0.;
00551         outputValues(15,i0, 7) =  0.;
00552         outputValues(15,i0, 8) =  8.*x;
00553         outputValues(15,i0, 9) =  0.;  
00554         
00555         outputValues(16,i0, 0) =  0.;
00556         outputValues(16,i0, 1) =  0.;
00557         outputValues(16,i0, 2) =  0.;
00558         outputValues(16,i0, 3) =  0.;
00559         outputValues(16,i0, 4) = -8.*z;
00560         outputValues(16,i0, 5) = -8.*y;
00561         outputValues(16,i0, 6) =  0.;
00562         outputValues(16,i0, 7) =  0.;
00563         outputValues(16,i0, 8) = -8.*x;
00564         outputValues(16,i0, 9) =  0.;  
00565         
00566         outputValues(17,i0, 0) =  0.;
00567         outputValues(17,i0, 1) =  0.;
00568         outputValues(17,i0, 2) =  0.;
00569         outputValues(17,i0, 3) =  0.;
00570         outputValues(17,i0, 4) =  8.*z;
00571         outputValues(17,i0, 5) =  8.*y;
00572         outputValues(17,i0, 6) =  0.;
00573         outputValues(17,i0, 7) =  16.*z;
00574         outputValues(17,i0, 8) =  8.*(-1 + x + 2*y);
00575         outputValues(17,i0, 9) =  0.;
00576         
00577       }
00578       break;
00579       
00580     case OPERATOR_D4:
00581       {
00582         // There are only few constant non-zero entries. Initialize by zero and then assign non-zero entries.
00583         int DkCardinality = Intrepid::getDkCardinality(operatorType, this -> basisCellTopology_.getDimension() );
00584         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00585           for (int i0 = 0; i0 < dim0; i0++) {
00586             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00587               outputValues(dofOrd, i0, dkOrd) = 0.0;
00588             }
00589           }
00590         }    
00591         
00592         for (int i0 = 0; i0 < dim0; i0++) {
00593           
00594           outputValues(0, i0, 5) = 4.;
00595           outputValues(0, i0, 8) = 4.;
00596           outputValues(0, i0,12) = 4.;
00597           
00598           outputValues(1, i0, 5) = 4.;
00599           
00600           outputValues(2, i0,12) = 4.;
00601           
00602           outputValues(3, i0, 5) = 4.;
00603           outputValues(3, i0, 8) = 4.;
00604           outputValues(3, i0,12) = 4.;
00605           
00606           outputValues(4, i0, 5) = 4.0;
00607           
00608           outputValues(5, i0,12) = 4.0;
00609           
00610           outputValues(6, i0, 5) =-8.;
00611           outputValues(6, i0, 8) =-4.;
00612           
00613           outputValues(7, i0, 8) = 4.;
00614           
00615           outputValues(8, i0, 8) =-4.;
00616           outputValues(8, i0,12) =-8.;
00617           
00618           outputValues(9, i0, 5) =-8.;
00619           outputValues(9, i0, 8) =-8.;
00620           outputValues(9, i0,12) =-8.;
00621           
00622           outputValues(10,i0, 5) =-8.;
00623           
00624           outputValues(11,i0,12) =-8.;
00625           
00626           outputValues(12,i0, 5) =-8.;
00627           outputValues(12,i0, 8) =-4.;
00628           
00629           outputValues(13,i0, 8) = 4.;
00630           
00631           outputValues(14,i0, 8) =-4;
00632           outputValues(14,i0,12) =-8.;
00633           
00634           outputValues(15,i0, 5) =16.;
00635           outputValues(15,i0, 8) = 8.;
00636           
00637           outputValues(16,i0, 8) =-8.;
00638           
00639           
00640           outputValues(17,i0, 8) = 8.;
00641           outputValues(17,i0,12) =16.;   
00642         }
00643       }
00644       break;
00645       
00646     case OPERATOR_D5:
00647     case OPERATOR_D6:
00648     case OPERATOR_D7:
00649     case OPERATOR_D8:
00650     case OPERATOR_D9:
00651     case OPERATOR_D10:
00652       {
00653         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
00654         int DkCardinality = Intrepid::getDkCardinality(operatorType, 
00655                                                        this -> basisCellTopology_.getDimension() );
00656         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00657           for (int i0 = 0; i0 < dim0; i0++) {
00658             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00659               outputValues(dofOrd, i0, dkOrd) = 0.0;
00660             }
00661           }
00662         }
00663       }
00664       break;
00665       
00666     default:
00667       TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00668                           ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): Invalid operator type");
00669   }
00670 }
00671 
00672 
00673   
00674 template<class Scalar, class ArrayScalar>
00675 void Basis_HGRAD_WEDGE_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&       outputValues,
00676                                                              const ArrayScalar &    inputPoints,
00677                                                              const ArrayScalar &    cellVertices,
00678                                                              const EOperator        operatorType) const {
00679   TEST_FOR_EXCEPTION( (true), std::logic_error,
00680                       ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): FEM Basis calling an FVD member function");
00681 }
00682 }// namespace Intrepid
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_WEDGE_C2_FEMDef.hpp
