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