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