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 
52 namespace 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 );
188  int idxp = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,0,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++) {
197  int idx0 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,0,0);
198  int idx1 = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,1,0);
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);
212  int idxpq = TabulatorTet<Scalar,ArrayScalar,0>::idx(p,q,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  }
281  FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np);
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
304 
305 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
306 #ifdef __GNUC__
307 #warning "The Intrepid package is deprecated"
308 #endif
309 #endif
310 
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...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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...