Intrepid
Intrepid_HGRAD_TET_Cn_FEMDef.hpp
1 #ifndef INTREPID_HGRAD_TET_CN_FEMDEF_HPP
2 #define INTREPID_HGRAD_TET_CN_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53  template<class Scalar, class ArrayScalar>
55  const EPointType pointType ):
56  Phis( n ),
57  V((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
58  Vinv((n+1)*(n+2)*(n+3)/6,(n+1)*(n+2)*(n+3)/6),
59  latticePts( (n+1)*(n+2)*(n+3)/6 , 3 )
60  {
61  const int N = (n+1)*(n+2)*(n+3)/6;
62  this -> basisCardinality_ = N;
63  this -> basisDegree_ = n;
64  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
65  this -> basisType_ = BASIS_FEM_FIAT;
66  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
67  this -> basisTagsAreSet_ = false;
68 
69  // construct lattice
70 
71  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
72  this->getBaseCellTopology() ,
73  n ,
74  0 ,
75  pointType );
76 
77 
78  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
79  // so we transpose on copy below.
80 
81  Phis.getValues( V , latticePts , OPERATOR_VALUE );
82 
83  // now I need to copy V into a Teuchos array to do the inversion
84  Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
85  for (int i=0;i<N;i++) {
86  for (int j=0;j<N;j++) {
87  Vsdm(i,j) = V(i,j);
88  }
89  }
90 
91  // invert the matrix
92  Teuchos::SerialDenseSolver<int,Scalar> solver;
93  solver.setMatrix( rcp( &Vsdm , false ) );
94  solver.invert( );
95 
96  // now I need to copy the inverse into Vinv
97  for (int i=0;i<N;i++) {
98  for (int j=0;j<N;j++) {
99  Vinv(i,j) = Vsdm(j,i);
100  }
101  }
102 
103  }
104 
105 
106  template<class Scalar, class ArrayScalar>
108 
109  // Basis-dependent initializations
110  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
111  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
112  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
113  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
114 
115  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
116 
117  int *tags = new int[ tagSize * this->getCardinality() ];
118  int *tag_cur = tags;
119  const int degree = this->getDegree();
120  const int numEdgeDof = degree - 1;
121  const int numFaceDof = PointTools::getLatticeSize( shards::CellTopology( shards::getCellTopologyData<shards::Triangle<3> >() ) ,
122  degree ,
123  1);
124  const int numCellDof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
125  degree ,
126  1 );
127  int edge_dof_cur[] = {0,0,0,0,0,0};
128  int face_dof_cur[] = {0,0,0,0};
129  int cell_dof_cur = 0;
130 
131  // this is the really big mess :(
132  // BEGIN DOF ON BOTTOM FACE
133  // first vertex: 0
134  tag_cur[0] = 0; tag_cur[1] = 0; tag_cur[2] = 0; tag_cur[3] = 1;
135  tag_cur += tagSize;
136  // end first vertex
137 
138  // internal points on line from vertex 0 to vertex 1. This is edge 0
139  for (int i=1;i<degree;i++) {
140  tag_cur[0] = 1; tag_cur[1] = 0; tag_cur[2] = edge_dof_cur[0]; tag_cur[3] = numEdgeDof;
141  edge_dof_cur[0]++;
142  tag_cur += tagSize;
143  }
144  // end line from vertex 0 to vertex 1
145 
146  // begin vertex 1
147  tag_cur[0] = 0; tag_cur[1] = 1; tag_cur[2] = 0; tag_cur[3] = 1;
148  tag_cur += tagSize;
149  // end vertex 1
150 
151  // internal lines on bottom face
152  for (int i=1;i<degree;i++) {
153  // first dof is on edge 2
154  tag_cur[0] = 1; tag_cur[1] = 2; tag_cur[2] = edge_dof_cur[2]; tag_cur[3] = numEdgeDof;
155  edge_dof_cur[2]++;
156  tag_cur += tagSize;
157  // end dof on edge 2
158 
159  // internal points are on bottom face, which is face 3
160  for (int j=1;j<degree-i;j++) {
161  tag_cur[0] = 2; tag_cur[1] = 3; tag_cur[2] = face_dof_cur[3]; tag_cur[3] = numFaceDof;
162  face_dof_cur[3]++;
163  tag_cur += tagSize;
164  }
165  // end internal points on face
166 
167  // last dof is on edge 1
168  tag_cur[0] = 1; tag_cur[1] = 1; tag_cur[2] = edge_dof_cur[1]; tag_cur[3] = numEdgeDof;
169  edge_dof_cur[1]++;
170  tag_cur += tagSize;
171  // end dof on edge 1
172  }
173  // end internal lines on bottom face
174 
175  // vertex 2 on bottom face
176  tag_cur[0] = 0; tag_cur[1] = 2; tag_cur[2] = 0; tag_cur[3] = 1;
177  tag_cur += tagSize;
178  // end vertex 2 on bottom face
179 
180  // END DOF ON BOTTOM FACE
181 
182  // BEGIN DOF ON INTERNAL FACE SLICES (ascending z)
183  for (int i=1;i<degree;i++) {
184  // bottom line of internal face
185  // first point is on edge 3 (from vertex 0 to 3)
186  tag_cur[0] = 1; tag_cur[1] = 3; tag_cur[2] = edge_dof_cur[3]; tag_cur[3] = numEdgeDof;
187  edge_dof_cur[3]++;
188  tag_cur += tagSize;
189  // end first point
190  // points internal to face of vertices (0,1,3), which is face 0
191  for (int j=1;j<degree-i;j++) {
192  tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = face_dof_cur[0]; tag_cur[3] = numFaceDof;
193  face_dof_cur[0]++;
194  tag_cur += tagSize;
195  }
196  // end points internal to face 0
197  // last point on bottom line is on edge 4
198  tag_cur[0] = 1; tag_cur[1] = 4; tag_cur[2] = edge_dof_cur[4]; tag_cur[3] = numEdgeDof;
199  edge_dof_cur[4]++;
200  tag_cur += tagSize;
201  // end last point on bottom edge
202  // end bottom line of internal face
203 
204  // begin internal lines of internal face
205  for (int j=1;j<degree-i;j++) {
206  // first point on line is on face of vertices (0,3,2), which is face 2
207  tag_cur[0] = 2; tag_cur[1] = 2; tag_cur[2] = face_dof_cur[2]; tag_cur[3] = numFaceDof;
208  face_dof_cur[2]++;
209  tag_cur += tagSize;
210  // end first point of line
211  // begin internal points on the cell
212  for (int k=1;k<degree-i-j;k++) {
213  tag_cur[0] = 3; tag_cur[1] = 0; tag_cur[2] = cell_dof_cur; tag_cur[3] = numCellDof;
214  cell_dof_cur++;
215  tag_cur += tagSize;
216  }
217  // end internal points on the cell
218  // last point on the line is on face with vertices (1,2,3) , which is face 1
219  tag_cur[0] = 2; tag_cur[1] = 1; tag_cur[2] = face_dof_cur[1]; tag_cur[3] = numFaceDof;
220  face_dof_cur[1]++;
221  tag_cur += tagSize;
222  // end last point of line
223  }
224  // end internal lines of internal face
225  // begin top point on current face slice: on edge 5
226  tag_cur[0] = 1; tag_cur[1] = 5; tag_cur[2] = edge_dof_cur[5]; tag_cur[3] = numEdgeDof;
227  edge_dof_cur[5]++;
228  tag_cur += 4;
229  // end top point on current face slice
230  }
231  // END DOF ON INTERNAL FACE SLICES
232 
233  // TOP VERTEX: 3
234  tag_cur[0] = 0; tag_cur[1] = 3; tag_cur[2] = 0; tag_cur[3] = 1;
235  // END TOP VERTEX:3
236 
237  // end of really big mess :)
238 
239 
240 
241  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
242  this -> ordinalToTag_,
243  tags,
244  this -> basisCardinality_,
245  tagSize,
246  posScDim,
247  posScOrd,
248  posDfOrd);
249 
250  delete []tags;
251 
252  }
253 
254 
255 
256  template<class Scalar, class ArrayScalar>
258  const ArrayScalar & inputPoints,
259  const EOperator operatorType) const {
260 
261  // Verify arguments
262 #ifdef HAVE_INTREPID_DEBUG
263  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
264  inputPoints,
265  operatorType,
266  this -> getBaseCellTopology(),
267  this -> getCardinality() );
268 #endif
269  const int numPts = inputPoints.dimension(0);
270  const int numBf = this->getCardinality();
271 
272  try {
273  switch (operatorType) {
274  case OPERATOR_VALUE:
275  {
276  FieldContainer<Scalar> phisCur( numBf , numPts );
277  Phis.getValues( phisCur , inputPoints , operatorType );
278  for (int i=0;i<outputValues.dimension(0);i++) {
279  for (int j=0;j<outputValues.dimension(1);j++) {
280  outputValues(i,j) = 0.0;
281  for (int k=0;k<this->getCardinality();k++) {
282  outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
283  }
284  }
285  }
286  }
287  break;
288  case OPERATOR_GRAD:
289  case OPERATOR_D1:
290  case OPERATOR_D2:
291  case OPERATOR_D3:
292  case OPERATOR_D4:
293  case OPERATOR_D5:
294  case OPERATOR_D6:
295  case OPERATOR_D7:
296  case OPERATOR_D8:
297  case OPERATOR_D9:
298  case OPERATOR_D10:
299  {
300  const int dkcard =
301  (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,3): getDkCardinality(operatorType,3);
302 
303  FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
304  Phis.getValues( phisCur , inputPoints , operatorType );
305 
306  for (int i=0;i<outputValues.dimension(0);i++) {
307  for (int j=0;j<outputValues.dimension(1);j++) {
308  for (int k=0;k<outputValues.dimension(2);k++) {
309  outputValues(i,j,k) = 0.0;
310  for (int l=0;l<this->getCardinality();l++) {
311  outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
312  }
313  }
314  }
315  }
316 
317  }
318  break;
319  default:
320  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
321  ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");
322  break;
323  }
324  }
325  catch (std::invalid_argument &exception){
326  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
327  ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): Operator type not implemented");
328  }
329 
330  }
331 
332 
333 
334  template<class Scalar, class ArrayScalar>
336  const ArrayScalar & inputPoints,
337  const ArrayScalar & cellVertices,
338  const EOperator operatorType) const {
339  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
340  ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): FEM Basis calling an FVD member function");
341  }
342 
343 
344 }// namespace Intrepid
345 #endif
346 
347 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
348 #ifdef __GNUC__
349 #warning "The Intrepid package is deprecated"
350 #endif
351 #endif
352 
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis
The orthogonal basis on triangles, out of which the nodal basis is constructed.
EBasis basisType_
Type of the basis.
Basis_HGRAD_TET_Cn_FEM(const int n, const EPointType pointType)
Constructor.
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
FieldContainer< Scalar > Vinv
The inverse of V. The columns of Vinv express the Lagrange basis in terms of the orthogonal basis...
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
FieldContainer< Scalar > latticePts
stores the points at which degrees of freedom are located.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
FieldContainer< Scalar > V
The Vandermonde matrix with V_{ij} = phi_i(x_j), where x_j is the j_th point in the lattice...
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...