Intrepid2
Intrepid2_HGRAD_TET_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_HGRAD_TET_CN_FEM_ORTH_DEF_HPP__
17#define __INTREPID2_HGRAD_TET_CN_FEM_ORTH_DEF_HPP__
18
19namespace Intrepid2 {
20// -------------------------------------------------------------------------------------
21
22namespace Impl {
23
24template<typename OutputViewType,
25typename inputViewType,
26typename workViewType,
27bool hasDeriv>
28KOKKOS_INLINE_FUNCTION
29void OrthPolynomialTet<OutputViewType,inputViewType,workViewType,hasDeriv,0>::generate(
30 OutputViewType output,
31 const inputViewType input,
32 workViewType /*work*/,
33 const ordinal_type order ) {
34
35 constexpr ordinal_type spaceDim = 3;
36 constexpr ordinal_type maxNumPts = Parameters::MaxNumPtsPerBasisEval;
37
38 typedef typename OutputViewType::value_type value_type;
39
40 auto output0 = (hasDeriv) ? Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL(),0) : Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL());
41
42 const ordinal_type
43 npts = input.extent(0);
44
45 const auto z = input;
46
47 // each point needs to be transformed from Pavel's element
48 // z(i,0) --> (2.0 * z(i,0) - 1.0)
49 // z(i,1) --> (2.0 * z(i,1) - 1.0)
50
51 // set D^{0,0,0} = 1.0
52 {
53 const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(0,0,0);
54 for (ordinal_type i=0;i<npts;++i) {
55 output0(loc, i) = 1.0;
56 if(hasDeriv) {
57 output.access(loc,i,1) = 0;
58 output.access(loc,i,2) = 0;
59 output.access(loc,i,3) = 0;
60 }
61 }
62 }
63
64 if (order > 0) {
65 value_type f1[maxNumPts]={},f2[maxNumPts]={}, f3[maxNumPts]={}, f4[maxNumPts]={},f5[maxNumPts]={},
66 df2_1[maxNumPts]={}, df2_2[maxNumPts]={}, df5_2[maxNumPts]={};
67 value_type df1_0, df1_1, df1_2, df3_1, df3_2, df4_2;
68
69 for (int i=0;i<npts;i++) {
70 f1[i] = 2.0*z(i,0) + z(i,1) + z(i,2)- 1.0;
71 f2[i] = pow(z(i,1) + z(i,2) -1.0, 2);
72 f3[i] = 2.0*z(i,1) + z(i,2) -1.0;
73 f4[i] = 1.0 - z(i,2);
74 f5[i] = f4[i] * f4[i];
75 if(hasDeriv) {
76 df2_1[i] = 2*(z(i,1) + z(i,2) -1.0);
77 df2_2[i] = df2_1[i];
78 df5_2[i] = -2*f4[i];
79 }
80 }
81
82 df1_0 = 2.0;
83 df1_1 = 1.0;
84 df1_2 = 1.0;
85 df3_1 = 2;
86 df3_2 = 1;
87 df4_2 = -1;
88
89 // set D^{1,0,0} = f1
90 {
91 const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(1,0,0);
92 for (ordinal_type i=0;i<npts;++i) {
93 output0(loc, i) = f1[i];
94 if(hasDeriv) {
95 output.access(loc,i,1) = df1_0;
96 output.access(loc,i,2) = df1_1;
97 output.access(loc,i,3) = df1_2;
98 }
99 }
100 }
101
102 // recurrence in p
103 // D^{p+1,0,0}
104 for (ordinal_type p=1;p<order;p++) {
105 const ordinal_type
106 loc = Intrepid2::getPnEnumeration<spaceDim>(p,0,0),
107 loc_p1 = Intrepid2::getPnEnumeration<spaceDim>(p+1,0,0),
108 loc_m1 = Intrepid2::getPnEnumeration<spaceDim>(p-1,0,0);
109
110 const value_type
111 a = (2.0*p+1.0)/(1.0+p),
112 b = p / (p+1.0);
113
114 for (ordinal_type i=0;i<npts;++i) {
115 output0(loc_p1,i) = ( a * f1[i] * output0(loc,i) -
116 b * f2[i] * output0(loc_m1,i) );
117 if(hasDeriv) {
118 output.access(loc_p1,i,1) = a * (f1[i] * output.access(loc,i,1) + df1_0 * output0(loc,i)) -
119 b * f2[i] * output.access(loc_m1,i,1) ;
120 output.access(loc_p1,i,2) = a * (f1[i] * output.access(loc,i,2) + df1_1 * output0(loc,i)) -
121 b * (df2_1[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,2)) ;
122 output.access(loc_p1,i,3) = a * (f1[i] * output.access(loc,i,3) + df1_2 * output0(loc,i)) -
123 b * (df2_2[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,3)) ;
124 }
125 }
126 }
127
128 // D^{p,1,0}
129 for (ordinal_type p=0;p<order;++p) {
130 const ordinal_type
131 loc_p_0 = Intrepid2::getPnEnumeration<spaceDim>(p,0,0),
132 loc_p_1 = Intrepid2::getPnEnumeration<spaceDim>(p,1,0);
133
134 for (ordinal_type i=0;i<npts;++i) {
135 const value_type coeff = (2.0 * p + 3.0) * z(i,1) + z(i,2)-1.0;
136 output0(loc_p_1,i) = output0(loc_p_0,i)*coeff;
137 if(hasDeriv) {
138 output.access(loc_p_1,i,1) = output.access(loc_p_0,i,1)*coeff;
139 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);
140 output.access(loc_p_1,i,3) = output.access(loc_p_0,i,3)*coeff + output0(loc_p_0,i);
141 }
142 }
143 }
144
145
146 // recurrence in q
147 // D^{p,q+1,0}
148 for (ordinal_type p=0;p<order-1;++p)
149 for (ordinal_type q=1;q<order-p;++q) {
150 const ordinal_type
151 loc_p_qp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q+1,0),
152 loc_p_q = Intrepid2::getPnEnumeration<spaceDim>(p,q,0),
153 loc_p_qm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q-1,0);
154
155 value_type a,b,c, coeff, dcoeff_1, dcoeff_2;
156 Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+1,0,q);
157 for (ordinal_type i=0;i<npts;++i) {
158 coeff = a * f3[i] + b * f4[i];
159 dcoeff_1 = a * df3_1;
160 dcoeff_2 = a * df3_2 + b * df4_2;
161 output0(loc_p_qp1,i) = coeff * output0(loc_p_q,i)
162 - c* f5[i] * output0(loc_p_qm1,i) ;
163 if(hasDeriv) {
164 output.access(loc_p_qp1,i,1) = coeff * output.access(loc_p_q,i,1) +
165 - c * f5[i] * output.access(loc_p_qm1,i,1) ;
166 output.access(loc_p_qp1,i,2) = coeff * output.access(loc_p_q,i,2) + dcoeff_1 * output0(loc_p_q,i)
167 - c * f5[i] * output.access(loc_p_qm1,i,2) ;
168 output.access(loc_p_qp1,i,3) = coeff * output.access(loc_p_q,i,3) + dcoeff_2 * output0(loc_p_q,i)
169 - c * f5[i] * output.access(loc_p_qm1,i,3) - c * df5_2[i] * output0(loc_p_qm1,i);
170 }
171 }
172 }
173
174
175 // D^{p,q,1}
176 for (ordinal_type p=0;p<order;++p)
177 for (ordinal_type q=0;q<order-p;++q) {
178
179 const ordinal_type
180 loc_p_q_0 = Intrepid2::getPnEnumeration<spaceDim>(p,q,0),
181 loc_p_q_1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,1);
182
183 for (ordinal_type i=0;i<npts;++i) {
184 const value_type coeff = 2.0 * ( 2.0 + q + p ) * z(i,2) - 1.0;
185 output0(loc_p_q_1,i) = output0(loc_p_q_0,i) * coeff;
186 if(hasDeriv) {
187 output.access(loc_p_q_1,i,1) = output.access(loc_p_q_0,i,1) * coeff;
188 output.access(loc_p_q_1,i,2) = output.access(loc_p_q_0,i,2) * coeff;
189 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);
190 }
191 }
192 }
193
194 // general r recurrence
195 // D^{p,q,r+1}
196 for (ordinal_type p=0;p<order-1;++p)
197 for (ordinal_type q=0;q<order-p-1;++q)
198 for (ordinal_type r=1;r<order-p-q;++r) {
199 const ordinal_type
200 loc_p_q_rp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,r+1),
201 loc_p_q_r = Intrepid2::getPnEnumeration<spaceDim>(p,q,r),
202 loc_p_q_rm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q,r-1);
203
204 value_type a,b,c, coeff;
205 Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+2*q+2,0,r);
206 for (ordinal_type i=0;i<npts;++i) {
207 coeff = 2.0 * a * z(i,2) - a + b;
208 output0(loc_p_q_rp1,i) = coeff * output0(loc_p_q_r,i) - c * output0(loc_p_q_rm1,i) ;
209 if(hasDeriv) {
210 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);
211 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);
212 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);
213 }
214 }
215 }
216
217 }
218
219 // orthogonalize
220 for (ordinal_type p=0;p<=order;++p)
221 for (ordinal_type q=0;q<=order-p;++q)
222 for (ordinal_type r=0;r<=order-p-q;++r) {
223 ordinal_type loc_p_q_r = Intrepid2::getPnEnumeration<spaceDim>(p,q,r);
224 value_type scal = std::sqrt( (p+0.5)*(p+q+1.0)*(p+q+r+1.5) );
225 for (ordinal_type i=0;i<npts;++i) {
226 output0(loc_p_q_r,i) *= scal;
227 if(hasDeriv) {
228 output.access(loc_p_q_r,i,1) *= scal;
229 output.access(loc_p_q_r,i,2) *= scal;
230 output.access(loc_p_q_r,i,3) *= scal;
231 }
232 }
233 }
234}
235
236template<typename OutputViewType,
237typename inputViewType,
238typename workViewType,
239bool hasDeriv>
240KOKKOS_INLINE_FUNCTION
241void OrthPolynomialTet<OutputViewType,inputViewType,workViewType,hasDeriv,1>::generate(
242 OutputViewType output,
243 const inputViewType input,
244 workViewType work,
245 const ordinal_type order ) {
246 constexpr ordinal_type spaceDim = 3;
247 const ordinal_type
248 npts = input.extent(0),
249 card = output.extent(0);
250
251 workViewType dummyView;
252 OrthPolynomialTet<workViewType,inputViewType,workViewType,hasDeriv,0>::generate(work, input, dummyView, order);
253 for (ordinal_type i=0;i<card;++i)
254 for (ordinal_type j=0;j<npts;++j)
255 for (ordinal_type k=0;k<spaceDim;++k)
256 output.access(i,j,k) = work(i,j,k+1);
257}
258
259
260// when n >= 2, use recursion
261template<typename OutputViewType,
262typename inputViewType,
263typename workViewType,
264bool hasDeriv,
265ordinal_type n>
266KOKKOS_INLINE_FUNCTION
267void OrthPolynomialTet<OutputViewType,inputViewType,workViewType,hasDeriv,n>::generate(
268 OutputViewType /* output */,
269 const inputViewType /* input */,
270 workViewType /* work */,
271 const ordinal_type /* order */ ) {
272#if 0 //#ifdef HAVE_INTREPID2_SACADO
273
274constexpr ordinal_type spaceDim = 3;
275constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
276
277typedef typename OutputViewType::value_type value_type;
278typedef Sacado::Fad::SFad<value_type,spaceDim> fad_type;
279
280const ordinal_type
281npts = input.extent(0),
282card = output.extent(0);
283
284// use stack buffer
285fad_type inBuf[Parameters::MaxNumPtsPerBasisEval][spaceDim];
286fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][n*(n+1)/2];
287
288typedef typename inputViewType::memory_space memory_space;
289typedef typename inputViewType::memory_space memory_space;
290typedef typename Kokkos::View<fad_type***, memory_space> outViewType;
291typedef typename Kokkos::View<fad_type**, memory_space> inViewType;
292auto vcprop = Kokkos::common_view_alloc_prop(input);
293
294inViewType in(Kokkos::view_wrap((value_type*)&inBuf[0][0], vcprop), npts, spaceDim);
295outViewType out(Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, n*(n+1)/2);
296
297for (ordinal_type i=0;i<npts;++i)
298 for (ordinal_type j=0;j<spaceDim;++j) {
299 in(i,j) = input(i,j);
300 in(i,j).diff(j,spaceDim);
301 }
302
303typedef typename Kokkos::DynRankView<fad_type, memory_space> outViewType_;
304outViewType_ workView;
305if (n==2) {
306 //char outBuf[bufSize*sizeof(typename inViewType::value_type)];
307 fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][spaceDim+1];
308 auto vcprop = Kokkos::common_view_alloc_prop(in);
309 workView = outViewType_( Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, spaceDim+1);
310}
311OrthPolynomialTet<outViewType,inViewType,outViewType_,hasDeriv,n-1>::generate(out, in, workView, order);
312
313for (ordinal_type i=0;i<card;++i)
314 for (ordinal_type j=0;j<npts;++j) {
315 for (ordinal_type i_dx = 0; i_dx <= n; ++i_dx)
316 for (ordinal_type i_dy = 0; i_dy <= n-i_dx; ++i_dy) {
317 ordinal_type i_dz = n - i_dx - i_dy;
318 ordinal_type i_Dn = Intrepid2::getDkEnumeration<spaceDim>(i_dx,i_dy,i_dz);
319 if(i_dx > 0) {
320 //n=2: (f_x)_x, (f_y)_x, (f_z)_x
321 //n=3: (f_xx)_x, (f_xy)_x, (f_xz)_x, (f_yy)_x, (f_yz)_x, (f_zz)_x
322 ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx-1, i_dy, i_dz);
323 output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(0);
324 }
325 else if (i_dy > 0) {
326 //n=2: (f_y)_y, (f_z)_y
327 //n=3: (f_yy)_y, (f_yz)_y, (f_zz)_y
328 ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx, i_dy-1, i_dz);
329 output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(1);
330 }
331 else {
332 //n=2: (f_z)_z;
333 //n=3: (f_zz)_z
334 ordinal_type i_Dnm1 = Intrepid2::getDkEnumeration<spaceDim>(i_dx, i_dy, i_dz-1);
335 output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(2);
336 }
337 }
338 }
339#else
340INTREPID2_TEST_FOR_ABORT( true,
341 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
342#endif
343}
344
345
346template<EOperator opType>
347template<typename OutputViewType,
348typename inputViewType,
349typename workViewType>
350KOKKOS_INLINE_FUNCTION
351void
352Basis_HGRAD_TET_Cn_FEM_ORTH::Serial<opType>::
353getValues( OutputViewType output,
354 const inputViewType input,
355 workViewType work,
356 const ordinal_type order) {
357 switch (opType) {
358 case OPERATOR_VALUE: {
359 OrthPolynomialTet<OutputViewType,inputViewType,workViewType,false,0>::generate( output, input, work, order );
360 break;
361 }
362 case OPERATOR_GRAD:
363 case OPERATOR_D1: {
364 OrthPolynomialTet<OutputViewType,inputViewType,workViewType,true,1>::generate( output, input, work, order );
365 break;
366 }
367 case OPERATOR_D2: {
368 OrthPolynomialTet<OutputViewType,inputViewType,workViewType,true,2>::generate( output, input, work, order );
369 break;
370 }
371 default: {
372 INTREPID2_TEST_FOR_ABORT( true,
373 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::Serial::getValues) operator is not supported");
374 }
375 }
376}
377
378// -------------------------------------------------------------------------------------
379
380template<typename DT, ordinal_type numPtsPerEval,
381typename outputValueValueType, class ...outputValueProperties,
382typename inputPointValueType, class ...inputPointProperties>
383void
384Basis_HGRAD_TET_Cn_FEM_ORTH::
385getValues(
386 const typename DT::execution_space& space,
387 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
388 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
389 const ordinal_type order,
390 const EOperator operatorType ) {
391 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
392 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
393 typedef typename DT::execution_space ExecSpaceType;
394
395 // loopSize corresponds to the # of points
396 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
397 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
398 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
399 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
400
401 typedef typename inputPointViewType::value_type inputPointType;
402 const ordinal_type cardinality = outputValues.extent(0);
403 const ordinal_type spaceDim = 3;
404
405 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
406 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
407
408 switch (operatorType) {
409 case OPERATOR_VALUE: {
410 workViewType dummyWorkView;
411 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
412 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, dummyWorkView, order) );
413 break;
414 }
415 case OPERATOR_GRAD:
416 case OPERATOR_D1: {
417 workViewType work(Kokkos::view_alloc(space, "Basis_HGRAD_TET_In_FEM_ORTH::getValues::work", vcprop), cardinality, inputPoints.extent(0), spaceDim+1);
418 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D1,numPtsPerEval> FunctorType;
419 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, work, order) );
420 break;
421 }
422 case OPERATOR_D2:{
423 workViewType dummyWorkView;
424 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D2,numPtsPerEval> FunctorType;
425 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints ,dummyWorkView, order) );
426 break;
427 }
428 default: {
429 INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
430 ">>> ERROR (Basis_HGRAD_TET_Cn_FEM_ORTH): invalid operator type");
431 }
432 }
433}
434}
435
436
437// -------------------------------------------------------------------------------------
438template<typename DT, typename OT, typename PT>
440Basis_HGRAD_TET_Cn_FEM_ORTH( const ordinal_type order ) {
441
442 constexpr ordinal_type spaceDim = 3;
443 this->basisCardinality_ = Intrepid2::getPnCardinality<spaceDim>(order);
444 this->basisDegree_ = order;
445 this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
446 this->basisType_ = BASIS_FEM_HIERARCHICAL;
447 this->basisCoordinates_ = COORDINATES_CARTESIAN;
448 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
449
450 // initialize tags
451 {
452 // Basis-dependent initializations
453 constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
454 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
455 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
456 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
457
458 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
459 ordinal_type tags[maxCard][tagSize];
460 const ordinal_type card = this->basisCardinality_;
461 for (ordinal_type i=0;i<card;++i) {
462 tags[i][0] = 2; // these are all "internal" i.e. "volume" DoFs
463 tags[i][1] = 0; // there is only one line
464 tags[i][2] = i; // local DoF id
465 tags[i][3] = card; // total number of DoFs
466 }
467
468 OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
469
470 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
471 // tags are constructed on host
472 this->setOrdinalTagData(this->tagToOrdinal_,
473 this->ordinalToTag_,
474 tagView,
475 this->basisCardinality_,
476 tagSize,
477 posScDim,
478 posScOrd,
479 posDfOrd);
480 }
481
482 // dof coords is not applicable to hierarchical functions
483}
484
485}
486#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.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.