Intrepid
Intrepid_HDIV_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+1)*(n+3)/2 )
57 {
58 const int N = n*(n+1)*(n+3)/2;
59 this -> basisCardinality_ = N;
60 this -> basisDegree_ = n;
61 this -> basisCellTopology_
62 = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
63 this -> basisType_ = BASIS_FEM_FIAT;
64 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
65 this -> basisTagsAreSet_ = false;
66
67
68 const int littleN = n*(n+1)*(n+2)/2; // dim of (P_{n-1})^3 -- smaller space
69 const int bigN = (n+1)*(n+2)*(n+3)/2; // dim of (P_{n})^2 -- larger space
70 const int start_PkH = (n-1)*n*(n+1)/6; // dim of P({n-2}), offset into
71 const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH;
72 const int scalarLittleN = littleN/3;
73 const int scalarBigN = bigN/3;
74
75 // first, need to project the basis for RT space onto the
76 // orthogonal basis of degree n
77 // get coefficients of PkHx
78
79 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
80
81 // basis for the space is
82 // { (phi_i,0,0) }_{i=0}^{scalarLittleN-1} ,
83 // { (0,phi_i,0) }_{i=0}^{scalarLittleN-1} ,
84 // { (0,0,phi_i) }_{i=0}^{scalarLittleN-1} ,
85 // { (x,y) . phi_i}_{i=startPKH}^{scalarLittleN-1}
86 // columns of V1 are expansion of this basis in terms of the orthogonal basis
87 // for P_{n}^3
88
89
90 // these two loops get the first three sets of basis functions
91 for (int i=0;i<scalarLittleN;i++) {
92 for (int k=0;k<3;k++) {
93 V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0;
94 }
95 }
96
97 // now I need to integrate { (x,y,z) phi } against the big basis
98 // first, get a cubature rule.
100 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 );
101 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
102 myCub.getCubature( cubPoints , cubWeights );
103
104 // tabulate the scalar orthonormal basis at cubature points
105 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
106 Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
107
108
109 // now do the integration
110 for (int i=0;i<dim_PkH;i++) {
111 for (int j=0;j<scalarBigN;j++) { // int (x,y,z) phi_i \cdot (phi_j,0,0)
112 V1(j,littleN+i) = 0.0;
113 for (int d=0;d<3;d++) {
114 for (int k=0;k<myCub.getNumPoints();k++) {
115 V1(j+d*scalarBigN,littleN+i) +=
116 cubWeights(k) * cubPoints(k,d)
117 * phisAtCubPoints(start_PkH+i,k)
118 * phisAtCubPoints(j,k);
119 }
120 }
121 }
122 }
123
124
125 // next, apply the RT nodes (rows) to the basis for (P_n)^3 (columns)
126 Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
127
128 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
129 const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
130 n+2 ,
131 1 );
132
133 FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 );
134 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
135 faceTop ,
136 n+2 ,
137 1 ,
138 pointType );
139
140 // holds the image of the triangle points on each face.
141 FieldContainer<Scalar> facePts( numPtsPerFace , 3 );
142 FieldContainer<Scalar> phisAtFacePoints( scalarBigN ,
143 numPtsPerFace );
144
145
146
147 // these are scaled by the appropriate face areas.
148 // area of faces 0,2,3 are 0.5
149 // area of face 1 is sqrt(3)/2
150
151 Scalar normal[][4] = { {0.0,0.5,-0.5,0.0},
152 {-0.5,0.5,0.0,0.0},
153 {0.0,0.5,0.0,-0.5} };
154
155 for (int i=0;i<4;i++) { // loop over faces
157 twoDPts ,
158 2 ,
159 i ,
160 this->basisCellTopology_ );
161
162 Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE );
163
164 // loop over points (rows of V2)
165 for (int j=0;j<numPtsPerFace;j++) {
166 // loop over orthonormal basis functions (columns of V2)
167 for (int k=0;k<scalarBigN;k++) {
168 for (int l=0;l<3;l++) {
169 V2(numPtsPerFace*i+j,k+l*scalarBigN) = normal[l][i] * phisAtFacePoints(k,j);
170 }
171 }
172 }
173 }
174
175 // remaining nodes point values of each vector component on interior
176 // points of a lattice of degree+2
177 // This way, RT0 --> degree = 1 and internal lattice has no points
178 // RT1 --> degree = 2, and internal lattice has one point (inside of quartic)
179 if (n > 1) {
180 const int numInternalPoints = PointTools::getLatticeSize( this->getBaseCellTopology() ,
181 n + 2 ,
182 1 );
183
184 FieldContainer<Scalar> internalPoints( numInternalPoints , 3 );
185 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
186 this->getBaseCellTopology() ,
187 n + 2 ,
188 1 ,
189 pointType );
190
191 FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
192 Phis_.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
193
194 // copy values into right positions of V2
195 for (int i=0;i<numInternalPoints;i++) {
196 for (int j=0;j<scalarBigN;j++) {
197 for (int k=0;k<3;k++) {
198 V2(4*numPtsPerFace+k*numInternalPoints+i,k*scalarBigN+j) = phisAtInternalPoints(j,i);
199 }
200 }
201 }
202 }
203
204 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
205
206 // multiply V2 * V1 --> V
207 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
208
209 // std::cout << "Vandermonde:\n";
210 // std::cout << Vsdm << "\n";
211 // std::cout << "End Vandermonde\n";
212
213 Teuchos::SerialDenseSolver<int,Scalar> solver;
214 solver.setMatrix( rcp( &Vsdm , false ) );
215 solver.invert( );
216
217
218 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
219 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
220
221 //std::cout << Csdm << "\n";
222
223 for (int i=0;i<bigN;i++) {
224 for (int j=0;j<N;j++) {
225 coeffs_(i,j) = Csdm(i,j);
226 }
227 }
228 }
229
230 template<class Scalar, class ArrayScalar>
232
233 // Basis-dependent initializations
234 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
235 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
236 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
237 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
238
239 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
240
241 int *tags = new int[ tagSize * this->getCardinality() ];
242 int *tag_cur = tags;
243 const int deg = this->getDegree();
244
245 const int numPtsPerFace = deg*(deg+1)/2;
246
247 // there are degree internal dofs on each edge -- normals. Let's do them
248 for (int f=0;f<4;f++) {
249 for (int i=0;i<numPtsPerFace;i++) {
250 tag_cur[0] = 2; tag_cur[1] = f; tag_cur[2] = i; tag_cur[3] = numPtsPerFace;
251 tag_cur += tagSize;
252 }
253 }
254 // end face dofs
255
256 // the rest of the dofs are internal
257 const int numInternalDof = this->getCardinality() - 4 * numPtsPerFace;
258 int internalDofCur=0;
259 for (int i=4*numPtsPerFace;i<this->getCardinality();i++) {
260 tag_cur[0] = 3; tag_cur[1] = 0; tag_cur[2] = internalDofCur; tag_cur[3] = numInternalDof;
261 tag_cur += tagSize;
262 internalDofCur++;
263 }
264
265
266 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
267 this -> ordinalToTag_,
268 tags,
269 this -> basisCardinality_,
270 tagSize,
271 posScDim,
272 posScOrd,
273 posDfOrd);
274
275 delete []tags;
276
277 }
278
279
280
281 template<class Scalar, class ArrayScalar>
283 const ArrayScalar & inputPoints,
284 const EOperator operatorType) const {
285
286 // Verify arguments
287#ifdef HAVE_INTREPID_DEBUG
288 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
289 inputPoints,
290 operatorType,
291 this -> getBaseCellTopology(),
292 this -> getCardinality() );
293#endif
294 const int numPts = inputPoints.dimension(0);
295 const int deg = this -> getDegree();
296 const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6;
297
298 try {
299 switch (operatorType) {
300 case OPERATOR_VALUE:
301 {
302 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
303 Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
304
305 for (int i=0;i<outputValues.dimension(0);i++) { // RT bf
306 for (int j=0;j<outputValues.dimension(1);j++) { // point
307 for (int l=0;l<3;l++) {
308 outputValues(i,j,l) = 0.0;
309 }
310 for (int k=0;k<scalarBigN;k++) { // Dubiner bf
311 for (int l=0;l<3;l++) { // vector components
312 outputValues(i,j,l) += coeffs_(k+l*scalarBigN,i) * phisCur(k,j);
313 }
314 }
315 }
316 }
317 }
318 break;
319 case OPERATOR_DIV:
320 {
321 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 );
322 Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
323 for (int i=0;i<outputValues.dimension(0);i++) { // bf loop
324 for (int j=0;j<outputValues.dimension(1);j++) { // point loop
325 outputValues(i,j) = 0.0;
326 for (int k=0;k<scalarBigN;k++) {
327 outputValues(i,j) += coeffs_(k,i) * phisCur(k,j,0);
328 outputValues(i,j) += coeffs_(k+scalarBigN,i) * phisCur(k,j,1);
329 outputValues(i,j) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,2);
330 }
331 }
332 }
333 }
334 break;
335 default:
336 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
337 ">>> ERROR (Basis_HDIV_TET_In_FEM): Operator type not implemented");
338 break;
339 }
340 }
341 catch (std::invalid_argument &exception){
342 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
343 ">>> ERROR (Basis_HDIV_TET_In_FEM): Operator type not implemented");
344 }
345
346 }
347
348
349
350 template<class Scalar, class ArrayScalar>
352 const ArrayScalar & inputPoints,
353 const ArrayScalar & cellVertices,
354 const EOperator operatorType) const {
355 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
356 ">>> ERROR (Basis_HDIV_TET_In_FEM): FEM Basis calling an FVD member function");
357 }
358
359
360}// 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_
expansion coefficients of the nodal basis in terms of the orthgonal one
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
Basis_HDIV_TET_In_FEM(const int n, const EPointType pointType)
Constructor.
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis_
Orthogonal basis out of which the nodal basis is constructed.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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 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.
Defines direct integration rules on a tetrahedron.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual int getNumPoints() const
Returns the number of cubature points.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
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...