|
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 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov) or 00026 // Robert Kirby (robert.c.kirby@ttu.edu) 00027 // 00028 // ************************************************************************ 00029 // @HEADER 00030 00036 namespace Intrepid { 00037 template<class Scalar, class ArrayScalar> 00038 Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int orderx , const int ordery, 00039 const ArrayScalar &pts_x , 00040 const ArrayScalar &pts_y ): 00041 xBasis_( orderx , pts_x ), 00042 yBasis_( ordery , pts_y ) 00043 { 00044 this->basisCardinality_ = (orderx+1)*(ordery+1); 00045 if (orderx > ordery) { 00046 this->basisDegree_ = orderx; 00047 } 00048 else { 00049 this->basisDegree_ = ordery; 00050 } 00051 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00052 this -> basisType_ = BASIS_FEM_FIAT; 00053 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00054 this -> basisTagsAreSet_ = false; 00055 } 00056 00057 template<class Scalar, class ArrayScalar> 00058 Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int order, 00059 const EPointType &pointType ): 00060 xBasis_( order , pointType ), 00061 yBasis_( order , pointType ) 00062 { 00063 this->basisCardinality_ = (order+1)*(order+1); 00064 this->basisDegree_ = order; 00065 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00066 this -> basisType_ = BASIS_FEM_FIAT; 00067 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00068 this -> basisTagsAreSet_ = false; 00069 } 00070 00071 template<class Scalar, class ArrayScalar> 00072 void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::initializeTags() 00073 { 00074 // Basis-dependent initializations 00075 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00076 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00077 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00078 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00079 00080 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00081 00082 std::vector<int> tags( tagSize * this->getCardinality() ); 00083 00084 // temporarily just put everything on the cell itself 00085 for (int i=0;i<this->getCardinality();i++) { 00086 tags[tagSize*i] = 2; 00087 tags[tagSize*i+1] = 0; 00088 tags[tagSize*i+2] = i; 00089 tags[tagSize*i+3] = this->getCardinality(); 00090 } 00091 00092 // now let's try to do it "right" 00093 // let's get the x and y bases and their dof 00094 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags(); 00095 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags(); 00096 00097 std::map<int,std::map<int,int> > total_dof_per_entity; 00098 std::map<int,std::map<int,int> > current_dof_per_entity; 00099 00100 for (int i=0;i<4;i++) { 00101 total_dof_per_entity[0][i] = 0; 00102 current_dof_per_entity[0][i] = 0; 00103 } 00104 for (int i=0;i<4;i++) { 00105 total_dof_per_entity[1][i] = 0; 00106 current_dof_per_entity[1][i] = 0; 00107 } 00108 total_dof_per_entity[2][0] = 0; 00109 current_dof_per_entity[2][0] = 0; 00110 00111 // let's tally the total degrees of freedom on each entity 00112 for (int j=0;j<yBasis_.getCardinality();j++) { 00113 const int ydim = ydoftags[j][0]; 00114 const int yent = ydoftags[j][1]; 00115 for (int i=0;i<xBasis_.getCardinality();i++) { 00116 const int xdim = xdoftags[i][0]; 00117 const int xent = xdoftags[i][1]; 00118 int dofdim; 00119 int dofent; 00120 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent ); 00121 total_dof_per_entity[dofdim][dofent] += 1; 00122 } 00123 } 00124 00125 int tagcur = 0; 00126 for (int j=0;j<yBasis_.getCardinality();j++) { 00127 const int ydim = ydoftags[j][0]; 00128 const int yent = ydoftags[j][1]; 00129 for (int i=0;i<xBasis_.getCardinality();i++) { 00130 const int xdim = xdoftags[i][0]; 00131 const int xent = xdoftags[i][1]; 00132 int dofdim; 00133 int dofent; 00134 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent ); 00135 tags[4*tagcur] = dofdim; 00136 tags[4*tagcur+1] = dofent; 00137 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00138 current_dof_per_entity[dofdim][dofent]++; 00139 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00140 tagcur++; 00141 } 00142 } 00143 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_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues , 00157 const ArrayScalar &inputPoints , 00158 const EOperator operatorType ) const 00159 { 00160 #ifdef HAVE_INTREPID_DEBUG 00161 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00162 inputPoints, 00163 operatorType, 00164 this -> getBaseCellTopology(), 00165 this -> getCardinality() ); 00166 #endif 00167 00168 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1); 00169 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1); 00170 00171 for (int i=0;i<inputPoints.dimension(0);i++) { 00172 xInputPoints(i,0) = inputPoints(i,0); 00173 yInputPoints(i,0) = inputPoints(i,1); 00174 } 00175 00176 switch (operatorType) { 00177 case OPERATOR_VALUE: 00178 { 00179 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00180 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00181 00182 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00183 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00184 00185 int bfcur = 0; 00186 for (int j=0;j<yBasis_.getCardinality();j++) { 00187 for (int i=0;i<xBasis_.getCardinality();i++) { 00188 for (int k=0;k<inputPoints.dimension(0);k++) { 00189 outputValues(bfcur,k) = xBasisValues(i,k) * yBasisValues(j,k); 00190 } 00191 bfcur++; 00192 } 00193 } 00194 } 00195 break; 00196 case OPERATOR_GRAD: 00197 case OPERATOR_D1: 00198 { 00199 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00200 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00201 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00202 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00203 00204 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00205 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00206 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00207 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00208 00209 // there are two multiindices: I need the (1,0) and (0,1) derivatives 00210 int bfcur = 0; 00211 00212 for (int j=0;j<yBasis_.getCardinality();j++) { 00213 for (int i=0;i<xBasis_.getCardinality();i++) { 00214 for (int k=0;k<inputPoints.dimension(0);k++) { 00215 outputValues(bfcur,k,0) = xBasisDerivs(i,k,0) * yBasisValues(j,k); 00216 outputValues(bfcur,k,1) = xBasisValues(i,k) * yBasisDerivs(j,k,0); 00217 } 00218 bfcur++; 00219 } 00220 } 00221 } 00222 break; 00223 case OPERATOR_D2: 00224 case OPERATOR_D3: 00225 case OPERATOR_D4: 00226 case OPERATOR_D5: 00227 case OPERATOR_D6: 00228 case OPERATOR_D7: 00229 case OPERATOR_D8: 00230 case OPERATOR_D9: 00231 case OPERATOR_D10: 00232 { 00233 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00234 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00235 00236 Teuchos::Array<int> partialMult; 00237 00238 for (int d=0;d<getDkCardinality(operatorType,2);d++) { 00239 getDkMultiplicities( partialMult , d , operatorType , 2 ); 00240 if (partialMult[0] == 0) { 00241 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00242 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE ); 00243 } 00244 else { 00245 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00246 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 ); 00247 xBasis_.getValues( xBasisValues , xInputPoints, xop ); 00248 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00249 } 00250 if (partialMult[1] == 0) { 00251 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00252 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE ); 00253 } 00254 else { 00255 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00256 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 ); 00257 yBasis_.getValues( yBasisValues , yInputPoints, yop ); 00258 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00259 } 00260 00261 00262 int bfcur = 0; 00263 for (int j=0;j<yBasis_.getCardinality();j++) { 00264 for (int i=0;i<xBasis_.getCardinality();i++) { 00265 for (int k=0;k<inputPoints.dimension(0);k++) { 00266 outputValues(bfcur,k,d) = xBasisValues(i,k) * yBasisValues(j,k); 00267 } 00268 bfcur++; 00269 } 00270 } 00271 } 00272 } 00273 break; 00274 case OPERATOR_CURL: 00275 { 00276 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00277 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00278 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00279 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00280 00281 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00282 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00283 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00284 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00285 00286 // there are two multiindices: I need the (1,0) and (0,1) derivatives 00287 int bfcur = 0; 00288 00289 for (int j=0;j<yBasis_.getCardinality();j++) { 00290 for (int i=0;i<xBasis_.getCardinality();i++) { 00291 for (int k=0;k<inputPoints.dimension(0);k++) { 00292 outputValues(bfcur,k,0) = xBasisValues(i,k) * yBasisDerivs(j,k,0); 00293 outputValues(bfcur,k,1) = -xBasisDerivs(i,k,0) * yBasisValues(j,k); 00294 } 00295 bfcur++; 00296 } 00297 } 00298 } 00299 break; 00300 default: 00301 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00302 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented"); 00303 break; 00304 } 00305 } 00306 00307 template<class Scalar,class ArrayScalar> 00308 void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00309 const ArrayScalar & inputPoints, 00310 const ArrayScalar & cellVertices, 00311 const EOperator operatorType) const { 00312 TEST_FOR_EXCEPTION( (true), std::logic_error, 00313 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): FEM Basis calling an FVD member function"); 00314 } 00315 00316 00317 }// namespace Intrepid 00318
1.7.4