Intrepid
Intrepid_HGRAD_TRI_Cn_FEMDef.hpp
1#ifndef INTREPID_HGRAD_TRI_CN_FEMDEF_HPP
2#define INTREPID_HGRAD_TRI_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)/2,(n+1)*(n+2)/2),
58 Vinv((n+1)*(n+2)/2,(n+1)*(n+2)/2),
59 latticePts( (n+1)*(n+2)/2 , 2 )
60 {
61 TEUCHOS_TEST_FOR_EXCEPTION( n <= 0, std::invalid_argument, "polynomial order must be >= 1");
62
63 const int N = (n+1)*(n+2)/2;
64 this -> basisCardinality_ = N;
65 this -> basisDegree_ = n;
66 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
67 this -> basisType_ = BASIS_FEM_FIAT;
68 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
69 this -> basisTagsAreSet_ = false;
70
71 // construct lattice
72
73 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
74
75 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts ,
76 myTri_3 ,
77 n ,
78 0 ,
79 pointType );
80
81
82 // form Vandermonde matrix. Actually, this is the transpose of the VDM,
83 // so we transpose on copy below.
84
85 Phis.getValues( V , latticePts , OPERATOR_VALUE );
86
87 // now I need to copy V into a Teuchos array to do the inversion
88 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
89 for (int i=0;i<N;i++) {
90 for (int j=0;j<N;j++) {
91 Vsdm(i,j) = V(i,j);
92 }
93 }
94
95 // invert the matrix
96 Teuchos::SerialDenseSolver<int,Scalar> solver;
97 solver.setMatrix( rcp( &Vsdm , false ) );
98 solver.invert( );
99
100 // now I need to copy the inverse into Vinv
101 for (int i=0;i<N;i++) {
102 for (int j=0;j<N;j++) {
103 Vinv(i,j) = Vsdm(j,i);
104 }
105 }
106
107 }
108
109
110 template<class Scalar, class ArrayScalar>
111 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::initializeTags() {
112
113 // Basis-dependent initializations
114 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
115 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
116 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
117 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
118
119 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
120
121 int *tags = new int[ tagSize * this->getCardinality() ];
122 int *tag_cur = tags;
123 const int degree = this->getDegree();
124
125 // BEGIN DOF ALONG BOTTOM EDGE
126
127 // the first dof is on vertex 0
128 tag_cur[0] = 0; tag_cur[1] = 0; tag_cur[2] = 0; tag_cur[3] = 1;
129 tag_cur += tagSize;
130
131 // next degree-1 dof are on edge 0
132 for (int i=1;i<degree;i++) {
133 tag_cur[0] = 1; tag_cur[1] = 0; tag_cur[2] = i-1; tag_cur[3] = degree-1;
134 tag_cur += tagSize;
135 }
136
137 // last dof is on vertex 1
138 tag_cur[0] = 0; tag_cur[1] = 1; tag_cur[2] = 0; tag_cur[3] = 1;
139 tag_cur += tagSize;
140
141 // END DOF ALONG BOTTOM EDGE
142
143 int num_internal_dof = PointTools::getLatticeSize( this->getBaseCellTopology() ,
144 this->getDegree() ,
145 1 );
146
147 int internal_dof_cur = 0;
148
149 // BEGIN DOF ALONG INTERNAL HORIZONTAL LINES
150 for (int i=1;i<degree;i++) {
151 // first dof along line is on edge #2
152 tag_cur[0] = 1; tag_cur[1] = 2; tag_cur[2] = i-1; tag_cur[3] = degree-1;
153 tag_cur += tagSize;
154
155 // next dof are internal
156 for (int j=1;j<degree-i;j++) {
157 tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = internal_dof_cur; tag_cur[3] = num_internal_dof;
158 internal_dof_cur++;
159 tag_cur += tagSize;
160 }
161
162 // last dof along line is on edge 1
163 tag_cur[0] = 1; tag_cur[1] = 1; tag_cur[2] = i-1; tag_cur[3] = degree-1;
164 tag_cur += tagSize;
165
166 }
167 // END DOF ALONG INTERNAL HORIZONTAL LINES
168
169 // LAST DOF IS ON VERTEX 2
170 tag_cur[0] = 0; tag_cur[1] = 2; tag_cur[2] = 0; tag_cur[3] = 1;
171 // END LAST DOF
172
173
174 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
175 this -> ordinalToTag_,
176 tags,
177 this -> basisCardinality_,
178 tagSize,
179 posScDim,
180 posScOrd,
181 posDfOrd);
182
183 delete []tags;
184
185 }
186
187
188
189 template<class Scalar, class ArrayScalar>
190 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
191 const ArrayScalar & inputPoints,
192 const EOperator operatorType) const {
193
194 // Verify arguments
195#ifdef HAVE_INTREPID_DEBUG
196 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
197 inputPoints,
198 operatorType,
199 this -> getBaseCellTopology(),
200 this -> getCardinality() );
201#endif
202 const int numPts = inputPoints.dimension(0);
203 const int numBf = this->getCardinality();
204
205 try {
206 switch (operatorType) {
207 case OPERATOR_VALUE:
208 {
209 FieldContainer<Scalar> phisCur( numBf , numPts );
210 Phis.getValues( phisCur , inputPoints , operatorType );
211 for (int i=0;i<outputValues.dimension(0);i++) {
212 for (int j=0;j<outputValues.dimension(1);j++) {
213 outputValues(i,j) = 0.0;
214 for (int k=0;k<this->getCardinality();k++) {
215 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j);
216 }
217 }
218 }
219 }
220 break;
221 case OPERATOR_GRAD:
222 case OPERATOR_D1:
223 case OPERATOR_D2:
224 case OPERATOR_D3:
225 case OPERATOR_D4:
226 case OPERATOR_D5:
227 case OPERATOR_D6:
228 case OPERATOR_D7:
229 case OPERATOR_D8:
230 case OPERATOR_D9:
231 case OPERATOR_D10:
232 {
233 const int dkcard =
234 (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,2): getDkCardinality(operatorType,2);
235
236 FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
237 Phis.getValues( phisCur , inputPoints , operatorType );
238
239 for (int i=0;i<outputValues.dimension(0);i++) {
240 for (int j=0;j<outputValues.dimension(1);j++) {
241 for (int k=0;k<outputValues.dimension(2);k++) {
242 outputValues(i,j,k) = 0.0;
243 for (int l=0;l<this->getCardinality();l++) {
244 outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k);
245 }
246 }
247 }
248 }
249 }
250 break;
251 case OPERATOR_CURL: // only works in 2d. first component is -d/dy, second is d/dx
252 {
253 FieldContainer<Scalar> phisCur( numBf , numPts , getDkCardinality( OPERATOR_D1 , 2 ) );
254 Phis.getValues( phisCur , inputPoints , OPERATOR_D1 );
255
256 for (int i=0;i<outputValues.dimension(0);i++) {
257 for (int j=0;j<outputValues.dimension(1);j++) {
258 outputValues(i,j,0) = 0.0;
259 outputValues(i,j,1) = 0.0;
260 for (int k=0;k<this->getCardinality();k++) {
261 outputValues(i,j,0) += this->Vinv(k,i) * phisCur(k,j,1);
262 }
263 for (int k=0;k<this->getCardinality();k++) {
264 outputValues(i,j,1) -= this->Vinv(k,i) * phisCur(k,j,0);
265 }
266 }
267 }
268 }
269 break;
270 default:
271 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
272 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");
273 break;
274 }
275 }
276 catch (std::invalid_argument &exception){
277 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
278 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented");
279 }
280
281 }
282
283
284
285 template<class Scalar, class ArrayScalar>
286 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
287 const ArrayScalar & inputPoints,
288 const ArrayScalar & cellVertices,
289 const EOperator operatorType) const {
290 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
291 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function");
292 }
293
294
295}// namespace Intrepid
296#endif
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
Basis_HGRAD_TRI_Cn_FEM(const int n, const EPointType pointType)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...