|
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 00068 template<class Scalar, class ArrayScalar> 00069 Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_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 0 , 00099 POINTTYPE_EQUISPACED ); 00100 00101 } 00102 00103 } 00104 00105 template<class Scalar, class ArrayScalar> 00106 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::initializeTags() { 00107 00108 // Basis-dependent intializations 00109 int tagSize = 4; // size of DoF tag 00110 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00111 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00112 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00113 00114 std::vector<int> tags( tagSize * this->getCardinality() ); 00115 00116 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags(); 00117 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags(); 00118 00119 std::map<int,std::map<int,int> > total_dof_per_entity; 00120 std::map<int,std::map<int,int> > current_dof_per_entity; 00121 00122 for (int i=0;i<4;i++) { 00123 total_dof_per_entity[0][i] = 0; 00124 current_dof_per_entity[0][i] = 0; 00125 } 00126 for (int i=0;i<4;i++) { 00127 total_dof_per_entity[1][i] = 0; 00128 current_dof_per_entity[1][i] = 0; 00129 } 00130 total_dof_per_entity[2][0] = 0; 00131 current_dof_per_entity[2][0] = 0; 00132 00133 // tally dof on each facet. none on vertex 00134 for (int i=0;i<4;i++) { 00135 total_dof_per_entity[1][i] = openBasis_.getCardinality(); 00136 } 00137 00138 total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality(); 00139 00140 int tagcur = 0; 00141 // loop over the x-component basis functions, which are (psi(x)phi(y),0) 00142 // for psi in the closed basis and phi in the open 00143 for (int j=0;j<openBasis_.getCardinality();j++) { 00144 const int odim = openDofTags[j][0]; 00145 const int oent = openDofTags[j][1]; 00146 for (int i=0;i<closedBasis_.getCardinality();i++) { 00147 const int cdim = closedDofTags[i][0]; 00148 const int cent = closedDofTags[i][1]; 00149 int dofdim; 00150 int dofent; 00151 ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent); 00152 tags[4*tagcur] = dofdim; 00153 tags[4*tagcur+1] = dofent; 00154 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00155 current_dof_per_entity[dofdim][dofent]++; 00156 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00157 tagcur++; 00158 } 00159 } 00160 // now we have to do it for the y-component basis functions, which are 00161 // (0,phi(x)psi(y)) for psi in the closed basis and phi in the open 00162 for (int j=0;j<closedBasis_.getCardinality();j++) { 00163 const int cdim = closedDofTags[j][0]; 00164 const int cent = closedDofTags[j][1]; 00165 for (int i=0;i<openBasis_.getCardinality();i++) { 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 // for (int i=0;i<this->getCardinality();i++) { 00181 // for (int j=0;j<4;j++) { 00182 // std::cout << tags[4*i+j] << " "; 00183 // } 00184 // std::cout << std::endl; 00185 // } 00186 00187 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00188 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00189 this -> ordinalToTag_, 00190 &(tags[0]), 00191 this -> basisCardinality_, 00192 tagSize, 00193 posScDim, 00194 posScOrd, 00195 posDfOrd); 00196 } 00197 00198 00199 template<class Scalar, class ArrayScalar> 00200 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00201 const ArrayScalar & inputPoints, 00202 const EOperator operatorType) const { 00203 00204 // Verify arguments 00205 #ifdef HAVE_INTREPID_DEBUG 00206 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues, 00207 inputPoints, 00208 operatorType, 00209 this -> getBaseCellTopology(), 00210 this -> getCardinality() ); 00211 #endif 00212 00213 // Number of evaluation points = dim 0 of inputPoints 00214 int dim0 = inputPoints.dimension(0); 00215 00216 // separate out points 00217 FieldContainer<Scalar> xPoints(dim0,1); 00218 FieldContainer<Scalar> yPoints(dim0,1); 00219 00220 for (int i=0;i<dim0;i++) { 00221 xPoints(i,0) = inputPoints(i,0); 00222 yPoints(i,0) = inputPoints(i,1); 00223 } 00224 00225 switch (operatorType) { 00226 case OPERATOR_VALUE: 00227 { 00228 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00229 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00230 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00231 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00232 00233 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00234 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00235 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00236 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00237 00238 int bfcur = 0; 00239 // x component bfs are (closed(x) open(y),0) 00240 for (int j=0;j<openBasis_.getCardinality();j++) { 00241 for (int i=0;i<closedBasis_.getCardinality();i++) { 00242 for (int l=0;l<dim0;l++) { 00243 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l); 00244 outputValues(bfcur,l,1) = 0.0; 00245 } 00246 bfcur++; 00247 } 00248 } 00249 00250 // y component bfs are (0,open(x) closed(y)) 00251 for (int j=0;j<closedBasis_.getCardinality();j++) { 00252 for (int i=0;i<openBasis_.getCardinality();i++) { 00253 for (int l=0;l<dim0;l++) { 00254 outputValues(bfcur,l,0) = 0.0; 00255 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l); 00256 } 00257 bfcur++; 00258 } 00259 } 00260 } 00261 break; 00262 case OPERATOR_DIV: 00263 { 00264 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00265 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00266 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00267 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00268 00269 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00270 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00271 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00272 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00273 00274 int bfcur = 0; 00275 00276 // x component basis functions first 00277 for (int j=0;j<openBasis_.getCardinality();j++) { 00278 for (int i=0;i<closedBasis_.getCardinality();i++) { 00279 for (int l=0;l<dim0;l++) { 00280 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l); 00281 } 00282 bfcur++; 00283 } 00284 } 00285 00286 // now y component basis functions 00287 for (int j=0;j<closedBasis_.getCardinality();j++) { 00288 for (int i=0;i<openBasis_.getCardinality();i++) { 00289 for (int l=0;l<dim0;l++) { 00290 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0); 00291 } 00292 bfcur++; 00293 } 00294 } 00295 } 00296 break; 00297 case OPERATOR_CURL: 00298 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00299 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): CURL is invalid operator for HDIV Basis Functions"); 00300 break; 00301 00302 case OPERATOR_GRAD: 00303 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00304 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): GRAD is invalid operator for HDIV Basis Functions"); 00305 break; 00306 00307 case OPERATOR_D1: 00308 case OPERATOR_D2: 00309 case OPERATOR_D3: 00310 case OPERATOR_D4: 00311 case OPERATOR_D5: 00312 case OPERATOR_D6: 00313 case OPERATOR_D7: 00314 case OPERATOR_D8: 00315 case OPERATOR_D9: 00316 case OPERATOR_D10: 00317 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00318 (operatorType == OPERATOR_D2) || 00319 (operatorType == OPERATOR_D3) || 00320 (operatorType == OPERATOR_D4) || 00321 (operatorType == OPERATOR_D5) || 00322 (operatorType == OPERATOR_D6) || 00323 (operatorType == OPERATOR_D7) || 00324 (operatorType == OPERATOR_D8) || 00325 (operatorType == OPERATOR_D9) || 00326 (operatorType == OPERATOR_D10) ), 00327 std::invalid_argument, 00328 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type"); 00329 break; 00330 00331 default: 00332 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00333 (operatorType != OPERATOR_GRAD) && 00334 (operatorType != OPERATOR_CURL) && 00335 (operatorType != OPERATOR_DIV) && 00336 (operatorType != OPERATOR_D1) && 00337 (operatorType != OPERATOR_D2) && 00338 (operatorType != OPERATOR_D3) && 00339 (operatorType != OPERATOR_D4) && 00340 (operatorType != OPERATOR_D5) && 00341 (operatorType != OPERATOR_D6) && 00342 (operatorType != OPERATOR_D7) && 00343 (operatorType != OPERATOR_D8) && 00344 (operatorType != OPERATOR_D9) && 00345 (operatorType != OPERATOR_D10) ), 00346 std::invalid_argument, 00347 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type"); 00348 } 00349 } 00350 00351 00352 00353 template<class Scalar, class ArrayScalar> 00354 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00355 const ArrayScalar & inputPoints, 00356 const ArrayScalar & cellVertices, 00357 const EOperator operatorType) const { 00358 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00359 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): FEM Basis calling an FVD member function"); 00360 } 00361 00362 template<class Scalar, class ArrayScalar> 00363 void Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const 00364 { 00365 // x-component basis functions 00366 int cur = 0; 00367 00368 for (int j=0;j<openPts_.dimension(0);j++) 00369 { 00370 for (int i=0;i<closedPts_.dimension(0);i++) 00371 { 00372 DofCoords(cur,0) = closedPts_(i,0); 00373 DofCoords(cur,1) = openPts_(j,0); 00374 cur++; 00375 } 00376 } 00377 00378 // y-component basis functions 00379 for (int j=0;j<closedPts_.dimension(0);j++) 00380 { 00381 for (int i=0;i<openPts_.dimension(0);i++) 00382 { 00383 DofCoords(cur,0) = openPts_(i,0); 00384 DofCoords(cur,1) = closedPts_(j,0); 00385 cur++; 00386 } 00387 } 00388 00389 return; 00390 } 00391 00392 00393 }// namespace Intrepid
1.7.6.1