|
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_HDIV_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_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() * openBasis_.getCardinality() * openBasis_.getCardinality(); 00049 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00050 this -> basisType_ = BASIS_FEM_FIAT; 00051 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00052 this -> basisTagsAreSet_ = false; 00053 } 00054 00055 template<class Scalar, class ArrayScalar> 00056 Basis_HDIV_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_HEX_In_FEM( int order , const EPointType &pointType ): 00057 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ), 00058 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ) 00059 { 00060 this -> basisDegree_ = order; 00061 this -> basisCardinality_ = 3 * closedBasis_.getCardinality() * openBasis_.getCardinality() * openBasis_.getCardinality(); 00062 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00063 this -> basisType_ = BASIS_FEM_FIAT; 00064 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00065 this -> basisTagsAreSet_ = false; 00066 } 00067 00068 00069 00070 template<class Scalar, class ArrayScalar> 00071 void Basis_HDIV_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 // four vertices 00088 for (int i=0;i<4;i++) { 00089 total_dof_per_entity[0][i] = 0; 00090 current_dof_per_entity[0][i] = 0; 00091 } 00092 // twelve edges 00093 for (int i=0;i<12;i++) { 00094 total_dof_per_entity[1][i] = 0; 00095 current_dof_per_entity[1][i] = 0; 00096 } 00097 // six faces 00098 for (int i=0;i<6;i++) { 00099 total_dof_per_entity[2][i] = 0; 00100 current_dof_per_entity[2][i] = 0; 00101 } 00102 total_dof_per_entity[3][0] = 0; 00103 current_dof_per_entity[3][0] = 0; 00104 00105 // tally dof on each facet. none on vertex or edge 00106 for (int i=0;i<6;i++) { 00107 total_dof_per_entity[2][i] = openBasis_.getCardinality() * openBasis_.getCardinality(); 00108 } 00109 00110 total_dof_per_entity[2][0] = this->getCardinality() - 6 * openBasis_.getCardinality()* openBasis_.getCardinality(); 00111 00112 int tagcur = 0; 00113 // loop over the x-component basis functions, which are (psi(x)phi(y)phi(z),0,0) 00114 // for psi in the closed basis and phi in the open 00115 for (int k=0;k<openBasis_.getCardinality();k++) { 00116 const int odimk = openDofTags[k][0]; 00117 const int oentk = openDofTags[k][1]; 00118 for (int j=0;j<openBasis_.getCardinality();j++) { 00119 const int odimj = openDofTags[j][0]; 00120 const int oentj = openDofTags[j][1]; 00121 for (int i=0;i<closedBasis_.getCardinality();i++) { 00122 const int cdim = closedDofTags[i][0]; 00123 const int cent = closedDofTags[i][1]; 00124 int dofdim; 00125 int dofent; 00126 ProductTopology::lineProduct3d(cdim,cent,odimj,oentj,odimk,oentk,dofdim,dofent); 00127 tags[4*tagcur] = dofdim; 00128 tags[4*tagcur+1] = dofent; 00129 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00130 current_dof_per_entity[dofdim][dofent]++; 00131 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00132 tagcur++; 00133 } 00134 } 00135 } 00136 00137 // now we have to do it for the y-component basis functions, which are 00138 // (0,phi(x)psi(y)phi(z),0) for psi in the closed basis and phi in the open 00139 for (int k=0;k<openBasis_.getCardinality();k++) { 00140 const int odimk = openDofTags[k][0]; 00141 const int oentk = openDofTags[k][1]; 00142 for (int j=0;j<closedBasis_.getCardinality();j++) { 00143 const int cdim = closedDofTags[j][0]; 00144 const int cent = closedDofTags[j][1]; 00145 for (int i=0;i<openBasis_.getCardinality();i++) { 00146 const int odimi = openDofTags[i][0]; 00147 const int oenti = openDofTags[i][1]; 00148 int dofdim; 00149 int dofent; 00150 ProductTopology::lineProduct3d(odimi,oenti,cdim,cent,odimk,oentk,dofdim,dofent); 00151 tags[4*tagcur] = dofdim; 00152 tags[4*tagcur+1] = dofent; 00153 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00154 current_dof_per_entity[dofdim][dofent]++; 00155 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00156 tagcur++; 00157 } 00158 } 00159 } 00160 00161 // now we have to do it for the z-component basis functions, which are 00162 // (0,0,phi(x)phi(y)psi(z)) for psi in the closed basis and phi in the open 00163 for (int k=0;k<closedBasis_.getCardinality();k++) { 00164 const int cdim = closedDofTags[k][0]; 00165 const int cent = closedDofTags[k][1]; 00166 for (int j=0;j<openBasis_.getCardinality();j++) { 00167 const int odimj = openDofTags[j][0]; 00168 const int oentj = openDofTags[j][1]; 00169 for (int i=0;i<openBasis_.getCardinality();i++) { 00170 const int odimi = openDofTags[i][0]; 00171 const int oenti = openDofTags[i][1]; 00172 int dofdim; 00173 int dofent; 00174 ProductTopology::lineProduct3d(odimi,oenti,odimj,oentj,cdim,cent,dofdim,dofent); 00175 tags[4*tagcur] = dofdim; 00176 tags[4*tagcur+1] = dofent; 00177 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00178 current_dof_per_entity[dofdim][dofent]++; 00179 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00180 tagcur++; 00181 } 00182 } 00183 } 00184 00185 // for (int i=0;i<this->getCardinality();i++) { 00186 // for (int j=0;j<4;j++) { 00187 // std::cout << tags[4*i+j] << " "; 00188 // } 00189 // std::cout << std::endl; 00190 // } 00191 00192 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00193 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00194 this -> ordinalToTag_, 00195 &(tags[0]), 00196 this -> basisCardinality_, 00197 tagSize, 00198 posScDim, 00199 posScOrd, 00200 posDfOrd); 00201 } 00202 00203 00204 template<class Scalar, class ArrayScalar> 00205 void Basis_HDIV_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00206 const ArrayScalar & inputPoints, 00207 const EOperator operatorType) const { 00208 00209 // Verify arguments 00210 #ifdef HAVE_INTREPID_DEBUG 00211 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues, 00212 inputPoints, 00213 operatorType, 00214 this -> getBaseCellTopology(), 00215 this -> getCardinality() ); 00216 #endif 00217 00218 // Number of evaluation points = dim 0 of inputPoints 00219 int dim0 = inputPoints.dimension(0); 00220 00221 // separate out points 00222 FieldContainer<Scalar> xPoints(dim0,1); 00223 FieldContainer<Scalar> yPoints(dim0,1); 00224 FieldContainer<Scalar> zPoints(dim0,1); 00225 00226 00227 for (int i=0;i<dim0;i++) { 00228 xPoints(i,0) = inputPoints(i,0); 00229 yPoints(i,0) = inputPoints(i,1); 00230 zPoints(i,0) = inputPoints(i,2); 00231 } 00232 00233 switch (operatorType) { 00234 case OPERATOR_VALUE: 00235 { 00236 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00237 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00238 FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 ); 00239 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00240 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00241 FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 ); 00242 00243 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00244 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00245 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE ); 00246 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00247 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00248 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00249 00250 int bfcur = 0; 00251 // x component bfs are (closed(x) open(y) open(z),0,0) 00252 for (int k=0;k<openBasis_.getCardinality();k++) { 00253 for (int j=0;j<openBasis_.getCardinality();j++) { 00254 for (int i=0;i<closedBasis_.getCardinality();i++) { 00255 for (int l=0;l<dim0;l++) { 00256 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00257 outputValues(bfcur,l,1) = 0.0; 00258 outputValues(bfcur,l,2) = 0.0; 00259 } 00260 bfcur++; 00261 } 00262 } 00263 } 00264 00265 // y component bfs are (0,open(x) closed(y) open(z),0) 00266 for (int k=0;k<openBasis_.getCardinality();k++) { 00267 for (int j=0;j<closedBasis_.getCardinality();j++) { 00268 for (int i=0;i<openBasis_.getCardinality();i++) { 00269 for (int l=0;l<dim0;l++) { 00270 outputValues(bfcur,l,0) = 0.0; 00271 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00272 outputValues(bfcur,l,2) = 0.0; 00273 } 00274 bfcur++; 00275 } 00276 } 00277 } 00278 00279 // z component bfs are (0,0,open(x) open(y) closed(z)) 00280 for (int k=0;k<closedBasis_.getCardinality();k++) { 00281 for (int j=0;j<openBasis_.getCardinality();j++) { 00282 for (int i=0;i<openBasis_.getCardinality();i++) { 00283 for (int l=0;l<dim0;l++) { 00284 outputValues(bfcur,l,0) = 0.0; 00285 outputValues(bfcur,l,1) = 0.0; 00286 outputValues(bfcur,l,2) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00287 } 00288 bfcur++; 00289 } 00290 } 00291 } 00292 00293 00294 } 00295 break; 00296 case OPERATOR_DIV: 00297 { 00298 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00299 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00300 FieldContainer<Scalar> closedBasisDerivsZPts( closedBasis_.getCardinality() , dim0 , 1 ); 00301 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00302 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00303 FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 ); 00304 00305 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00306 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00307 closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 ); 00308 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00309 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00310 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00311 00312 int bfcur = 0; 00313 00314 // x component basis functions first 00315 for (int k=0;k<openBasis_.getCardinality();k++) { 00316 for (int j=0;j<openBasis_.getCardinality();j++) { 00317 for (int i=0;i<closedBasis_.getCardinality();i++) { 00318 for (int l=0;l<dim0;l++) { 00319 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00320 } 00321 bfcur++; 00322 } 00323 } 00324 } 00325 00326 // now y component basis functions 00327 for (int k=0;k<openBasis_.getCardinality();k++) { 00328 for (int j=0;j<closedBasis_.getCardinality();j++) { 00329 for (int i=0;i<openBasis_.getCardinality();i++) { 00330 for (int l=0;l<dim0;l++) { 00331 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l); 00332 } 00333 bfcur++; 00334 } 00335 } 00336 } 00337 00338 // now z component basis functions 00339 for (int k=0;k<closedBasis_.getCardinality();k++) { 00340 for (int j=0;j<openBasis_.getCardinality();j++) { 00341 for (int i=0;i<openBasis_.getCardinality();i++) { 00342 for (int l=0;l<dim0;l++) { 00343 outputValues(bfcur,l) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0); 00344 } 00345 bfcur++; 00346 } 00347 } 00348 } 00349 } 00350 break; 00351 case OPERATOR_CURL: 00352 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00353 ">>> ERROR (Basis_HDIV_HEX_In_FEM): CURL is invalid operator for HDIV Basis Functions"); 00354 break; 00355 00356 case OPERATOR_GRAD: 00357 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00358 ">>> ERROR (Basis_HDIV_HEX_In_FEM): GRAD is invalid operator for HDIV Basis Functions"); 00359 break; 00360 00361 case OPERATOR_D1: 00362 case OPERATOR_D2: 00363 case OPERATOR_D3: 00364 case OPERATOR_D4: 00365 case OPERATOR_D5: 00366 case OPERATOR_D6: 00367 case OPERATOR_D7: 00368 case OPERATOR_D8: 00369 case OPERATOR_D9: 00370 case OPERATOR_D10: 00371 TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00372 (operatorType == OPERATOR_D2) || 00373 (operatorType == OPERATOR_D3) || 00374 (operatorType == OPERATOR_D4) || 00375 (operatorType == OPERATOR_D5) || 00376 (operatorType == OPERATOR_D6) || 00377 (operatorType == OPERATOR_D7) || 00378 (operatorType == OPERATOR_D8) || 00379 (operatorType == OPERATOR_D9) || 00380 (operatorType == OPERATOR_D10) ), 00381 std::invalid_argument, 00382 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type"); 00383 break; 00384 00385 default: 00386 TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00387 (operatorType != OPERATOR_GRAD) && 00388 (operatorType != OPERATOR_CURL) && 00389 (operatorType != OPERATOR_DIV) && 00390 (operatorType != OPERATOR_D1) && 00391 (operatorType != OPERATOR_D2) && 00392 (operatorType != OPERATOR_D3) && 00393 (operatorType != OPERATOR_D4) && 00394 (operatorType != OPERATOR_D5) && 00395 (operatorType != OPERATOR_D6) && 00396 (operatorType != OPERATOR_D7) && 00397 (operatorType != OPERATOR_D8) && 00398 (operatorType != OPERATOR_D9) && 00399 (operatorType != OPERATOR_D10) ), 00400 std::invalid_argument, 00401 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type"); 00402 } 00403 } 00404 00405 00406 00407 template<class Scalar, class ArrayScalar> 00408 void Basis_HDIV_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00409 const ArrayScalar & inputPoints, 00410 const ArrayScalar & cellVertices, 00411 const EOperator operatorType) const { 00412 TEST_FOR_EXCEPTION( (true), std::logic_error, 00413 ">>> ERROR (Basis_HDIV_HEX_In_FEM): FEM Basis calling an FVD member function"); 00414 } 00415 00416 }// namespace Intrepid
1.7.4