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 
50 namespace 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 getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
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.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
FieldContainer< Scalar > coeffs_
expansion coefficients of the nodal basis in terms of the orthgonal one
EBasis basisType_
Type of the basis.
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.
Defines direct integration rules on a tetrahedron.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Basis_HDIV_TET_In_FEM(const int n, const EPointType pointType)
Constructor.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis_
Orthogonal basis out of which the nodal basis is constructed.
virtual int getNumPoints() const
Returns the number of cubature points.
EPointType
Enumeration of types of point distributions in Intrepid.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
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...