|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP 00002 #define INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP 00003 // @HEADER 00004 // ************************************************************************ 00005 // 00006 // Intrepid Package 00007 // Copyright (2007) Sandia Corporation 00008 // 00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00010 // license for use of this work by or on behalf of the U.S. Government. 00011 // 00012 // Redistribution and use in source and binary forms, with or without 00013 // modification, are permitted provided that the following conditions are 00014 // met: 00015 // 00016 // 1. Redistributions of source code must retain the above copyright 00017 // notice, this list of conditions and the following disclaimer. 00018 // 00019 // 2. Redistributions in binary form must reproduce the above copyright 00020 // notice, this list of conditions and the following disclaimer in the 00021 // documentation and/or other materials provided with the distribution. 00022 // 00023 // 3. Neither the name of the Corporation nor the names of the 00024 // contributors may be used to endorse or promote products derived from 00025 // this software without specific prior written permission. 00026 // 00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00038 // 00039 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00040 // Denis Ridzal (dridzal@sandia.gov), or 00041 // Kara Peterson (kjpeter@sandia.gov) 00042 // 00043 // ************************************************************************ 00044 // @HEADER 00045 00051 using Teuchos::Array; 00052 using Teuchos::RCP; 00053 00054 namespace Intrepid { 00055 template<class Scalar, class ArrayScalar> 00056 Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int orderx , const int ordery, 00057 const ArrayScalar &pts_x , 00058 const ArrayScalar &pts_y ): 00059 ptsx_( pts_x.dimension(0) , 1 ) , 00060 ptsy_( pts_y.dimension(0) , 1 ) 00061 { 00062 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1); 00063 bases[0].resize(2); 00064 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( orderx , pts_x ) ); 00065 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( ordery , pts_y ) ); 00066 this->setBases( bases ); 00067 00068 this->basisCardinality_ = (orderx+1)*(ordery+1); 00069 if (orderx > ordery) { 00070 this->basisDegree_ = orderx; 00071 } 00072 else { 00073 this->basisDegree_ = ordery; 00074 } 00075 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00076 this -> basisType_ = BASIS_FEM_FIAT; 00077 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00078 this -> basisTagsAreSet_ = false; 00079 00080 for (int i=0;i<pts_x.dimension(0);i++) 00081 { 00082 ptsx_(i,0) = pts_x(i,0); 00083 } 00084 00085 for (int i=0;i<pts_y.dimension(0);i++) 00086 { 00087 ptsy_(i,0) = pts_y(i,0); 00088 } 00089 00090 } 00091 00092 template<class Scalar, class ArrayScalar> 00093 Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int order, 00094 const EPointType &pointType ): 00095 ptsx_( order+1 , 1 ) , 00096 ptsy_( order+1 , 1 ) 00097 { 00098 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1); 00099 bases[0].resize(2); 00100 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) ); 00101 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) ); 00102 this->setBases( bases ); 00103 00104 this->basisCardinality_ = (order+1)*(order+1); 00105 this->basisDegree_ = order; 00106 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00107 this -> basisType_ = BASIS_FEM_FIAT; 00108 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00109 this -> basisTagsAreSet_ = false; 00110 00111 // fill up the pt arrays with calls to the lattice 00112 EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND; 00113 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ , 00114 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) , 00115 order , 00116 0 , 00117 pt ); 00118 00119 for (int i=0;i<order+1;i++) 00120 { 00121 ptsy_(i,0) = ptsx_(i,0); 00122 } 00123 } 00124 00125 template<class Scalar, class ArrayScalar> 00126 void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::initializeTags() 00127 { 00128 // Basis-dependent initializations 00129 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00130 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00131 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00132 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00133 00134 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00135 00136 std::vector<int> tags( tagSize * this->getCardinality() ); 00137 00138 // temporarily just put everything on the cell itself 00139 for (int i=0;i<this->getCardinality();i++) { 00140 tags[tagSize*i] = 2; 00141 tags[tagSize*i+1] = 0; 00142 tags[tagSize*i+2] = i; 00143 tags[tagSize*i+3] = this->getCardinality(); 00144 } 00145 00146 Basis<Scalar,ArrayScalar> &xBasis_ = *this->getBases()[0][0]; 00147 Basis<Scalar,ArrayScalar> &yBasis_ = *this->getBases()[0][1]; 00148 00149 // now let's try to do it "right" 00150 // let's get the x and y bases and their dof 00151 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags(); 00152 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags(); 00153 00154 std::map<int,std::map<int,int> > total_dof_per_entity; 00155 std::map<int,std::map<int,int> > current_dof_per_entity; 00156 00157 for (int i=0;i<4;i++) { 00158 total_dof_per_entity[0][i] = 0; 00159 current_dof_per_entity[0][i] = 0; 00160 } 00161 for (int i=0;i<4;i++) { 00162 total_dof_per_entity[1][i] = 0; 00163 current_dof_per_entity[1][i] = 0; 00164 } 00165 total_dof_per_entity[2][0] = 0; 00166 current_dof_per_entity[2][0] = 0; 00167 00168 // let's tally the total degrees of freedom on each entity 00169 for (int j=0;j<yBasis_.getCardinality();j++) { 00170 const int ydim = ydoftags[j][0]; 00171 const int yent = ydoftags[j][1]; 00172 for (int i=0;i<xBasis_.getCardinality();i++) { 00173 const int xdim = xdoftags[i][0]; 00174 const int xent = xdoftags[i][1]; 00175 int dofdim; 00176 int dofent; 00177 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent ); 00178 total_dof_per_entity[dofdim][dofent] += 1; 00179 } 00180 } 00181 00182 int tagcur = 0; 00183 for (int j=0;j<yBasis_.getCardinality();j++) { 00184 const int ydim = ydoftags[j][0]; 00185 const int yent = ydoftags[j][1]; 00186 for (int i=0;i<xBasis_.getCardinality();i++) { 00187 const int xdim = xdoftags[i][0]; 00188 const int xent = xdoftags[i][1]; 00189 int dofdim; 00190 int dofent; 00191 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent ); 00192 tags[4*tagcur] = dofdim; 00193 tags[4*tagcur+1] = dofent; 00194 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00195 current_dof_per_entity[dofdim][dofent]++; 00196 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00197 tagcur++; 00198 } 00199 } 00200 00201 setOrdinalTagData(this -> tagToOrdinal_, 00202 this -> ordinalToTag_, 00203 &(tags[0]), 00204 this -> basisCardinality_, 00205 tagSize, 00206 posScDim, 00207 posScOrd, 00208 posDfOrd); 00209 00210 } 00211 00212 template<class Scalar, class ArrayScalar> 00213 void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues , 00214 const ArrayScalar &inputPoints , 00215 const EOperator operatorType ) const 00216 { 00217 #ifdef HAVE_INTREPID_DEBUG 00218 getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00219 inputPoints, 00220 operatorType, 00221 this -> getBaseCellTopology(), 00222 this -> getCardinality() ); 00223 #endif 00224 00225 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1); 00226 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1); 00227 00228 const Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0]; 00229 const Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1]; 00230 00231 for (int i=0;i<inputPoints.dimension(0);i++) { 00232 xInputPoints(i,0) = inputPoints(i,0); 00233 yInputPoints(i,0) = inputPoints(i,1); 00234 } 00235 00236 switch (operatorType) { 00237 case OPERATOR_VALUE: 00238 { 00239 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00240 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00241 00242 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00243 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00244 00245 int bfcur = 0; 00246 for (int j=0;j<yBasis_.getCardinality();j++) { 00247 for (int i=0;i<xBasis_.getCardinality();i++) { 00248 for (int k=0;k<inputPoints.dimension(0);k++) { 00249 outputValues(bfcur,k) = xBasisValues(i,k) * yBasisValues(j,k); 00250 } 00251 bfcur++; 00252 } 00253 } 00254 } 00255 break; 00256 case OPERATOR_GRAD: 00257 case OPERATOR_D1: 00258 { 00259 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00260 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00261 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00262 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00263 00264 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00265 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00266 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00267 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00268 00269 // there are two multiindices: I need the (1,0) and (0,1) derivatives 00270 int bfcur = 0; 00271 00272 for (int j=0;j<yBasis_.getCardinality();j++) { 00273 for (int i=0;i<xBasis_.getCardinality();i++) { 00274 for (int k=0;k<inputPoints.dimension(0);k++) { 00275 outputValues(bfcur,k,0) = xBasisDerivs(i,k,0) * yBasisValues(j,k); 00276 outputValues(bfcur,k,1) = xBasisValues(i,k) * yBasisDerivs(j,k,0); 00277 } 00278 bfcur++; 00279 } 00280 } 00281 } 00282 break; 00283 case OPERATOR_D2: 00284 case OPERATOR_D3: 00285 case OPERATOR_D4: 00286 case OPERATOR_D5: 00287 case OPERATOR_D6: 00288 case OPERATOR_D7: 00289 case OPERATOR_D8: 00290 case OPERATOR_D9: 00291 case OPERATOR_D10: 00292 { 00293 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00294 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00295 00296 Teuchos::Array<int> partialMult; 00297 00298 for (int d=0;d<getDkCardinality(operatorType,2);d++) { 00299 getDkMultiplicities( partialMult , d , operatorType , 2 ); 00300 if (partialMult[0] == 0) { 00301 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00302 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE ); 00303 } 00304 else { 00305 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00306 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 ); 00307 xBasis_.getValues( xBasisValues , xInputPoints, xop ); 00308 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00309 } 00310 if (partialMult[1] == 0) { 00311 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00312 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE ); 00313 } 00314 else { 00315 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00316 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 ); 00317 yBasis_.getValues( yBasisValues , yInputPoints, yop ); 00318 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00319 } 00320 00321 00322 int bfcur = 0; 00323 for (int j=0;j<yBasis_.getCardinality();j++) { 00324 for (int i=0;i<xBasis_.getCardinality();i++) { 00325 for (int k=0;k<inputPoints.dimension(0);k++) { 00326 outputValues(bfcur,k,d) = xBasisValues(i,k) * yBasisValues(j,k); 00327 } 00328 bfcur++; 00329 } 00330 } 00331 } 00332 } 00333 break; 00334 case OPERATOR_CURL: 00335 { 00336 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00337 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00338 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00339 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00340 00341 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00342 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00343 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00344 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00345 00346 // there are two multiindices: I need the (1,0) and (0,1) derivatives 00347 int bfcur = 0; 00348 00349 for (int j=0;j<yBasis_.getCardinality();j++) { 00350 for (int i=0;i<xBasis_.getCardinality();i++) { 00351 for (int k=0;k<inputPoints.dimension(0);k++) { 00352 outputValues(bfcur,k,0) = xBasisValues(i,k) * yBasisDerivs(j,k,0); 00353 outputValues(bfcur,k,1) = -xBasisDerivs(i,k,0) * yBasisValues(j,k); 00354 } 00355 bfcur++; 00356 } 00357 } 00358 } 00359 break; 00360 default: 00361 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument, 00362 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented"); 00363 break; 00364 } 00365 } 00366 00367 template<class Scalar,class ArrayScalar> 00368 void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00369 const ArrayScalar & inputPoints, 00370 const ArrayScalar & cellVertices, 00371 const EOperator operatorType) const { 00372 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00373 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): FEM Basis calling an FVD member function"); 00374 } 00375 00376 template<class Scalar,class ArrayScalar> 00377 void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const 00378 { 00379 int cur = 0; 00380 for (int j=0;j<ptsy_.dimension(0);j++) 00381 { 00382 for (int i=0;i<ptsx_.dimension(0);i++) 00383 { 00384 dofCoords(cur,0) = ptsx_(i,0); 00385 dofCoords(cur,1) = ptsy_(j,0); 00386 cur++; 00387 } 00388 } 00389 } 00390 00391 00392 }// namespace Intrepid 00393 00394 #endif
1.7.6.1