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
51namespace 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
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.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Basis_HGRAD_TET_Cn_FEM(const int n, const EPointType pointType)
Constructor.
FieldContainer< Scalar > latticePts
stores the points at which degrees of freedom are located.
FieldContainer< Scalar > Vinv
The inverse of V. The columns of Vinv express the Lagrange basis in terms of the orthogonal basis.
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis
The orthogonal basis on triangles, out of which the nodal basis is constructed.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
FieldContainer< Scalar > V
The Vandermonde matrix with V_{ij} = phi_i(x_j), where x_j is the j_th point in the lattice.
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 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...