|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_HEX_CN_FEMDEF_HPP 00002 #define INTREPID_HGRAD_HEX_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 namespace Intrepid { 00052 template<class Scalar, class ArrayScalar> 00053 Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int orderx , 00054 const int ordery , 00055 const int orderz , 00056 const ArrayScalar &pts_x , 00057 const ArrayScalar &pts_y , 00058 const ArrayScalar &pts_z ): 00059 ptsx_( pts_x.dimension(0) , 1 ), 00060 ptsy_( pts_y.dimension(0) , 1 ), 00061 ptsz_( pts_z.dimension(0) , 1 ) 00062 { 00063 for (int i=0;i<pts_x.dimension(0);i++) 00064 { 00065 ptsx_(i,0) = pts_x(i,0); 00066 } 00067 for (int i=0;i<pts_y.dimension(0);i++) 00068 { 00069 ptsy_(i,0) = pts_y(i,0); 00070 } 00071 for (int i=0;i<pts_z.dimension(0);i++) 00072 { 00073 ptsz_(i,0) = pts_z(i,0); 00074 } 00075 00076 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1); 00077 bases[0].resize(3); 00078 00079 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderx , pts_x ) ); 00080 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( ordery , pts_y ) ); 00081 bases[0][2] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderz , pts_z ) ); 00082 00083 this->setBases( bases ); 00084 00085 this->basisCardinality_ = (orderx+1)*(ordery+1)*(orderz+1); 00086 if (orderx >= ordery && orderx >= orderz ) { 00087 this->basisDegree_ = orderx; 00088 } 00089 else if (ordery >= orderx && ordery >= orderz) { 00090 this->basisDegree_ = ordery; 00091 } 00092 else { 00093 this->basisDegree_ = orderz; 00094 } 00095 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00096 this -> basisType_ = BASIS_FEM_FIAT; 00097 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00098 this -> basisTagsAreSet_ = false; 00099 } 00100 00101 template<class Scalar, class ArrayScalar> 00102 Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_Cn_FEM( const int order , 00103 const EPointType & pointType ): 00104 ptsx_( order+1 , 1 ), 00105 ptsy_( order+1 , 1 ), 00106 ptsz_( order+1 , 1 ) 00107 { 00108 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1); 00109 bases[0].resize(3); 00110 00111 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( order , pointType ) ); 00112 // basis is same in each direction, so I only need to instantiate it once! 00113 bases[0][1] = bases[0][0]; 00114 bases[0][2] = bases[0][0]; 00115 00116 this->setBases( bases ); 00117 00118 this->basisCardinality_ = (order+1)*(order+1)*(order+1); 00119 this->basisDegree_ = order; 00120 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00121 this -> basisType_ = BASIS_FEM_FIAT; 00122 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00123 this -> basisTagsAreSet_ = false; 00124 00125 // get points 00126 EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND; 00127 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ , 00128 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) , 00129 order , 00130 0 , 00131 pt ); 00132 for (int i=0;i<order+1;i++) 00133 { 00134 ptsy_(i,0) = ptsx_(i,0); 00135 ptsz_(i,0) = ptsx_(i,0); 00136 } 00137 00138 } 00139 00140 template<class Scalar, class ArrayScalar> 00141 void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::initializeTags() 00142 { 00143 // Basis-dependent initializations 00144 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00145 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00146 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00147 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00148 00149 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00150 00151 std::vector<int> tags( tagSize * this->getCardinality() ); 00152 00153 // temporarily just put everything on the cell itself 00154 for (int i=0;i<this->getCardinality();i++) { 00155 tags[tagSize*i] = 3; 00156 tags[tagSize*i+1] = 0; 00157 tags[tagSize*i+2] = i; 00158 tags[tagSize*i+3] = this->getCardinality(); 00159 } 00160 00161 Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0]; 00162 Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1]; 00163 Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2]; 00164 00165 00166 // now let's try to do it "right" 00167 // let's get the x, y, z bases and their dof 00168 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags(); 00169 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags(); 00170 const std::vector<std::vector<int> >& zdoftags = zBasis_.getAllDofTags(); 00171 00172 std::map<int,std::map<int,int> > total_dof_per_entity; 00173 std::map<int,std::map<int,int> > current_dof_per_entity; 00174 00175 // vertex dof 00176 for (int i=0;i<8;i++) { 00177 total_dof_per_entity[0][i] = 0; 00178 current_dof_per_entity[0][i] = 0; 00179 } 00180 // edge dof (12 edges) 00181 for (int i=0;i<12;i++) { 00182 total_dof_per_entity[1][i] = 0; 00183 current_dof_per_entity[1][i] = 0; 00184 } 00185 // face dof (6 faces) 00186 for (int i=0;i<6;i++) { 00187 total_dof_per_entity[2][i] = 0; 00188 current_dof_per_entity[2][i] = 0; 00189 } 00190 // internal dof 00191 total_dof_per_entity[3][0] = 0; 00192 //total_dof_per_entity[3][1] = 0; 00193 00194 // let's tally the total degrees of freedom on each entity 00195 for (int k=0;k<zBasis_.getCardinality();k++) { 00196 const int zdim = zdoftags[k][0]; 00197 const int zent = zdoftags[k][1]; 00198 for (int j=0;j<yBasis_.getCardinality();j++) { 00199 const int ydim = ydoftags[j][0]; 00200 const int yent = ydoftags[j][1]; 00201 for (int i=0;i<xBasis_.getCardinality();i++) { 00202 const int xdim = xdoftags[i][0]; 00203 const int xent = xdoftags[i][1]; 00204 int dofdim; 00205 int dofent; 00206 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent ); 00207 std::cout << i << " " << j << " " << k << "\t" << dofdim << " " << dofent << std::endl; 00208 total_dof_per_entity[dofdim][dofent] += 1; 00209 00210 } 00211 } 00212 } 00213 00214 int tagcur = 0; 00215 for (int k=0;k<zBasis_.getCardinality();k++) { 00216 const int zdim = zdoftags[k][0]; 00217 const int zent = zdoftags[k][1]; 00218 for (int j=0;j<yBasis_.getCardinality();j++) { 00219 const int ydim = ydoftags[j][0]; 00220 const int yent = ydoftags[j][1]; 00221 for (int i=0;i<xBasis_.getCardinality();i++) { 00222 const int xdim = xdoftags[i][0]; 00223 const int xent = xdoftags[i][1]; 00224 int dofdim; 00225 int dofent; 00226 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent ); 00227 tags[4*tagcur] = dofdim; 00228 tags[4*tagcur+1] = dofent; 00229 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent]; 00230 current_dof_per_entity[dofdim][dofent]++; 00231 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent]; 00232 tagcur++; 00233 } 00234 } 00235 } 00236 00237 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00238 this -> ordinalToTag_, 00239 &(tags[0]), 00240 this -> basisCardinality_, 00241 tagSize, 00242 posScDim, 00243 posScOrd, 00244 posDfOrd); 00245 00246 } 00247 00248 template<class Scalar, class ArrayScalar> 00249 void Basis_HGRAD_HEX_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues , 00250 const ArrayScalar &inputPoints , 00251 const EOperator operatorType ) const 00252 { 00253 #ifdef HAVE_INTREPID_DEBUG 00254 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00255 inputPoints, 00256 operatorType, 00257 this -> getBaseCellTopology(), 00258 this -> getCardinality() ); 00259 #endif 00260 00261 Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0]; 00262 Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1]; 00263 Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2]; 00264 00265 00266 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1); 00267 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1); 00268 FieldContainer<Scalar> zInputPoints(inputPoints.dimension(0),1); 00269 00270 for (int i=0;i<inputPoints.dimension(0);i++) { 00271 xInputPoints(i,0) = inputPoints(i,0); 00272 yInputPoints(i,0) = inputPoints(i,1); 00273 zInputPoints(i,0) = inputPoints(i,2); 00274 } 00275 00276 switch (operatorType) { 00277 case OPERATOR_VALUE: 00278 { 00279 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00280 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00281 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00282 00283 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00284 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00285 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE); 00286 00287 int bfcur = 0; 00288 for (int k=0;k<zBasis_.getCardinality();k++) { 00289 for (int j=0;j<yBasis_.getCardinality();j++) { 00290 for (int i=0;i<xBasis_.getCardinality();i++) { 00291 for (int l=0;l<inputPoints.dimension(0);l++) { 00292 outputValues(bfcur,l) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l); 00293 } 00294 bfcur++; 00295 } 00296 } 00297 } 00298 } 00299 break; 00300 case OPERATOR_D1: 00301 case OPERATOR_GRAD: 00302 { 00303 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00304 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00305 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00306 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00307 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00308 FieldContainer<Scalar> zBasisDerivs(zBasis_.getCardinality(),zInputPoints.dimension(0),1); 00309 00310 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE); 00311 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE); 00312 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE); 00313 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1); 00314 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1); 00315 zBasis_.getValues(zBasisDerivs,zInputPoints,OPERATOR_D1); 00316 00317 int bfcur = 0; 00318 for (int k=0;k<zBasis_.getCardinality();k++) { 00319 for (int j=0;j<yBasis_.getCardinality();j++) { 00320 for (int i=0;i<xBasis_.getCardinality();i++) { 00321 for (int l=0;l<inputPoints.dimension(0);l++) { 00322 outputValues(bfcur,l,0) = xBasisDerivs(i,l,0) * yBasisValues(j,l) * zBasisValues(k,l); 00323 outputValues(bfcur,l,1) = xBasisValues(i,l) * yBasisDerivs(j,l,0) * zBasisValues(k,l); 00324 outputValues(bfcur,l,2) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisDerivs(k,l,0); 00325 } 00326 bfcur++; 00327 } 00328 } 00329 } 00330 } 00331 break; 00332 case OPERATOR_D2: 00333 case OPERATOR_D3: 00334 case OPERATOR_D4: 00335 case OPERATOR_D5: 00336 case OPERATOR_D6: 00337 case OPERATOR_D7: 00338 case OPERATOR_D8: 00339 case OPERATOR_D9: 00340 case OPERATOR_D10: 00341 { 00342 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00343 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00344 FieldContainer<Scalar> zBasisValues(yBasis_.getCardinality(),zInputPoints.dimension(0)); 00345 00346 Teuchos::Array<int> partialMult; 00347 00348 for (int d=0;d<getDkCardinality(operatorType,3);d++) { 00349 getDkMultiplicities( partialMult , d , operatorType , 3 ); 00350 if (partialMult[0] == 0) { 00351 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00352 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE ); 00353 } 00354 else { 00355 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1); 00356 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 ); 00357 xBasis_.getValues( xBasisValues , xInputPoints, xop ); 00358 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0)); 00359 } 00360 if (partialMult[1] == 0) { 00361 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00362 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE ); 00363 } 00364 else { 00365 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1); 00366 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 ); 00367 yBasis_.getValues( yBasisValues , yInputPoints, yop ); 00368 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0)); 00369 } 00370 if (partialMult[2] == 0) { 00371 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00372 zBasis_.getValues( zBasisValues , zInputPoints, OPERATOR_VALUE ); 00373 } 00374 else { 00375 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0),1); 00376 EOperator zop = (EOperator) ( (int) OPERATOR_D1 + partialMult[2] - 1 ); 00377 zBasis_.getValues( zBasisValues , zInputPoints, zop ); 00378 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0)); 00379 } 00380 00381 00382 int bfcur = 0; 00383 for (int k=0;k<zBasis_.getCardinality();k++) { 00384 for (int j=0;j<yBasis_.getCardinality();j++) { 00385 for (int i=0;i<xBasis_.getCardinality();i++) { 00386 for (int l=0;l<inputPoints.dimension(0);l++) { 00387 outputValues(bfcur,l,d) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l); 00388 } 00389 bfcur++; 00390 } 00391 } 00392 } 00393 } 00394 } 00395 break; 00396 default: 00397 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument, 00398 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented"); 00399 break; 00400 } 00401 } 00402 00403 template<class Scalar,class ArrayScalar> 00404 void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00405 const ArrayScalar & inputPoints, 00406 const ArrayScalar & cellVertices, 00407 const EOperator operatorType) const { 00408 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00409 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): FEM Basis calling an FVD member function"); 00410 } 00411 00412 template<class Scalar,class ArrayScalar> 00413 void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const 00414 { 00415 int cur = 0; 00416 for (int k=0;k<ptsz_.dimension(0);k++) 00417 { 00418 for (int j=0;j<ptsy_.dimension(0);j++) 00419 { 00420 for (int i=0;i<ptsx_.dimension(0);i++) 00421 { 00422 dofCoords(cur,0) = ptsx_(i); 00423 dofCoords(cur,1) = ptsy_(j); 00424 dofCoords(cur,2) = ptsz_(k); 00425 cur++; 00426 } 00427 } 00428 } 00429 } 00430 00431 }// namespace Intrepid 00432 00433 #endif
1.7.6.1