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