Intrepid
Intrepid_HGRAD_TET_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
1#ifndef INTREPID_HGRAD_TET_CN_FEM_ORTHDEF_HPP
2#define INTREPID_HGRAD_TET_CN_FEM_ORTHDEF_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
52namespace Intrepid {
53
54 template<class Scalar, class ArrayScalar>
56 {
57 this -> basisCardinality_ = (degree+1)*(degree+2)*(degree+3)/6;
58 this -> basisDegree_ = degree;
59 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
60 this -> basisType_ = BASIS_FEM_HIERARCHICAL;
61 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62 this -> basisTagsAreSet_ = false;
63 }
64
65
66
67 template<class Scalar, class ArrayScalar>
69
70 // Basis-dependent initializations
71 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
72 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
73 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
74 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
75
76 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
77 int *tags = new int[tagSize * this->getCardinality()];
78 for (int i=0;i<this->getCardinality();i++) {
79 tags[4*i] = 2;
80 tags[4*i+1] = 0;
81 tags[4*i+2] = i;
82 tags[4*i+3] = this->getCardinality();
83 }
84
85 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
86 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
87 this -> ordinalToTag_,
88 tags,
89 this -> basisCardinality_,
90 tagSize,
91 posScDim,
92 posScOrd,
93 posDfOrd);
94 }
95
96
97
98 template<class Scalar, class ArrayScalar>
100 const ArrayScalar & inputPoints,
101 const EOperator operatorType) const {
102
103 // Verify arguments
104#ifdef HAVE_INTREPID_DEBUG
105 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
106 inputPoints,
107 operatorType,
108 this -> getBaseCellTopology(),
109 this -> getCardinality() );
110#endif
111 const int deg = this->getDegree();
112
113 switch (operatorType) {
114 case OPERATOR_VALUE:
115 {
117 deg ,
118 inputPoints );
119 }
120 break;
121 case OPERATOR_GRAD:
122 case OPERATOR_D1:
123 {
125 deg ,
126 inputPoints );
127 }
128 break;
129 default:
130 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
131 ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): invalid or unsupported operator" );
132 }
133
134 return;
135 }
136
137 template<class Scalar, class ArrayScalar>
139 const ArrayScalar & inputPoints,
140 const ArrayScalar & cellVertices,
141 const EOperator operatorType) const {
142 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
143 ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): FEM Basis calling an FVD member function");
144 }
145
146 template<class Scalar, class ArrayScalar>
147 void TabulatorTet<Scalar,ArrayScalar,0>::tabulate( ArrayScalar &outputValues ,
148 const int deg ,
149 const ArrayScalar &z )
150 {
151 const int np = z.dimension( 0 );
152 int idxcur;
153
154 // each point needs to be transformed from Pavel's element
155 // z(i,0) --> (2.0 * z(i,0) - 1.0)
156 // z(i,1) --> (2.0 * z(i,1) - 1.0)
157 // z(i,2) --> (2.0 * z(i,2) - 1.0)
158
159 Teuchos::Array<Scalar> f1(np),f2(np),f3(np),f4(np),f5(np);
160
161 for (int i=0;i<np;i++) {
162 f1[i] = 0.5 * ( 2.0 + 2.0*(2.0*z(i,0)-1.0) + (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
163 Scalar foo = 0.5 * ( (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
164 f2[i] = foo * foo;
165 f3[i] = 0.5 * ( 1.0 + 2.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
166 f4[i] = 0.5 * ( 1.0 - (2.0*z(i,2)-1.0) );
167 f5[i] = f4[i] * f4[i];
168 }
169
170 // constant term
172 for (int i=0;i<np;i++) {
173 outputValues(idxcur,i) = 1.0 + z(i,0) - z(i,0) + z(i,1) - z(i,1) + z(i,2) - z(i,2);
174 }
175
176 if (deg > 0) {
177
178 // D^{1,0,0}
180 for (int i=0;i<np;i++) {
181 outputValues(idxcur,i) = f1[i];
182 }
183
184 // p recurrence
185 for (int p=1;p<deg;p++) {
186 Scalar a1 = (2.0 * p + 1.0) / ( p + 1.0);
187 Scalar a2 = p / ( p + 1.0 );
189 int idxpp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p+1,0,0);
190 int idxpm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p-1,0,0);
191 for (int i=0;i<np;i++) {
192 outputValues(idxpp1,i) = a1 * f1[i] * outputValues(idxp,i) - a2 * f2[i] * outputValues(idxpm1,i);
193 }
194 }
195 // q = 1
196 for (int p=0;p<deg;p++) {
199 for (int i=0;i<np;i++) {
200 outputValues(idx1,i) = outputValues(idx0,i) * ( p * ( 1.0 + (2.0*z(i,1)-1.0) ) +
201 0.5 * ( 2.0 + 3.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) );
202 }
203 }
204
205 // q recurrence
206 for (int p=0;p<deg-1;p++) {
207 for (int q=1;q<deg-p;q++) {
208 Scalar aq,bq,cq;
209
210 TabulatorTet<Scalar,ArrayScalar,0>::jrc(2.0*p+1.0 ,0 ,q, aq, bq, cq);
211 int idxpqp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q+1,0);
213 int idxpqm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q-1,0);
214 for (int i=0;i<np;i++) {
215 outputValues(idxpqp1,i) = ( aq * f3[i] + bq * f4[i] ) * outputValues(idxpq,i)
216 - ( cq * f5[i] ) * outputValues(idxpqm1,i);
217 }
218 }
219 }
220
221 // r = 1
222 for (int p=0;p<deg;p++) {
223 for (int q=0;q<deg-p;q++) {
224 int idxpq1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,1);
225 int idxpq0 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,0);
226 for (int i=0;i<np;i++) {
227 outputValues(idxpq1,i) = outputValues(idxpq0,i) * ( 1.0 + p + q + ( 2.0 + q +
228 p ) * (2.0*z(i,2)-1.0) );
229 }
230 }
231 }
232 // general r recurrence
233 for (int p=0;p<deg-1;p++) {
234 for (int q=0;q<deg-p-1;q++) {
235 for (int r=1;r<deg-p-q;r++) {
236 Scalar ar,br,cr;
237 int idxpqrp1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r+1);
238 int idxpqr = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r);
239 int idxpqrm1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,r-1);
240 jrc(2.0*p+2.0*q+2.0, 0.0, r, ar, br, cr);
241 for (int i=0;i<np;i++) {
242 outputValues(idxpqrp1,i) = (ar * (2.0*z(i,2)-1.0) + br) * outputValues( idxpqr , i ) - cr * outputValues(idxpqrm1,i);
243 }
244 }
245 }
246 }
247
248 }
249 // normalize
250 for (int p=0;p<=deg;p++) {
251 for (int q=0;q<=deg-p;q++) {
252 for (int r=0;r<=deg-p-q;r++) {
254 Scalar scal = sqrt( (p+0.5)*(p+q+1.0)*(p+q+r+1.5) );
255 for (int i=0;i<np;i++) {
256 outputValues(idxcur,i) *= scal;
257 }
258 }
259 }
260 }
261
262 return;
263
264 }
265
266
267 template<typename Scalar, typename ArrayScalar>
268 void TabulatorTet<Scalar,ArrayScalar,1>::tabulate( ArrayScalar &outputValues ,
269 const int deg ,
270 const ArrayScalar &z )
271 {
272 const int np = z.dimension(0);
273 const int card = outputValues.dimension(0);
274 FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
275 for (int i=0;i<np;i++) {
276 for (int j=0;j<3;j++) {
277 dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
278 dZ(i,j).diff(j,3);
279 }
280 }
282
284 deg ,
285 dZ );
286
287 for (int i=0;i<card;i++) {
288 for (int j=0;j<np;j++) {
289 for (int k=0;k<3;k++) {
290 outputValues(i,j,k) = dResult(i,j).dx(k);
291 }
292 }
293 }
294
295 return;
296
297 }
298
299
300}// namespace Intrepid
301
302
303#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.
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 Tetrahedron cell.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
This is an internal class with a static member function for tabulating derivatives of orthogonal expa...
static void tabulate(ArrayScalar &outputValues, const int deg, const ArrayScalar &inputPoints)
basic tabulate mathod evaluates the derivOrder^th derivatives of the basis functions at inputPoints i...