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