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
52namespace Intrepid2 {
53// -------------------------------------------------------------------------------------
54
55namespace Impl {
56
57template<typename OutputViewType,
58typename inputViewType,
59typename workViewType,
60bool hasDeriv>
61KOKKOS_INLINE_FUNCTION
62void 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
269template<typename OutputViewType,
270typename inputViewType,
271typename workViewType,
272bool hasDeriv>
273KOKKOS_INLINE_FUNCTION
274void 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
294template<typename OutputViewType,
295typename inputViewType,
296typename workViewType,
297bool hasDeriv,
298ordinal_type n>
299KOKKOS_INLINE_FUNCTION
300void 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
307constexpr ordinal_type spaceDim = 3;
308constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
309
310typedef typename OutputViewType::value_type value_type;
311typedef Sacado::Fad::SFad<value_type,spaceDim> fad_type;
312
313const ordinal_type
314npts = input.extent(0),
315card = output.extent(0);
316
317// use stack buffer
318fad_type inBuf[Parameters::MaxNumPtsPerBasisEval][spaceDim];
319fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][n*(n+1)/2];
320
321typedef typename inputViewType::memory_space memory_space;
322typedef typename inputViewType::memory_space memory_space;
323typedef typename Kokkos::View<fad_type***, memory_space> outViewType;
324typedef typename Kokkos::View<fad_type**, memory_space> inViewType;
325auto vcprop = Kokkos::common_view_alloc_prop(input);
326
327inViewType in(Kokkos::view_wrap((value_type*)&inBuf[0][0], vcprop), npts, spaceDim);
328outViewType out(Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, n*(n+1)/2);
329
330for (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
336typedef typename Kokkos::DynRankView<fad_type, memory_space> outViewType_;
337outViewType_ workView;
338if (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}
344OrthPolynomialTet<outViewType,inViewType,outViewType_,hasDeriv,n-1>::generate(out, in, workView, order);
345
346for (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
373INTREPID2_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
379template<EOperator opType>
380template<typename OutputViewType,
381typename inputViewType,
382typename workViewType>
383KOKKOS_INLINE_FUNCTION
384void
385Basis_HGRAD_TET_Cn_FEM_ORTH::Serial<opType>::
386getValues( 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
413template<typename DT, ordinal_type numPtsPerEval,
414typename outputValueValueType, class ...outputValueProperties,
415typename inputPointValueType, class ...inputPointProperties>
416void
417Basis_HGRAD_TET_Cn_FEM_ORTH::
418getValues( 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 DT::execution_space 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// -------------------------------------------------------------------------------------
469template<typename DT, typename OT, typename PT>
471Basis_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 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
480
481 // initialize tags
482 {
483 // Basis-dependent initializations
484 constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
485 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
486 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
487 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
488
489 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
490 ordinal_type tags[maxCard][tagSize];
491 const ordinal_type card = this->basisCardinality_;
492 for (ordinal_type i=0;i<card;++i) {
493 tags[i][0] = 2; // these are all "internal" i.e. "volume" DoFs
494 tags[i][1] = 0; // there is only one line
495 tags[i][2] = i; // local DoF id
496 tags[i][3] = card; // total number of DoFs
497 }
498
499 OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
500
501 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
502 // tags are constructed on host
503 this->setOrdinalTagData(this->tagToOrdinal_,
504 this->ordinalToTag_,
505 tagView,
506 this->basisCardinality_,
507 tagSize,
508 posScDim,
509 posScOrd,
510 posDfOrd);
511 }
512
513 // dof coords is not applicable to hierarchical functions
514}
515
516}
517#endif
KOKKOS_INLINE_FUNCTION void getJacobyRecurrenceCoeffs(value_type &an, value_type &bn, value_type &cn, const ordinal_type alpha, const ordinal_type beta, const ordinal_type n)
function for computing the Jacobi recurrence coefficients so that
Basis_HGRAD_TET_Cn_FEM_ORTH(const ordinal_type order)
Constructor.
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.