Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_WEDGE_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_WEDGE_I1_FEM<Scalar,ArrayScalar>::Basis_HCURL_WEDGE_I1_FEM()
00040   {
00041     this -> basisCardinality_  = 9;
00042     this -> basisDegree_       = 1;
00043     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >() );
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_WEDGE_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   
00070   // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00071   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00072                               this -> ordinalToTag_,
00073                               tags,
00074                               this -> basisCardinality_,
00075                               tagSize,
00076                               posScDim,
00077                               posScOrd,
00078                               posDfOrd);
00079 }
00080 
00081 
00082 
00083 template<class Scalar, class ArrayScalar>
00084 void Basis_HCURL_WEDGE_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00085                                                             const ArrayScalar &  inputPoints,
00086                                                             const EOperator      operatorType) const {
00087 
00088 // Verify arguments
00089 #ifdef HAVE_INTREPID_DEBUG
00090   Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
00091                                                       inputPoints,
00092                                                       operatorType,
00093                                                       this -> getBaseCellTopology(),
00094                                                       this -> getCardinality() );
00095 #endif
00096   
00097  // Number of evaluation points = dim 0 of inputPoints
00098   int dim0 = inputPoints.dimension(0);
00099 
00100   // Temporaries: (x,y,z) coordinates of the evaluation point
00101   Scalar x = 0.0;                                    
00102   Scalar y = 0.0;                                    
00103   Scalar z = 0.0;                                    
00104   
00105   switch (operatorType) {
00106     case OPERATOR_VALUE:
00107       for (int i0 = 0; i0 < dim0; i0++) {
00108         x = inputPoints(i0, 0);
00109         y = inputPoints(i0, 1);
00110         z = inputPoints(i0, 2);
00111         
00112         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00113         outputValues(0, i0, 0) = (1.0 - z)*(1.0 - y)/2.0;
00114         outputValues(0, i0, 1) = x*(1.0 - z)/2.0;
00115         outputValues(0, i0, 2) = 0.0;
00116 
00117         outputValues(1, i0, 0) = y*(z - 1.0)/2.0;
00118         outputValues(1, i0, 1) = x*(1.0 - z)/2.0;
00119         outputValues(1, i0, 2) = 0.0;
00120 
00121         outputValues(2, i0, 0) = y*(z - 1.0)/2.0;
00122         outputValues(2, i0, 1) = (1.0 - x)*(z - 1.0)/2.0;
00123         outputValues(2, i0, 2) = 0.0;
00124 
00125         outputValues(3, i0, 0) = (1.0 - y)*(1.0 + z)/2.0;
00126         outputValues(3, i0, 1) = x*(1.0 + z)/2.0;
00127         outputValues(3, i0, 2) = 0.0;
00128 
00129         outputValues(4, i0, 0) =-y*(1.0 + z)/2.0;
00130         outputValues(4, i0, 1) = x*(1.0 + z)/2.0;
00131         outputValues(4, i0, 2) = 0.0;
00132 
00133         outputValues(5, i0, 0) = -y*(1.0 + z)/2.0;
00134         outputValues(5, i0, 1) = (x - 1.0)*(1.0 + z)/2.0;
00135         outputValues(5, i0, 2) = 0.0;
00136 
00137         outputValues(6, i0, 0) = 0.0;
00138         outputValues(6, i0, 1) = 0.0;
00139         outputValues(6, i0, 2) = (1.0 - x - y)/2.0;
00140 
00141         outputValues(7, i0, 0) = 0.0;
00142         outputValues(7, i0, 1) = 0.0;
00143         outputValues(7, i0, 2) = x/2.0;
00144 
00145         outputValues(8, i0, 0) = 0.0;
00146         outputValues(8, i0, 1) = 0.0;
00147         outputValues(8, i0, 2) = y/2.0;
00148 
00149       }
00150       break;
00151       
00152     case OPERATOR_CURL:
00153       for (int i0 = 0; i0 < dim0; i0++) {
00154         x = inputPoints(i0, 0);
00155         y = inputPoints(i0, 1);
00156         z = inputPoints(i0, 2);
00157         
00158         // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
00159         outputValues(0, i0, 0) = x/2.0;
00160         outputValues(0, i0, 1) = (y - 1.0)/2.0;
00161         outputValues(0, i0, 2) = 1.0 - z;
00162 
00163         outputValues(1, i0, 0) = x/2.0;
00164         outputValues(1, i0, 1) = y/2.0;
00165         outputValues(1, i0, 2) = 1.0 - z;
00166 
00167         outputValues(2, i0, 0) = (x - 1.0)/2.0;
00168         outputValues(2, i0, 1) = y/2.0; 
00169         outputValues(2, i0, 2) = 1.0 - z;
00170 
00171         outputValues(3, i0, 0) = -x/2.0;
00172         outputValues(3, i0, 1) = (1.0 - y)/2.0;
00173         outputValues(3, i0, 2) = 1.0 + z;
00174 
00175         outputValues(4, i0, 0) = -x/2.0;
00176         outputValues(4, i0, 1) = -y/2.0;
00177         outputValues(4, i0, 2) = 1.0 + z;
00178 
00179         outputValues(5, i0, 0) = (1.0 - x)/2.0;
00180         outputValues(5, i0, 1) = -y/2.0;
00181         outputValues(5, i0, 2) = 1.0 + z;
00182 
00183         outputValues(6, i0, 0) =-0.5;
00184         outputValues(6, i0, 1) = 0.5;
00185         outputValues(6, i0, 2) = 0.0;
00186 
00187         outputValues(7, i0, 0) = 0.0;
00188         outputValues(7, i0, 1) =-0.5;
00189         outputValues(7, i0, 2) = 0.0;
00190 
00191         outputValues(8, i0, 0) = 0.5;
00192         outputValues(8, i0, 1) = 0.0;
00193         outputValues(8, i0, 2) = 0.0;
00194       }
00195       break;
00196       
00197     case OPERATOR_DIV:
00198        TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
00199                           ">>> ERROR (Basis_HCURL_WEDGE_I1_FEM): DIV is invalid operator for HCURL Basis Functions");
00200       break;
00201       
00202     case OPERATOR_GRAD:
00203        TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
00204                           ">>> ERROR (Basis_HCURL_WEDGE_I1_FEM): GRAD is invalid operator for HCURL Basis Functions");
00205       break;
00206 
00207     case OPERATOR_D1:
00208     case OPERATOR_D2:
00209     case OPERATOR_D3:
00210     case OPERATOR_D4:
00211     case OPERATOR_D5:
00212     case OPERATOR_D6:
00213     case OPERATOR_D7:
00214     case OPERATOR_D8:
00215     case OPERATOR_D9:
00216     case OPERATOR_D10:
00217       TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1)    ||
00218                             (operatorType == OPERATOR_D2)    ||
00219                             (operatorType == OPERATOR_D3)    ||
00220                             (operatorType == OPERATOR_D4)    ||
00221                             (operatorType == OPERATOR_D5)    ||
00222                             (operatorType == OPERATOR_D6)    ||
00223                             (operatorType == OPERATOR_D7)    ||
00224                             (operatorType == OPERATOR_D8)    ||
00225                             (operatorType == OPERATOR_D9)    ||
00226                             (operatorType == OPERATOR_D10) ),
00227                           std::invalid_argument,
00228                           ">>> ERROR (Basis_HCURL_WEDGE_I1_FEM): Invalid operator type");
00229       break;
00230 
00231     default:
00232       TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
00233                             (operatorType != OPERATOR_GRAD)  &&
00234                             (operatorType != OPERATOR_CURL)  &&
00235                             (operatorType != OPERATOR_DIV)   &&
00236                             (operatorType != OPERATOR_D1)    &&
00237                             (operatorType != OPERATOR_D2)    &&
00238                             (operatorType != OPERATOR_D3)    &&
00239                             (operatorType != OPERATOR_D4)    &&
00240                             (operatorType != OPERATOR_D5)    &&
00241                             (operatorType != OPERATOR_D6)    &&
00242                             (operatorType != OPERATOR_D7)    &&
00243                             (operatorType != OPERATOR_D8)    &&
00244                             (operatorType != OPERATOR_D9)    &&
00245                             (operatorType != OPERATOR_D10) ),
00246                           std::invalid_argument,
00247                           ">>> ERROR (Basis_HCURL_WEDGE_I1_FEM): Invalid operator type");
00248   }
00249 }
00250 
00251 
00252   
00253 template<class Scalar, class ArrayScalar>
00254 void Basis_HCURL_WEDGE_I1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00255                                                             const ArrayScalar &    inputPoints,
00256                                                             const ArrayScalar &    cellVertices,
00257                                                             const EOperator        operatorType) const {
00258   TEST_FOR_EXCEPTION( (true), std::logic_error,
00259                       ">>> ERROR (Basis_HCURL_WEDGE_I1_FEM): FEM Basis calling an FVD member function");
00260                                                              }
00261 
00262 }// namespace Intrepid
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/Basis/Intrepid_HCURL_WEDGE_I1_FEMDef.hpp
