|
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 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 00017 00018 // WITHOUT ANY WARRANTY; without even the implied warranty of 00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00020 // Lesser General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU Lesser General Public 00023 // License along with this library; if not, write to the Free Software 00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00025 // USA 00026 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00027 // Denis Ridzal (dridzal@sandia.gov) 00028 // Kara Peterson (kjpeter@sandia.gov). 00029 // 00030 // ************************************************************************ 00031 // @HEADER 00032 00038 namespace Intrepid { 00039 00040 template<class Scalar, class ArrayScalar> 00041 Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order , 00042 const ArrayScalar & ptsClosed , 00043 const ArrayScalar & ptsOpen): 00044 closedBasis_( order , ptsClosed ), 00045 openBasis_( order-1 , ptsOpen ) 00046 { 00047 this -> basisDegree_ = order; 00048 this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality(); 00049 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00050 this -> basisType_ = BASIS_FEM_FIAT; 00051 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00052 this -> basisTagsAreSet_ = false; 00053 } 00054 00055 template<class Scalar, class ArrayScalar> 00056 Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order , const EPointType &pointType ): 00057 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ), 00058 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ) 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 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::initializeTags() { 00070 00071 // Basis-dependent intializations 00072 int tagSize = 4; // size of DoF tag 00073 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00074 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00075 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00076 00077 std::vector<int> tags( tagSize * this->getCardinality() ); 00078 00079 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags(); 00080 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags(); 00081 00082 std::map<int,std::map<int,int> > total_dof_per_entity; 00083 std::map<int,std::map<int,int> > current_dof_per_entity; 00084 00085 for (int i=0;i<4;i++) { 00086 total_dof_per_entity[0][i] = 0; 00087 current_dof_per_entity[0][i] = 0; 00088 } 00089 for (int i=0;i<4;i++) { 00090 total_dof_per_entity[1][i] = 0; 00091 current_dof_per_entity[1][i] = 0; 00092 } 00093 total_dof_per_entity[2][0] = 0; 00094 current_dof_per_entity[2][0] = 0; 00095 00096 // tally dof on each facet. none on vertex 00097 for (int i=0;i<4;i++) { 00098 total_dof_per_entity[1][i] = openBasis_.getCardinality(); 00099 } 00100 00101 total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality(); 00102 00103 int tagcur = 0; 00104 // loop over the x-component basis functions, which are (psi(x)phi(y),0) 00105 // for psi in the closed basis and phi in the open 00106 for (int j=0;j<openBasis_.getCardinality();j++) { 00107 const int odim = openDofTags[j][0]; 00108 const int oent = openDofTags[j][1]; 00109 for (int i=0;i<closedBasis_.getCardinality();i++) { 00110 const int cdim = closedDofTags[i][0]; 00111 const int cent = closedDofTags[i][1]; 00112 int dofdim; 00113 int dofent; 00114 ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent); 00115 tags[4*tagcur] = dofdim; 00116 tags[4*tagcur+1] = dofent; 00117 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00118 current_dof_per_entity[dofdim][dofent]++; 00119 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00120 tagcur++; 00121 } 00122 } 00123 // now we have to do it for the y-component basis functions, which are 00124 // (0,phi(x)psi(y)) for psi in the closed basis and phi in the open 00125 for (int j=0;j<closedBasis_.getCardinality();j++) { 00126 const int cdim = closedDofTags[j][0]; 00127 const int cent = closedDofTags[j][1]; 00128 for (int i=0;i<openBasis_.getCardinality();i++) { 00129 const int odim = openDofTags[i][0]; 00130 const int oent = openDofTags[i][1]; 00131 int dofdim; 00132 int dofent; 00133 ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent); 00134 tags[4*tagcur] = dofdim; 00135 tags[4*tagcur+1] = dofent; 00136 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00137 current_dof_per_entity[dofdim][dofent]++; 00138 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00139 tagcur++; 00140 } 00141 } 00142 00143 // for (int i=0;i<this->getCardinality();i++) { 00144 // for (int j=0;j<4;j++) { 00145 // std::cout << tags[4*i+j] << " "; 00146 // } 00147 // std::cout << std::endl; 00148 // } 00149 00150 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00151 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00152 this -> ordinalToTag_, 00153 &(tags[0]), 00154 this -> basisCardinality_, 00155 tagSize, 00156 posScDim, 00157 posScOrd, 00158 posDfOrd); 00159 } 00160 00161 00162 template<class Scalar, class ArrayScalar> 00163 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00164 const ArrayScalar & inputPoints, 00165 const EOperator operatorType) const { 00166 00167 // Verify arguments 00168 #ifdef HAVE_INTREPID_DEBUG 00169 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues, 00170 inputPoints, 00171 operatorType, 00172 this -> getBaseCellTopology(), 00173 this -> getCardinality() ); 00174 #endif 00175 00176 // Number of evaluation points = dim 0 of inputPoints 00177 int dim0 = inputPoints.dimension(0); 00178 00179 // separate out points 00180 FieldContainer<Scalar> xPoints(dim0,1); 00181 FieldContainer<Scalar> yPoints(dim0,1); 00182 00183 for (int i=0;i<dim0;i++) { 00184 xPoints(i,0) = inputPoints(i,0); 00185 yPoints(i,0) = inputPoints(i,1); 00186 } 00187 00188 switch (operatorType) { 00189 case OPERATOR_VALUE: 00190 { 00191 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00192 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00193 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00194 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00195 00196 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00197 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00198 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00199 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00200 00201 int bfcur = 0; 00202 // x component bfs are (closed(x) open(y),0) 00203 for (int j=0;j<openBasis_.getCardinality();j++) { 00204 for (int i=0;i<closedBasis_.getCardinality();i++) { 00205 for (int l=0;l<dim0;l++) { 00206 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l); 00207 outputValues(bfcur,l,1) = 0.0; 00208 } 00209 bfcur++; 00210 } 00211 } 00212 00213 // y component bfs are (0,open(x) closed(y)) 00214 for (int j=0;j<closedBasis_.getCardinality();j++) { 00215 for (int i=0;i<openBasis_.getCardinality();i++) { 00216 for (int l=0;l<dim0;l++) { 00217 outputValues(bfcur,l,0) = 0.0; 00218 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l); 00219 } 00220 bfcur++; 00221 } 00222 } 00223 } 00224 break; 00225 case OPERATOR_DIV: 00226 { 00227 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00228 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00229 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00230 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00231 00232 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00233 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00234 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00235 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00236 00237 int bfcur = 0; 00238 00239 // x component basis functions first 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) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l); 00244 } 00245 bfcur++; 00246 } 00247 } 00248 00249 // now y component basis functions 00250 for (int j=0;j<closedBasis_.getCardinality();j++) { 00251 for (int i=0;i<openBasis_.getCardinality();i++) { 00252 for (int l=0;l<dim0;l++) { 00253 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0); 00254 } 00255 bfcur++; 00256 } 00257 } 00258 } 00259 break; 00260 case OPERATOR_CURL: 00261 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00262 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): CURL is invalid operator for HDIV Basis Functions"); 00263 break; 00264 00265 case OPERATOR_GRAD: 00266 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00267 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): GRAD is invalid operator for HDIV Basis Functions"); 00268 break; 00269 00270 case OPERATOR_D1: 00271 case OPERATOR_D2: 00272 case OPERATOR_D3: 00273 case OPERATOR_D4: 00274 case OPERATOR_D5: 00275 case OPERATOR_D6: 00276 case OPERATOR_D7: 00277 case OPERATOR_D8: 00278 case OPERATOR_D9: 00279 case OPERATOR_D10: 00280 TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00281 (operatorType == OPERATOR_D2) || 00282 (operatorType == OPERATOR_D3) || 00283 (operatorType == OPERATOR_D4) || 00284 (operatorType == OPERATOR_D5) || 00285 (operatorType == OPERATOR_D6) || 00286 (operatorType == OPERATOR_D7) || 00287 (operatorType == OPERATOR_D8) || 00288 (operatorType == OPERATOR_D9) || 00289 (operatorType == OPERATOR_D10) ), 00290 std::invalid_argument, 00291 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type"); 00292 break; 00293 00294 default: 00295 TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00296 (operatorType != OPERATOR_GRAD) && 00297 (operatorType != OPERATOR_CURL) && 00298 (operatorType != OPERATOR_DIV) && 00299 (operatorType != OPERATOR_D1) && 00300 (operatorType != OPERATOR_D2) && 00301 (operatorType != OPERATOR_D3) && 00302 (operatorType != OPERATOR_D4) && 00303 (operatorType != OPERATOR_D5) && 00304 (operatorType != OPERATOR_D6) && 00305 (operatorType != OPERATOR_D7) && 00306 (operatorType != OPERATOR_D8) && 00307 (operatorType != OPERATOR_D9) && 00308 (operatorType != OPERATOR_D10) ), 00309 std::invalid_argument, 00310 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type"); 00311 } 00312 } 00313 00314 00315 00316 template<class Scalar, class ArrayScalar> 00317 void Basis_HDIV_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00318 const ArrayScalar & inputPoints, 00319 const ArrayScalar & cellVertices, 00320 const EOperator operatorType) const { 00321 TEST_FOR_EXCEPTION( (true), std::logic_error, 00322 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): FEM Basis calling an FVD member function"); 00323 } 00324 00325 }// namespace Intrepid
1.7.4