|
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_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int orderx , 00039 const int ordery , 00040 const int orderz , 00041 const ArrayScalar &pts_x , 00042 const ArrayScalar &pts_y , 00043 const ArrayScalar &pts_z ): 00044 xBasis_( orderx , pts_x ), 00045 yBasis_( ordery , pts_y ), 00046 zBasis_( orderz , pts_z ) 00047 { 00048 this->basisCardinality_ = (orderx+1)*(ordery+1)*(orderz+1); 00049 if (orderx >= ordery && orderx >= orderz ) { 00050 this->basisDegree_ = orderx; 00051 } 00052 else if (ordery >= orderx && ordery >= orderz) { 00053 this->basisDegree_ = ordery; 00054 } 00055 else { 00056 this->basisDegree_ = orderz; 00057 } 00058 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00059 this -> basisType_ = BASIS_FEM_FIAT; 00060 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00061 this -> basisTagsAreSet_ = false; 00062 } 00063 00064 template<class Scalar, class ArrayScalar> 00065 Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int order , 00066 const EPointType & pointType ): 00067 xBasis_( order , pointType ), 00068 yBasis_( order , pointType ), 00069 zBasis_( order , pointType ) 00070 { 00071 this->basisCardinality_ = (order+1)*(order+1)*(order+1); 00072 this->basisDegree_ = order; 00073 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00074 this -> basisType_ = BASIS_FEM_FIAT; 00075 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00076 this -> basisTagsAreSet_ = false; 00077 } 00078 00079 template<class Scalar, class ArrayScalar> 00080 void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::initializeTags() 00081 { 00082 // Basis-dependent initializations 00083 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00084 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00085 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00086 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00087 00088 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00089 00090 std::vector<int> tags( tagSize * this->getCardinality() ); 00091 00092 // temporarily just put everything on the cell itself 00093 for (int i=0;i<this->getCardinality();i++) { 00094 tags[tagSize*i] = 3; 00095 tags[tagSize*i+1] = 0; 00096 tags[tagSize*i+2] = i; 00097 tags[tagSize*i+3] = this->getCardinality(); 00098 } 00099 00100 // now let's try to do it "right" 00101 // let's get the x, y, z bases and their dof 00102 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags(); 00103 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags(); 00104 const std::vector<std::vector<int> >& zdoftags = zBasis_.getAllDofTags(); 00105 00106 std::map<int,std::map<int,int> > total_dof_per_entity; 00107 std::map<int,std::map<int,int> > current_dof_per_entity; 00108 00109 // vertex dof 00110 for (int i=0;i<8;i++) { 00111 total_dof_per_entity[0][i] = 0; 00112 current_dof_per_entity[0][i] = 0; 00113 } 00114 // edge dof (12 edges) 00115 for (int i=0;i<12;i++) { 00116 total_dof_per_entity[1][i] = 0; 00117 current_dof_per_entity[1][i] = 0; 00118 } 00119 // face dof (6 faces) 00120 for (int i=0;i<6;i++) { 00121 total_dof_per_entity[2][i] = 0; 00122 current_dof_per_entity[2][i] = 0; 00123 } 00124 // internal dof 00125 total_dof_per_entity[3][0] = 0; 00126 //total_dof_per_entity[3][1] = 0; 00127 00128 // let's tally the total degrees of freedom on each entity 00129 for (int k=0;k<zBasis_.getCardinality();k++) { 00130 const int zdim = zdoftags[k][0]; 00131 const int zent = zdoftags[k][1]; 00132 for (int j=0;j<yBasis_.getCardinality();j++) { 00133 const int ydim = ydoftags[j][0]; 00134 const int yent = ydoftags[j][1]; 00135 for (int i=0;i<xBasis_.getCardinality();i++) { 00136 const int xdim = xdoftags[i][0]; 00137 const int xent = xdoftags[i][1]; 00138 int dofdim; 00139 int dofent; 00140 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent ); 00141 std::cout << i << " " << j << " " << k << "\t" << dofdim << " " << dofent << std::endl; 00142 total_dof_per_entity[dofdim][dofent] += 1; 00143 00144 } 00145 } 00146 } 00147 00148 int tagcur = 0; 00149 for (int k=0;k<zBasis_.getCardinality();k++) { 00150 const int zdim = zdoftags[k][0]; 00151 const int zent = zdoftags[k][1]; 00152 for (int j=0;j<yBasis_.getCardinality();j++) { 00153 const int ydim = ydoftags[j][0]; 00154 const int yent = ydoftags[j][1]; 00155 for (int i=0;i<xBasis_.getCardinality();i++) { 00156 const int xdim = xdoftags[i][0]; 00157 const int xent = xdoftags[i][1]; 00158 int dofdim; 00159 int dofent; 00160 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent ); 00161 tags[4*tagcur] = dofdim; 00162 tags[4*tagcur+1] = dofent; 00163 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00164 current_dof_per_entity[dofdim][dofent]++; 00165 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00166 tagcur++; 00167 } 00168 } 00169 } 00170 00171 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00172 this -> ordinalToTag_, 00173 &(tags[0]), 00174 this -> basisCardinality_, 00175 tagSize, 00176 posScDim, 00177 posScOrd, 00178 posDfOrd); 00179 00180 } 00181 00182 template<class Scalar, class ArrayScalar> 00183 void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues , 00184 const ArrayScalar &inputPoints , 00185 const EOperator operatorType ) const 00186 { 00187 #ifdef HAVE_INTREPID_DEBUG 00188 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00189 inputPoints, 00190 operatorType, 00191 this -> getBaseCellTopology(), 00192 this -> getCardinality() ); 00193 #endif 00194 00195 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1); 00196 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1); 00197 FieldContainer<Scalar> zInputPoints(inputPoints.dimension(0),1); 00198 00199 for (int i=0;i<inputPoints.dimension(0);i++) { 00200 xInputPoints(i,0) = inputPoints(i,0); 00201 yInputPoints(i,0) = inputPoints(i,1); 00202 zInputPoints(i,0) = inputPoints(i,2); 00203 } 00204 00205 switch (operatorType) { 00206 case OPERATOR_VALUE: 00207 { 00208 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00209 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00210 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00211 00212 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00213 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00214 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE); 00215 00216 int bfcur = 0; 00217 for (int k=0;k<zBasis_.getCardinality();k++) { 00218 for (int j=0;j<yBasis_.getCardinality();j++) { 00219 for (int i=0;i<xBasis_.getCardinality();i++) { 00220 for (int l=0;l<inputPoints.dimension(0);l++) { 00221 outputValues(bfcur,l) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l); 00222 } 00223 bfcur++; 00224 } 00225 } 00226 } 00227 } 00228 break; 00229 case OPERATOR_D1: 00230 case OPERATOR_GRAD: 00231 { 00232 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00233 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00234 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00235 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00236 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00237 FieldContainer<Scalar> zBasisDerivs(zBasis_.getCardinality(),zInputPoints.dimension(0),1); 00238 00239 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00240 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00241 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE); 00242 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00243 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00244 zBasis_.getValues(zBasisDerivs,zInputPoints,OPERATOR_D1); 00245 00246 int bfcur = 0; 00247 for (int k=0;k<zBasis_.getCardinality();k++) { 00248 for (int j=0;j<yBasis_.getCardinality();j++) { 00249 for (int i=0;i<xBasis_.getCardinality();i++) { 00250 for (int l=0;l<inputPoints.dimension(0);l++) { 00251 outputValues(bfcur,l,0) = xBasisDerivs(i,l,0) * yBasisValues(j,l) * zBasisValues(k,l); 00252 outputValues(bfcur,l,1) = xBasisValues(i,l) * yBasisDerivs(j,l,0) * zBasisValues(k,l); 00253 outputValues(bfcur,l,2) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisDerivs(k,l,0); 00254 } 00255 bfcur++; 00256 } 00257 } 00258 } 00259 } 00260 break; 00261 case OPERATOR_D2: 00262 case OPERATOR_D3: 00263 case OPERATOR_D4: 00264 case OPERATOR_D5: 00265 case OPERATOR_D6: 00266 case OPERATOR_D7: 00267 case OPERATOR_D8: 00268 case OPERATOR_D9: 00269 case OPERATOR_D10: 00270 { 00271 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00272 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00273 FieldContainer<Scalar> zBasisValues(yBasis_.getCardinality(),zInputPoints.dimension(0)); 00274 00275 Teuchos::Array<int> partialMult; 00276 00277 for (int d=0;d<getDkCardinality(operatorType,3);d++) { 00278 getDkMultiplicities( partialMult , d , operatorType , 3 ); 00279 if (partialMult[0] == 0) { 00280 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00281 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE ); 00282 } 00283 else { 00284 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00285 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 ); 00286 xBasis_.getValues( xBasisValues , xInputPoints, xop ); 00287 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00288 } 00289 if (partialMult[1] == 0) { 00290 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00291 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE ); 00292 } 00293 else { 00294 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00295 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 ); 00296 yBasis_.getValues( yBasisValues , yInputPoints, yop ); 00297 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00298 } 00299 if (partialMult[2] == 0) { 00300 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00301 zBasis_.getValues( zBasisValues , zInputPoints, OPERATOR_VALUE ); 00302 } 00303 else { 00304 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0),1); 00305 EOperator zop = (EOperator) ( (int) OPERATOR_D1 + partialMult[2] - 1 ); 00306 zBasis_.getValues( zBasisValues , zInputPoints, zop ); 00307 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00308 } 00309 00310 00311 int bfcur = 0; 00312 for (int k=0;k<zBasis_.getCardinality();k++) { 00313 for (int j=0;j<yBasis_.getCardinality();j++) { 00314 for (int i=0;i<xBasis_.getCardinality();i++) { 00315 for (int l=0;l<inputPoints.dimension(0);l++) { 00316 outputValues(bfcur,l,d) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l); 00317 } 00318 bfcur++; 00319 } 00320 } 00321 } 00322 } 00323 } 00324 break; 00325 default: 00326 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00327 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented"); 00328 break; 00329 } 00330 } 00331 00332 template<class Scalar,class ArrayScalar> 00333 void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00334 const ArrayScalar & inputPoints, 00335 const ArrayScalar & cellVertices, 00336 const EOperator operatorType) const { 00337 TEST_FOR_EXCEPTION( (true), std::logic_error, 00338 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): FEM Basis calling an FVD member function"); 00339 } 00340 00341 00342 }// namespace Intrepid 00343
1.7.4