|
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). 00026 // Kara Peterson (kjpeter@sandia.gov). 00027 // 00028 // ************************************************************************ 00029 // @HEADER 00030 00037 namespace Intrepid { 00038 00039 template<class Scalar, class ArrayScalar> 00040 Basis_HCURL_TET_In_FEM<Scalar,ArrayScalar>::Basis_HCURL_TET_In_FEM( const int n , 00041 const EPointType pointType ): 00042 Phis_( n ), 00043 coeffs_( (n+1)*(n+2)*(n+3)/2 , n*(n+2)*(n+3)/2 ) 00044 { 00045 const int N = n*(n+2)*(n+3)/2; 00046 this -> basisCardinality_ = N; 00047 this -> basisDegree_ = n; 00048 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() ); 00049 this -> basisType_ = BASIS_FEM_FIAT; 00050 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00051 this -> basisTagsAreSet_ = false; 00052 00053 const int littleN = n*(n+1)*(n+2)/2; // dim of (P_{n-1})^3 -- smaller space 00054 const int bigN = (n+1)*(n+2)*(n+3)/2; // dim of (P_{n})^3 -- larger space 00055 const int start_PkH = (n-1)*n*(n+1)/6; // dim of P({n-2}), offset into 00056 const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH; 00057 const int scalarLittleN = littleN/3; 00058 const int scalarBigN = bigN/3; 00059 00060 // first, need to project the basis for Nedelec space onto the 00061 // orthogonal basis of degree n 00062 // get coefficients of PkHx 00063 00064 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, littleN + 3 * dim_PkH); 00065 00066 // these two loops get the first three sets of basis functions 00067 for (int i=0;i<scalarLittleN;i++) { 00068 for (int k=0;k<3;k++) { 00069 V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0; 00070 } 00071 } 00072 00073 // first 3*scalarLittleN columns are (P_{n-1})^3 space 00074 00075 00076 // now I need to integrate { (x,y,z) \times } against the big basis 00077 // first, get a cubature rule. 00078 CubatureDirectTetDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n ); 00079 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 ); 00080 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() ); 00081 myCub.getCubature( cubPoints , cubWeights ); 00082 00083 // tabulate the scalar orthonormal basis at cubature points 00084 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() ); 00085 Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE ); 00086 00087 00088 00089 // first set of these functions will write into the first dimPkH columns of remainder 00090 00091 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials 00092 // write into second spatial component, where 00093 // I integrate z phi_j phi_i 00094 for (int i=0;i<scalarBigN;i++) { 00095 V1(scalarBigN+i,littleN+j) = 0.0; 00096 for (int k=0;k<myCub.getNumPoints();k++) { 00097 V1(scalarBigN+i,littleN+j) -= cubWeights(k) * cubPoints(k,2) 00098 * phisAtCubPoints(start_PkH+j,k) 00099 * phisAtCubPoints(i,k); 00100 } 00101 } 00102 // write into third spatial component (-y phi_j, phi_i) 00103 for (int i=0;i<scalarBigN;i++) { 00104 V1(2*scalarBigN+i,littleN+j) = 0.0; 00105 for (int k=0;k<myCub.getNumPoints();k++) { 00106 V1(2*scalarBigN+i,littleN+j) += cubWeights(k) * cubPoints(k,1) 00107 * phisAtCubPoints(start_PkH+j,k) 00108 * phisAtCubPoints(i,k); 00109 } 00110 } 00111 } 00112 00113 // second set of basis functions, write into second set of dimPkH columns 00114 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials 00115 // write into first spatial component, where 00116 // I integrate -z phi_j phi_i 00117 for (int i=0;i<scalarBigN;i++) { 00118 V1(i,littleN+dim_PkH+j) = 0.0; 00119 for (int k=0;k<myCub.getNumPoints();k++) { 00120 V1(i,littleN+dim_PkH+j) += cubWeights(k) * cubPoints(k,2) 00121 * phisAtCubPoints(start_PkH+j,k) 00122 * phisAtCubPoints(i,k); 00123 } 00124 } 00125 00126 // third spatial component, x phi_j phi_i 00127 for (int i=0;i<scalarBigN;i++) { 00128 V1(2*scalarBigN+i,littleN+dim_PkH+j) = 0.0; 00129 for (int k=0;k<myCub.getNumPoints();k++) { 00130 V1(2*scalarBigN+i,littleN+dim_PkH+j) -= cubWeights(k) * cubPoints(k,0) 00131 * phisAtCubPoints(start_PkH+j,k) 00132 * phisAtCubPoints(i,k); 00133 } 00134 } 00135 } 00136 00137 // third clump of dimPkH columns 00138 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials 00139 // write into first spatial component, where 00140 // I integrate y phi_j phi_i 00141 for (int i=0;i<scalarBigN;i++) { 00142 V1(i,littleN+2*dim_PkH+j) = 0.0; 00143 for (int k=0;k<myCub.getNumPoints();k++) { 00144 V1(i,littleN+2*dim_PkH+j) -= cubWeights(k) * cubPoints(k,1) 00145 * phisAtCubPoints(start_PkH+j,k) 00146 * phisAtCubPoints(i,k); 00147 } 00148 } 00149 // second spatial component, -x phi_j phi_i 00150 for (int i=0;i<scalarBigN;i++) { 00151 V1(scalarBigN+i,littleN+2*dim_PkH+j) = 0.0; 00152 for (int k=0;k<myCub.getNumPoints();k++) { 00153 V1(scalarBigN+i,littleN+2*dim_PkH+j) += cubWeights(k) * cubPoints(k,0) 00154 * phisAtCubPoints(start_PkH+j,k) 00155 * phisAtCubPoints(i,k); 00156 } 00157 } 00158 } 00159 00160 // now I need to set up an SVD to get a basis for the space 00161 Teuchos::SerialDenseMatrix<int,Scalar> S(bigN,1); 00162 Teuchos::SerialDenseMatrix<int,Scalar> U(bigN, bigN); 00163 Teuchos::SerialDenseMatrix<int,Scalar> Vt(bigN,bigN); 00164 Teuchos::SerialDenseMatrix<int,Scalar> work(5*bigN,1); 00165 Teuchos::SerialDenseMatrix<int,Scalar> rWork(1,1); 00166 int info; 00167 00168 Teuchos::LAPACK<int,Scalar> lala; 00169 00170 lala.GESVD( 'A', 00171 'N', 00172 V1.numRows() , 00173 V1.numCols() , 00174 V1.values() , 00175 V1.stride() , 00176 S.values() , 00177 U.values() , 00178 U.stride() , 00179 Vt.values() , 00180 Vt.stride() , 00181 work.values() , 00182 5*bigN , 00183 rWork.values() , 00184 &info ); 00185 00186 int num_nonzero_sv = 0; 00187 for (int i=0;i<bigN;i++) { 00188 if (S(i,0) > INTREPID_TOL) { 00189 num_nonzero_sv++; 00190 } 00191 } 00192 00193 Teuchos::SerialDenseMatrix<int,Scalar> Uslender(bigN, num_nonzero_sv); 00194 for (int j=0;j<num_nonzero_sv;j++) { 00195 for (int i=0;i<bigN;i++) { 00196 Uslender(i,j) = U(i,j); 00197 } 00198 } 00199 00200 // apply nodes to big space 00201 Teuchos::SerialDenseMatrix<int,Scalar> V2(N, bigN); 00202 00203 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() ); 00204 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() ); 00205 00206 00207 const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop , 00208 n+1 , 00209 1 ); 00210 00211 const int numPtsPerFace = PointTools::getLatticeSize( faceTop , 00212 n+1 , 00213 1 ); 00214 00215 const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ , 00216 n+1 , 00217 1 ); 00218 00219 // these hold the reference domain points that will be mapped to each edge or face 00220 FieldContainer<Scalar> oneDPts( numPtsPerEdge , 1 ); 00221 FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 ); 00222 00223 if (pointType == POINTTYPE_WARPBLEND) { 00224 CubatureDirectLineGauss<Scalar> edgeRule( numPtsPerEdge ); 00225 FieldContainer<Scalar> edgeCubWts( numPtsPerEdge ); 00226 edgeRule.getCubature( oneDPts , edgeCubWts ); 00227 } 00228 else if (pointType == POINTTYPE_EQUISPACED ) { 00229 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( oneDPts , 00230 edgeTop , 00231 n+1 , 00232 1 , 00233 pointType ); 00234 } 00235 00236 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts , 00237 faceTop , 00238 n+1 , 00239 1 , 00240 pointType ); 00241 00242 FieldContainer<Scalar> edgePts( numPtsPerEdge , 3 ); 00243 FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , numPtsPerEdge ); 00244 00245 FieldContainer<Scalar> facePts( numPtsPerFace , 3 ); 00246 FieldContainer<Scalar> phisAtFacePoints( scalarBigN , 00247 numPtsPerFace ); 00248 00249 FieldContainer<Scalar> edgeTan( 3 ); 00250 00251 // loop over the edges 00252 for (int edge=0;edge<6;edge++) { 00253 CellTools<Scalar>::getReferenceEdgeTangent( edgeTan , 00254 edge , 00255 this->basisCellTopology_ ); 00256 /* multiply by 2.0 to account for a scaling in Pavel's definition */ 00257 for (int j=0;j<3;j++) { 00258 edgeTan(j) *= 2.0; 00259 } 00260 00261 CellTools<Scalar>::mapToReferenceSubcell( edgePts , 00262 oneDPts , 00263 1 , 00264 edge , 00265 this->basisCellTopology_ ); 00266 00267 Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE ); 00268 00269 // loop over points (rows of V2) 00270 for (int j=0;j<numPtsPerEdge;j++) { 00271 // loop over orthonormal basis functions (columns of V2) 00272 for (int k=0;k<scalarBigN;k++) { 00273 for (int d=0;d<3;d++) { 00274 V2(edge*numPtsPerEdge+j,k+scalarBigN*d) = edgeTan(d) * phisAtEdgePoints(k,j); 00275 } 00276 } 00277 } 00278 } 00279 00280 // handle the faces, if needed 00281 if (n > 1) { 00282 FieldContainer<Scalar> refFaceTanU(3); 00283 FieldContainer<Scalar> refFaceTanV(3); 00284 for (int face=0;face<4;face++) { 00285 CellTools<Scalar>::getReferenceFaceTangents( refFaceTanU , 00286 refFaceTanV , 00287 face , 00288 this->basisCellTopology_ ); 00289 CellTools<Scalar>::mapToReferenceSubcell( facePts , 00290 twoDPts , 00291 2 , 00292 face , 00293 this->basisCellTopology_ ); 00294 Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE ); 00295 for (int j=0;j<numPtsPerFace;j++) { 00296 for (int k=0;k<scalarBigN;k++) { 00297 for (int d=0;d<3;d++) { 00298 V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j,k+scalarBigN*d) = 00299 refFaceTanU(d) * phisAtFacePoints(k,j); 00300 V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j+1,k+scalarBigN*d) = 00301 refFaceTanV(d) * phisAtFacePoints(k,j); 00302 } 00303 } 00304 } 00305 } 00306 } 00307 00308 // internal dof, if needed 00309 if (n > 2) { 00310 FieldContainer<Scalar> cellPoints( numPtsPerCell , 3 ); 00311 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( cellPoints , 00312 this->getBaseCellTopology() , 00313 n + 1 , 00314 1 , 00315 pointType ); 00316 FieldContainer<Scalar> phisAtCellPoints( scalarBigN , numPtsPerCell ); 00317 Phis_.getValues( phisAtCellPoints , cellPoints , OPERATOR_VALUE ); 00318 for (int i=0;i<numPtsPerCell;i++) { 00319 for (int j=0;j<scalarBigN;j++) { 00320 for (int k=0;k<3;k++) { 00321 V2(6*numPtsPerEdge+8*numPtsPerFace+k*numPtsPerCell+i,k*scalarBigN+j) = phisAtCellPoints(j,i); 00322 } 00323 } 00324 } 00325 } 00326 00327 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N ); 00328 00329 // multiply V2 * U --> V 00330 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , Uslender , 0.0 ); 00331 00332 Teuchos::SerialDenseSolver<int,Scalar> solver; 00333 solver.setMatrix( rcp( &Vsdm , false ) ); 00334 00335 solver.invert( ); 00336 00337 00338 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N ); 00339 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , Uslender , Vsdm , 0.0 ); 00340 00341 //std::cout << Csdm << "\n"; 00342 00343 for (int i=0;i<bigN;i++) { 00344 for (int j=0;j<N;j++) { 00345 coeffs_(i,j) = Csdm(i,j); 00346 } 00347 } 00348 00349 //std::cout << coeffs_ << std::endl; 00350 00351 } 00352 00353 template<class Scalar, class ArrayScalar> 00354 void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::initializeTags() { 00355 // Basis-dependent initializations 00356 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00357 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00358 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00359 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00360 00361 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00362 00363 int *tags = new int[ tagSize * this->getCardinality() ]; 00364 int *tag_cur = tags; 00365 const int deg = this->getDegree(); 00366 00367 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() ); 00368 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() ); 00369 00370 00371 const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop , 00372 deg+1 , 00373 1 ); 00374 00375 const int numPtsPerFace = PointTools::getLatticeSize( faceTop , 00376 deg+1 , 00377 1 ); 00378 00379 const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ , 00380 deg+1 , 00381 1 ); 00382 00383 // edge dof first 00384 for (int e=0;e<6;e++) { 00385 for (int i=0;i<numPtsPerEdge;i++) { 00386 tag_cur[0] = 1; tag_cur[1] = e; tag_cur[2] = i; tag_cur[3] = numPtsPerEdge; 00387 tag_cur += tagSize; 00388 } 00389 } 00390 00391 // face dof, 2 * numPtsPerFace dof per face 00392 for (int f=0;f<4;f++) { 00393 for (int i=0;i<2*numPtsPerFace;i++) { 00394 tag_cur[0] = 2; tag_cur[1] = f; tag_cur[2] = i; tag_cur[3] = 2*numPtsPerFace; 00395 tag_cur+= tagSize; 00396 } 00397 } 00398 00399 // internal dof, 3 * numPtsPerCell 00400 for (int i=0;i<3*numPtsPerCell;i++) { 00401 tag_cur[0] = 3; tag_cur[1] = 0; tag_cur[2] = i; tag_cur[3] = 3*numPtsPerCell; 00402 tag_cur += tagSize; 00403 } 00404 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00405 this -> ordinalToTag_, 00406 tags, 00407 this -> basisCardinality_, 00408 tagSize, 00409 posScDim, 00410 posScOrd, 00411 posDfOrd); 00412 00413 delete []tags; 00414 00415 } 00416 00417 00418 00419 template<class Scalar, class ArrayScalar> 00420 void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00421 const ArrayScalar & inputPoints, 00422 const EOperator operatorType) const { 00423 00424 // Verify arguments 00425 #ifdef HAVE_INTREPID_DEBUG 00426 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues, 00427 inputPoints, 00428 operatorType, 00429 this -> getBaseCellTopology(), 00430 this -> getCardinality() ); 00431 #endif 00432 const int numPts = inputPoints.dimension(0); 00433 const int deg = this -> getDegree(); 00434 const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6; 00435 00436 try { 00437 switch (operatorType) { 00438 case OPERATOR_VALUE: 00439 { 00440 FieldContainer<Scalar> phisCur( scalarBigN , numPts ); 00441 Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE ); 00442 00443 for (int i=0;i<outputValues.dimension(0);i++) { // RT bf 00444 for (int j=0;j<outputValues.dimension(1);j++) { // point 00445 for (int d=0;d<3;d++) { 00446 outputValues(i,j,d) = 0.0; 00447 } 00448 for (int k=0;k<scalarBigN;k++) { // Dubiner bf 00449 for (int d=0;d<3;d++) { 00450 outputValues(i,j,d) += coeffs_(k+d*scalarBigN,i) * phisCur(k,j); 00451 } 00452 } 00453 } 00454 } 00455 } 00456 break; 00457 case OPERATOR_CURL: 00458 { 00459 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 ); 00460 Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD ); 00461 for (int i=0;i<outputValues.dimension(0);i++) { // bf loop 00462 for (int j=0;j<outputValues.dimension(1);j++) { // point loop 00463 outputValues(i,j,0) = 0.0; 00464 for (int k=0;k<scalarBigN;k++) { 00465 outputValues(i,j,0) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,1); 00466 outputValues(i,j,0) -= coeffs_(k+scalarBigN,i) * phisCur(k,j,2); 00467 } 00468 00469 outputValues(i,j,1) = 0.0; 00470 for (int k=0;k<scalarBigN;k++) { 00471 outputValues(i,j,1) += coeffs_(k,i) * phisCur(k,j,2); 00472 outputValues(i,j,1) -= coeffs_(k+2*scalarBigN,i) * phisCur(k,j,0); 00473 } 00474 00475 outputValues(i,j,2) = 0.0; 00476 for (int k=0;k<scalarBigN;k++) { 00477 outputValues(i,j,2) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0); 00478 outputValues(i,j,2) -= coeffs_(k,i) * phisCur(k,j,1); 00479 } 00480 } 00481 } 00482 } 00483 break; 00484 default: 00485 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00486 ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented"); 00487 break; 00488 } 00489 } 00490 catch (std::invalid_argument &exception){ 00491 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00492 ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented"); 00493 } 00494 00495 } 00496 00497 00498 00499 template<class Scalar, class ArrayScalar> 00500 void Basis_HCURL_TET_In_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00501 const ArrayScalar & inputPoints, 00502 const ArrayScalar & cellVertices, 00503 const EOperator operatorType) const { 00504 TEST_FOR_EXCEPTION( (true), std::logic_error, 00505 ">>> ERROR (Basis_HCURL_TET_In_FEM): FEM Basis calling an FVD member function"); 00506 } 00507 00508 00509 }// namespace Intrepid
1.7.4