|
Intrepid
|
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
1.7.4