Intrepid2
Intrepid2_HGRAD_TRI_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
17#ifndef __INTREPID2_HGRAD_TRI_CN_FEM_ORTH_DEF_HPP__
18#define __INTREPID2_HGRAD_TRI_CN_FEM_ORTH_DEF_HPP__
19
20namespace Intrepid2 {
21// -------------------------------------------------------------------------------------
22
23namespace Impl {
24
25template<typename OutputViewType,
26typename inputViewType,
27typename workViewType,
28bool hasDeriv>
29KOKKOS_INLINE_FUNCTION
30void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,0>::generate(
31 OutputViewType output,
32 const inputViewType input,
33 workViewType /*work*/,
34 const ordinal_type order ) {
35
36 constexpr ordinal_type spaceDim = 2;
37 constexpr ordinal_type maxNumPts = Parameters::MaxNumPtsPerBasisEval;
38
39 typedef typename OutputViewType::value_type value_type;
40
41 auto output0 = (hasDeriv) ? Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL(),0) : Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL());
42
43 const ordinal_type
44 npts = input.extent(0);
45
46 const auto z = input;
47
48 // each point needs to be transformed from Pavel's element
49 // z(i,0) --> (2.0 * z(i,0) - 1.0)
50 // z(i,1) --> (2.0 * z(i,1) - 1.0)
51
52
53 // set D^{0,0} = 1.0
54 {
55 const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(0,0);
56 for (ordinal_type i=0;i<npts;++i) {
57 output0(loc, i) = 1.0;
58 if(hasDeriv) {
59 output.access(loc,i,1) = 0;
60 output.access(loc,i,2) = 0;
61 }
62 }
63 }
64
65 if (order > 0) {
66 value_type f1[maxNumPts]={},f2[maxNumPts]={}, df2_1[maxNumPts]={};
67 value_type df1_0, df1_1;
68
69 for (ordinal_type i=0;i<npts;++i) {
70 f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0)); // \eta_1 * (1 - \eta_2)/2
71 f2[i] = pow(z(i,1)-1,2); //( (1 - \eta_2)/2 )^2
72 if(hasDeriv) {
73 df1_0 = 2.0;
74 df1_1 = 1.0;
75 df2_1[i] = 2.0*(z(i,1)-1);
76 }
77 }
78
79 // set D^{1,0} = f1
80 {
81 const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(1,0);
82 for (ordinal_type i=0;i<npts;++i) {
83 output0(loc, i) = f1[i];
84 if(hasDeriv) {
85 output.access(loc,i,1) = df1_0;
86 output.access(loc,i,2) = df1_1;
87 }
88 }
89 }
90
91 // recurrence in p
92 for (ordinal_type p=1;p<order;p++) {
93 const ordinal_type
94 loc = Intrepid2::getPnEnumeration<spaceDim>(p,0),
95 loc_p1 = Intrepid2::getPnEnumeration<spaceDim>(p+1,0),
96 loc_m1 = Intrepid2::getPnEnumeration<spaceDim>(p-1,0);
97
98 const value_type
99 a = (2.0*p+1.0)/(1.0+p),
100 b = p / (p+1.0);
101
102 for (ordinal_type i=0;i<npts;++i) {
103 output0(loc_p1,i) = ( a * f1[i] * output0(loc,i) -
104 b * f2[i] * output0(loc_m1,i) );
105 if(hasDeriv) {
106 output.access(loc_p1,i,1) = a * (f1[i] * output.access(loc,i,1) + df1_0 * output0(loc,i)) -
107 b * f2[i] * output.access(loc_m1,i,1) ;
108 output.access(loc_p1,i,2) = a * (f1[i] * output.access(loc,i,2) + df1_1 * output0(loc,i)) -
109 b * (df2_1[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,2)) ;
110 }
111 }
112 }
113
114 // D^{p,1}
115 for (ordinal_type p=0;p<order;++p) {
116 const ordinal_type
117 loc_p_0 = Intrepid2::getPnEnumeration<spaceDim>(p,0),
118 loc_p_1 = Intrepid2::getPnEnumeration<spaceDim>(p,1);
119
120 for (ordinal_type i=0;i<npts;++i) {
121 output0(loc_p_1,i) = output0(loc_p_0,i)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
122 if(hasDeriv) {
123 output.access(loc_p_1,i,1) = output.access(loc_p_0,i,1)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
124 output.access(loc_p_1,i,2) = output.access(loc_p_0,i,2)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0)) + output0(loc_p_0,i)*(3.0+2.0*p);
125 }
126 }
127 }
128
129
130 // recurrence in q
131 for (ordinal_type p=0;p<order-1;++p)
132 for (ordinal_type q=1;q<order-p;++q) {
133 const ordinal_type
134 loc_p_qp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q+1),
135 loc_p_q = Intrepid2::getPnEnumeration<spaceDim>(p,q),
136 loc_p_qm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q-1);
137
138 value_type a,b,c;
139 Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+1,0,q);
140 for (ordinal_type i=0;i<npts;++i) {
141 output0(loc_p_qp1,i) = (a*(2.0*z(i,1)-1.0)+b)*output0(loc_p_q,i)
142 - c*output0(loc_p_qm1,i) ;
143 if(hasDeriv) {
144 output.access(loc_p_qp1,i,1) = (a*(2.0*z(i,1)-1.0)+b)*output.access(loc_p_q,i,1)
145 - c*output.access(loc_p_qm1,i,1) ;
146 output.access(loc_p_qp1,i,2) = (a*(2.0*z(i,1)-1.0)+b)*output.access(loc_p_q,i,2) +2*a*output0(loc_p_q,i)
147 - c*output.access(loc_p_qm1,i,2) ;
148 }
149 }
150 }
151 }
152
153 // orthogonalize
154 for (ordinal_type p=0;p<=order;++p)
155 for (ordinal_type q=0;q<=order-p;++q)
156 for (ordinal_type i=0;i<npts;++i) {
157 output0(Intrepid2::getPnEnumeration<spaceDim>(p,q),i) *= std::sqrt( (p+0.5)*(p+q+1.0));
158 if(hasDeriv) {
159 output.access(Intrepid2::getPnEnumeration<spaceDim>(p,q),i,1) *= std::sqrt( (p+0.5)*(p+q+1.0));
160 output.access(Intrepid2::getPnEnumeration<spaceDim>(p,q),i,2) *= std::sqrt( (p+0.5)*(p+q+1.0));
161 }
162 }
163}
164
165template<typename OutputViewType,
166typename inputViewType,
167typename workViewType,
168bool hasDeriv>
169KOKKOS_INLINE_FUNCTION
170void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,1>::generate(
171 OutputViewType output,
172 const inputViewType input,
173 workViewType work,
174 const ordinal_type order ) {
175 constexpr ordinal_type spaceDim = 2;
176 const ordinal_type
177 npts = input.extent(0),
178 card = output.extent(0);
179
180 workViewType dummyView;
181 OrthPolynomialTri<workViewType,inputViewType,workViewType,hasDeriv,0>::generate(work, input, dummyView, order);
182 for (ordinal_type i=0;i<card;++i)
183 for (ordinal_type j=0;j<npts;++j)
184 for (ordinal_type k=0;k<spaceDim;++k)
185 output.access(i,j,k) = work(i,j,k+1);
186}
187
188
189// when n >= 2, use recursion
190template<typename OutputViewType,
191typename inputViewType,
192typename workViewType,
193bool hasDeriv,
194ordinal_type n>
195KOKKOS_INLINE_FUNCTION
196void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,n>::generate(
197 OutputViewType /* output */,
198 const inputViewType /* input */,
199 workViewType /* work */,
200 const ordinal_type /* order */ ) {
201INTREPID2_TEST_FOR_ABORT( true,
202 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
203}
204
205
206
207template<EOperator opType>
208template<typename OutputViewType,
209typename inputViewType,
210typename workViewType>
211KOKKOS_INLINE_FUNCTION
212void
213Basis_HGRAD_TRI_Cn_FEM_ORTH::Serial<opType>::
214getValues( OutputViewType output,
215 const inputViewType input,
216 workViewType work,
217 const ordinal_type order) {
218 switch (opType) {
219 case OPERATOR_VALUE: {
220 OrthPolynomialTri<OutputViewType,inputViewType,workViewType,false,0>::generate( output, input, work, order );
221 break;
222 }
223 case OPERATOR_GRAD:
224 case OPERATOR_D1: {
225 OrthPolynomialTri<OutputViewType,inputViewType,workViewType,true,1>::generate( output, input, work, order );
226 break;
227 }
228 case OPERATOR_D2: {
229 OrthPolynomialTri<OutputViewType,inputViewType,workViewType,true,2>::generate( output, input, work, order );
230 break;
231 }
232 default: {
233 INTREPID2_TEST_FOR_ABORT( true,
234 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::Serial::getValues) operator is not supported");
235 }
236 }
237}
238
239template<typename DT, ordinal_type numPtsPerEval,
240typename outputValueValueType, class ...outputValueProperties,
241typename inputPointValueType, class ...inputPointProperties>
242void
243Basis_HGRAD_TRI_Cn_FEM_ORTH::
244getValues(
245 const typename DT::execution_space& space,
246 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
247 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
248 const ordinal_type order,
249 const EOperator operatorType ) {
250 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
251 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
252 typedef typename DT::execution_space ExecSpaceType;
253
254 // loopSize corresponds to the # of points
255 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
256 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
257 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
258 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
259
260 typedef typename inputPointViewType::value_type inputPointType;
261 const ordinal_type cardinality = outputValues.extent(0);
262 const ordinal_type spaceDim = 2;
263
264 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
265 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
266
267 switch (operatorType) {
268 case OPERATOR_VALUE: {
269 workViewType dummyWorkView;
270 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
271 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, dummyWorkView, order) );
272 break;
273 }
274 case OPERATOR_GRAD:
275 case OPERATOR_D1: {
276 workViewType work(Kokkos::view_alloc(space, "Basis_HGRAD_TRI_In_FEM_ORTH::getValues::work", vcprop), cardinality, inputPoints.extent(0), spaceDim+1);
277 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D1,numPtsPerEval> FunctorType;
278 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, work, order) );
279 break;
280 }
281 case OPERATOR_D2:{
282 workViewType dummyWorkView;
283 typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D2,numPtsPerEval> FunctorType;
284 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints ,dummyWorkView, order) );
285 break;
286 }
287 default: {
288 INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
289 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid operator type");
290 }
291 }
292}
293}
294
295
296// -------------------------------------------------------------------------------------
297
298template<typename DT, typename OT, typename PT>
299Basis_HGRAD_TRI_Cn_FEM_ORTH<DT,OT,PT>::
300Basis_HGRAD_TRI_Cn_FEM_ORTH( const ordinal_type order ) {
301
302 constexpr ordinal_type spaceDim = 2;
303 this->basisCardinality_ = Intrepid2::getPnCardinality<spaceDim>(order);
304 this->basisDegree_ = order;
305 this->basisCellTopologyKey_ = shards::Triangle<3>::key;
306 this->basisType_ = BASIS_FEM_HIERARCHICAL;
307 this->basisCoordinates_ = COORDINATES_CARTESIAN;
308 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
309
310 // initialize tags
311 {
312 // Basis-dependent initializations
313 constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
314 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
315 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
316 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
317
318 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
319 ordinal_type tags[maxCard][tagSize];
320 const ordinal_type card = this->basisCardinality_;
321 for (ordinal_type i=0;i<card;++i) {
322 tags[i][0] = 2; // these are all "internal" i.e. "volume" DoFs
323 tags[i][1] = 0; // there is only one line
324 tags[i][2] = i; // local DoF id
325 tags[i][3] = card; // total number of DoFs
326 }
327
328 OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
329
330 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
331 // tags are constructed on host
332 this->setOrdinalTagData(this->tagToOrdinal_,
333 this->ordinalToTag_,
334 tagView,
335 this->basisCardinality_,
336 tagSize,
337 posScDim,
338 posScOrd,
339 posDfOrd);
340 }
341
342 // dof coords is not applicable to hierarchical functions
343}
344
345}
346#endif
Intrepid2::getJacobyRecurrenceCoeffs
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
Definition Intrepid2_BasisDef.hpp:237
Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::Basis_HGRAD_TRI_Cn_FEM_ORTH
Basis_HGRAD_TRI_Cn_FEM_ORTH(const ordinal_type order)
Constructor.
Definition Intrepid2_HGRAD_TRI_Cn_FEM_ORTHDef.hpp:300
Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125
Intrepid2::Parameters::MaxNumPtsPerBasisEval
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.
Definition Intrepid2_Types.hpp:100
Generated by doxygen 1.9.8