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