Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_HEX_C2_FEMDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_HEX_C2_FEMDEF_HPP
00002 #define INTREPID_HGRAD_HEX_C2_FEMDEF_HPP
00003 // @HEADER
00004 // ************************************************************************
00005 //
00006 //                           Intrepid Package
00007 //                 Copyright (2007) Sandia Corporation
00008 //
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 //
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00040 //                    Denis Ridzal  (dridzal@sandia.gov), or
00041 //                    Kara Peterson (kjpeter@sandia.gov)
00042 //
00043 // ************************************************************************
00044 // @HEADER
00045 
00051 namespace Intrepid {
00052 
00053   
00054 template<class Scalar, class ArrayScalar>
00055 Basis_HGRAD_HEX_C2_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_C2_FEM()
00056   {
00057     this -> basisCardinality_  = 27;
00058     this -> basisDegree_       = 2;    
00059     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00060     this -> basisType_         = BASIS_FEM_DEFAULT;
00061     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00062     this -> basisTagsAreSet_   = false;
00063   }
00064   
00065   
00066   
00067 template<class Scalar, class ArrayScalar>
00068 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
00069   
00070   // Basis-dependent intializations
00071   int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00072   int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00073   int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00074   int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00075 
00076   // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 
00077   int tags[]  = { 0, 0, 0, 1,     // Nodes 0 to 7 follow vertex order of the topology
00078                   0, 1, 0, 1,
00079                   0, 2, 0, 1,
00080                   0, 3, 0, 1,
00081                   0, 4, 0, 1,
00082                   0, 5, 0, 1,
00083                   0, 6, 0, 1,
00084                   0, 7, 0, 1,
00085                   1, 0, 0, 1,      // Node 8  -> edge 0
00086                   1, 1, 0, 1,      // Node 9  -> edge 1
00087                   1, 2, 0, 1,      // Node 10 -> edge 2
00088                   1, 3, 0, 1,      // Node 11 -> edge 3
00089                   1, 8, 0, 1,      // Node 12 -> edge 8
00090                   1, 9, 0, 1,      // Node 13 -> edge 9
00091                   1,10, 0, 1,      // Node 14 -> edge 10
00092                   1,11, 0, 1,      // Node 15 -> edge 11
00093                   1, 4, 0, 1,      // Node 16 -> edge 4
00094                   1, 5, 0, 1,      // Node 17 -> edge 5
00095                   1, 6, 0, 1,      // Node 18 -> edge 6
00096                   1, 7, 0, 1,      // Node 19 -> edge 7
00097                   3, 0, 0, 1,      // Node 20 -> Hexahedron
00098                   2, 4, 0, 1,      // Node 21 -> face 4
00099                   2, 5, 0, 1,      // Node 22 -> face 5
00100                   2, 3, 0, 1,      // Node 23 -> face 3
00101                   2, 1, 0, 1,      // Node 24 -> face 1
00102                   2, 0, 0, 1,      // Node 25 -> face 0
00103                   2, 2, 0, 1,      // Node 26 -> face 2
00104   };
00105   
00106   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00107   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00108                               this -> ordinalToTag_,
00109                               tags,
00110                               this -> basisCardinality_,
00111                               tagSize,
00112                               posScDim,
00113                               posScOrd,
00114                               posDfOrd);
00115 }
00116 
00117 
00118 
00119 template<class Scalar, class ArrayScalar>
00120 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00121                                                             const ArrayScalar &  inputPoints,
00122                                                             const EOperator      operatorType) const {
00123   
00124   // Verify arguments
00125 #ifdef HAVE_INTREPID_DEBUG
00126   Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00127                                                       inputPoints,
00128                                                       operatorType,
00129                                                       this -> getBaseCellTopology(),
00130                                                       this -> getCardinality() );
00131 #endif
00132   
00133   // Number of evaluation points = dim 0 of inputPoints
00134   int dim0 = inputPoints.dimension(0);  
00135   
00136   // Temporaries: (x,y,z) coordinates of the evaluation point
00137   Scalar x = 0.0;                                    
00138   Scalar y = 0.0;                                    
00139   Scalar z = 0.0;                                    
00140   
00141   switch (operatorType) {
00142     
00143     case OPERATOR_VALUE:
00144       for (int i0 = 0; i0 < dim0; i0++) {
00145         x = inputPoints(i0, 0);
00146         y = inputPoints(i0, 1);
00147         z = inputPoints(i0, 2);
00148         
00149         // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
00150         outputValues( 0, i0) = 0.125*(-1. + x)*x*(-1. + y)*y*(-1. + z)*z;
00151         outputValues( 1, i0) = 0.125*x*(1.+ x)*(-1. + y)*y*(-1. + z)*z;
00152         outputValues( 2, i0) = 0.125*x*(1.+ x)*y*(1.+ y)*(-1. + z)*z;
00153         outputValues( 3, i0) = 0.125*(-1. + x)*x*y*(1.+ y)*(-1. + z)*z;
00154         outputValues( 4, i0) = 0.125*(-1. + x)*x*(-1. + y)*y*z*(1.+ z);
00155         outputValues( 5, i0) = 0.125*x*(1.+ x)*(-1. + y)*y*z*(1.+ z);
00156         outputValues( 6, i0) = 0.125*x*(1.+ x)*y*(1.+ y)*z*(1.+ z);
00157         outputValues( 7, i0) = 0.125*(-1. + x)*x*y*(1.+ y)*z*(1.+ z);
00158         outputValues( 8, i0) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*(-1. + z)*z;
00159         outputValues( 9, i0) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*(-1. + z)*z;
00160         outputValues(10, i0) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*(-1. + z)*z;
00161         outputValues(11, i0) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*(-1. + z)*z;
00162         outputValues(12, i0) = 0.25*(-1. + x)*x*(-1. + y)*y*(1. - z)*(1. + z);
00163         outputValues(13, i0) = 0.25*x*(1.+ x)*(-1. + y)*y*(1. - z)*(1. + z);
00164         outputValues(14, i0) = 0.25*x*(1.+ x)*y*(1.+ y)*(1. - z)*(1. + z);
00165         outputValues(15, i0) = 0.25*(-1. + x)*x*y*(1.+ y)*(1. - z)*(1. + z);
00166         outputValues(16, i0) = 0.25*(1. - x)*(1. + x)*(-1. + y)*y*z*(1.+ z);
00167         outputValues(17, i0) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*z*(1.+ z);
00168         outputValues(18, i0) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*z*(1.+ z);
00169         outputValues(19, i0) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*z*(1.+ z);
00170         outputValues(20, i0) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00171         outputValues(21, i0) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z;
00172         outputValues(22, i0) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z);
00173         outputValues(23, i0) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00174         outputValues(24, i0) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00175         outputValues(25, i0) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z);
00176         outputValues(26, i0) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z);
00177       }
00178       break;
00179       
00180     case OPERATOR_GRAD:
00181     case OPERATOR_D1:
00182       for (int i0 = 0; i0 < dim0; i0++) {
00183         x = inputPoints(i0,0);
00184         y = inputPoints(i0,1);
00185         z = inputPoints(i0,2);
00186 
00187         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00188         outputValues(0, i0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
00189         outputValues(0, i0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*(-1. + z)*z;
00190         outputValues(0, i0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.125 + 0.25*z);
00191           
00192         outputValues(1, i0, 0) = (0.125 + 0.25*x)*(-1. + y)*y*(-1. + z)*z;
00193         outputValues(1, i0, 1) = x*(1. + x)*(-0.125 + 0.25*y)*(-1. + z)*z;
00194         outputValues(1, i0, 2) = x*(1. + x)*(-1. + y)*y*(-0.125 + 0.25*z);
00195           
00196         outputValues(2, i0, 0) = (0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
00197         outputValues(2, i0, 1) = x*(1. + x)*(0.125 + 0.25*y)*(-1. + z)*z;
00198         outputValues(2, i0, 2) = x*(1. + x)*y*(1. + y)*(-0.125 + 0.25*z);
00199           
00200         outputValues(3, i0, 0) = (-0.125 + 0.25*x)*y*(1. + y)*(-1. + z)*z;
00201         outputValues(3, i0, 1) = (-1. + x)*x*(0.125 + 0.25*y)*(-1. + z)*z;
00202         outputValues(3, i0, 2) = (-1. + x)*x*y*(1. + y)*(-0.125 + 0.25*z);
00203           
00204         outputValues(4, i0, 0) = (-0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
00205         outputValues(4, i0, 1) = (-1. + x)*x*(-0.125 + 0.25*y)*z*(1. + z);
00206         outputValues(4, i0, 2) = (-1. + x)*x*(-1. + y)*y*(0.125 + 0.25*z);
00207           
00208         outputValues(5, i0, 0) = (0.125 + 0.25*x)*(-1. + y)*y*z*(1. + z);
00209         outputValues(5, i0, 1) = x*(1. + x)*(-0.125 + 0.25*y)*z*(1. + z);
00210         outputValues(5, i0, 2) = x*(1. + x)*(-1. + y)*y*(0.125 + 0.25*z);
00211           
00212         outputValues(6, i0, 0) = (0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
00213         outputValues(6, i0, 1) = x*(1. + x)*(0.125 + 0.25*y)*z*(1. + z);
00214         outputValues(6, i0, 2) = x*(1. + x)*y*(1. + y)*(0.125 + 0.25*z);
00215           
00216         outputValues(7, i0, 0) = (-0.125 + 0.25*x)*y*(1. + y)*z*(1. + z);
00217         outputValues(7, i0, 1) = (-1. + x)*x*(0.125 + 0.25*y)*z*(1. + z);
00218         outputValues(7, i0, 2) = (-1. + x)*x*y*(1. + y)*(0.125 + 0.25*z);
00219           
00220         outputValues(8, i0, 0) = -0.5*x*(-1. + y)*y*(-1. + z)*z;
00221         outputValues(8, i0, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*(-1. + z)*z;
00222         outputValues(8, i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-0.25 + 0.5*z);
00223         
00224         outputValues(9, i0, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
00225         outputValues(9, i0, 1) = x*(1. + x)*(-0.5*y)*(-1. + z)*z;
00226         outputValues(9, i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
00227           
00228         outputValues(10,i0, 0) = -0.5*x*y*(1. + y)*(-1. + z)*z;
00229         outputValues(10,i0, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*(-1. + z)*z;
00230         outputValues(10,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-0.25 + 0.5*z);
00231           
00232         outputValues(11,i0, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*(-1. + z)*z;
00233         outputValues(11,i0, 1) = (-1. + x)*x*(-0.5*y)*(-1. + z)*z;
00234         outputValues(11,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-0.25 + 0.5*z);
00235           
00236         outputValues(12,i0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
00237         outputValues(12,i0, 1) = (-1. + x)*x*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
00238         outputValues(12,i0, 2) = (-1. + x)*x*(-1. + y)*y*(-0.5*z);
00239           
00240         outputValues(13,i0, 0) = (0.25 + 0.5*x)*(-1. + y)*y*(1. - z)*(1. + z);
00241         outputValues(13,i0, 1) = x*(1. + x)*(-0.25 + 0.5*y)*(1. - z)*(1. + z);
00242         outputValues(13,i0, 2) = x*(1. + x)*(-1. + y)*y*(-0.5*z);
00243           
00244         outputValues(14,i0, 0) = (0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
00245         outputValues(14,i0, 1) = x*(1. + x)*(0.25 + 0.5*y)*(1. - z)*(1. + z);
00246         outputValues(14,i0, 2) = x*(1. + x)*y*(1. + y)*(-0.5*z);
00247           
00248         outputValues(15,i0, 0) = (-0.25 + 0.5*x)*y*(1. + y)*(1. - z)*(1. + z);
00249         outputValues(15,i0, 1) = (-1. + x)*x*(0.25 + 0.5*y)*(1. - z)*(1. + z);
00250         outputValues(15,i0, 2) = (-1. + x)*x*y*(1. + y)*(-0.5*z);
00251           
00252         outputValues(16,i0, 0) = -0.5*x*(-1. + y)*y*z*(1. + z);
00253         outputValues(16,i0, 1) = (1. - x)*(1. + x)*(-0.25 + 0.5*y)*z*(1. + z);
00254         outputValues(16,i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(0.25 + 0.5*z);
00255           
00256         outputValues(17,i0, 0) = (0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
00257         outputValues(17,i0, 1) = x*(1. + x)*(-0.5*y)*z*(1. + z);
00258         outputValues(17,i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(0.25 + 0.5*z);
00259           
00260         outputValues(18,i0, 0) = -0.5*x*y*(1. + y)*z*(1. + z);
00261         outputValues(18,i0, 1) = (1. - x)*(1. + x)*(0.25 + 0.5*y)*z*(1. + z);
00262         outputValues(18,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(0.25 + 0.5*z);
00263           
00264         outputValues(19,i0, 0) = (-0.25 + 0.5*x)*(1. - y)*(1. + y)*z*(1. + z);
00265         outputValues(19,i0, 1) = (-1. + x)*x*(-0.5*y)*z*(1. + z);
00266         outputValues(19,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(0.25 + 0.5*z);
00267           
00268         outputValues(20,i0, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00269         outputValues(20,i0, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z);
00270         outputValues(20,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z);
00271           
00272         outputValues(21,i0, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z;
00273         outputValues(21,i0, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z;
00274         outputValues(21,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z);
00275           
00276         outputValues(22,i0, 0) = -x*(1. - y)*(1. + y)*z*(1. + z);
00277         outputValues(22,i0, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z);
00278         outputValues(22,i0, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z);
00279           
00280         outputValues(23,i0, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00281         outputValues(23,i0, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z);
00282         outputValues(23,i0, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z);
00283           
00284         outputValues(24,i0, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z);
00285         outputValues(24,i0, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z);
00286         outputValues(24,i0, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z);
00287           
00288         outputValues(25,i0, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z);
00289         outputValues(25,i0, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z);
00290         outputValues(25,i0, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z);
00291           
00292         outputValues(26,i0, 0) = -x*y*(1. + y)*(1. - z)*(1. + z);
00293         outputValues(26,i0, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z);
00294         outputValues(26,i0, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z);
00295       }
00296       break;
00297       
00298     case OPERATOR_CURL:
00299       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
00300                           ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
00301       break;
00302       
00303     case OPERATOR_DIV:
00304       TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00305                           ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
00306       break;
00307       
00308     case OPERATOR_D2:
00309       for (int i0 = 0; i0 < dim0; i0++) {
00310         x = inputPoints(i0,0);
00311         y = inputPoints(i0,1);
00312         z = inputPoints(i0,2);
00313         
00314         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6) 
00315         outputValues(0, i0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z; 
00316         outputValues(0, i0, 1) = (-0.125 + y*(0.25 - 0.25*z) + x*(0.25 + y*(-0.5 + 0.5*z) - 0.25*z) + 0.125*z)*z; 
00317         outputValues(0, i0, 2) = y*(-0.125 + x*(0.25 + y*(-0.25 + 0.5*z) - 0.5*z) + y*(0.125 - 0.25*z) + 0.25*z); 
00318         outputValues(0, i0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z; 
00319         outputValues(0, i0, 4) = x*(-0.125 + y*(0.25 - 0.5*z) + x*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z) + 0.25*z); 
00320         outputValues(0, i0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y; 
00321         
00322         outputValues(1, i0, 0) = 0.25*(-1. + y)*y*(-1. + z)*z; 
00323         outputValues(1, i0, 1) = (0.125 + x*(0.25 + y*(-0.5 + 0.5*z) - 0.25*z) + y*(-0.25 + 0.25*z) - 0.125*z)*z; 
00324         outputValues(1, i0, 2) = y*(0.125 + x*(0.25 + y*(-0.25 + 0.5*z) - 0.5*z) + y*(-0.125 + 0.25*z) - 0.25*z); 
00325         outputValues(1, i0, 3) = 0.25*x*(1 + x)*(-1. + z)*z; 
00326         outputValues(1, i0, 4) = x*(1. + x)*(0.125 + y*(-0.25 + 0.5*z) - 0.25*z); 
00327         outputValues(1, i0, 5) = 0.25*x*(1 + x)*(-1. + y)*y; 
00328         
00329         outputValues(2, i0, 0) = 0.25*y*(1 + y)*(-1. + z)*z; 
00330         outputValues(2, i0, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*(-1. + z)*z; 
00331         outputValues(2, i0, 2) = y*(1. + y)*(-0.125 + x*(-0.25 + 0.5*z) + 0.25*z); 
00332         outputValues(2, i0, 3) = 0.25*x*(1 + x)*(-1. + z)*z; 
00333         outputValues(2, i0, 4) = x*(1. + x)*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z); 
00334         outputValues(2, i0, 5) = 0.25*x*(1 + x)*y*(1 + y); 
00335         
00336         outputValues(3, i0, 0) = 0.25*y*(1 + y)*(-1. + z)*z; 
00337         outputValues(3, i0, 1) = (0.125 + y*(0.25 - 0.25*z) + x*(-0.25 + y*(-0.5 + 0.5*z) + 0.25*z) - 0.125*z)*z; 
00338         outputValues(3, i0, 2) = y*(1. + y)*(0.125 + x*(-0.25 + 0.5*z) - 0.25*z); 
00339         outputValues(3, i0, 3) = 0.25*(-1. + x)*x*(-1. + z)*z; 
00340         outputValues(3, i0, 4) = x*(0.125 + y*(0.25 - 0.5*z) + x*(-0.125 + y*(-0.25 + 0.5*z) + 0.25*z) - 0.25*z); 
00341         outputValues(3, i0, 5) = 0.25*(-1. + x)*x*y*(1 + y); 
00342         
00343         outputValues(4, i0, 0) = 0.25*(-1. + y)*y*z*(1 + z); 
00344         outputValues(4, i0, 1) = (0.125 + x*(-0.25 + 0.5*y) - 0.25*y)*z*(1. + z); 
00345         outputValues(4, i0, 2) = y*(0.125 + x*(-0.25 + y*(0.25 + 0.5*z) - 0.5*z) + y*(-0.125 - 0.25*z) + 0.25*z); 
00346         outputValues(4, i0, 3) = 0.25*(-1. + x)*x*z*(1 + z); 
00347         outputValues(4, i0, 4) = x*(0.125 + y*(-0.25 - 0.5*z) + x*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z) + 0.25*z); 
00348         outputValues(4, i0, 5) = 0.25*(-1. + x)*x*(-1. + y)*y; 
00349         
00350         outputValues(5, i0, 0) = 0.25*(-1. + y)*y*z*(1 + z); 
00351         outputValues(5, i0, 1) = (-0.125 + x*(-0.25 + 0.5*y) + 0.25*y)*z*(1. + z); 
00352         outputValues(5, i0, 2) = (-1. + y)*y*(0.125 + x*(0.25 + 0.5*z) + 0.25*z); 
00353         outputValues(5, i0, 3) = 0.25*x*(1 + x)*z*(1 + z); 
00354         outputValues(5, i0, 4) = x*(1. + x)*(-0.125 + y*(0.25 + 0.5*z) - 0.25*z); 
00355         outputValues(5, i0, 5) = 0.25*x*(1 + x)*(-1. + y)*y; 
00356         
00357         outputValues(6, i0, 0) = 0.25*y*(1 + y)*z*(1 + z); 
00358         outputValues(6, i0, 1) = (0.125 + x*(0.25 + 0.5*y) + 0.25*y)*z*(1. + z); 
00359         outputValues(6, i0, 2) = y*(1. + y)*(0.125 + x*(0.25 + 0.5*z) + 0.25*z); 
00360         outputValues(6, i0, 3) = 0.25*x*(1 + x)*z*(1 + z); 
00361         outputValues(6, i0, 4) = x*(1. + x)*(0.125 + y*(0.25 + 0.5*z) + 0.25*z); 
00362         outputValues(6, i0, 5) = 0.25*x*(1 + x)*y*(1 + y); 
00363         
00364         outputValues(7, i0, 0) = 0.25*y*(1 + y)*z*(1 + z); 
00365         outputValues(7, i0, 1) = (-0.125 + x*(0.25 + 0.5*y) - 0.25*y)*z*(1. + z); 
00366         outputValues(7, i0, 2) = y*(1. + y)*(-0.125 + x*(0.25 + 0.5*z) - 0.25*z); 
00367         outputValues(7, i0, 3) = 0.25*(-1. + x)*x*z*(1 + z); 
00368         outputValues(7, i0, 4) = (-1. + x)*x*(0.125 + y*(0.25 + 0.5*z) + 0.25*z); 
00369         outputValues(7, i0, 5) = 0.25*(-1. + x)*x*y*(1 + y); 
00370         
00371         outputValues(8, i0, 0) = -0.5*(-1. + y)*y*(-1. + z)*z; 
00372         outputValues(8, i0, 1) = (0.  + x*(-0.5 + y))*z + (x*(0.5 - y) )*(z*z); 
00373         outputValues(8, i0, 2) = (y*y)*(x*(0.5 - z) ) + y*(x*(-0.5 + z)); 
00374         outputValues(8, i0, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z; 
00375         outputValues(8, i0, 4) = 0.25 + (x*x)*(-0.25 + y*(0.5 - z) + 0.5*z) - 0.5*z + y*(-0.5 + z); 
00376         outputValues(8, i0, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y; 
00377         
00378         outputValues(9, i0, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z; 
00379         outputValues(9, i0, 1) = (0.5*y + x*(y))*z + (x*(-y) - 0.5*y)*(z*z); 
00380         outputValues(9, i0, 2) = -0.25 + (y*y)*(0.25 - 0.5*z) + 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z); 
00381         outputValues(9, i0, 3) = -0.5*x*(1 + x)*(-1. + z)*z; 
00382         outputValues(9, i0, 4) = x*(y*(0.5 - z) ) + (x*x)*(y*(0.5 - z) ); 
00383         outputValues(9, i0, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y); 
00384         
00385         outputValues(10,i0, 0) = -0.5*y*(1 + y)*(-1. + z)*z; 
00386         outputValues(10,i0, 1) = (0.  + x*(0.5 + y))*z + (x*(-0.5 - y) )*(z*z); 
00387         outputValues(10,i0, 2) = y*(x*(0.5 - z) ) + (y*y)*(x*(0.5 - z) ); 
00388         outputValues(10,i0, 3) = 0.5*(1. - x)*(1. + x)*(-1. + z)*z; 
00389         outputValues(10,i0, 4) = -0.25 + (x*x)*(0.25 + y*(0.5 - z) - 0.5*z) + 0.5*z + y*(-0.5 + z); 
00390         outputValues(10,i0, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y); 
00391         
00392         outputValues(11,i0, 0) = 0.5*(1. - y)*(1. + y)*(-1. + z)*z; 
00393         outputValues(11,i0, 1) = (-0.5*y + x*(y))*z + (x*(-y) + 0.5*y)*(z*z); 
00394         outputValues(11,i0, 2) = 0.25 + (y*y)*(-0.25 + 0.5*z) - 0.5*z + x*(-0.5 + (y*y)*(0.5 - z) + z); 
00395         outputValues(11,i0, 3) = -0.5*(-1. + x)*x*(-1. + z)*z; 
00396         outputValues(11,i0, 4) = (x*x)*(y*(0.5 - z) ) + x*(y*(-0.5 + z)); 
00397         outputValues(11,i0, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y); 
00398         
00399         outputValues(12,i0, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z); 
00400         outputValues(12,i0, 1) = 0.25  - 0.25*(z*z) + y*(-0.5  + 0.5*(z*z)) + x*(-0.5  + 0.5*(z*z) + y*(1.  - (z*z))); 
00401         outputValues(12,i0, 2) = (y*y)*(x*(-z) + 0.5*z) + y*(-0.5*z + x*(z)); 
00402         outputValues(12,i0, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z); 
00403         outputValues(12,i0, 4) = (x*x)*(y*(-z) + 0.5*z) + x*(-0.5*z + y*(z)); 
00404         outputValues(12,i0, 5) = -0.5*(-1. + x)*x*(-1. + y)*y; 
00405         
00406         outputValues(13,i0, 0) = 0.5*(-1. + y)*y*(1. - z)*(1. + z); 
00407         outputValues(13,i0, 1) = -0.25  + 0.25*(z*z) + y*(0.5  - 0.5*(z*z)) + x*(-0.5  + 0.5*(z*z) + y*(1.  - (z*z))); 
00408         outputValues(13,i0, 2) = (y*y)*(x*(-z) - 0.5*z) + y*(0.5*z + x*(z)); 
00409         outputValues(13,i0, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z); 
00410         outputValues(13,i0, 4) = x*(y*(-z) + 0.5*z) + (x*x)*(y*(-z) + 0.5*z); 
00411         outputValues(13,i0, 5) = -0.5*x*(1 + x)*(-1. + y)*y; 
00412 
00413         outputValues(14,i0, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z); 
00414         outputValues(14,i0, 1) = 0.25  - 0.25*(z*z) + y*(0.5  - 0.5*(z*z)) + x*(0.5  - 0.5*(z*z) + y*(1.  - (z*z))); 
00415         outputValues(14,i0, 2) = y*(x*(-z) - 0.5*z) + (y*y)*(x*(-z) - 0.5*z); 
00416         outputValues(14,i0, 3) = 0.5*x*(1 + x)*(1. - z)*(1. + z); 
00417         outputValues(14,i0, 4) = x*(y*(-z) - 0.5*z) + (x*x)*(y*(-z) - 0.5*z); 
00418         outputValues(14,i0, 5) = -0.5*x*(1 + x)*y*(1 + y); 
00419         
00420         outputValues(15,i0, 0) = 0.5*y*(1 + y)*(1. - z)*(1. + z); 
00421         outputValues(15,i0, 1) = -0.25  + 0.25*(z*z) + y*(-0.5  + 0.5*(z*z)) + x*(0.5  - 0.5*(z*z) + y*(1.  - (z*z))); 
00422         outputValues(15,i0, 2) = y*(x*(-z) + 0.5*z) + (y*y)*(x*(-z) + 0.5*z); 
00423         outputValues(15,i0, 3) = 0.5*(-1. + x)*x*(1. - z)*(1. + z); 
00424         outputValues(15,i0, 4) = (x*x)*(y*(-z) - 0.5*z) + x*(0.5*z + y*(z)); 
00425         outputValues(15,i0, 5) = -0.5*(-1. + x)*x*y*(1 + y); 
00426         
00427         outputValues(16,i0, 0) = -0.5*(-1. + y)*y*z*(1 + z); 
00428         outputValues(16,i0, 1) = (x*(0.5 - y) )*z + (x*(0.5 - y) )*(z*z); 
00429         outputValues(16,i0, 2) = (y*y)*(x*(-0.5 - z) ) + y*(x*(0.5 + z)); 
00430         outputValues(16,i0, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z); 
00431         outputValues(16,i0, 4) = -0.25 + (x*x)*(0.25 + y*(-0.5 - z) + 0.5*z) - 0.5*z + y*(0.5 + z); 
00432         outputValues(16,i0, 5) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y; 
00433         
00434         outputValues(17,i0, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z); 
00435         outputValues(17,i0, 1) = (x*(-y) - 0.5*y)*z + (x*(-y) - 0.5*y)*(z*z); 
00436         outputValues(17,i0, 2) = 0.25 + (y*y)*(-0.25 - 0.5*z) + 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z); 
00437         outputValues(17,i0, 3) = -0.5*x*(1 + x)*z*(1 + z); 
00438         outputValues(17,i0, 4) = x*(y*(-0.5 - z) ) + (x*x)*(y*(-0.5 - z) ); 
00439         outputValues(17,i0, 5) = 0.5*x*(1 + x)*(1. - y)*(1. + y); 
00440         
00441         outputValues(18,i0, 0) = -0.5*y*(1 + y)*z*(1 + z); 
00442         outputValues(18,i0, 1) = (x*(-0.5 - y) )*z + (x*(-0.5 - y) )*(z*z); 
00443         outputValues(18,i0, 2) = y*(x*(-0.5 - z) ) + (y*y)*(x*(-0.5 - z) ); 
00444         outputValues(18,i0, 3) = 0.5*(1. - x)*(1. + x)*z*(1 + z); 
00445         outputValues(18,i0, 4) = 0.25 + (x*x)*(-0.25 + y*(-0.5 - z) - 0.5*z) + 0.5*z + y*(0.5 + z); 
00446         outputValues(18,i0, 5) = 0.5*(1. - x)*(1. + x)*y*(1 + y); 
00447         
00448         outputValues(19,i0, 0) = 0.5*(1. - y)*(1. + y)*z*(1 + z); 
00449         outputValues(19,i0, 1) = (x*(-y) + 0.5*y)*z + (x*(-y) + 0.5*y)*(z*z); 
00450         outputValues(19,i0, 2) = -0.25 + (y*y)*(0.25 + 0.5*z) - 0.5*z + x*(0.5 + (y*y)*(-0.5 - z) + z); 
00451         outputValues(19,i0, 3) = -0.5*(-1. + x)*x*z*(1 + z); 
00452         outputValues(19,i0, 4) = (x*x)*(y*(-0.5 - z) ) + x*(y*(0.5 + z)); 
00453         outputValues(19,i0, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y); 
00454         
00455         outputValues(20,i0, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z); 
00456         outputValues(20,i0, 1) = -4.*x*y*(-1. + z*z);         
00457         outputValues(20,i0, 2) = x*((y*y)*(-4.*z) + 4.*z); 
00458         outputValues(20,i0, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z); 
00459         outputValues(20,i0, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z); 
00460         outputValues(20,i0, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y); 
00461         
00462         outputValues(21,i0, 0) = -(1. - y)*(1. + y)*(-1. + z)*z; 
00463         outputValues(21,i0, 1) = (x*(-2.*y) )*z + (0.  + x*(2.*y))*(z*z); 
00464         outputValues(21,i0, 2) =  x*(1. - 2.*z + (y*y)*(-1. + 2.*z)); 
00465         outputValues(21,i0, 3) = -(1. - x)*(1. + x)*(-1. + z)*z; 
00466         outputValues(21,i0, 4) = y*(1. - 2.*z)  + (x*x)*(y*(-1. + 2.*z)); 
00467         outputValues(21,i0, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); 
00468         
00469         outputValues(22,i0, 0) = -(1. - y)*(1. + y)*z*(1 + z); 
00470         outputValues(22,i0, 1) = (0.  + x*(2.*y))*z + (0.  + x*(2.*y))*(z*z); 
00471         outputValues(22,i0, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z)); 
00472         outputValues(22,i0, 3) = -(1. - x)*(1. + x)*z*(1 + z); 
00473         outputValues(22,i0, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z)); 
00474         outputValues(22,i0, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); 
00475         
00476         outputValues(23,i0, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); 
00477         outputValues(23,i0, 1) = (-1. + 2.*x)*y*(-1. + z*z);
00478         outputValues(23,i0, 2) = (-1. + 2.*x)*(-1. + y*y)*z; 
00479         outputValues(23,i0, 3) =-(-1. + x)*x*(1. - z)*(1. + z); 
00480         outputValues(23,i0, 4) =  2.*(-1. + x)*x*y*z; 
00481         outputValues(23,i0, 5) =-(-1. + x)*x*(1. - y)*(1. + y); 
00482         
00483         outputValues(24,i0, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); 
00484         outputValues(24,i0, 1) = (1. + 2.*x)*y*(-1. + z*z); 
00485         outputValues(24,i0, 2) = (1. + 2.*x)*(-1. + y*y)*z; 
00486         outputValues(24,i0, 3) = x*(1. + x)*(-1. + z)*(1. + z); 
00487         outputValues(24,i0, 4) = 2.*x*(1. + x)*y*z; 
00488         outputValues(24,i0, 5) = x*(1. + x)*(-1. + y)*(1. + y); 
00489         
00490         outputValues(25,i0, 0) = -(-1. + y)*y*(1. - z)*(1. + z); 
00491         outputValues(25,i0, 1) = x*(-1. + 2.*y)*(-1. + z*z); 
00492         outputValues(25,i0, 2) = 2.*x*(-1. + y)*y*z; 
00493         outputValues(25,i0, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z); 
00494         outputValues(25,i0, 4) = (-1. + x*x)*(-1. + 2.*y)*z; 
00495         outputValues(25,i0, 5) =-(1. - x)*(1. + x)*(-1. + y)*y; 
00496 
00497         outputValues(26,i0, 0) =  y*(1. + y)*(-1. + z)*(1. + z);
00498         outputValues(26,i0, 1) =  x*(1. + 2.*y)*(-1. + z*z);
00499         outputValues(26,i0, 2) =  2.*x*y*(1. + y)*z;
00500         outputValues(26,i0, 3) =  (-1. + x)*(1. + x)*(-1. + z)*(1. + z);
00501         outputValues(26,i0, 4) =  (-1. + x*x)*(1. + 2.*y)*z;
00502         outputValues(26,i0, 5) =  (-1. + x)*(1. + x)*y*(1. + y);
00503       }
00504       break;
00505       
00506     case OPERATOR_D3:
00507       for (int i0 = 0; i0 < dim0; i0++) {
00508         x = inputPoints(i0,0);
00509         y = inputPoints(i0,1);
00510         z = inputPoints(i0,2);
00511                 
00512         outputValues(0,i0, 0) = 0.;
00513         outputValues(0,i0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
00514         outputValues(0,i0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
00515         outputValues(0,i0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
00516         outputValues(0,i0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
00517         outputValues(0,i0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
00518         outputValues(0,i0, 6) = 0.;    
00519         outputValues(0,i0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
00520         outputValues(0,i0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
00521         outputValues(0,i0, 9) = 0.;
00522                 
00523         outputValues(1, i0, 0) = 0.;
00524         outputValues(1, i0, 1) = ((-1.+ 2.*y)*(-1.+ z)*z)/4.;
00525         outputValues(1, i0, 2) = ((-1.+ y)*y*(-1.+ 2.*z))/4.;
00526         outputValues(1, i0, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
00527         outputValues(1, i0, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(-1.+ 2.*z))/8.;
00528         outputValues(1, i0, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
00529         outputValues(1, i0, 6) = 0.;
00530         outputValues(1, i0, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
00531         outputValues(1, i0, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
00532         outputValues(1, i0, 9) = 0.;
00533         
00534         outputValues(2, i0, 0) = 0.;
00535         outputValues(2, i0, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
00536         outputValues(2, i0, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
00537         outputValues(2, i0, 3) = ((1.+ 2.*x)*(-1.+ z)*z)/4.;
00538         outputValues(2, i0, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
00539         outputValues(2, i0, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
00540         outputValues(2, i0, 6) = 0.;
00541         outputValues(2, i0, 7) = (x*(1.+ x)*(-1.+ 2.*z))/4.;
00542         outputValues(2, i0, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
00543         outputValues(2, i0, 9) = 0.;
00544         
00545         outputValues(3, i0, 0) = 0.;
00546         outputValues(3, i0, 1) = ((1.+ 2.*y)*(-1.+ z)*z)/4.;
00547         outputValues(3, i0, 2) = (y*(1.+ y)*(-1.+ 2.*z))/4.;
00548         outputValues(3, i0, 3) = ((-1.+ 2.*x)*(-1.+ z)*z)/4.;
00549         outputValues(3, i0, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(-1.+ 2.*z))/8.;
00550         outputValues(3, i0, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
00551         outputValues(3, i0, 6) = 0.;
00552         outputValues(3, i0, 7) = ((-1.+ x)*x*(-1.+ 2.*z))/4.;
00553         outputValues(3, i0, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
00554         outputValues(3, i0, 9) = 0.;
00555         
00556         outputValues(4, i0, 0) = 0.;
00557         outputValues(4, i0, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
00558         outputValues(4, i0, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
00559         outputValues(4, i0, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
00560         outputValues(4, i0, 4) = ((-1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
00561         outputValues(4, i0, 5) = ((-1.+ 2.*x)*(-1.+ y)*y)/4.;
00562         outputValues(4, i0, 6) = 0.;
00563         outputValues(4, i0, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
00564         outputValues(4, i0, 8) = ((-1.+ x)*x*(-1.+ 2.*y))/4.;
00565         outputValues(4, i0, 9) = 0.;
00566         
00567         outputValues(5, i0, 0) = 0.;
00568         outputValues(5, i0, 1) = ((-1.+ 2.*y)*z*(1.+ z))/4.;
00569         outputValues(5, i0, 2) = ((-1.+ y)*y*(1.+ 2.*z))/4.;
00570         outputValues(5, i0, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
00571         outputValues(5, i0, 4) = ((1.+ 2.*x)*(-1.+ 2.*y)*(1.+ 2.*z))/8.;
00572         outputValues(5, i0, 5) = ((1.+ 2.*x)*(-1.+ y)*y)/4.;
00573         outputValues(5, i0, 6) = 0.;
00574         outputValues(5, i0, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
00575         outputValues(5, i0, 8) = (x*(1.+ x)*(-1.+ 2.*y))/4.;
00576         outputValues(5, i0, 9) = 0.;
00577         
00578         outputValues(6, i0, 0) = 0.;
00579         outputValues(6, i0, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
00580         outputValues(6, i0, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
00581         outputValues(6, i0, 3) = ((1.+ 2.*x)*z*(1.+ z))/4.;
00582         outputValues(6, i0, 4) = ((1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
00583         outputValues(6, i0, 5) = ((1.+ 2.*x)*y*(1.+ y))/4.;
00584         outputValues(6, i0, 6) = 0.;
00585         outputValues(6, i0, 7) = (x*(1.+ x)*(1.+ 2.*z))/4.;
00586         outputValues(6, i0, 8) = (x*(1.+ x)*(1.+ 2.*y))/4.;
00587         outputValues(6, i0, 9) = 0.;
00588         
00589         outputValues(7, i0, 0) = 0.;
00590         outputValues(7, i0, 1) = ((1.+ 2.*y)*z*(1.+ z))/4.;
00591         outputValues(7, i0, 2) = (y*(1.+ y)*(1.+ 2.*z))/4.;
00592         outputValues(7, i0, 3) = ((-1.+ 2.*x)*z*(1.+ z))/4.;
00593         outputValues(7, i0, 4) = ((-1.+ 2.*x)*(1.+ 2.*y)*(1.+ 2.*z))/8.;
00594         outputValues(7, i0, 5) = ((-1.+ 2.*x)*y*(1.+ y))/4.;
00595         outputValues(7, i0, 6) = 0.;
00596         outputValues(7, i0, 7) = ((-1.+ x)*x*(1.+ 2.*z))/4.;
00597         outputValues(7, i0, 8) = ((-1.+ x)*x*(1.+ 2.*y))/4.;
00598         outputValues(7, i0, 9) = 0.; 
00599         
00600         outputValues(8, i0, 0) = 0.;
00601         outputValues(8, i0, 1) = -((-1.+ 2.*y)*(-1.+ z)*z)/2.;
00602         outputValues(8, i0, 2) = -((-1.+ y)*y*(-1.+ 2.*z))/2.;
00603         outputValues(8, i0, 3) = -(x*(-1.+ z)*z);
00604         outputValues(8, i0, 4) = -(x*(-1.+ 2.*y)*(-1.+ 2.*z))/2.;
00605         outputValues(8, i0, 5) = -(x*(-1.+ y)*y);
00606         outputValues(8, i0, 6) = 0.;
00607         outputValues(8, i0, 7) = -((-1.+ (x*x))*(-1.+ 2.*z))/2.;    
00608         outputValues(8, i0, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
00609         outputValues(8, i0, 9) = 0.;
00610         
00611         outputValues(9, i0, 0) = 0.;
00612         outputValues(9, i0, 1) = -(y*(-1.+ z)*z);
00613         outputValues(9, i0, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
00614         outputValues(9, i0, 3) = -((1.+ 2.*x)*(-1.+ z)*z)/2.;
00615         outputValues(9, i0, 4) = -((1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
00616         outputValues(9, i0, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
00617         outputValues(9, i0, 6) = 0.;
00618         outputValues(9, i0, 7) = -(x*(1.+ x)*(-1.+ 2.*z))/2.;
00619         outputValues(9, i0, 8) = -(x*(1.+ x)*y);
00620         outputValues(9, i0, 9) = 0.; 
00621         
00622         outputValues(10,i0, 0) = 0.;
00623         outputValues(10,i0, 1) = -((1.+ 2.*y)*(-1.+ z)*z)/2.;
00624         outputValues(10,i0, 2) = -(y*(1.+ y)*(-1.+ 2.*z))/2.;
00625         outputValues(10,i0, 3) = -(x*(-1.+ z)*z);
00626         outputValues(10,i0, 4) = -(x*(1.+ 2.*y)*(-1.+ 2.*z))/2.;
00627         outputValues(10,i0, 5) = -(x*y*(1.+ y));
00628         outputValues(10,i0, 6) = 0.;
00629         outputValues(10,i0, 7) =  -((-1.+ (x*x))*(-1.+ 2.*z))/2.;    
00630         outputValues(10,i0, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
00631         outputValues(10,i0, 9) = 0.;
00632         
00633         outputValues(11,i0, 0) = 0.;
00634         outputValues(11,i0, 1) = -(y*(-1.+ z)*z);
00635         outputValues(11,i0, 2) = -((-1.+ (y*y))*(-1.+ 2.*z))/2.;
00636         outputValues(11,i0, 3) = -((-1.+ 2.*x)*(-1.+ z)*z)/2.;
00637         outputValues(11,i0, 4) = -((-1.+ 2.*x)*y*(-1.+ 2.*z))/2.;
00638         outputValues(11,i0, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
00639         outputValues(11,i0, 6) = 0.;
00640         outputValues(11,i0, 7) = -((-1.+ x)*x*(-1.+ 2.*z))/2.;
00641         outputValues(11,i0, 8) = -((-1.+ x)*x*y);
00642         outputValues(11,i0, 9) = 0.;   
00643         
00644         outputValues(12,i0, 0) = 0.;
00645         outputValues(12,i0, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
00646         outputValues(12,i0, 2) = -((-1.+ y)*y*z);
00647         outputValues(12,i0, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
00648         outputValues(12,i0, 4) = -((-1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
00649         outputValues(12,i0, 5) = -((-1.+ 2.*x)*(-1.+ y)*y)/2.;
00650         outputValues(12,i0, 6) = 0.;
00651         outputValues(12,i0, 7) = -((-1.+ x)*x*z);    
00652         outputValues(12,i0, 8) = -((-1.+ x)*x*(-1.+ 2.*y))/2.;
00653         outputValues(12,i0, 9) = 0.;
00654         
00655         outputValues(13,i0, 0) = 0.;
00656         outputValues(13,i0, 1) = -((-1.+ 2.*y)*(-1.+ (z*z)))/2.;
00657         outputValues(13,i0, 2) = -((-1.+ y)*y*z);
00658         outputValues(13,i0, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
00659         outputValues(13,i0, 4) = -((1.+ 2.*x)*(-1.+ 2.*y)*z)/2.;
00660         outputValues(13,i0, 5) = -((1.+ 2.*x)*(-1.+ y)*y)/2.;
00661         outputValues(13,i0, 6) = 0.;
00662         outputValues(13,i0, 7) = -(x*(1.+ x)*z);
00663         outputValues(13,i0, 8) = -(x*(1.+ x)*(-1.+ 2.*y))/2.;
00664         outputValues(13,i0, 9) = 0.; 
00665         
00666         outputValues(14,i0, 0) = 0.;
00667         outputValues(14,i0, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
00668         outputValues(14,i0, 2) = -(y*(1.+ y)*z);
00669         outputValues(14,i0, 3) = -((1.+ 2.*x)*(-1.+ (z*z)))/2.;
00670         outputValues(14,i0, 4) = -((1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
00671         outputValues(14,i0, 5) = -((1.+ 2.*x)*y*(1.+ y))/2.;
00672         outputValues(14,i0, 6) = 0.;
00673         outputValues(14,i0, 7) = -(x*(1.+ x)*z);    
00674         outputValues(14,i0, 8) = -(x*(1.+ x)*(1.+ 2.*y))/2.;
00675         outputValues(14,i0, 9) = 0.;
00676         
00677         outputValues(15,i0, 0) = 0.;
00678         outputValues(15,i0, 1) = -((1.+ 2.*y)*(-1.+ (z*z)))/2.;
00679         outputValues(15,i0, 2) = -(y*(1.+ y)*z);
00680         outputValues(15,i0, 3) = -((-1.+ 2.*x)*(-1.+ (z*z)))/2.;
00681         outputValues(15,i0, 4) = -((-1.+ 2.*x)*(1.+ 2.*y)*z)/2.;
00682         outputValues(15,i0, 5) = -((-1.+ 2.*x)*y*(1.+ y))/2.;
00683         outputValues(15,i0, 6) = 0.;
00684         outputValues(15,i0, 7) = -((-1.+ x)*x*z);
00685         outputValues(15,i0, 8) = -((-1.+ x)*x*(1.+ 2.*y))/2.;
00686         outputValues(15,i0, 9) = 0.; 
00687         
00688         outputValues(16,i0, 0) = 0.;
00689         outputValues(16,i0, 1) = -((-1.+ 2.*y)*z*(1.+ z))/2.;
00690         outputValues(16,i0, 2) = -((-1.+ y)*y*(1.+ 2.*z))/2.;
00691         outputValues(16,i0, 3) = -(x*z*(1.+ z));
00692         outputValues(16,i0, 4) = -(x*(-1.+ 2.*y)*(1.+ 2.*z))/2.;
00693         outputValues(16,i0, 5) = -(x*(-1.+ y)*y);
00694         outputValues(16,i0, 6) = 0.;
00695         outputValues(16,i0, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;    
00696         outputValues(16,i0, 8) = -((-1.+ (x*x))*(-1.+ 2.*y))/2.;
00697         outputValues(16,i0, 9) = 0.;
00698         
00699         outputValues(17,i0, 0) = 0.;
00700         outputValues(17,i0, 1) = -(y*z*(1.+ z));
00701         outputValues(17,i0, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
00702         outputValues(17,i0, 3) = -((1.+ 2.*x)*z*(1.+ z))/2.;
00703         outputValues(17,i0, 4) = -((1.+ 2.*x)*y*(1.+ 2.*z))/2.;
00704         outputValues(17,i0, 5) = -((1.+ 2.*x)*(-1.+ (y*y)))/2.;
00705         outputValues(17,i0, 6) = 0.;
00706         outputValues(17,i0, 7) = -(x*(1.+ x)*(1.+ 2.*z))/2.;
00707         outputValues(17,i0, 8) = -(x*(1.+ x)*y);
00708         outputValues(17,i0, 9) = 0.;
00709         
00710         outputValues(18,i0, 0) = 0.;
00711         outputValues(18,i0, 1) = -((1.+ 2.*y)*z*(1.+ z))/2.;
00712         outputValues(18,i0, 2) = -(y*(1.+ y)*(1.+ 2.*z))/2.;
00713         outputValues(18,i0, 3) = -(x*z*(1.+ z));
00714         outputValues(18,i0, 4) = -(x*(1.+ 2.*y)*(1.+ 2.*z))/2.;
00715         outputValues(18,i0, 5) = -(x*y*(1.+ y));
00716         outputValues(18,i0, 6) = 0.;
00717         outputValues(18,i0, 7) = -((-1.+ (x*x))*(1.+ 2.*z))/2.;    
00718         outputValues(18,i0, 8) = -((-1.+ (x*x))*(1.+ 2.*y))/2.;
00719         outputValues(18,i0, 9) = 0.;
00720         
00721         outputValues(19,i0, 0) = 0.;
00722         outputValues(19,i0, 1) = -(y*z*(1.+ z));
00723         outputValues(19,i0, 2) = -((-1.+ (y*y))*(1.+ 2.*z))/2.;
00724         outputValues(19,i0, 3) = -((-1.+ 2.*x)*z*(1.+ z))/2.;
00725         outputValues(19,i0, 4) = -((-1.+ 2.*x)*y*(1.+ 2.*z))/2.;
00726         outputValues(19,i0, 5) = -((-1.+ 2.*x)*(-1.+ (y*y)))/2.;
00727         outputValues(19,i0, 6) = 0.;
00728         outputValues(19,i0, 7) = -((-1.+ x)*x*(1.+ 2.*z))/2.;
00729         outputValues(19,i0, 8) = -((-1.+ x)*x*y);
00730         outputValues(19,i0, 9) = 0.;
00731         
00732         outputValues(20,i0, 0) = 0.;
00733         outputValues(20,i0, 1) = -4*y*(-1.+ (z*z));
00734         outputValues(20,i0, 2) = -4*(-1.+ (y*y))*z;
00735         outputValues(20,i0, 3) = -4*x*(-1.+ (z*z));
00736         outputValues(20,i0, 4) = -8*x*y*z;
00737         outputValues(20,i0, 5) = -4*x*(-1.+ (y*y));
00738         outputValues(20,i0, 6) = 0.;
00739         outputValues(20,i0, 7) = -4*(-1.+ (x*x))*z;
00740         outputValues(20,i0, 8) = -4*(-1.+ (x*x))*y;
00741         outputValues(20,i0, 9) = 0.;    
00742         
00743         outputValues(21,i0, 0) = 0.;
00744         outputValues(21,i0, 1) = 2.*y*(-1.+ z)*z;
00745         outputValues(21,i0, 2) = (-1.+ (y*y))*(-1.+ 2.*z);
00746         outputValues(21,i0, 3) = 2.*x*(-1.+ z)*z;
00747         outputValues(21,i0, 4) = 2.*x*y*(-1.+ 2.*z);
00748         outputValues(21,i0, 5) = 2.*x*(-1.+ (y*y));
00749         outputValues(21,i0, 6) = 0.;
00750         outputValues(21,i0, 7) = (-1.+ (x*x))*(-1.+ 2.*z);
00751         outputValues(21,i0, 8) = 2.*(-1.+ (x*x))*y;
00752         outputValues(21,i0, 9) = 0.;    
00753         
00754         outputValues(22,i0, 0) = 0.;
00755         outputValues(22,i0, 1) = 2.*y*z*(1.+ z);
00756         outputValues(22,i0, 2) = (-1.+ (y*y))*(1.+ 2.*z);
00757         outputValues(22,i0, 3) = 2.*x*z*(1.+ z);
00758         outputValues(22,i0, 4) = 2.*x*y*(1.+ 2.*z);
00759         outputValues(22,i0, 5) = 2.*x*(-1.+ (y*y));
00760         outputValues(22,i0, 6) = 0.;
00761         outputValues(22,i0, 7) = (-1.+ (x*x))*(1.+ 2.*z);
00762         outputValues(22,i0, 8) = 2.*(-1.+ (x*x))*y;
00763         outputValues(22,i0, 9) = 0.;    
00764         
00765         outputValues(23,i0, 0) = 0.;
00766         outputValues(23,i0, 1) = 2.*y*(-1.+ (z*z));
00767         outputValues(23,i0, 2) = 2.*(-1.+ (y*y))*z;
00768         outputValues(23,i0, 3) = (-1.+ 2.*x)*(-1.+ (z*z));
00769         outputValues(23,i0, 4) = 2.*(-1.+ 2.*x)*y*z;
00770         outputValues(23,i0, 5) = (-1.+ 2.*x)*(-1.+ (y*y));
00771         outputValues(23,i0, 6) = 0.;
00772         outputValues(23,i0, 7) = 2.*(-1.+ x)*x*z;
00773         outputValues(23,i0, 8) = 2.*(-1.+ x)*x*y;
00774         outputValues(23,i0, 9) = 0.;
00775         
00776         outputValues(24,i0, 0) = 0.;
00777         outputValues(24,i0, 1) = 2.*y*(-1.+ (z*z));
00778         outputValues(24,i0, 2) = 2.*(-1.+ (y*y))*z;
00779         outputValues(24,i0, 3) = (1.+ 2.*x)*(-1.+ (z*z));
00780         outputValues(24,i0, 4) = 2.*(1.+ 2.*x)*y*z;
00781         outputValues(24,i0, 5) = (1.+ 2.*x)*(-1.+ (y*y));
00782         outputValues(24,i0, 6) = 0.;
00783         outputValues(24,i0, 7) = 2.*x*(1.+ x)*z;
00784         outputValues(24,i0, 8) = 2.*x*(1.+ x)*y;
00785         outputValues(24,i0, 9) = 0.;   
00786         
00787         outputValues(25,i0, 0) = 0.;
00788         outputValues(25,i0, 1) = (-1.+ 2.*y)*(-1.+ (z*z));
00789         outputValues(25,i0, 2) = 2.*(-1.+ y)*y*z;
00790         outputValues(25,i0, 3) = 2.*x*(-1.+ (z*z));
00791         outputValues(25,i0, 4) = 2.*x*(-1.+ 2.*y)*z;
00792         outputValues(25,i0, 5) = 2.*x*(-1.+ y)*y;
00793         outputValues(25,i0, 6) = 0.;
00794         outputValues(25,i0, 7) = 2.*(-1.+ (x*x))*z;
00795         outputValues(25,i0, 8) = (-1.+ (x*x))*(-1.+ 2.*y);
00796         outputValues(25,i0, 9) = 0.;  
00797         
00798         outputValues(26,i0, 0) = 0.;
00799         outputValues(26,i0, 1) = (1.+ 2.*y)*(-1.+ (z*z));    
00800         outputValues(26,i0, 2) = 2.*y*(1.+ y)*z;
00801         outputValues(26,i0, 3) = 2.*x*(-1.+ (z*z));
00802         outputValues(26,i0, 4) = 2.*x*(1.+ 2.*y)*z;
00803         outputValues(26,i0, 5) = 2.*x*y*(1.+ y);
00804         outputValues(26,i0, 6) = 0.;
00805         outputValues(26,i0, 7) = 2.*(-1.+ (x*x))*z;
00806         outputValues(26,i0, 8) = (-1.+ (x*x))*(1.+ 2.*y);
00807         outputValues(26,i0, 9) = 0.;
00808       }
00809       break;
00810       
00811     case OPERATOR_D4:
00812       {
00813         // Non-zero entries have Dk (derivative cardinality) indices {3,4,5,7,8,12}, all other entries are 0.
00814         // Intitialize array by zero and then fill only non-zero entries.
00815         int DkCardinality = Intrepid::getDkCardinality(operatorType, this -> basisCellTopology_.getDimension() );
00816         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
00817           for (int i0 = 0; i0 < dim0; i0++) {
00818             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
00819               outputValues(dofOrd, i0, dkOrd) = 0.0;
00820             }
00821           }
00822         }
00823         
00824         for (int i0 = 0; i0 < dim0; i0++) {
00825           x = inputPoints(i0,0);
00826           y = inputPoints(i0,1);
00827           z = inputPoints(i0,2);
00828           
00829           outputValues(0, i0, 3) = ((-1.+ z)*z)/2.;
00830           outputValues(0, i0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.; 
00831           outputValues(0, i0, 5) = ((-1.+ y)*y)/2.;
00832           outputValues(0, i0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
00833           outputValues(0, i0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
00834           outputValues(0, i0, 12)= ((-1.+ x)*x)/2.;
00835           
00836           outputValues(1, i0, 3) = ((-1.+ z)*z)/2.;
00837           outputValues(1, i0, 4) = ((-1.+ 2.*y)*(-1.+ 2.*z))/4.;
00838           outputValues(1, i0, 5) = ((-1.+ y)*y)/2.;
00839           outputValues(1, i0, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
00840           outputValues(1, i0, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
00841           outputValues(1, i0, 12)= (x*(1. + x))/2.;
00842           
00843           outputValues(2, i0, 3) = ((-1.+ z)*z)/2.;
00844           outputValues(2, i0, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
00845           outputValues(2, i0, 5) = (y*(1. + y))/2.;
00846           outputValues(2, i0, 7) = ((1. + 2.*x)*(-1.+ 2.*z))/4.;
00847           outputValues(2, i0, 8) =  ((1. + 2.*x)*(1. + 2.*y))/4.;
00848           outputValues(2, i0, 12)= (x*(1. + x))/2.;
00849           
00850           outputValues(3, i0, 3) = ((-1.+ z)*z)/2.;
00851           outputValues(3, i0, 4) = ((1. + 2.*y)*(-1.+ 2.*z))/4.;
00852           outputValues(3, i0, 5) = (y*(1. + y))/2.;
00853           outputValues(3, i0, 7) = ((-1.+ 2.*x)*(-1.+ 2.*z))/4.;
00854           outputValues(3, i0, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.;
00855           outputValues(3, i0, 12)= ((-1.+ x)*x)/2.;
00856           
00857           outputValues(4, i0, 3) = (z*(1. + z))/2.;
00858           outputValues(4, i0, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
00859           outputValues(4, i0, 5) = ((-1.+ y)*y)/2.;
00860           outputValues(4, i0, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
00861           outputValues(4, i0, 8) = ((-1.+ 2.*x)*(-1.+ 2.*y))/4.;
00862           outputValues(4, i0, 12)= ((-1.+ x)*x)/2.;
00863           
00864           outputValues(5, i0, 3) = (z*(1. + z))/2.;
00865           outputValues(5, i0, 4) = ((-1.+ 2.*y)*(1. + 2.*z))/4.;
00866           outputValues(5, i0, 5) = ((-1.+ y)*y)/2.;
00867           outputValues(5, i0, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.; 
00868           outputValues(5, i0, 8) = ((1. + 2.*x)*(-1.+ 2.*y))/4.;
00869           outputValues(5, i0, 12)= (x*(1. + x))/2.;
00870           
00871           outputValues(6, i0, 3) = (z*(1. + z))/2.;
00872           outputValues(6, i0, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
00873           outputValues(6, i0, 5) = (y*(1. + y))/2.;
00874           outputValues(6, i0, 7) = ((1. + 2.*x)*(1. + 2.*z))/4.;
00875           outputValues(6, i0, 8) = ((1. + 2.*x)*(1. + 2.*y))/4.;
00876           outputValues(6, i0, 12)=  (x*(1. + x))/2.; 
00877           
00878           outputValues(7, i0, 3) = (z*(1. + z))/2.;
00879           outputValues(7, i0, 4) = ((1. + 2.*y)*(1. + 2.*z))/4.;
00880           outputValues(7, i0, 5) = (y*(1. + y))/2.;
00881           outputValues(7, i0, 7) = ((-1.+ 2.*x)*(1. + 2.*z))/4.;
00882           outputValues(7, i0, 8) = ((-1.+ 2.*x)*(1. + 2.*y))/4.; 
00883           outputValues(7, i0, 12)= ((-1.+ x)*x)/2.;
00884           
00885           outputValues(8, i0, 3) = -((-1.+ z)*z);
00886           outputValues(8, i0, 4) = -0.5 + y + z - 2.*y*z;
00887           outputValues(8, i0, 5) = -((-1.+ y)*y);
00888           outputValues(8, i0, 7) = x - 2.*x*z;  
00889           outputValues(8, i0, 8) = x - 2.*x*y;
00890           outputValues(8, i0, 12)= 1. - x*x;  
00891           
00892           outputValues(9, i0, 3) = -((-1.+ z)*z);
00893           outputValues(9, i0, 4) = y - 2.*y*z;   
00894           outputValues(9, i0, 5) = 1 - y*y;
00895           outputValues(9, i0, 7) = 0.5 + x - z - 2.*x*z;
00896           outputValues(9, i0, 8) = -((1. + 2.*x)*y);
00897           outputValues(9, i0, 12)= -(x*(1. + x)); 
00898           
00899           outputValues(10,i0, 3) = -((-1.+ z)*z);
00900           outputValues(10,i0, 4) = 0.5 + y - z - 2.*y*z; 
00901           outputValues(10,i0, 5) = -(y*(1. + y));
00902           outputValues(10,i0, 7) = x - 2.*x*z; 
00903           outputValues(10,i0, 8) = -(x*(1. + 2.*y)); 
00904           outputValues(10,i0, 12)=  1. - x*x;   
00905           
00906           outputValues(11,i0, 3) = -((-1.+ z)*z);
00907           outputValues(11,i0, 4) =  y - 2.*y*z;   
00908           outputValues(11,i0, 5) =  1. - y*y;   
00909           outputValues(11,i0, 7) = -0.5 + x + z - 2.*x*z;
00910           outputValues(11,i0, 8) =  y - 2.*x*y;
00911           outputValues(11,i0, 12)= -((-1.+ x)*x);
00912           
00913           outputValues(12,i0, 3) = 1. - z*z;
00914           outputValues(12,i0, 4) = z - 2.*y*z;
00915           outputValues(12,i0, 5) = -((-1.+ y)*y);
00916           outputValues(12,i0, 7) =  z - 2.*x*z;
00917           outputValues(12,i0, 8) = -0.5 + x + y - 2.*x*y;
00918           outputValues(12,i0, 12)= -((-1.+ x)*x); 
00919           
00920           outputValues(13,i0, 3) =  1. - z*z;  
00921           outputValues(13,i0, 4) = z - 2.*y*z;
00922           outputValues(13,i0, 5) = -((-1.+ y)*y);
00923           outputValues(13,i0, 7) =  -((1. + 2.*x)*z); 
00924           outputValues(13,i0, 8) = 0.5 + x - y - 2.*x*y; 
00925           outputValues(13,i0, 12)= -(x*(1. + x));
00926           
00927           outputValues(14,i0, 3) = 1. - z*z;
00928           outputValues(14,i0, 4) = -((1. + 2.*y)*z);
00929           outputValues(14,i0, 5) = -(y*(1. + y));
00930           outputValues(14,i0, 7) = -((1. + 2.*x)*z);
00931           outputValues(14,i0, 8) = -((1. + 2.*x)*(1. + 2.*y))/2.;
00932           outputValues(14,i0, 12)= -(x*(1. + x)); 
00933           
00934           outputValues(15,i0, 3) =  1. - z*z;
00935           outputValues(15,i0, 4) = -((1. + 2.*y)*z);
00936           outputValues(15,i0, 5) = -(y*(1. + y));
00937           outputValues(15,i0, 7) = z - 2.*x*z; 
00938           outputValues(15,i0, 8) = 0.5 + y - x*(1. + 2.*y); 
00939           outputValues(15,i0, 12)= -((-1.+ x)*x);
00940           
00941           outputValues(16,i0, 3) = -(z*(1. + z)); 
00942           outputValues(16,i0, 4) = 0.5 + z - y*(1. + 2.*z);
00943           outputValues(16,i0, 5) = -((-1.+ y)*y);
00944           outputValues(16,i0, 7) = -(x*(1. + 2.*z));
00945           outputValues(16,i0, 8) = x - 2.*x*y;
00946           outputValues(16,i0, 12)= 1. - x*x;  
00947           
00948           outputValues(17,i0, 3) = -(z*(1. + z));
00949           outputValues(17,i0, 4) = -(y*(1. + 2.*z));
00950           outputValues(17,i0, 5) = 1. - y*y;
00951           outputValues(17,i0, 7) = -((1. + 2.*x)*(1. + 2.*z))/2.;
00952           outputValues(17,i0, 8) = -((1. + 2.*x)*y); 
00953           outputValues(17,i0, 12)= -(x*(1. + x));
00954           
00955           outputValues(18,i0, 3) = -(z*(1. + z));
00956           outputValues(18,i0, 4) = -((1. + 2.*y)*(1. + 2.*z))/2.;
00957           outputValues(18,i0, 5) = -(y*(1. + y));
00958           outputValues(18,i0, 7) =  -(x*(1. + 2.*z)); 
00959           outputValues(18,i0, 8) =  -(x*(1. + 2.*y)); 
00960           outputValues(18,i0, 12)= 1. - x*x; 
00961           
00962           outputValues(19,i0, 3) = -(z*(1. + z));
00963           outputValues(19,i0, 4) = -(y*(1. + 2.*z));
00964           outputValues(19,i0, 5) = 1. - y*y; 
00965           outputValues(19,i0, 7) = 0.5 + z - x*(1. + 2.*z);
00966           outputValues(19,i0, 8) = y - 2.*x*y;
00967           outputValues(19,i0, 12)= -((-1.+ x)*x); 
00968           
00969           outputValues(20,i0, 3) = 4. - 4.*z*z;
00970           outputValues(20,i0, 4) = -8.*y*z;
00971           outputValues(20,i0, 5) = 4. - 4.*y*y; 
00972           outputValues(20,i0, 7) = -8.*x*z;
00973           outputValues(20,i0, 8) = -8.*x*y; 
00974           outputValues(20,i0, 12)= 4. - 4.*x*x; 
00975           
00976           outputValues(21,i0, 3) = 2.*(-1.+ z)*z;
00977           outputValues(21,i0, 4) = 2.*y*(-1.+ 2.*z);  
00978           outputValues(21,i0, 5) = 2.*(-1.+ y*y);
00979           outputValues(21,i0, 7) = 2.*x*(-1.+ 2.*z);  
00980           outputValues(21,i0, 8) = 4.*x*y;  
00981           outputValues(21,i0, 12)= 2.*(-1.+ x*x); 
00982           
00983           outputValues(22,i0, 3) = 2.*z*(1. + z);
00984           outputValues(22,i0, 4) = 2.*(y + 2.*y*z);
00985           outputValues(22,i0, 5) = 2.*(-1.+ y*y);
00986           outputValues(22,i0, 7) = 2.*(x + 2.*x*z); 
00987           outputValues(22,i0, 8) = 4.*x*y;
00988           outputValues(22,i0, 12)= 2.*(-1.+ x*x);
00989           
00990           outputValues(23,i0, 3) = 2.*(-1.+ z*z); 
00991           outputValues(23,i0, 4) = 4.*y*z; 
00992           outputValues(23,i0, 5) = 2.*(-1.+ y*y);
00993           outputValues(23,i0, 7) = 2.*(-1.+ 2.*x)*z;
00994           outputValues(23,i0, 8) = 2.*(-1.+ 2.*x)*y; 
00995           outputValues(23,i0, 12)= 2.*(-1.+ x)*x; 
00996           
00997           outputValues(24,i0, 3) = 2.*(-1.+ z*z);
00998           outputValues(24,i0, 4) = 4.*y*z; 
00999           outputValues(24,i0, 5) = 2.*(-1.+ y*y);
01000           outputValues(24,i0, 7) = 2.*(z + 2.*x*z); 
01001           outputValues(24,i0, 8) = 2.*(y + 2.*x*y);
01002           outputValues(24,i0, 12)= 2.*x*(1. + x);
01003           
01004           outputValues(25,i0, 3) =  2.*(-1.+ z*z); 
01005           outputValues(25,i0, 4) = 2.*(-1.+ 2.*y)*z;
01006           outputValues(25,i0, 5) = 2.*(-1.+ y)*y;
01007           outputValues(25,i0, 7) = 4.*x*z;  
01008           outputValues(25,i0, 8) = 2.*x*(-1.+ 2.*y);  
01009           outputValues(25,i0, 12)= 2.*(-1.+ x*x);
01010           
01011           outputValues(26,i0, 3) = 2.*(-1.+ z*z);
01012           outputValues(26,i0, 4) = 2.*(z + 2.*y*z);
01013           outputValues(26,i0, 5) = 2.*y*(1. + y);
01014           outputValues(26,i0, 7) =  4.*x*z;
01015           outputValues(26,i0, 8) = 2.*(x + 2.*x*y); 
01016           outputValues(26,i0, 12)= 2.*(-1.+ x*x);   
01017         }
01018       }
01019       break;
01020       
01021     case OPERATOR_D5:
01022     case OPERATOR_D6:
01023       TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
01024                           ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): operator not supported");
01025       break;
01026       
01027     case OPERATOR_D7:
01028     case OPERATOR_D8:
01029     case OPERATOR_D9:
01030     case OPERATOR_D10:
01031       {
01032         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
01033         int DkCardinality = Intrepid::getDkCardinality(operatorType, 
01034                                                        this -> basisCellTopology_.getDimension() );
01035         for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
01036           for (int i0 = 0; i0 < dim0; i0++) {
01037             for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
01038               outputValues(dofOrd, i0, dkOrd) = 0.0;
01039             }
01040           }
01041         }
01042       }
01043       break;
01044       
01045     default:
01046       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
01047                           ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): Invalid operator type");
01048   }
01049 }
01050 
01051 
01052 
01053 template<class Scalar, class ArrayScalar>
01054 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
01055                                                             const ArrayScalar &    inputPoints,
01056                                                             const ArrayScalar &    cellVertices,
01057                                                             const EOperator        operatorType) const {
01058   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
01059                       ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): FEM Basis calling an FVD member function");
01060                                                             }
01061 
01062 template<class Scalar, class ArrayScalar>
01063 void Basis_HGRAD_HEX_C2_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const {
01064 #ifdef HAVE_INTREPID_DEBUG
01065   // Verify rank of output array.
01066   TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
01067                       ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
01068   // Verify 0th dimension of output array.
01069   TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
01070                       ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
01071   // Verify 1st dimension of output array.
01072   TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
01073                       ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
01074 #endif
01075 
01076   DofCoords(0,0) = -1.0;   DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0;  
01077   DofCoords(1,0) =  1.0;   DofCoords(1,1) = -1.0; DofCoords(1,2) = -1.0;  
01078   DofCoords(2,0) =  1.0;   DofCoords(2,1) =  1.0; DofCoords(2,2) = -1.0;  
01079   DofCoords(3,0) = -1.0;   DofCoords(3,1) =  1.0; DofCoords(3,2) = -1.0;  
01080   DofCoords(4,0) = -1.0;   DofCoords(4,1) = -1.0; DofCoords(4,2) =  1.0;  
01081   DofCoords(5,0) =  1.0;   DofCoords(5,1) = -1.0; DofCoords(5,2) =  1.0;  
01082   DofCoords(6,0) =  1.0;   DofCoords(6,1) =  1.0; DofCoords(6,2) =  1.0;  
01083   DofCoords(7,0) = -1.0;   DofCoords(7,1) =  1.0; DofCoords(7,2) =  1.0;  
01084 
01085   DofCoords(8,0) =   0.0;   DofCoords(8,1) =  -1.0; DofCoords(8,2) =  -1.0;  
01086   DofCoords(9,0) =   1.0;   DofCoords(9,1) =   0.0; DofCoords(9,2) =  -1.0;  
01087   DofCoords(10,0) =  0.0;   DofCoords(10,1) =  1.0; DofCoords(10,2) = -1.0;  
01088   DofCoords(11,0) = -1.0;   DofCoords(11,1) =  0.0; DofCoords(11,2) = -1.0;  
01089   DofCoords(12,0) = -1.0;   DofCoords(12,1) = -1.0; DofCoords(12,2) =  0.0;  
01090   DofCoords(13,0) =  1.0;   DofCoords(13,1) = -1.0; DofCoords(13,2) =  0.0;  
01091   DofCoords(14,0) =  1.0;   DofCoords(14,1) =  1.0; DofCoords(14,2) =  0.0;  
01092   DofCoords(15,0) = -1.0;   DofCoords(15,1) =  1.0; DofCoords(15,2) =  0.0;  
01093   DofCoords(16,0) =  0.0;   DofCoords(16,1) = -1.0; DofCoords(16,2) =  1.0;  
01094   DofCoords(17,0) =  1.0;   DofCoords(17,1) =  0.0; DofCoords(17,2) =  1.0;  
01095   DofCoords(18,0) =  0.0;   DofCoords(18,1) =  1.0; DofCoords(18,2) =  1.0;  
01096   DofCoords(19,0) = -1.0;   DofCoords(19,1) =  0.0; DofCoords(19,2) =  1.0;  
01097 
01098   DofCoords(20,0) =  0.0;   DofCoords(20,1) =  0.0; DofCoords(20,2) =  0.0;  
01099 
01100   DofCoords(21,0) =  0.0;   DofCoords(21,1) =  0.0; DofCoords(21,2) = -1.0;  
01101   DofCoords(22,0) =  0.0;   DofCoords(22,1) =  0.0; DofCoords(22,2) =  1.0;  
01102   DofCoords(23,0) = -1.0;   DofCoords(23,1) =  0.0; DofCoords(23,2) =  0.0;  
01103   DofCoords(24,0) =  1.0;   DofCoords(24,1) =  0.0; DofCoords(24,2) =  0.0;  
01104   DofCoords(25,0) =  0.0;   DofCoords(25,1) = -1.0; DofCoords(25,2) =  0.0;  
01105   DofCoords(26,0) =  0.0;   DofCoords(26,1) =  1.0; DofCoords(26,2) =  0.0;  
01106 }
01107 
01108 }// namespace Intrepid
01109 #endif