Intrepid
Intrepid_HCURL_TET_In_FEMDef.hpp
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50namespace Intrepid {
51
52 template<class Scalar, class ArrayScalar>
54 const EPointType pointType ):
55 Phis_( n ),
56 coeffs_( (n+1)*(n+2)*(n+3)/2 , n*(n+2)*(n+3)/2 )
57 {
58 const int N = n*(n+2)*(n+3)/2;
59 this -> basisCardinality_ = N;
60 this -> basisDegree_ = n;
61 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
62 this -> basisType_ = BASIS_FEM_FIAT;
63 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
64 this -> basisTagsAreSet_ = false;
65
66 const int littleN = n*(n+1)*(n+2)/2; // dim of (P_{n-1})^3 -- smaller space
67 const int bigN = (n+1)*(n+2)*(n+3)/2; // dim of (P_{n})^3 -- larger space
68 const int start_PkH = (n-1)*n*(n+1)/6; // dim of P({n-2}), offset into
69 const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH;
70 const int scalarLittleN = littleN/3;
71 const int scalarBigN = bigN/3;
72
73 // first, need to project the basis for Nedelec space onto the
74 // orthogonal basis of degree n
75 // get coefficients of PkHx
76
77 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, littleN + 3 * dim_PkH);
78
79 // these two loops get the first three sets of basis functions
80 for (int i=0;i<scalarLittleN;i++) {
81 for (int k=0;k<3;k++) {
82 V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0;
83 }
84 }
85
86 // first 3*scalarLittleN columns are (P_{n-1})^3 space
87
88
89 // now I need to integrate { (x,y,z) \times } against the big basis
90 // first, get a cubature rule.
91 CubatureDirectTetDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
92 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 );
93 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
94 myCub.getCubature( cubPoints , cubWeights );
95
96 // tabulate the scalar orthonormal basis at cubature points
97 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
98 Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
99
100
101
102 // first set of these functions will write into the first dimPkH columns of remainder
103
104 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
105 // write into second spatial component, where
106 // I integrate z phi_j phi_i
107 for (int i=0;i<scalarBigN;i++) {
108 V1(scalarBigN+i,littleN+j) = 0.0;
109 for (int k=0;k<myCub.getNumPoints();k++) {
110 V1(scalarBigN+i,littleN+j) -= cubWeights(k) * cubPoints(k,2)
111 * phisAtCubPoints(start_PkH+j,k)
112 * phisAtCubPoints(i,k);
113 }
114 }
115 // write into third spatial component (-y phi_j, phi_i)
116 for (int i=0;i<scalarBigN;i++) {
117 V1(2*scalarBigN+i,littleN+j) = 0.0;
118 for (int k=0;k<myCub.getNumPoints();k++) {
119 V1(2*scalarBigN+i,littleN+j) += cubWeights(k) * cubPoints(k,1)
120 * phisAtCubPoints(start_PkH+j,k)
121 * phisAtCubPoints(i,k);
122 }
123 }
124 }
125
126 // second set of basis functions, write into second set of dimPkH columns
127 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
128 // write into first spatial component, where
129 // I integrate -z phi_j phi_i
130 for (int i=0;i<scalarBigN;i++) {
131 V1(i,littleN+dim_PkH+j) = 0.0;
132 for (int k=0;k<myCub.getNumPoints();k++) {
133 V1(i,littleN+dim_PkH+j) += cubWeights(k) * cubPoints(k,2)
134 * phisAtCubPoints(start_PkH+j,k)
135 * phisAtCubPoints(i,k);
136 }
137 }
138
139 // third spatial component, x phi_j phi_i
140 for (int i=0;i<scalarBigN;i++) {
141 V1(2*scalarBigN+i,littleN+dim_PkH+j) = 0.0;
142 for (int k=0;k<myCub.getNumPoints();k++) {
143 V1(2*scalarBigN+i,littleN+dim_PkH+j) -= cubWeights(k) * cubPoints(k,0)
144 * phisAtCubPoints(start_PkH+j,k)
145 * phisAtCubPoints(i,k);
146 }
147 }
148 }
149
150 // third clump of dimPkH columns
151 for (int j=0;j<dim_PkH;j++) { // loop over homogeneous polynomials
152 // write into first spatial component, where
153 // I integrate y phi_j phi_i
154 for (int i=0;i<scalarBigN;i++) {
155 V1(i,littleN+2*dim_PkH+j) = 0.0;
156 for (int k=0;k<myCub.getNumPoints();k++) {
157 V1(i,littleN+2*dim_PkH+j) -= cubWeights(k) * cubPoints(k,1)
158 * phisAtCubPoints(start_PkH+j,k)
159 * phisAtCubPoints(i,k);
160 }
161 }
162 // second spatial component, -x phi_j phi_i
163 for (int i=0;i<scalarBigN;i++) {
164 V1(scalarBigN+i,littleN+2*dim_PkH+j) = 0.0;
165 for (int k=0;k<myCub.getNumPoints();k++) {
166 V1(scalarBigN+i,littleN+2*dim_PkH+j) += cubWeights(k) * cubPoints(k,0)
167 * phisAtCubPoints(start_PkH+j,k)
168 * phisAtCubPoints(i,k);
169 }
170 }
171 }
172
173 // now I need to set up an SVD to get a basis for the space
174 Teuchos::SerialDenseMatrix<int,Scalar> S(bigN,1);
175 Teuchos::SerialDenseMatrix<int,Scalar> U(bigN, bigN);
176 Teuchos::SerialDenseMatrix<int,Scalar> Vt(bigN,bigN);
177 Teuchos::SerialDenseMatrix<int,Scalar> work(5*bigN,1);
178 Teuchos::SerialDenseMatrix<int,Scalar> rWork(1,1);
179 int info;
180
181 Teuchos::LAPACK<int,Scalar> lala;
182
183 lala.GESVD( 'A',
184 'N',
185 V1.numRows() ,
186 V1.numCols() ,
187 V1.values() ,
188 V1.stride() ,
189 S.values() ,
190 U.values() ,
191 U.stride() ,
192 Vt.values() ,
193 Vt.stride() ,
194 work.values() ,
195 5*bigN ,
196 rWork.values() ,
197 &info );
198
199 int num_nonzero_sv = 0;
200 for (int i=0;i<bigN;i++) {
201 if (S(i,0) > INTREPID_TOL) {
202 num_nonzero_sv++;
203 }
204 }
205
206 Teuchos::SerialDenseMatrix<int,Scalar> Uslender(bigN, num_nonzero_sv);
207 for (int j=0;j<num_nonzero_sv;j++) {
208 for (int i=0;i<bigN;i++) {
209 Uslender(i,j) = U(i,j);
210 }
211 }
212
213 // apply nodes to big space
214 Teuchos::SerialDenseMatrix<int,Scalar> V2(N, bigN);
215
216 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
217 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
218
219
220 const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
221 n+1 ,
222 1 );
223
224 const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
225 n+1 ,
226 1 );
227
228 const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
229 n+1 ,
230 1 );
231
232 // these hold the reference domain points that will be mapped to each edge or face
233 FieldContainer<Scalar> oneDPts( numPtsPerEdge , 1 );
234 FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 );
235
236 if (pointType == POINTTYPE_WARPBLEND) {
237 CubatureDirectLineGauss<Scalar> edgeRule( numPtsPerEdge );
238 FieldContainer<Scalar> edgeCubWts( numPtsPerEdge );
239 edgeRule.getCubature( oneDPts , edgeCubWts );
240 }
241 else if (pointType == POINTTYPE_EQUISPACED ) {
242 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( oneDPts ,
243 edgeTop ,
244 n+1 ,
245 1 ,
246 pointType );
247 }
248
249 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
250 faceTop ,
251 n+1 ,
252 1 ,
253 pointType );
254
255 FieldContainer<Scalar> edgePts( numPtsPerEdge , 3 );
256 FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , numPtsPerEdge );
257
258 FieldContainer<Scalar> facePts( numPtsPerFace , 3 );
259 FieldContainer<Scalar> phisAtFacePoints( scalarBigN ,
260 numPtsPerFace );
261
262 FieldContainer<Scalar> edgeTan( 3 );
263
264 // loop over the edges
265 for (int edge=0;edge<6;edge++) {
267 edge ,
268 this->basisCellTopology_ );
269 /* multiply by 2.0 to account for a scaling in Pavel's definition */
270 for (int j=0;j<3;j++) {
271 edgeTan(j) *= 2.0;
272 }
273
275 oneDPts ,
276 1 ,
277 edge ,
278 this->basisCellTopology_ );
279
280 Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
281
282 // loop over points (rows of V2)
283 for (int j=0;j<numPtsPerEdge;j++) {
284 // loop over orthonormal basis functions (columns of V2)
285 for (int k=0;k<scalarBigN;k++) {
286 for (int d=0;d<3;d++) {
287 V2(edge*numPtsPerEdge+j,k+scalarBigN*d) = edgeTan(d) * phisAtEdgePoints(k,j);
288 }
289 }
290 }
291 }
292
293 // handle the faces, if needed
294 if (n > 1) {
295 FieldContainer<Scalar> refFaceTanU(3);
296 FieldContainer<Scalar> refFaceTanV(3);
297 for (int face=0;face<4;face++) {
299 refFaceTanV ,
300 face ,
301 this->basisCellTopology_ );
303 twoDPts ,
304 2 ,
305 face ,
306 this->basisCellTopology_ );
307 Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE );
308 for (int j=0;j<numPtsPerFace;j++) {
309 for (int k=0;k<scalarBigN;k++) {
310 for (int d=0;d<3;d++) {
311 V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j,k+scalarBigN*d) =
312 refFaceTanU(d) * phisAtFacePoints(k,j);
313 V2(6*numPtsPerEdge+2*face*numPtsPerFace+2*j+1,k+scalarBigN*d) =
314 refFaceTanV(d) * phisAtFacePoints(k,j);
315 }
316 }
317 }
318 }
319 }
320
321 // internal dof, if needed
322 if (n > 2) {
323 FieldContainer<Scalar> cellPoints( numPtsPerCell , 3 );
324 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( cellPoints ,
325 this->getBaseCellTopology() ,
326 n + 1 ,
327 1 ,
328 pointType );
329 FieldContainer<Scalar> phisAtCellPoints( scalarBigN , numPtsPerCell );
330 Phis_.getValues( phisAtCellPoints , cellPoints , OPERATOR_VALUE );
331 for (int i=0;i<numPtsPerCell;i++) {
332 for (int j=0;j<scalarBigN;j++) {
333 for (int k=0;k<3;k++) {
334 V2(6*numPtsPerEdge+8*numPtsPerFace+k*numPtsPerCell+i,k*scalarBigN+j) = phisAtCellPoints(j,i);
335 }
336 }
337 }
338 }
339
340 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
341
342 // multiply V2 * U --> V
343 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , Uslender , 0.0 );
344
345 Teuchos::SerialDenseSolver<int,Scalar> solver;
346 solver.setMatrix( rcp( &Vsdm , false ) );
347
348 solver.invert( );
349
350
351 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
352 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , Uslender , Vsdm , 0.0 );
353
354 //std::cout << Csdm << "\n";
355
356 for (int i=0;i<bigN;i++) {
357 for (int j=0;j<N;j++) {
358 coeffs_(i,j) = Csdm(i,j);
359 }
360 }
361
362 //std::cout << coeffs_ << std::endl;
363
364 }
365
366 template<class Scalar, class ArrayScalar>
368 // Basis-dependent initializations
369 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
370 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
371 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
372 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
373
374 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
375
376 int *tags = new int[ tagSize * this->getCardinality() ];
377 int *tag_cur = tags;
378 const int deg = this->getDegree();
379
380 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
381 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
382
383
384 const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
385 deg+1 ,
386 1 );
387
388 const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
389 deg+1 ,
390 1 );
391
392 const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
393 deg+1 ,
394 1 );
395
396 // edge dof first
397 for (int e=0;e<6;e++) {
398 for (int i=0;i<numPtsPerEdge;i++) {
399 tag_cur[0] = 1; tag_cur[1] = e; tag_cur[2] = i; tag_cur[3] = numPtsPerEdge;
400 tag_cur += tagSize;
401 }
402 }
403
404 // face dof, 2 * numPtsPerFace dof per face
405 for (int f=0;f<4;f++) {
406 for (int i=0;i<2*numPtsPerFace;i++) {
407 tag_cur[0] = 2; tag_cur[1] = f; tag_cur[2] = i; tag_cur[3] = 2*numPtsPerFace;
408 tag_cur+= tagSize;
409 }
410 }
411
412 // internal dof, 3 * numPtsPerCell
413 for (int i=0;i<3*numPtsPerCell;i++) {
414 tag_cur[0] = 3; tag_cur[1] = 0; tag_cur[2] = i; tag_cur[3] = 3*numPtsPerCell;
415 tag_cur += tagSize;
416 }
417 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
418 this -> ordinalToTag_,
419 tags,
420 this -> basisCardinality_,
421 tagSize,
422 posScDim,
423 posScOrd,
424 posDfOrd);
425
426 delete []tags;
427
428 }
429
430
431
432 template<class Scalar, class ArrayScalar>
434 const ArrayScalar & inputPoints,
435 const EOperator operatorType) const {
436
437 // Verify arguments
438#ifdef HAVE_INTREPID_DEBUG
439 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
440 inputPoints,
441 operatorType,
442 this -> getBaseCellTopology(),
443 this -> getCardinality() );
444#endif
445 const int numPts = inputPoints.dimension(0);
446 const int deg = this -> getDegree();
447 const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6;
448
449 try {
450 switch (operatorType) {
451 case OPERATOR_VALUE:
452 {
453 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
454 Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
455
456 for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
457 for (int j=0;j<outputValues.dimension(1);j++) { // point
458 for (int d=0;d<3;d++) {
459 outputValues(i,j,d) = 0.0;
460 }
461 for (int k=0;k<scalarBigN;k++) { // Dubiner bf
462 for (int d=0;d<3;d++) {
463 outputValues(i,j,d) += coeffs_(k+d*scalarBigN,i) * phisCur(k,j);
464 }
465 }
466 }
467 }
468 }
469 break;
470 case OPERATOR_CURL:
471 {
472 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 );
473 Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
474 for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
475 for (int j=0;j<outputValues.dimension(1);j++) { // point loop
476 outputValues(i,j,0) = 0.0;
477 for (int k=0;k<scalarBigN;k++) {
478 outputValues(i,j,0) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,1);
479 outputValues(i,j,0) -= coeffs_(k+scalarBigN,i) * phisCur(k,j,2);
480 }
481
482 outputValues(i,j,1) = 0.0;
483 for (int k=0;k<scalarBigN;k++) {
484 outputValues(i,j,1) += coeffs_(k,i) * phisCur(k,j,2);
485 outputValues(i,j,1) -= coeffs_(k+2*scalarBigN,i) * phisCur(k,j,0);
486 }
487
488 outputValues(i,j,2) = 0.0;
489 for (int k=0;k<scalarBigN;k++) {
490 outputValues(i,j,2) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0);
491 outputValues(i,j,2) -= coeffs_(k,i) * phisCur(k,j,1);
492 }
493 }
494 }
495 }
496 break;
497 default:
498 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
499 ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
500 break;
501 }
502 }
503 catch (std::invalid_argument &exception){
504 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
505 ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
506 }
507
508 }
509
510
511
512 template<class Scalar, class ArrayScalar>
514 const ArrayScalar & inputPoints,
515 const ArrayScalar & cellVertices,
516 const EOperator operatorType) const {
517 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
518 ">>> ERROR (Basis_HCURL_TET_In_FEM): FEM Basis calling an FVD member function");
519 }
520
521
522}// namespace Intrepid
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
FieldContainer< Scalar > coeffs_
Array holding the expansion coefficients of the nodal basis in terms of Phis_.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis_
Orthogonal basis of ofder n, in terms of which the H(curl) basis functions are expressed.
Basis_HCURL_TET_In_FEM(const int n, const EPointType pointType)
Constructor.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....
static void getReferenceEdgeTangent(ArrayEdgeTangent &refEdgeTangent, const int &edgeOrd, const shards::CellTopology &parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getReferenceFaceTangents(ArrayFaceTangentU &refFaceTanU, ArrayFaceTangentV &refFaceTanV, const int &faceOrd, const shards::CellTopology &parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...