Intrepid2
Intrepid2_HGRAD_TET_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
49 #ifndef __INTREPID2_HGRAD_TET_CN_FEM_ORTH_DEF_HPP__
50 #define __INTREPID2_HGRAD_TET_CN_FEM_ORTH_DEF_HPP__
51 
52 namespace Intrepid2 {
53 // -------------------------------------------------------------------------------------
54 
55 namespace Impl {
56 
57 template<typename outputViewType,
58 typename inputViewType,
59 typename workViewType,
60 bool hasDeriv>
61 KOKKOS_INLINE_FUNCTION
62 void OrthPolynomialTet<outputViewType,inputViewType,workViewType,hasDeriv,0>::generate(
63  outputViewType output,
64  const inputViewType input,
65  workViewType /*work*/,
66  const ordinal_type order ) {
67 
68  constexpr ordinal_type spaceDim = 3;
69  constexpr ordinal_type maxNumPts = Parameters::MaxNumPtsPerBasisEval;
70 
71  typedef typename outputViewType::value_type value_type;
72 
73  auto output0 = (hasDeriv) ? Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL(),0) : Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL());
74 
75  const ordinal_type
76  npts = input.extent(0);
77 
78  const auto z = input;
79 
80  // each point needs to be transformed from Pavel's element
81  // z(i,0) --> (2.0 * z(i,0) - 1.0)
82  // z(i,1) --> (2.0 * z(i,1) - 1.0)
83 
84  // set D^{0,0,0} = 1.0
85  {
86  const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(0,0,0);
87  for (ordinal_type i=0;i<npts;++i) {
88  output0(loc, i) = 1.0;
89  if(hasDeriv) {
90  output.access(loc,i,1) = 0;
91  output.access(loc,i,2) = 0;
92  output.access(loc,i,3) = 0;
93  }
94  }
95  }
96 
97  if (order > 0) {
98  value_type f1[maxNumPts]={},f2[maxNumPts]={}, f3[maxNumPts]={}, f4[maxNumPts]={},f5[maxNumPts]={},
99  df2_1[maxNumPts]={}, df2_2[maxNumPts]={}, df5_2[maxNumPts]={};
100  value_type df1_0, df1_1, df1_2, df3_1, df3_2, df4_2;
101 
102  for (int i=0;i<npts;i++) {
103  f1[i] = 2.0*z(i,0) + z(i,1) + z(i,2)- 1.0;
104  f2[i] = pow(z(i,1) + z(i,2) -1.0, 2);
105  f3[i] = 2.0*z(i,1) + z(i,2) -1.0;
106  f4[i] = 1.0 - z(i,2);
107  f5[i] = f4[i] * f4[i];
108  if(hasDeriv) {
109  df2_1[i] = 2*(z(i,1) + z(i,2) -1.0);
110  df2_2[i] = df2_1[i];
111  df5_2[i] = -2*f4[i];
112  }
113  }
114 
115  df1_0 = 2.0;
116  df1_1 = 1.0;
117  df1_2 = 1.0;
118  df3_1 = 2;
119  df3_2 = 1;
120  df4_2 = -1;
121 
122  // set D^{1,0,0} = f1
123  {
124  const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(1,0,0);
125  for (ordinal_type i=0;i<npts;++i) {
126  output0(loc, i) = f1[i];
127  if(hasDeriv) {
128  output.access(loc,i,1) = df1_0;
129  output.access(loc,i,2) = df1_1;
130  output.access(loc,i,3) = df1_2;
131  }
132  }
133  }
134 
135  // recurrence in p
136  // D^{p+1,0,0}
137  for (ordinal_type p=1;p<order;p++) {
138  const ordinal_type
139  loc = Intrepid2::getPnEnumeration<spaceDim>(p,0,0),
140  loc_p1 = Intrepid2::getPnEnumeration<spaceDim>(p+1,0,0),
141  loc_m1 = Intrepid2::getPnEnumeration<spaceDim>(p-1,0,0);
142 
143  const value_type
144  a = (2.0*p+1.0)/(1.0+p),
145  b = p / (p+1.0);
146 
147  for (ordinal_type i=0;i<npts;++i) {
148  output0(loc_p1,i) = ( a * f1[i] * output0(loc,i) -
149  b * f2[i] * output0(loc_m1,i) );
150  if(hasDeriv) {
151  output.access(loc_p1,i,1) = a * (f1[i] * output.access(loc,i,1) + df1_0 * output0(loc,i)) -
152  b * f2[i] * output.access(loc_m1,i,1) ;
153  output.access(loc_p1,i,2) = a * (f1[i] * output.access(loc,i,2) + df1_1 * output0(loc,i)) -
154  b * (df2_1[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,2)) ;
155  output.access(loc_p1,i,3) = a * (f1[i] * output.access(loc,i,3) + df1_2 * output0(loc,i)) -
156  b * (df2_2[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,3)) ;
157  }
158  }
159  }
160 
161  // D^{p,1,0}
162  for (ordinal_type p=0;p<order;++p) {
163  const ordinal_type
164  loc_p_0 = Intrepid2::getPnEnumeration<spaceDim>(p,0,0),
165  loc_p_1 = Intrepid2::getPnEnumeration<spaceDim>(p,1,0);
166 
167  for (ordinal_type i=0;i<npts;++i) {
168  const value_type coeff = (2.0 * p + 3.0) * z(i,1) + z(i,2)-1.0;
169  output0(loc_p_1,i) = output0(loc_p_0,i)*coeff;
170  if(hasDeriv) {
171  output.access(loc_p_1,i,1) = output.access(loc_p_0,i,1)*coeff;
172  output.access(loc_p_1,i,2) = output.access(loc_p_0,i,2)*coeff + output0(loc_p_0,i)*(2.0 * p + 3.0);
173  output.access(loc_p_1,i,3) = output.access(loc_p_0,i,3)*coeff + output0(loc_p_0,i);
174  }
175  }
176  }
177 
178 
179  // recurrence in q
180  // D^{p,q+1,0}
181  for (ordinal_type p=0;p<order-1;++p)
182  for (ordinal_type q=1;q<order-p;++q) {
183  const ordinal_type
184  loc_p_qp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q+1,0),
185  loc_p_q = Intrepid2::getPnEnumeration<spaceDim>(p,q,0),
186  loc_p_qm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q-1,0);
187 
188  value_type a,b,c, coeff, dcoeff_1, dcoeff_2;
189  Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+1,0,q);
190  for (ordinal_type i=0;i<npts;++i) {
191  coeff = a * f3[i] + b * f4[i];
192  dcoeff_1 = a * df3_1;
193  dcoeff_2 = a * df3_2 + b * df4_2;
194  output0(loc_p_qp1,i) = coeff * output0(loc_p_q,i)
195  - c* f5[i] * output0(loc_p_qm1,i) ;
196  if(hasDeriv) {
197  output.access(loc_p_qp1,i,1) = coeff * output.access(loc_p_q,i,1) +
198  - c * f5[i] * output.access(loc_p_qm1,i,1) ;
199  output.access(loc_p_qp1,i,2) = coeff * output.access(loc_p_q,i,2) + dcoeff_1 * output0(loc_p_q,i)
200  - c * f5[i] * output.access(loc_p_qm1,i,2) ;
201  output.access(loc_p_qp1,i,3) = coeff * output.access(loc_p_q,i,3) + dcoeff_2 * output0(loc_p_q,i)
202  - c * f5[i] * output.access(loc_p_qm1,i,3) - c * df5_2[i] * output0(loc_p_qm1,i);
203  }
204  }
205  }
206 
207 
208  // D^{p,q,1}
209  for (ordinal_type p=0;p<order;++p)
210  for (ordinal_type q=0;q<order-p;++q) {
211 
212  const ordinal_type
213  loc_p_q_0 = Intrepid2::getPnEnumeration<spaceDim>(p,q,0),
214  loc_p_q_1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,1);
215 
216  for (ordinal_type i=0;i<npts;++i) {
217  const value_type coeff = 2.0 * ( 2.0 + q + p ) * z(i,2) - 1.0;
218  output0(loc_p_q_1,i) = output0(loc_p_q_0,i) * coeff;
219  if(hasDeriv) {
220  output.access(loc_p_q_1,i,1) = output.access(loc_p_q_0,i,1) * coeff;
221  output.access(loc_p_q_1,i,2) = output.access(loc_p_q_0,i,2) * coeff;
222  output.access(loc_p_q_1,i,3) = output.access(loc_p_q_0,i,3) * coeff + 2.0 * ( 2.0 + q + p ) * output0(loc_p_q_0,i);
223  }
224  }
225  }
226 
227  // general r recurrence
228  // D^{p,q,r+1}
229  for (ordinal_type p=0;p<order-1;++p)
230  for (ordinal_type q=0;q<order-p-1;++q)
231  for (ordinal_type r=1;r<order-p-q;++r) {
232  const ordinal_type
233  loc_p_q_rp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,r+1),
234  loc_p_q_r = Intrepid2::getPnEnumeration<spaceDim>(p,q,r),
235  loc_p_q_rm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,r-1);
236 
237  value_type a,b,c, coeff;
238  Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+2*q+2,0,r);
239  for (ordinal_type i=0;i<npts;++i) {
240  coeff = 2.0 * a * z(i,2) - a + b;
241  output0(loc_p_q_rp1,i) = coeff * output0(loc_p_q_r,i) - c * output0(loc_p_q_rm1,i) ;
242  if(hasDeriv) {
243  output.access(loc_p_q_rp1,i,1) = coeff * output.access(loc_p_q_r,i,1) - c * output.access(loc_p_q_rm1,i,1);
244  output.access(loc_p_q_rp1,i,2) = coeff * output.access(loc_p_q_r,i,2) - c * output.access(loc_p_q_rm1,i,2);
245  output.access(loc_p_q_rp1,i,3) = coeff * output.access(loc_p_q_r,i,3) + 2 * a * output0(loc_p_q_r,i) - c * output.access(loc_p_q_rm1,i,3);
246  }
247  }
248  }
249 
250  }
251 
252  // orthogonalize
253  for (ordinal_type p=0;p<=order;++p)
254  for (ordinal_type q=0;q<=order-p;++q)
255  for (ordinal_type r=0;r<=order-p-q;++r) {
256  ordinal_type loc_p_q_r = Intrepid2::getPnEnumeration<spaceDim>(p,q,r);
257  value_type scal = std::sqrt( (p+0.5)*(p+q+1.0)*(p+q+r+1.5) );
258  for (ordinal_type i=0;i<npts;++i) {
259  output0(loc_p_q_r,i) *= scal;
260  if(hasDeriv) {
261  output.access(loc_p_q_r,i,1) *= scal;
262  output.access(loc_p_q_r,i,2) *= scal;
263  output.access(loc_p_q_r,i,3) *= scal;
264  }
265  }
266  }
267 }
268 
269 template<typename outputViewType,
270 typename inputViewType,
271 typename workViewType,
272 bool hasDeriv>
273 KOKKOS_INLINE_FUNCTION
274 void OrthPolynomialTet<outputViewType,inputViewType,workViewType,hasDeriv,1>::generate(
275  outputViewType output,
276  const inputViewType input,
277  workViewType work,
278  const ordinal_type order ) {
279  constexpr ordinal_type spaceDim = 3;
280  const ordinal_type
281  npts = input.extent(0),
282  card = output.extent(0);
283 
284  workViewType dummyView;
285  OrthPolynomialTet<workViewType,inputViewType,workViewType,hasDeriv,0>::generate(work, input, dummyView, order);
286  for (ordinal_type i=0;i<card;++i)
287  for (ordinal_type j=0;j<npts;++j)
288  for (ordinal_type k=0;k<spaceDim;++k)
289  output.access(i,j,k) = work(i,j,k+1);
290 }
291 
292 
293 // when n >= 2, use recursion
294 template<typename outputViewType,
295 typename inputViewType,
296 typename workViewType,
297 bool hasDeriv,
298 ordinal_type n>
299 KOKKOS_INLINE_FUNCTION
300 void OrthPolynomialTet<outputViewType,inputViewType,workViewType,hasDeriv,n>::generate(
301  outputViewType /* output */,
302  const inputViewType /* input */,
303  workViewType /* work */,
304  const ordinal_type /* order */ ) {
305 #if 0 //#ifdef HAVE_INTREPID2_SACADO
306 
307 constexpr ordinal_type spaceDim = 3;
308 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
309 
310 typedef typename outputViewType::value_type value_type;
311 typedef Sacado::Fad::SFad<value_type,spaceDim> fad_type;
312 
313 const ordinal_type
314 npts = input.extent(0),
315 card = output.extent(0);
316 
317 // use stack buffer
318 fad_type inBuf[Parameters::MaxNumPtsPerBasisEval][spaceDim];
319 fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][n*(n+1)/2];
320 
321 typedef typename inputViewType::memory_space memory_space;
322 typedef typename inputViewType::memory_space memory_space;
323 typedef typename Kokkos::View<fad_type***, memory_space> outViewType;
324 typedef typename Kokkos::View<fad_type**, memory_space> inViewType;
325 auto vcprop = Kokkos::common_view_alloc_prop(input);
326 
327 inViewType in(Kokkos::view_wrap((value_type*)&inBuf[0][0], vcprop), npts, spaceDim);
328 outViewType out(Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, n*(n+1)/2);
329 
330 for (ordinal_type i=0;i<npts;++i)
331  for (ordinal_type j=0;j<spaceDim;++j) {
332  in(i,j) = input(i,j);
333  in(i,j).diff(j,spaceDim);
334  }
335 
336 typedef typename Kokkos::DynRankView<fad_type, memory_space> outViewType_;
337 outViewType_ workView;
338 if (n==2) {
339  //char outBuf[bufSize*sizeof(typename inViewType::value_type)];
340  fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][spaceDim+1];
341  auto vcprop = Kokkos::common_view_alloc_prop(in);
342  workView = outViewType_( Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, spaceDim+1);
343 }
344 OrthPolynomialTet<outViewType,inViewType,outViewType_,hasDeriv,n-1>::generate(out, in, workView, order);
345 
346 for (ordinal_type i=0;i<card;++i)
347  for (ordinal_type j=0;j<npts;++j) {
348  for (ordinal_type i_dx = 0; i_dx <= n; ++i_dx)
349  for (ordinal_type i_dy = 0; i_dy <= n-i_dx; ++i_dy) {
350  ordinal_type i_dz = n - i_dx - i_dy;
351  ordinal_type i_Dn = Intrepid2::getDkEnumeration<spaceDim>(i_dx,i_dy,i_dz);
352  if(i_dx > 0) {
353  //n=2: (f_x)_x, (f_y)_x, (f_z)_x
354  //n=3: (f_xx)_x, (f_xy)_x, (f_xz)_x, (f_yy)_x, (f_yz)_x, (f_zz)_x
355  ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx-1, i_dy, i_dz);
356  output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(0);
357  }
358  else if (i_dy > 0) {
359  //n=2: (f_y)_y, (f_z)_y
360  //n=3: (f_yy)_y, (f_yz)_y, (f_zz)_y
361  ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx, i_dy-1, i_dz);
362  output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(1);
363  }
364  else {
365  //n=2: (f_z)_z;
366  //n=3: (f_zz)_z
367  ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx, i_dy, i_dz-1);
368  output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(2);
369  }
370  }
371  }
372 #else
373 INTREPID2_TEST_FOR_ABORT( true,
374  ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
375 #endif
376 }
377 
378 
379 template<EOperator opType>
380 template<typename outputViewType,
381 typename inputViewType,
382 typename workViewType>
383 KOKKOS_INLINE_FUNCTION
384 void
385 Basis_HGRAD_TET_Cn_FEM_ORTH::Serial<opType>::
386 getValues( outputViewType output,
387  const inputViewType input,
388  workViewType work,
389  const ordinal_type order) {
390  switch (opType) {
391  case OPERATOR_VALUE: {
392  OrthPolynomialTet<outputViewType,inputViewType,workViewType,false,0>::generate( output, input, work, order );
393  break;
394  }
395  case OPERATOR_GRAD:
396  case OPERATOR_D1: {
397  OrthPolynomialTet<outputViewType,inputViewType,workViewType,true,1>::generate( output, input, work, order );
398  break;
399  }
400  case OPERATOR_D2: {
401  OrthPolynomialTet<outputViewType,inputViewType,workViewType,true,2>::generate( output, input, work, order );
402  break;
403  }
404  default: {
405  INTREPID2_TEST_FOR_ABORT( true,
406  ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::Serial::getValues) operator is not supported");
407  }
408  }
409 }
410 
411 // -------------------------------------------------------------------------------------
412 
413 template<typename SpT, ordinal_type numPtsPerEval,
414 typename outputValueValueType, class ...outputValueProperties,
415 typename inputPointValueType, class ...inputPointProperties>
416 void
417 Basis_HGRAD_TET_Cn_FEM_ORTH::
418 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
419  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
420  const ordinal_type order,
421  const EOperator operatorType ) {
422  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
423  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
424  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
425 
426  // loopSize corresponds to the # of points
427  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
428  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
429  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
430  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
431 
432  typedef typename inputPointViewType::value_type inputPointType;
433  const ordinal_type cardinality = outputValues.extent(0);
434  const ordinal_type spaceDim = 3;
435 
436  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
437  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
438 
439  switch (operatorType) {
440  case OPERATOR_VALUE: {
441  workViewType dummyWorkView;
442  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
443  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, dummyWorkView, order) );
444  break;
445  }
446  case OPERATOR_GRAD:
447  case OPERATOR_D1: {
448  workViewType work(Kokkos::view_alloc("Basis_HGRAD_TET_In_FEM_ORTH::getValues::work", vcprop), cardinality, inputPoints.extent(0), spaceDim+1);
449  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D1,numPtsPerEval> FunctorType;
450  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, work, order) );
451  break;
452  }
453  case OPERATOR_D2:{
454  workViewType dummyWorkView;
455  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D2,numPtsPerEval> FunctorType;
456  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints ,dummyWorkView, order) );
457  break;
458  }
459  default: {
460  INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
461  ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): invalid operator type");
462  }
463  }
464 }
465 }
466 
467 
468 // -------------------------------------------------------------------------------------
469 template<typename SpT, typename OT, typename PT>
470 Basis_HGRAD_TET_Cn_FEM_ORTH<SpT,OT,PT>::
471 Basis_HGRAD_TET_Cn_FEM_ORTH( const ordinal_type order ) {
472 
473  constexpr ordinal_type spaceDim = 3;
474  this->basisCardinality_ = Intrepid2::getPnCardinality<spaceDim>(order);
475  this->basisDegree_ = order;
476  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
477  this->basisType_ = BASIS_FEM_HIERARCHICAL;
478  this->basisCoordinates_ = COORDINATES_CARTESIAN;
479 
480  // initialize tags
481  {
482  // Basis-dependent initializations
483  constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
484  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
485  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
486  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
487 
488  constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
489  ordinal_type tags[maxCard][tagSize];
490  const ordinal_type card = this->basisCardinality_;
491  for (ordinal_type i=0;i<card;++i) {
492  tags[i][0] = 2; // these are all "internal" i.e. "volume" DoFs
493  tags[i][1] = 0; // there is only one line
494  tags[i][2] = i; // local DoF id
495  tags[i][3] = card; // total number of DoFs
496  }
497 
498  ordinal_type_array_1d_host tagView(&tags[0][0], card*tagSize);
499 
500  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
501  // tags are constructed on host
502  this->setOrdinalTagData(this->tagToOrdinal_,
503  this->ordinalToTag_,
504  tagView,
505  this->basisCardinality_,
506  tagSize,
507  posScDim,
508  posScOrd,
509  posDfOrd);
510  }
511 
512  // dof coords is not applicable to hierarchical functions
513 }
514 
515 }
516 #endif
Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::Basis_HGRAD_TET_Cn_FEM_ORTH
Basis_HGRAD_TET_Cn_FEM_ORTH(const ordinal_type order)
Constructor.
Definition: Intrepid2_HGRAD_TET_Cn_FEM_ORTHDef.hpp:471
Intrepid2::Parameters::MaxNumPtsPerBasisEval
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.
Definition: Intrepid2_Types.hpp:121
Intrepid2::Basis::ordinal_type_array_1d_host
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:113
Generated by   doxygen 1.8.5