|
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 00017 00018 // WITHOUT ANY WARRANTY; without even the implied warranty of 00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00020 // Lesser General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU Lesser General Public 00023 // License along with this library; if not, write to the Free Software 00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00025 // USA 00026 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00027 // Denis Ridzal (dridzal@sandia.gov) 00028 // Kara Peterson (kjpeter@sandia.gov). 00029 // 00030 // ************************************************************************ 00031 // @HEADER 00032 00038 namespace Intrepid { 00039 00040 template<class Scalar, class ArrayScalar> 00041 Basis_HCURL_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_In_FEM( int order , 00042 const ArrayScalar & ptsClosed , 00043 const ArrayScalar & ptsOpen): 00044 closedBasis_( order , ptsClosed ), 00045 openBasis_( order-1 , ptsOpen ) 00046 { 00047 this -> basisDegree_ = order; 00048 this -> basisCardinality_ = 3 * closedBasis_.getCardinality() 00049 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 00050 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00051 this -> basisType_ = BASIS_FEM_FIAT; 00052 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00053 this -> basisTagsAreSet_ = false; 00054 } 00055 00056 template<class Scalar, class ArrayScalar> 00057 Basis_HCURL_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_In_FEM( int order , const EPointType &pointType ): 00058 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ), 00059 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ) 00060 { 00061 this -> basisDegree_ = order; 00062 this -> basisCardinality_ = 3 * closedBasis_.getCardinality() 00063 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 00064 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00065 this -> basisType_ = BASIS_FEM_FIAT; 00066 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00067 this -> basisTagsAreSet_ = false; 00068 } 00069 00070 template<class Scalar, class ArrayScalar> 00071 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::initializeTags() { 00072 00073 // Basis-dependent intializations 00074 int tagSize = 4; // size of DoF tag 00075 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00076 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00077 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00078 00079 std::vector<int> tags( tagSize * this->getCardinality() ); 00080 00081 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags(); 00082 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags(); 00083 00084 std::map<int,std::map<int,int> > total_dof_per_entity; 00085 std::map<int,std::map<int,int> > current_dof_per_entity; 00086 00087 const int order = this->basisDegree_; 00088 00089 // vertices 00090 for (int i=0;i<4;i++) { 00091 total_dof_per_entity[0][i] = 0; 00092 current_dof_per_entity[0][i] = 0; 00093 } 00094 // edges 00095 for (int i=0;i<12;i++) { 00096 total_dof_per_entity[1][i] = 0; 00097 current_dof_per_entity[1][i] = 0; 00098 } 00099 // faces 00100 for (int i=0;i<6;i++) { 00101 total_dof_per_entity[2][i] = 0; 00102 current_dof_per_entity[2][i] = 0; 00103 } 00104 total_dof_per_entity[3][0] = 0; 00105 current_dof_per_entity[3][0] = 0; 00106 00107 // tally dof on each facet. none on vertex 00108 // edge dof 00109 for (int i=0;i<12;i++) { 00110 total_dof_per_entity[1][i] = order; 00111 } 00112 // face dof 00113 for (int i=0;i<6;i++) { 00114 total_dof_per_entity[2][i] = 2 * (order - 1) * order; 00115 } 00116 // internal dof 00117 total_dof_per_entity[3][0] = this->basisCardinality_ - 12 * order - 6 * total_dof_per_entity[2][0]; 00118 00119 int tagcur = 0; 00120 for (int k=0;k<closedBasis_.getCardinality();k++) { 00121 const int dimk = closedDofTags[k][0]; 00122 const int entk = closedDofTags[k][1]; 00123 for (int j=0;j<closedBasis_.getCardinality();j++) { 00124 const int dimj = closedDofTags[j][0]; 00125 const int entj = closedDofTags[j][1]; 00126 for (int i=0;i<openBasis_.getCardinality();i++) { 00127 const int dimi = openDofTags[i][0]; 00128 const int enti = openDofTags[i][1]; 00129 int dofdim; 00130 int dofent; 00131 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00132 tags[4*tagcur] = dofdim; 00133 tags[4*tagcur+1] = dofent; 00134 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00135 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00136 current_dof_per_entity[dofdim][dofent]++; 00137 tagcur++; 00138 } 00139 } 00140 } 00141 for (int k=0;k<closedBasis_.getCardinality();k++) { 00142 const int dimk = closedDofTags[k][0]; 00143 const int entk = closedDofTags[k][1]; 00144 for (int j=0;j<openBasis_.getCardinality();j++) { 00145 const int dimj = openDofTags[j][0]; 00146 const int entj = openDofTags[j][1]; 00147 for (int i=0;i<closedBasis_.getCardinality();i++) { 00148 const int dimi = closedDofTags[i][0]; 00149 const int enti = closedDofTags[i][1]; 00150 int dofdim; 00151 int dofent; 00152 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00153 tags[4*tagcur] = dofdim; 00154 tags[4*tagcur+1] = dofent; 00155 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00156 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00157 current_dof_per_entity[dofdim][dofent]++; 00158 tagcur++; 00159 } 00160 } 00161 } 00162 for (int k=0;k<openBasis_.getCardinality();k++) { 00163 const int dimk = openDofTags[k][0]; 00164 const int entk = openDofTags[k][1]; 00165 for (int j=0;j<closedBasis_.getCardinality();j++) { 00166 const int dimj = closedDofTags[j][0]; 00167 const int entj = closedDofTags[j][1]; 00168 for (int i=0;i<closedBasis_.getCardinality();i++) { 00169 const int dimi = closedDofTags[i][0]; 00170 const int enti = closedDofTags[i][1]; 00171 int dofdim; 00172 int dofent; 00173 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00174 tags[4*tagcur] = dofdim; 00175 tags[4*tagcur+1] = dofent; 00176 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00177 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00178 current_dof_per_entity[dofdim][dofent]++; 00179 tagcur++; 00180 } 00181 } 00182 } 00183 00184 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00185 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00186 this -> ordinalToTag_, 00187 &(tags[0]), 00188 this -> basisCardinality_, 00189 tagSize, 00190 posScDim, 00191 posScOrd, 00192 posDfOrd); 00193 } 00194 00195 00196 template<class Scalar, class ArrayScalar> 00197 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00198 const ArrayScalar & inputPoints, 00199 const EOperator operatorType) const { 00200 00201 // Verify arguments 00202 #ifdef HAVE_INTREPID_DEBUG 00203 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues, 00204 inputPoints, 00205 operatorType, 00206 this -> getBaseCellTopology(), 00207 this -> getCardinality() ); 00208 #endif 00209 00210 // Number of evaluation points = dim 0 of inputPoints 00211 int dim0 = inputPoints.dimension(0); 00212 00213 // separate out points 00214 FieldContainer<Scalar> xPoints(dim0,1); 00215 FieldContainer<Scalar> yPoints(dim0,1); 00216 FieldContainer<Scalar> zPoints(dim0,1); 00217 00218 for (int i=0;i<dim0;i++) { 00219 xPoints(i,0) = inputPoints(i,0); 00220 yPoints(i,0) = inputPoints(i,1); 00221 zPoints(i,0) = inputPoints(i,2); 00222 } 00223 00224 switch (operatorType) { 00225 case OPERATOR_VALUE: 00226 { 00227 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00228 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00229 FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 ); 00230 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00231 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00232 FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 ); 00233 00234 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00235 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00236 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE ); 00237 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00238 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00239 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00240 00241 // first we get the "horizontal basis functions that are tangent to the x-varying edges 00242 int bfcur = 0; 00243 for (int k=0;k<closedBasis_.getCardinality();k++) { 00244 for (int j=0;j<closedBasis_.getCardinality();j++) { 00245 for (int i=0;i<openBasis_.getCardinality();i++) { 00246 for (int l=0;l<dim0;l++) { 00247 outputValues(bfcur,l,0) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00248 outputValues(bfcur,l,1) = 0.0; 00249 outputValues(bfcur,l,2) = 0.0; 00250 } 00251 bfcur++; 00252 } 00253 } 00254 } 00255 00256 // second we get the basis functions in the direction of the y-varying edges 00257 for (int k=0;k<closedBasis_.getCardinality();k++) { 00258 for (int j=0;j<openBasis_.getCardinality();j++) { 00259 for (int i=0;i<closedBasis_.getCardinality();i++) { 00260 for (int l=0;l<dim0;l++) { 00261 outputValues(bfcur,l,0) = 0.0; 00262 outputValues(bfcur,l,1) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00263 outputValues(bfcur,l,2) = 0.0; 00264 } 00265 bfcur++; 00266 } 00267 } 00268 } 00269 00270 // third we get the basis functions in the direction of the y-varying edges 00271 for (int k=0;k<openBasis_.getCardinality();k++) { 00272 for (int j=0;j<closedBasis_.getCardinality();j++) { 00273 for (int i=0;i<closedBasis_.getCardinality();i++) { 00274 for (int l=0;l<dim0;l++) { 00275 outputValues(bfcur,l,0) = 0.0; 00276 outputValues(bfcur,l,1) = 0.0; 00277 outputValues(bfcur,l,2) = closedBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00278 } 00279 bfcur++; 00280 } 00281 } 00282 } 00283 } 00284 break; 00285 case OPERATOR_CURL: 00286 { 00287 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00288 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00289 FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 ); 00290 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00291 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00292 FieldContainer<Scalar> closedBasisDerivsZPts( closedBasis_.getCardinality() , dim0 , 1 ); 00293 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00294 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00295 FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 ); 00296 00297 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00298 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00299 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE ); 00300 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00301 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00302 closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 ); 00303 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00304 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00305 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00306 00307 int bfcur = 0; 00308 00309 // first we get the basis functions that are tangent to the x-varying edges 00310 for (int k=0;k<closedBasis_.getCardinality();k++) { 00311 for (int j=0;j<closedBasis_.getCardinality();j++) { 00312 for (int i=0;i<openBasis_.getCardinality();i++) { 00313 for (int l=0;l<dim0;l++) { 00314 outputValues(bfcur,l,0) = 0.0; 00315 outputValues(bfcur,l,1) = -openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0); 00316 outputValues(bfcur,l,2) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * closedBasisValsZPts(k,l); 00317 } 00318 bfcur++; 00319 } 00320 } 00321 } 00322 00323 // second we get the basis functions in the direction of the y-varying edges 00324 for (int k=0;k<closedBasis_.getCardinality();k++) { 00325 for (int j=0;j<openBasis_.getCardinality();j++) { 00326 for (int i=0;i<closedBasis_.getCardinality();i++) { 00327 for (int l=0;l<dim0;l++) { 00328 outputValues(bfcur,l,0) = -closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0); 00329 outputValues(bfcur,l,1) = 0.0; 00330 outputValues(bfcur,l,2) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00331 } 00332 bfcur++; 00333 } 00334 } 00335 } 00336 00337 // third we get the basis functions in the direction of the y-varying edges 00338 for (int k=0;k<openBasis_.getCardinality();k++) { 00339 for (int j=0;j<closedBasis_.getCardinality();j++) { 00340 for (int i=0;i<closedBasis_.getCardinality();i++) { 00341 for (int l=0;l<dim0;l++) { 00342 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l); 00343 outputValues(bfcur,l,1) = -closedBasisDerivsXPts(i,l,0) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00344 outputValues(bfcur,l,2) = 0.0; 00345 } 00346 bfcur++; 00347 } 00348 } 00349 } 00350 } 00351 break; 00352 case OPERATOR_DIV: 00353 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00354 ">>> ERROR (Basis_HCURL_HEX_In_FEM): DIV is invalid operator for HCURL Basis Functions"); 00355 break; 00356 00357 case OPERATOR_GRAD: 00358 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00359 ">>> ERROR (Basis_HCURL_HEX_In_FEM): GRAD is invalid operator for HCURL Basis Functions"); 00360 break; 00361 00362 case OPERATOR_D1: 00363 case OPERATOR_D2: 00364 case OPERATOR_D3: 00365 case OPERATOR_D4: 00366 case OPERATOR_D5: 00367 case OPERATOR_D6: 00368 case OPERATOR_D7: 00369 case OPERATOR_D8: 00370 case OPERATOR_D9: 00371 case OPERATOR_D10: 00372 TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00373 (operatorType == OPERATOR_D2) || 00374 (operatorType == OPERATOR_D3) || 00375 (operatorType == OPERATOR_D4) || 00376 (operatorType == OPERATOR_D5) || 00377 (operatorType == OPERATOR_D6) || 00378 (operatorType == OPERATOR_D7) || 00379 (operatorType == OPERATOR_D8) || 00380 (operatorType == OPERATOR_D9) || 00381 (operatorType == OPERATOR_D10) ), 00382 std::invalid_argument, 00383 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type"); 00384 break; 00385 00386 default: 00387 TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00388 (operatorType != OPERATOR_GRAD) && 00389 (operatorType != OPERATOR_CURL) && 00390 (operatorType != OPERATOR_CURL) && 00391 (operatorType != OPERATOR_D1) && 00392 (operatorType != OPERATOR_D2) && 00393 (operatorType != OPERATOR_D3) && 00394 (operatorType != OPERATOR_D4) && 00395 (operatorType != OPERATOR_D5) && 00396 (operatorType != OPERATOR_D6) && 00397 (operatorType != OPERATOR_D7) && 00398 (operatorType != OPERATOR_D8) && 00399 (operatorType != OPERATOR_D9) && 00400 (operatorType != OPERATOR_D10) ), 00401 std::invalid_argument, 00402 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type"); 00403 } 00404 } 00405 00406 00407 00408 template<class Scalar, class ArrayScalar> 00409 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00410 const ArrayScalar & inputPoints, 00411 const ArrayScalar & cellVertices, 00412 const EOperator operatorType) const { 00413 TEST_FOR_EXCEPTION( (true), std::logic_error, 00414 ">>> ERROR (Basis_HCURL_HEX_In_FEM): FEM Basis calling an FVD member function"); 00415 } 00416 00417 }// namespace Intrepid
1.7.4