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