|
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 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00049 namespace Intrepid { 00050 00051 template<class Scalar, class ArrayScalar> 00052 Basis_HCURL_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_In_FEM( int order , 00053 const ArrayScalar & ptsClosed , 00054 const ArrayScalar & ptsOpen): 00055 closedBasis_( order , ptsClosed ), 00056 openBasis_( order-1 , ptsOpen ) , 00057 closedPts_( ptsClosed ), 00058 openPts_( ptsOpen ) 00059 { 00060 this -> basisDegree_ = order; 00061 this -> basisCardinality_ = 3 * closedBasis_.getCardinality() 00062 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 00063 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00064 this -> basisType_ = BASIS_FEM_FIAT; 00065 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00066 this -> basisTagsAreSet_ = false; 00067 00068 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3); 00069 bases[0].resize(3); bases[1].resize(3); bases[2].resize(3); 00070 bases[0][0] = rcp( &openBasis_ , false ); 00071 bases[0][1] = rcp( &closedBasis_ , false ); 00072 bases[0][2] = rcp( &closedBasis_ , false ); 00073 bases[1][0] = rcp( &closedBasis_ , false ); 00074 bases[1][1] = rcp( &openBasis_ , false ); 00075 bases[1][2] = rcp( &closedBasis_ , false ); 00076 bases[2][0] = rcp( &closedBasis_ , false ); 00077 bases[2][1] = rcp( &closedBasis_ , false ); 00078 bases[2][2] = rcp( &openBasis_ , false ); 00079 this->setBases( bases ); 00080 00081 } 00082 00083 template<class Scalar, class ArrayScalar> 00084 Basis_HCURL_HEX_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_HEX_In_FEM( int order , const EPointType &pointType ): 00085 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ), 00086 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ), 00087 closedPts_( order+1 , 1 ), 00088 openPts_( order , 1 ) 00089 { 00090 this -> basisDegree_ = order; 00091 this -> basisCardinality_ = 3 * closedBasis_.getCardinality() 00092 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 00093 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00094 this -> basisType_ = BASIS_FEM_FIAT; 00095 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00096 this -> basisTagsAreSet_ = false; 00097 00098 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ , 00099 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) , 00100 order , 00101 0 , 00102 pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED ); 00103 00104 if (pointType == POINTTYPE_SPECTRAL) 00105 { 00106 PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ , 00107 order - 1 ); 00108 } 00109 else 00110 { 00111 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ , 00112 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) , 00113 order - 1, 00114 0 , 00115 POINTTYPE_EQUISPACED ); 00116 } 00117 00118 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3); 00119 bases[0].resize(3); bases[1].resize(3); bases[2].resize(3); 00120 bases[0][0] = rcp( &openBasis_ , false ); 00121 bases[0][1] = rcp( &closedBasis_ , false ); 00122 bases[0][2] = rcp( &closedBasis_ , false ); 00123 bases[1][0] = rcp( &closedBasis_ , false ); 00124 bases[1][1] = rcp( &openBasis_ , false ); 00125 bases[1][2] = rcp( &closedBasis_ , false ); 00126 bases[2][0] = rcp( &closedBasis_ , false ); 00127 bases[2][1] = rcp( &closedBasis_ , false ); 00128 bases[2][2] = rcp( &openBasis_ , false ); 00129 this->setBases( bases ); 00130 00131 } 00132 00133 template<class Scalar, class ArrayScalar> 00134 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::initializeTags() { 00135 00136 // Basis-dependent intializations 00137 int tagSize = 4; // size of DoF tag 00138 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00139 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00140 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00141 00142 std::vector<int> tags( tagSize * this->getCardinality() ); 00143 00144 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags(); 00145 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags(); 00146 00147 std::map<int,std::map<int,int> > total_dof_per_entity; 00148 std::map<int,std::map<int,int> > current_dof_per_entity; 00149 00150 const int order = this->basisDegree_; 00151 00152 // vertices 00153 for (int i=0;i<4;i++) { 00154 total_dof_per_entity[0][i] = 0; 00155 current_dof_per_entity[0][i] = 0; 00156 } 00157 // edges 00158 for (int i=0;i<12;i++) { 00159 total_dof_per_entity[1][i] = 0; 00160 current_dof_per_entity[1][i] = 0; 00161 } 00162 // faces 00163 for (int i=0;i<6;i++) { 00164 total_dof_per_entity[2][i] = 0; 00165 current_dof_per_entity[2][i] = 0; 00166 } 00167 total_dof_per_entity[3][0] = 0; 00168 current_dof_per_entity[3][0] = 0; 00169 00170 // tally dof on each facet. none on vertex 00171 // edge dof 00172 for (int i=0;i<12;i++) { 00173 total_dof_per_entity[1][i] = order; 00174 } 00175 // face dof 00176 for (int i=0;i<6;i++) { 00177 total_dof_per_entity[2][i] = 2 * (order - 1) * order; 00178 } 00179 // internal dof 00180 total_dof_per_entity[3][0] = this->basisCardinality_ - 12 * order - 6 * total_dof_per_entity[2][0]; 00181 00182 int tagcur = 0; 00183 for (int k=0;k<closedBasis_.getCardinality();k++) { 00184 const int dimk = closedDofTags[k][0]; 00185 const int entk = closedDofTags[k][1]; 00186 for (int j=0;j<closedBasis_.getCardinality();j++) { 00187 const int dimj = closedDofTags[j][0]; 00188 const int entj = closedDofTags[j][1]; 00189 for (int i=0;i<openBasis_.getCardinality();i++) { 00190 const int dimi = openDofTags[i][0]; 00191 const int enti = openDofTags[i][1]; 00192 int dofdim; 00193 int dofent; 00194 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00195 tags[4*tagcur] = dofdim; 00196 tags[4*tagcur+1] = dofent; 00197 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00198 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00199 current_dof_per_entity[dofdim][dofent]++; 00200 tagcur++; 00201 } 00202 } 00203 } 00204 for (int k=0;k<closedBasis_.getCardinality();k++) { 00205 const int dimk = closedDofTags[k][0]; 00206 const int entk = closedDofTags[k][1]; 00207 for (int j=0;j<openBasis_.getCardinality();j++) { 00208 const int dimj = openDofTags[j][0]; 00209 const int entj = openDofTags[j][1]; 00210 for (int i=0;i<closedBasis_.getCardinality();i++) { 00211 const int dimi = closedDofTags[i][0]; 00212 const int enti = closedDofTags[i][1]; 00213 int dofdim; 00214 int dofent; 00215 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00216 tags[4*tagcur] = dofdim; 00217 tags[4*tagcur+1] = dofent; 00218 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00219 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00220 current_dof_per_entity[dofdim][dofent]++; 00221 tagcur++; 00222 } 00223 } 00224 } 00225 for (int k=0;k<openBasis_.getCardinality();k++) { 00226 const int dimk = openDofTags[k][0]; 00227 const int entk = openDofTags[k][1]; 00228 for (int j=0;j<closedBasis_.getCardinality();j++) { 00229 const int dimj = closedDofTags[j][0]; 00230 const int entj = closedDofTags[j][1]; 00231 for (int i=0;i<closedBasis_.getCardinality();i++) { 00232 const int dimi = closedDofTags[i][0]; 00233 const int enti = closedDofTags[i][1]; 00234 int dofdim; 00235 int dofent; 00236 ProductTopology::lineProduct3d(dimi,enti,dimj,entj,dimk,entk,dofdim,dofent); 00237 tags[4*tagcur] = dofdim; 00238 tags[4*tagcur+1] = dofent; 00239 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00240 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00241 current_dof_per_entity[dofdim][dofent]++; 00242 tagcur++; 00243 } 00244 } 00245 } 00246 00247 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00248 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00249 this -> ordinalToTag_, 00250 &(tags[0]), 00251 this -> basisCardinality_, 00252 tagSize, 00253 posScDim, 00254 posScOrd, 00255 posDfOrd); 00256 } 00257 00258 00259 template<class Scalar, class ArrayScalar> 00260 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00261 const ArrayScalar & inputPoints, 00262 const EOperator operatorType) const { 00263 00264 // Verify arguments 00265 #ifdef HAVE_INTREPID_DEBUG 00266 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues, 00267 inputPoints, 00268 operatorType, 00269 this -> getBaseCellTopology(), 00270 this -> getCardinality() ); 00271 #endif 00272 00273 // Number of evaluation points = dim 0 of inputPoints 00274 int dim0 = inputPoints.dimension(0); 00275 00276 // separate out points 00277 FieldContainer<Scalar> xPoints(dim0,1); 00278 FieldContainer<Scalar> yPoints(dim0,1); 00279 FieldContainer<Scalar> zPoints(dim0,1); 00280 00281 for (int i=0;i<dim0;i++) { 00282 xPoints(i,0) = inputPoints(i,0); 00283 yPoints(i,0) = inputPoints(i,1); 00284 zPoints(i,0) = inputPoints(i,2); 00285 } 00286 00287 switch (operatorType) { 00288 case OPERATOR_VALUE: 00289 { 00290 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00291 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00292 FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 ); 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 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00301 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00302 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00303 00304 // first we get the "horizontal basis functions that are tangent to the x-varying edges 00305 int bfcur = 0; 00306 for (int k=0;k<closedBasis_.getCardinality();k++) { 00307 for (int j=0;j<closedBasis_.getCardinality();j++) { 00308 for (int i=0;i<openBasis_.getCardinality();i++) { 00309 for (int l=0;l<dim0;l++) { 00310 outputValues(bfcur,l,0) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00311 outputValues(bfcur,l,1) = 0.0; 00312 outputValues(bfcur,l,2) = 0.0; 00313 } 00314 bfcur++; 00315 } 00316 } 00317 } 00318 00319 // second we get the basis functions in the direction of the y-varying edges 00320 for (int k=0;k<closedBasis_.getCardinality();k++) { 00321 for (int j=0;j<openBasis_.getCardinality();j++) { 00322 for (int i=0;i<closedBasis_.getCardinality();i++) { 00323 for (int l=0;l<dim0;l++) { 00324 outputValues(bfcur,l,0) = 0.0; 00325 outputValues(bfcur,l,1) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00326 outputValues(bfcur,l,2) = 0.0; 00327 } 00328 bfcur++; 00329 } 00330 } 00331 } 00332 00333 // third we get the basis functions in the direction of the y-varying edges 00334 for (int k=0;k<openBasis_.getCardinality();k++) { 00335 for (int j=0;j<closedBasis_.getCardinality();j++) { 00336 for (int i=0;i<closedBasis_.getCardinality();i++) { 00337 for (int l=0;l<dim0;l++) { 00338 outputValues(bfcur,l,0) = 0.0; 00339 outputValues(bfcur,l,1) = 0.0; 00340 outputValues(bfcur,l,2) = closedBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00341 } 00342 bfcur++; 00343 } 00344 } 00345 } 00346 } 00347 break; 00348 case OPERATOR_CURL: 00349 { 00350 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00351 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00352 FieldContainer<Scalar> closedBasisValsZPts( closedBasis_.getCardinality() , dim0 ); 00353 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00354 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00355 FieldContainer<Scalar> closedBasisDerivsZPts( closedBasis_.getCardinality() , dim0 , 1 ); 00356 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00357 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00358 FieldContainer<Scalar> openBasisValsZPts( openBasis_.getCardinality() , dim0 ); 00359 00360 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00361 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00362 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE ); 00363 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00364 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00365 closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 ); 00366 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00367 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00368 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE ); 00369 00370 int bfcur = 0; 00371 00372 // first we get the basis functions that are tangent to the x-varying edges 00373 for (int k=0;k<closedBasis_.getCardinality();k++) { 00374 for (int j=0;j<closedBasis_.getCardinality();j++) { 00375 for (int i=0;i<openBasis_.getCardinality();i++) { 00376 for (int l=0;l<dim0;l++) { 00377 outputValues(bfcur,l,0) = 0.0; 00378 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0); 00379 outputValues(bfcur,l,2) = -openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * closedBasisValsZPts(k,l); 00380 } 00381 bfcur++; 00382 } 00383 } 00384 } 00385 00386 // second we get the basis functions in the direction of the y-varying edges 00387 for (int k=0;k<closedBasis_.getCardinality();k++) { 00388 for (int j=0;j<openBasis_.getCardinality();j++) { 00389 for (int i=0;i<closedBasis_.getCardinality();i++) { 00390 for (int l=0;l<dim0;l++) { 00391 outputValues(bfcur,l,0) = -closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0); 00392 outputValues(bfcur,l,1) = 0.0; 00393 outputValues(bfcur,l,2) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l); 00394 } 00395 bfcur++; 00396 } 00397 } 00398 } 00399 00400 // third we get the basis functions in the direction of the y-varying edges 00401 for (int k=0;k<openBasis_.getCardinality();k++) { 00402 for (int j=0;j<closedBasis_.getCardinality();j++) { 00403 for (int i=0;i<closedBasis_.getCardinality();i++) { 00404 for (int l=0;l<dim0;l++) { 00405 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l); 00406 outputValues(bfcur,l,1) = -closedBasisDerivsXPts(i,l,0) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l); 00407 outputValues(bfcur,l,2) = 0.0; 00408 } 00409 bfcur++; 00410 } 00411 } 00412 } 00413 } 00414 break; 00415 case OPERATOR_DIV: 00416 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00417 ">>> ERROR (Basis_HCURL_HEX_In_FEM): DIV is invalid operator for HCURL Basis Functions"); 00418 break; 00419 00420 case OPERATOR_GRAD: 00421 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00422 ">>> ERROR (Basis_HCURL_HEX_In_FEM): GRAD is invalid operator for HCURL Basis Functions"); 00423 break; 00424 00425 case OPERATOR_D1: 00426 case OPERATOR_D2: 00427 case OPERATOR_D3: 00428 case OPERATOR_D4: 00429 case OPERATOR_D5: 00430 case OPERATOR_D6: 00431 case OPERATOR_D7: 00432 case OPERATOR_D8: 00433 case OPERATOR_D9: 00434 case OPERATOR_D10: 00435 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00436 (operatorType == OPERATOR_D2) || 00437 (operatorType == OPERATOR_D3) || 00438 (operatorType == OPERATOR_D4) || 00439 (operatorType == OPERATOR_D5) || 00440 (operatorType == OPERATOR_D6) || 00441 (operatorType == OPERATOR_D7) || 00442 (operatorType == OPERATOR_D8) || 00443 (operatorType == OPERATOR_D9) || 00444 (operatorType == OPERATOR_D10) ), 00445 std::invalid_argument, 00446 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type"); 00447 break; 00448 00449 default: 00450 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00451 (operatorType != OPERATOR_GRAD) && 00452 (operatorType != OPERATOR_CURL) && 00453 (operatorType != OPERATOR_CURL) && 00454 (operatorType != OPERATOR_D1) && 00455 (operatorType != OPERATOR_D2) && 00456 (operatorType != OPERATOR_D3) && 00457 (operatorType != OPERATOR_D4) && 00458 (operatorType != OPERATOR_D5) && 00459 (operatorType != OPERATOR_D6) && 00460 (operatorType != OPERATOR_D7) && 00461 (operatorType != OPERATOR_D8) && 00462 (operatorType != OPERATOR_D9) && 00463 (operatorType != OPERATOR_D10) ), 00464 std::invalid_argument, 00465 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type"); 00466 } 00467 } 00468 00469 00470 00471 template<class Scalar, class ArrayScalar> 00472 void Basis_HCURL_HEX_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00473 const ArrayScalar & inputPoints, 00474 const ArrayScalar & cellVertices, 00475 const EOperator operatorType) const { 00476 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00477 ">>> ERROR (Basis_HCURL_HEX_In_FEM): FEM Basis calling an FVD member function"); 00478 } 00479 00480 template<class Scalar, class ArrayScalar> 00481 void Basis_HCURL_HEX_In_FEM<Scalar,ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const 00482 { 00483 // x-component basis functions 00484 int cur = 0; 00485 00486 for (int k=0;k<closedPts_.dimension(0);k++) 00487 { 00488 for (int j=0;j<closedPts_.dimension(0);j++) 00489 { 00490 for (int i=0;i<openPts_.dimension(0);i++) 00491 { 00492 DofCoords(cur,0) = openPts_(i,0); 00493 DofCoords(cur,1) = closedPts_(j,0); 00494 DofCoords(cur,2) = closedPts_(k,0); 00495 cur++; 00496 } 00497 } 00498 } 00499 // y-component basis functions 00500 for (int k=0;k<closedPts_.dimension(0);k++) 00501 { 00502 for (int j=0;j<openPts_.dimension(0);j++) 00503 { 00504 for (int i=0;i<closedPts_.dimension(0);i++) 00505 { 00506 DofCoords(cur,0) = closedPts_(i,0); 00507 DofCoords(cur,1) = openPts_(j,0); 00508 DofCoords(cur,2) = closedPts_(k,0); 00509 cur++; 00510 } 00511 } 00512 } 00513 00514 // z-component basis functions 00515 for (int k=0;k<openPts_.dimension(0);k++) 00516 { 00517 for (int j=0;j<closedPts_.dimension(0);j++) 00518 { 00519 for (int i=0;i<closedPts_.dimension(0);i++) 00520 { 00521 DofCoords(cur,0) = closedPts_(i,0); 00522 DofCoords(cur,1) = closedPts_(j,0); 00523 DofCoords(cur,2) = openPts_(k,0); 00524 cur++; 00525 } 00526 } 00527 } 00528 00529 return; 00530 } 00531 00532 }// namespace Intrepid
1.7.6.1