Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_HEX_I1

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 //                    Kara Peterson (kjpeter@sandia.gov).
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00036 namespace Intrepid {
00037 
00038 template<class Scalar, class ArrayScalar>
00039 Basis_HCURL_HEX_I1_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_I1_FEM()
00040   {
00041     this -> basisCardinality_  = 12;
00042     this -> basisDegree_       = 1;
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 template<class Scalar, class ArrayScalar>
00050 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::initializeTags() {
00051   
00052   // Basis-dependent intializations
00053   int tagSize  = 4;        // size of DoF tag
00054   int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00055   int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00056   int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00057 
00058   // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 
00059   int tags[]  = {
00060                   1, 0, 0, 1,
00061                   1, 1, 0, 1,
00062                   1, 2, 0, 1,
00063                   1, 3, 0, 1,
00064                   1, 4, 0, 1,
00065                   1, 5, 0, 1,
00066                   1, 6, 0, 1,
00067                   1, 7, 0, 1,
00068                   1, 8, 0, 1,
00069                   1, 9, 0, 1,
00070                   1, 10, 0, 1,
00071                   1, 11, 0, 1 };
00072   
00073   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00074   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00075                               this -> ordinalToTag_,
00076                               tags,
00077                               this -> basisCardinality_,
00078                               tagSize,
00079                               posScDim,
00080                               posScOrd,
00081                               posDfOrd);
00082 }
00083 
00084 
00085 
00086 template<class Scalar, class ArrayScalar>
00087 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00088                                                             const ArrayScalar &  inputPoints,
00089                                                             const EOperator      operatorType) const {
00090   
00091 // Verify arguments
00092 #ifdef HAVE_INTREPID_DEBUG
00093   Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
00094                                                       inputPoints,
00095                                                       operatorType,
00096                                                       this -> getBaseCellTopology(),
00097                                                       this -> getCardinality() );
00098 #endif
00099   
00100  // Number of evaluation points = dim 0 of inputPoints
00101   int dim0 = inputPoints.dimension(0);
00102 
00103   // Temporaries: (x,y,z) coordinates of the evaluation point
00104   Scalar x = 0.0;                                    
00105   Scalar y = 0.0;                                    
00106   Scalar z = 0.0;                                    
00107   
00108   switch (operatorType) {
00109     case OPERATOR_VALUE:
00110       for (int i0 = 0; i0 < dim0; i0++) {
00111         x = inputPoints(i0, 0);
00112         y = inputPoints(i0, 1);
00113         z = inputPoints(i0, 2);
00114         
00115         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00116         outputValues(0, i0, 0) = (1.0 - y)*(1.0 - z)/8.0;
00117         outputValues(0, i0, 1) = 0.0;
00118         outputValues(0, i0, 2) = 0.0;
00119 
00120         outputValues(1, i0, 0) = 0.0;
00121         outputValues(1, i0, 1) = (1.0 + x)*(1.0 - z)/8.0;
00122         outputValues(1, i0, 2) = 0.0;
00123 
00124         outputValues(2, i0, 0) = -(1.0 + y)*(1.0 - z)/8.0;
00125         outputValues(2, i0, 1) = 0.0;
00126         outputValues(2, i0, 2) = 0.0;
00127 
00128         outputValues(3, i0, 0) = 0.0;
00129         outputValues(3, i0, 1) = -(1.0 - x)*(1.0 - z)/8.0;
00130         outputValues(3, i0, 2) = 0.0;
00131 
00132         outputValues(4, i0, 0) = (1.0 - y)*(1.0 + z)/8.0;
00133         outputValues(4, i0, 1) = 0.0;
00134         outputValues(4, i0, 2) = 0.0;
00135 
00136         outputValues(5, i0, 0) = 0.0;
00137         outputValues(5, i0, 1) = (1.0 + x)*(1.0 + z)/8.0;
00138         outputValues(5, i0, 2) = 0.0;
00139 
00140         outputValues(6, i0, 0) = -(1.0 + y)*(1.0 + z)/8.0;
00141         outputValues(6, i0, 1) = 0.0;
00142         outputValues(6, i0, 2) = 0.0;
00143 
00144         outputValues(7, i0, 0) = 0.0;
00145         outputValues(7, i0, 1) = -(1.0 - x)*(1.0 + z)/8.0;
00146         outputValues(7, i0, 2) = 0.0;
00147 
00148         outputValues(8, i0, 0) = 0.0;
00149         outputValues(8, i0, 1) = 0.0;
00150         outputValues(8, i0, 2) = (1.0 - x)*(1.0 - y)/8.0;
00151 
00152         outputValues(9, i0, 0) = 0.0;
00153         outputValues(9, i0, 1) = 0.0;
00154         outputValues(9, i0, 2) = (1.0 + x)*(1.0 - y)/8.0;
00155 
00156         outputValues(10, i0, 0) = 0.0;
00157         outputValues(10, i0, 1) = 0.0;
00158         outputValues(10, i0, 2) = (1.0 + x)*(1.0 + y)/8.0;
00159 
00160         outputValues(11, i0, 0) = 0.0;
00161         outputValues(11, i0, 1) = 0.0;
00162         outputValues(11, i0, 2) = (1.0 - x)*(1.0 + y)/8.0;
00163       }
00164       break;
00165       
00166     case OPERATOR_CURL:
00167       for (int i0 = 0; i0 < dim0; i0++) {
00168         x = inputPoints(i0, 0);
00169         y = inputPoints(i0, 1);
00170         z = inputPoints(i0, 2);
00171         
00172         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00173         outputValues(0, i0, 0) = 0.0;
00174         outputValues(0, i0, 1) = -(1.0 - y)/8.0;
00175         outputValues(0, i0, 2) = (1.0 - z)/8.0;
00176 
00177         outputValues(1, i0, 0) = (1.0 + x)/8.0;
00178         outputValues(1, i0, 1) = 0.0;
00179         outputValues(1, i0, 2) = (1.0 - z)/8.0;
00180 
00181         outputValues(2, i0, 0) = 0.0;
00182         outputValues(2, i0, 1) = (1.0 + y)/8.0;
00183         outputValues(2, i0, 2) = (1.0 - z)/8.0;
00184 
00185         outputValues(3, i0, 0) = -(1.0 - x)/8.0;
00186         outputValues(3, i0, 1) = 0.0;
00187         outputValues(3, i0, 2) = (1.0 - z)/8.0;
00188 
00189         outputValues(4, i0, 0) = 0.0;
00190         outputValues(4, i0, 1) = (1.0 - y)/8.0;
00191         outputValues(4, i0, 2) = (1.0 + z)/8.0;
00192 
00193         outputValues(5, i0, 0) = -(1.0 + x)/8.0;
00194         outputValues(5, i0, 1) = 0.0;
00195         outputValues(5, i0, 2) = (1.0 + z)/8.0;
00196 
00197         outputValues(6, i0, 0) = 0.0;
00198         outputValues(6, i0, 1) = -(1.0 + y)/8.0;
00199         outputValues(6, i0, 2) = (1.0 + z)/8.0;
00200 
00201         outputValues(7, i0, 0) = (1.0 - x)/8.0;
00202         outputValues(7, i0, 1) = 0.0;
00203         outputValues(7, i0, 2) = (1.0 + z)/8.0;
00204 
00205         outputValues(8, i0, 0) = -(1.0 - x)/8.0;
00206         outputValues(8, i0, 1) = (1.0 - y)/8.0;
00207         outputValues(8, i0, 2) = 0.0;
00208 
00209         outputValues(9, i0, 0) = -(1.0 + x)/8.0;
00210         outputValues(9, i0, 1) = -(1.0 - y)/8.0;
00211         outputValues(9, i0, 2) = 0.0;
00212 
00213         outputValues(10, i0, 0) = (1.0 + x)/8.0;
00214         outputValues(10, i0, 1) = -(1.0 + y)/8.0;
00215         outputValues(10, i0, 2) = 0.0;
00216 
00217         outputValues(11, i0, 0) = (1.0 - x)/8.0;
00218         outputValues(11, i0, 1) = (1.0 + y)/8.0;
00219         outputValues(11, i0, 2) = 0.0;
00220       }
00221       break;
00222       
00223     case OPERATOR_DIV:
00224        TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00225                           ">>> ERROR (Basis_HCURL_HEX_I1_FEM): DIV is invalid operator for HCURL Basis Functions");
00226       break;
00227       
00228     case OPERATOR_GRAD:
00229        TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
00230                           ">>> ERROR (Basis_HCURL_HEX_I1_FEM): GRAD is invalid operator for HCURL Basis Functions");
00231       break;
00232 
00233     case OPERATOR_D1:
00234     case OPERATOR_D2:
00235     case OPERATOR_D3:
00236     case OPERATOR_D4:
00237     case OPERATOR_D5:
00238     case OPERATOR_D6:
00239     case OPERATOR_D7:
00240     case OPERATOR_D8:
00241     case OPERATOR_D9:
00242     case OPERATOR_D10:
00243       TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1)    ||
00244                             (operatorType == OPERATOR_D2)    ||
00245                             (operatorType == OPERATOR_D3)    ||
00246                             (operatorType == OPERATOR_D4)    ||
00247                             (operatorType == OPERATOR_D5)    ||
00248                             (operatorType == OPERATOR_D6)    ||
00249                             (operatorType == OPERATOR_D7)    ||
00250                             (operatorType == OPERATOR_D8)    ||
00251                             (operatorType == OPERATOR_D9)    ||
00252                             (operatorType == OPERATOR_D10) ),
00253                           std::invalid_argument,
00254                           ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type");
00255       break;
00256 
00257     default:
00258       TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
00259                             (operatorType != OPERATOR_GRAD)  &&
00260                             (operatorType != OPERATOR_CURL)  &&
00261                             (operatorType != OPERATOR_DIV)   &&
00262                             (operatorType != OPERATOR_D1)    &&
00263                             (operatorType != OPERATOR_D2)    &&
00264                             (operatorType != OPERATOR_D3)    &&
00265                             (operatorType != OPERATOR_D4)    &&
00266                             (operatorType != OPERATOR_D5)    &&
00267                             (operatorType != OPERATOR_D6)    &&
00268                             (operatorType != OPERATOR_D7)    &&
00269                             (operatorType != OPERATOR_D8)    &&
00270                             (operatorType != OPERATOR_D9)    &&
00271                             (operatorType != OPERATOR_D10) ),
00272                           std::invalid_argument,
00273                           ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type");
00274   }
00275 }
00276 
00277 
00278   
00279 template<class Scalar, class ArrayScalar>
00280 void Basis_HCURL_HEX_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00281                                                             const ArrayScalar &    inputPoints,
00282                                                             const ArrayScalar &    cellVertices,
00283                                                             const EOperator        operatorType) const {
00284   TEST_FOR_EXCEPTION( (true), std::logic_error,
00285                       ">>> ERROR (Basis_HCURL_HEX_I1_FEM): FEM Basis calling an FVD member function");
00286                                                              }
00287 
00288 }// namespace Intrepid
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_HEX_I1_FEMDef.hpp
