|
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_HCURL_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_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_HCURL_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_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_HCURL_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-face tangent basis functions, which are (0,phi(x)psi(y)) 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 other basis functions, which are 00124 // (psi(x)phi(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 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00144 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00145 this -> ordinalToTag_, 00146 &(tags[0]), 00147 this -> basisCardinality_, 00148 tagSize, 00149 posScDim, 00150 posScOrd, 00151 posDfOrd); 00152 } 00153 00154 00155 template<class Scalar, class ArrayScalar> 00156 void Basis_HCURL_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00157 const ArrayScalar & inputPoints, 00158 const EOperator operatorType) const { 00159 00160 // Verify arguments 00161 #ifdef HAVE_INTREPID_DEBUG 00162 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues, 00163 inputPoints, 00164 operatorType, 00165 this -> getBaseCellTopology(), 00166 this -> getCardinality() ); 00167 #endif 00168 00169 // Number of evaluation points = dim 0 of inputPoints 00170 int dim0 = inputPoints.dimension(0); 00171 00172 // separate out points 00173 FieldContainer<Scalar> xPoints(dim0,1); 00174 FieldContainer<Scalar> yPoints(dim0,1); 00175 00176 for (int i=0;i<dim0;i++) { 00177 xPoints(i,0) = inputPoints(i,0); 00178 yPoints(i,0) = inputPoints(i,1); 00179 } 00180 00181 switch (operatorType) { 00182 case OPERATOR_VALUE: 00183 { 00184 FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 ); 00185 FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 ); 00186 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00187 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00188 00189 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE ); 00190 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE ); 00191 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00192 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00193 00194 int bfcur = 0; 00195 // x component bfs are (closed(x) open(y),0) 00196 for (int j=0;j<openBasis_.getCardinality();j++) { 00197 for (int i=0;i<closedBasis_.getCardinality();i++) { 00198 for (int l=0;l<dim0;l++) { 00199 outputValues(bfcur,l,0) = 0.0; 00200 outputValues(bfcur,l,1) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l); 00201 } 00202 bfcur++; 00203 } 00204 } 00205 00206 // y component bfs are (0,open(x) closed(y)) 00207 for (int j=0;j<closedBasis_.getCardinality();j++) { 00208 for (int i=0;i<openBasis_.getCardinality();i++) { 00209 for (int l=0;l<dim0;l++) { 00210 outputValues(bfcur,l,0) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l); 00211 outputValues(bfcur,l,1) = 0.0; 00212 } 00213 bfcur++; 00214 } 00215 } 00216 } 00217 break; 00218 case OPERATOR_CURL: 00219 { 00220 FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 ); 00221 FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 ); 00222 FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 ); 00223 FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 ); 00224 00225 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 ); 00226 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 ); 00227 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE ); 00228 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE ); 00229 00230 int bfcur = 0; 00231 00232 // x component basis functions first 00233 for (int j=0;j<openBasis_.getCardinality();j++) { 00234 for (int i=0;i<closedBasis_.getCardinality();i++) { 00235 for (int l=0;l<dim0;l++) { 00236 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l); 00237 } 00238 bfcur++; 00239 } 00240 } 00241 00242 // now y component basis functions 00243 for (int j=0;j<closedBasis_.getCardinality();j++) { 00244 for (int i=0;i<openBasis_.getCardinality();i++) { 00245 for (int l=0;l<dim0;l++) { 00246 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0); 00247 } 00248 bfcur++; 00249 } 00250 } 00251 } 00252 break; 00253 case OPERATOR_DIV: 00254 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00255 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): DIV is invalid operator for HCURL Basis Functions"); 00256 break; 00257 00258 case OPERATOR_GRAD: 00259 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument, 00260 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): GRAD is invalid operator for HCURL Basis Functions"); 00261 break; 00262 00263 case OPERATOR_D1: 00264 case OPERATOR_D2: 00265 case OPERATOR_D3: 00266 case OPERATOR_D4: 00267 case OPERATOR_D5: 00268 case OPERATOR_D6: 00269 case OPERATOR_D7: 00270 case OPERATOR_D8: 00271 case OPERATOR_D9: 00272 case OPERATOR_D10: 00273 TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) || 00274 (operatorType == OPERATOR_D2) || 00275 (operatorType == OPERATOR_D3) || 00276 (operatorType == OPERATOR_D4) || 00277 (operatorType == OPERATOR_D5) || 00278 (operatorType == OPERATOR_D6) || 00279 (operatorType == OPERATOR_D7) || 00280 (operatorType == OPERATOR_D8) || 00281 (operatorType == OPERATOR_D9) || 00282 (operatorType == OPERATOR_D10) ), 00283 std::invalid_argument, 00284 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type"); 00285 break; 00286 00287 default: 00288 TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) && 00289 (operatorType != OPERATOR_GRAD) && 00290 (operatorType != OPERATOR_CURL) && 00291 (operatorType != OPERATOR_CURL) && 00292 (operatorType != OPERATOR_D1) && 00293 (operatorType != OPERATOR_D2) && 00294 (operatorType != OPERATOR_D3) && 00295 (operatorType != OPERATOR_D4) && 00296 (operatorType != OPERATOR_D5) && 00297 (operatorType != OPERATOR_D6) && 00298 (operatorType != OPERATOR_D7) && 00299 (operatorType != OPERATOR_D8) && 00300 (operatorType != OPERATOR_D9) && 00301 (operatorType != OPERATOR_D10) ), 00302 std::invalid_argument, 00303 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type"); 00304 } 00305 } 00306 00307 00308 00309 template<class Scalar, class ArrayScalar> 00310 void Basis_HCURL_QUAD_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00311 const ArrayScalar & inputPoints, 00312 const ArrayScalar & cellVertices, 00313 const EOperator operatorType) const { 00314 TEST_FOR_EXCEPTION( (true), std::logic_error, 00315 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): FEM Basis calling an FVD member function"); 00316 } 00317 00318 }// namespace Intrepid
1.7.4