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