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