Intrepid2
Intrepid2_HGRAD_QUAD_C2_FEMDef.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_QUAD_C2_FEM_DEF_HPP__
17#define __INTREPID2_HGRAD_QUAD_C2_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22
23 namespace Impl {
24
25 template<bool serendipity>
26 template<EOperator opType>
27 template<typename OutputViewType,
28 typename inputViewType>
29 KOKKOS_INLINE_FUNCTION
30 void
31 Basis_HGRAD_QUAD_DEG2_FEM<serendipity>::Serial<opType>::
32 getValues( OutputViewType output,
33 const inputViewType input ) {
34 switch (opType) {
35 case OPERATOR_VALUE : {
36 const auto x = input(0);
37 const auto y = input(1);
38
39 // output is a rank-1 array with dimensions (basisCardinality_)
40 if constexpr (!serendipity) {
41 output.access(0) = x*(x - 1.0)*y*(y - 1.0)/4.0;
42 output.access(1) = x*(x + 1.0)*y*(y - 1.0)/4.0;
43 output.access(2) = x*(x + 1.0)*y*(y + 1.0)/4.0;
44 output.access(3) = x*(x - 1.0)*y*(y + 1.0)/4.0;
45 // edge midpoints basis functions
46 output.access(4) = (1.0 - x)*(1.0 + x)*y*(y - 1.0)/2.0;
47 output.access(5) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0;
48 output.access(6) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0;
49 output.access(7) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0;
50
51 // quad bubble basis function
52 output.access(8) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y);
53
54 } else { //serendipity
55
56 output.access(0) = 0.25*(1.0 - x)*(1.0 - y)*(-x - y - 1.0);
57 output.access(1) = 0.25*(1.0 + x)*(1.0 - y)*( x - y - 1.0);
58 output.access(2) = 0.25*(1.0 + x)*(1.0 + y)*( x + y - 1.0);
59 output.access(3) = 0.25*(1.0 - x)*(1.0 + y)*(-x + y - 1.0);
60
61 output.access(4) = 0.5*(1.0 - x*x)*(1.0 - y);
62 output.access(5) = 0.5*(1.0 + x)*(1.0 - y*y);
63 output.access(6) = 0.5*(1.0 - x*x)*(1.0 + y);
64 output.access(7) = 0.5*(1.0 - x)*(1.0 - y*y);
65 }
66
67 break;
68 }
69 case OPERATOR_D1 :
70 case OPERATOR_GRAD : {
71 const auto x = input(0);
72 const auto y = input(1);
73
74 if constexpr (!serendipity) {
75
76 output.access(0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y;
77 output.access(0, 1) = (-1.0 + x)*x*(-0.25 + 0.5*y);
78
79 output.access(1, 0) = (0.25 + 0.5*x)*(-1. + y)*y;
80 output.access(1, 1) = x*(1. + x)*(-0.25 + 0.5*y);
81
82 output.access(2, 0) = (0.25 + 0.5*x)*y*(1. + y);
83 output.access(2, 1) = x*(1. + x)*(0.25 + 0.5*y);
84
85 output.access(3, 0) = (-0.25 + 0.5*x)*y*(1. + y);
86 output.access(3, 1) = (-1. + x)*x*(0.25 + 0.5*y);
87
88 output.access(4, 0) = x*(1.0 - y)*y;
89 output.access(4, 1) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
90
91 output.access(5, 0) = 0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
92 output.access(5, 1) =-x*(1.0 + x)*y;
93
94 output.access(6, 0) =-y*(1.0 + y)*x;
95 output.access(6, 1) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
96
97 output.access(7, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x);
98 output.access(7, 1) = (1.0 - x)*x*y;
99
100 output.access(8, 0) =-2.0*(1.0 - y)*(1.0 + y)*x;
101 output.access(8, 1) =-2.0*(1.0 - x)*(1.0 + x)*y;
102
103 } else { //serendipity
104
105 output.access(0, 0) = -0.25*(1.0-y)*(-x-y-1.0) - 0.25*(1.0-x)*(1.0-y);
106 output.access(0, 1) = -0.25*(1.0-x)*(-x-y-1.0) - 0.25*(1.0-x)*(1.0-y);
107
108 output.access(1, 0) = 0.25*(1.0-y)*( x-y-1.0) + 0.25*(1.0+x)*(1.0-y);
109 output.access(1, 1) = -0.25*(1.0+x)*( x-y-1.0) - 0.25*(1.0+x)*(1.0-y);
110
111 output.access(2, 0) = 0.25*(1.0+y)*( x+y-1.0) + 0.25*(1.0+x)*(1.0+y);
112 output.access(2, 1) = 0.25*(1.0+x)*( x+y-1.0) + 0.25*(1.0+x)*(1.0+y);
113
114 output.access(3, 0) = -0.25*(1.0+y)*(-x+y-1.0) - 0.25*(1.0-x)*(1.0+y);
115 output.access(3, 1) = 0.25*(1.0-x)*(-x+y-1.0) + 0.25*(1.0-x)*(1.0+y);
116
117 output.access(4, 0) = -x*(1.0-y);
118 output.access(4, 1) = -0.5*(1.0-x*x);
119
120 output.access(5, 0) = 0.5*(1.0-y*y);
121 output.access(5, 1) = -y*(1.0+x);
122
123 output.access(6, 0) = -x*(1.0+y);
124 output.access(6, 1) = 0.5*(1.0-x*x);
125
126 output.access(7, 0) = -0.5*(1.0-y*y);
127 output.access(7, 1) = -y*(1.0-x);
128 }
129 break;
130 }
131 case OPERATOR_CURL : {
132 const auto x = input(0);
133 const auto y = input(1);
134
135 // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
136 // CURL(u) = (u_y, -u_x), is rotated GRAD
137
138 if constexpr (!serendipity) {
139 output.access(0, 1) =-(-0.25 + 0.5*x)*(-1. + y)*y;
140 output.access(0, 0) = (-1.0 + x)*x*(-0.25 + 0.5*y);
141
142 output.access(1, 1) =-(0.25 + 0.5*x)*(-1. + y)*y;
143 output.access(1, 0) = x*(1. + x)*(-0.25 + 0.5*y);
144
145 output.access(2, 1) =-(0.25 + 0.5*x)*y*(1. + y);
146 output.access(2, 0) = x*(1. + x)*(0.25 + 0.5*y);
147
148 output.access(3, 1) =-(-0.25 + 0.5*x)*y*(1. + y);
149 output.access(3, 0) = (-1. + x)*x*(0.25 + 0.5*y);
150
151 output.access(4, 1) =-x*(1.0 - y)*y;
152 output.access(4, 0) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
153
154 output.access(5, 1) =-0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
155 output.access(5, 0) =-x*(1.0 + x)*y;
156
157 output.access(6, 1) = y*(1.0 + y)*x;
158 output.access(6, 0) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
159
160 output.access(7, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x);
161 output.access(7, 0) = (1.0 - x)*x*y;
162
163 output.access(8, 1) = 2.0*(1.0 - y)*(1.0 + y)*x;
164 output.access(8, 0) =-2.0*(1.0 - x)*(1.0 + x)*y;
165
166 } else { //serendipity
167 output.access(0, 1) = 0.25*(1.0-y)*(-x-y-1.0) + 0.25*(1.0-x)*(1.0-y);
168 output.access(0, 0) = -0.25*(1.0-x)*(-x-y-1.0) - 0.25*(1.0-x)*(1.0-y);
169
170 output.access(1, 1) = -0.25*(1.0-y)*( x-y-1.0) - 0.25*(1.0+x)*(1.0-y);
171 output.access(1, 0) = -0.25*(1.0+x)*( x-y-1.0) - 0.25*(1.0+x)*(1.0-y);
172
173 output.access(2, 1) = -0.25*(1.0+y)*( x+y-1.0) - 0.25*(1.0+x)*(1.0+y);
174 output.access(2, 0) = 0.25*(1.0+x)*( x+y-1.0) + 0.25*(1.0+x)*(1.0+y);
175
176 output.access(3, 1) = 0.25*(1.0+y)*(-x+y-1.0) + 0.25*(1.0-x)*(1.0+y);
177 output.access(3, 0) = 0.25*(1.0-x)*(-x+y-1.0) + 0.25*(1.0-x)*(1.0+y);
178
179 output.access(4, 1) = x*(1.0-y);
180 output.access(4, 0) = -0.5*(1.0-x*x);
181
182 output.access(5, 1) = -0.5*(1.0-y*y);
183 output.access(5, 0) = -y*(1.0+x);
184
185 output.access(6, 1) = x*(1.0+y);
186 output.access(6, 0) = 0.5*(1.0-x*x);
187
188 output.access(7, 1) = 0.5*(1.0-y*y);
189 output.access(7, 0) = -y*(1.0-x);
190 }
191 break;
192 }
193 case OPERATOR_D2 : {
194 const auto x = input(0);
195 const auto y = input(1);
196 // output.access is a rank-3 array with dimensions (basisCardinality_, D2Cardinality=3)
197
198 if constexpr (!serendipity) {
199
200 output.access(0, 0) = 0.5*(-1.0 + y)*y;
201 output.access(0, 1) = 0.25 - 0.5*y + x*(-0.5 + 1.*y);
202 output.access(0, 2) = 0.5*(-1.0 + x)*x;
203
204 output.access(1, 0) = 0.5*(-1.0 + y)*y;
205 output.access(1, 1) =-0.25 + 0.5*y + x*(-0.5 + 1.*y);
206 output.access(1, 2) = 0.5*x*(1.0 + x);
207
208 output.access(2, 0) = 0.5*y*(1.0 + y);
209 output.access(2, 1) = 0.25 + 0.5*y + x*(0.5 + 1.*y);
210 output.access(2, 2) = 0.5*x*(1.0 + x);
211
212 output.access(3, 0) = 0.5*y*(1.0 + y);
213 output.access(3, 1) =-0.25 - 0.5*y + x*(0.5 + 1.*y);
214 output.access(3, 2) = 0.5*(-1.0 + x)*x;
215
216 output.access(4, 0) = (1.0 - y)*y;
217 output.access(4, 1) = x*(1. - 2.*y);
218 output.access(4, 2) = (1.0 - x)*(1.0 + x);
219
220 output.access(5, 0) = (1.0 - y)*(1.0 + y);
221 output.access(5, 1) = x*(0. - 2.*y) - 1.*y;
222 output.access(5, 2) =-x*(1.0 + x);
223
224 output.access(6, 0) =-y*(1.0 + y);
225 output.access(6, 1) = x*(-1. - 2.*y);
226 output.access(6, 2) = (1.0 - x)*(1.0 + x);
227
228 output.access(7, 0) = (1.0 - y)*(1.0 + y);
229 output.access(7, 1) = x*(0. - 2.*y) + 1.*y;
230 output.access(7, 2) = (1.0 - x)*x;
231
232 output.access(8, 0) =-2.0 + 2.0*y*y;
233 output.access(8, 1) = 4*x*y;
234 output.access(8, 2) =-2.0 + 2.0*x*x;
235
236 } else { //serendipity
237
238 output.access(0, 0) = 0.5*(1.0 - y);
239 output.access(0, 1) = 0.25*(1.0 - 2.0*x - 2.0*y);
240 output.access(0, 2) = 0.5*(1.0 - x);
241
242 output.access(1, 0) = 0.5*(1.0 - y);
243 output.access(1, 1) = -0.25*(1.0 + 2.0*x - 2.0*y);
244 output.access(1, 2) = 0.5*(1.0 + x);
245
246 output.access(2, 0) = 0.5*(1.0 + y);
247 output.access(2, 1) = 0.25*(1.0 + 2.0*x + 2.0*y);
248 output.access(2, 2) = 0.5*(1.0 + x);
249
250 output.access(3, 0) = 0.5*(1.0 + y);
251 output.access(3, 1) = -0.25*(1.0 - 2.0*x + 2.0*y);
252 output.access(3, 2) = 0.5*(1.0 - x);
253
254 output.access(4, 0) = -(1.0 - y);
255 output.access(4, 1) = x;
256 output.access(4, 2) = 0.0;
257
258 output.access(5, 0) = 0.0;
259 output.access(5, 1) = -y;
260 output.access(5, 2) = -(1.0 + x);
261
262 output.access(6, 0) = -(1.0 + y);
263 output.access(6, 1) = -x;
264 output.access(6, 2) = 0.0;
265
266 output.access(7, 0) = 0.0;
267 output.access(7, 1) = y;
268 output.access(7, 2) = -(1.0 - x);
269
270 }
271 break;
272 }
273 case OPERATOR_D3 : {
274 if constexpr (!serendipity) {
275 const auto x = input(0);
276 const auto y = input(1);
277 output.access(0, 0) = 0.0;
278 output.access(0, 1) =-0.5 + y;
279 output.access(0, 2) =-0.5 + x;
280 output.access(0, 3) = 0.0;
281
282 output.access(1, 0) = 0.0;
283 output.access(1, 1) =-0.5 + y;
284 output.access(1, 2) = 0.5 + x;
285 output.access(1, 3) = 0.0;
286
287 output.access(2, 0) = 0.0;
288 output.access(2, 1) = 0.5 + y;
289 output.access(2, 2) = 0.5 + x;
290 output.access(2, 3) = 0.0;
291
292 output.access(3, 0) = 0.0;
293 output.access(3, 1) = 0.5 + y;
294 output.access(3, 2) =-0.5 + x;
295 output.access(3, 3) = 0.0;
296
297 output.access(4, 0) = 0.0;
298 output.access(4, 1) = 1.0 - 2.0*y;
299 output.access(4, 2) =-2.0*x;
300 output.access(4, 3) = 0.0;
301
302 output.access(5, 0) = 0.0;
303 output.access(5, 1) =-2.0*y;
304 output.access(5, 2) =-1.0 - 2.0*x;
305 output.access(5, 3) = 0.0;
306
307 output.access(6, 0) = 0.0;
308 output.access(6, 1) =-1.0 - 2.0*y;
309 output.access(6, 2) =-2.0*x;
310 output.access(6, 3) = 0.0;
311
312 output.access(7, 0) = 0.0;
313 output.access(7, 1) =-2.0*y;
314 output.access(7, 2) = 1.0 - 2.0*x;
315 output.access(7, 3) = 0.0;
316
317 output.access(8, 0) = 0.0;
318 output.access(8, 1) = 4.0*y;
319 output.access(8, 2) = 4.0*x;
320 output.access(8, 3) = 0.0;
321
322 } else { //serendipity
323
324 output.access(0, 0) = 0.0;
325 output.access(0, 1) =-0.5;
326 output.access(0, 2) =-0.5;
327 output.access(0, 3) = 0.0;
328
329 output.access(1, 0) = 0.0;
330 output.access(1, 1) =-0.5;
331 output.access(1, 2) = 0.5;
332 output.access(1, 3) = 0.0;
333
334 output.access(2, 0) = 0.0;
335 output.access(2, 1) = 0.5;
336 output.access(2, 2) = 0.5;
337 output.access(2, 3) = 0.0;
338
339 output.access(3, 0) = 0.0;
340 output.access(3, 1) = 0.5;
341 output.access(3, 2) =-0.5;
342 output.access(3, 3) = 0.0;
343
344 output.access(4, 0) = 0.0;
345 output.access(4, 1) = 1.0;
346 output.access(4, 2) = 0.0;
347 output.access(4, 3) = 0.0;
348
349 output.access(5, 0) = 0.0;
350 output.access(5, 1) = 0.0;
351 output.access(5, 2) =-1.0;
352 output.access(5, 3) = 0.0;
353
354 output.access(6, 0) = 0.0;
355 output.access(6, 1) =-1.0;
356 output.access(6, 2) = 0.0;
357 output.access(6, 3) = 0.0;
358
359 output.access(7, 0) = 0.0;
360 output.access(7, 1) = 0.0;
361 output.access(7, 2) = 1.0;
362 output.access(7, 3) = 0.0;
363 }
364 break;
365 }
366 case OPERATOR_D4 : {
367
368 const ordinal_type jend = output.extent(1);
369 const ordinal_type iend = output.extent(0);
370
371 for (ordinal_type j=0;j<jend;++j)
372 for (ordinal_type i=0;i<iend;++i)
373 output.access(i, j) = 0.0;
374
375 if constexpr (!serendipity) {
376 output.access(0, 2) = 1.0;
377 output.access(1, 2) = 1.0;
378 output.access(2, 2) = 1.0;
379 output.access(3, 2) = 1.0;
380
381 output.access(4, 2) =-2.0;
382 output.access(5, 2) =-2.0;
383 output.access(6, 2) =-2.0;
384 output.access(7, 2) =-2.0;
385
386 output.access(8, 2) = 4.0;
387 }
388 break;
389 }
390 case OPERATOR_MAX : {
391 const ordinal_type jend = output.extent(1);
392 const ordinal_type iend = output.extent(0);
393
394 for (ordinal_type j=0;j<jend;++j)
395 for (ordinal_type i=0;i<iend;++i)
396 output.access(i, j) = 0.0;
397 break;
398 }
399 default: {
400 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
401 opType != OPERATOR_GRAD &&
402 opType != OPERATOR_CURL &&
403 opType != OPERATOR_D1 &&
404 opType != OPERATOR_D2 &&
405 opType != OPERATOR_D3 &&
406 opType != OPERATOR_D4 &&
407 opType != OPERATOR_MAX,
408 ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::Serial::getValues) operator is not supported");
409
410 }
411 }
412 }
413
414 template<bool serendipity>
415 template<typename DT,
416 typename outputValueValueType, class ...outputValueProperties,
417 typename inputPointValueType, class ...inputPointProperties>
418 void
419 Basis_HGRAD_QUAD_DEG2_FEM<serendipity>::
420 getValues( const typename DT::execution_space& space,
421 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
422 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
423 const EOperator operatorType ) {
424 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
425 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
426 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
427
428 // Number of evaluation points = dim 0 of inputPoints
429 const auto loopSize = inputPoints.extent(0);
430 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
431
432 switch (operatorType) {
433
434 case OPERATOR_VALUE: {
435 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
436 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
437 break;
438 }
439 case OPERATOR_GRAD:
440 case OPERATOR_D1: {
441 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
442 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
443 break;
444 }
445 case OPERATOR_CURL: {
446 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
447 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
448 break;
449 }
450 case OPERATOR_DIV: {
451 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
452 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
453 break;
454 }
455 case OPERATOR_D2: {
456 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
457 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
458 break;
459 }
460 case OPERATOR_D3: {
461 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D3Cardinality=4)
462 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D3> FunctorType;
463 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
464 break;
465 }
466 case OPERATOR_D4: {
467 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D4Cardinality=5)
468 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D4> FunctorType;
469 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
470 break;
471 }
472 case OPERATOR_D5:
473 case OPERATOR_D6:
474 case OPERATOR_D7:
475 case OPERATOR_D8:
476 case OPERATOR_D9:
477 case OPERATOR_D10: {
478 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
479 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
480 break;
481 }
482 default: {
483 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
484 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
485 }
486 }
487 }
488
489 }
490 // -------------------------------------------------------------------------------------
491
492
493 template<bool serendipity, typename DT, typename OT, typename PT>
494 Basis_HGRAD_QUAD_DEG2_FEM<serendipity, DT,OT,PT>::
495 Basis_HGRAD_QUAD_DEG2_FEM() {
496 const ordinal_type spaceDim = 2;
497 this->basisCardinality_ = serendipity ? 8 : 9;
498 this->basisDegree_ = 2;
499 this->basisCellTopologyKey_ = shards::Quadrilateral<4>::key;
500 this->basisType_ = BASIS_FEM_DEFAULT;
501 this->basisCoordinates_ = COORDINATES_CARTESIAN;
502 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
503
504 {
505 // Basis-dependent intializations
506 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
507 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
508 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
509 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
510
511 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
512 ordinal_type tags[36] = { 0, 0, 0, 1,
513 0, 1, 0, 1,
514 0, 2, 0, 1,
515 0, 3, 0, 1,
516 // edge midpoints
517 1, 0, 0, 1,
518 1, 1, 0, 1,
519 1, 2, 0, 1,
520 1, 3, 0, 1,
521 // quad center, not used for serendipity elements
522 2, 0, 0, 1};
523
524 //host view
525 OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 32 : 36);
526
527 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
528 this->setOrdinalTagData(this->tagToOrdinal_,
529 this->ordinalToTag_,
530 tagView,
531 this->basisCardinality_,
532 tagSize,
533 posScDim,
534 posScOrd,
535 posDfOrd);
536 }
537
538 // dofCoords on host and create its mirror view to device
539 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
540 dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
541
542 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0;
543 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0;
544 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0;
545 dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0;
546
547 dofCoords(4,0) = 0.0; dofCoords(4,1) = -1.0;
548 dofCoords(5,0) = 1.0; dofCoords(5,1) = 0.0;
549 dofCoords(6,0) = 0.0; dofCoords(6,1) = 1.0;
550 dofCoords(7,0) = -1.0; dofCoords(7,1) = 0.0;
551
552 if constexpr (!serendipity) {
553 dofCoords(8,0) = 0.0; dofCoords(8,1) = 0.0;
554 }
555
556 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
557 Kokkos::deep_copy(this->dofCoords_, dofCoords);
558 }
559
560 template<bool serendipity, typename DT, typename OT, typename PT>
561 void
562 Basis_HGRAD_QUAD_DEG2_FEM<serendipity,DT,OT,PT>::getScratchSpaceSize(
563 ordinal_type& perTeamSpaceSize,
564 ordinal_type& perThreadSpaceSize,
565 const PointViewType inputPoints,
566 const EOperator operatorType) const {
567 perTeamSpaceSize = 0;
568 perThreadSpaceSize = 0;
569 }
570
571 template<bool serendipity, typename DT, typename OT, typename PT>
572 KOKKOS_INLINE_FUNCTION
573 void
574 Basis_HGRAD_QUAD_DEG2_FEM<serendipity,DT,OT,PT>::getValues(
575 OutputViewType outputValues,
576 const PointViewType inputPoints,
577 const EOperator operatorType,
578 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
579 const typename DT::execution_space::scratch_memory_space & scratchStorage,
580 const ordinal_type subcellDim,
581 const ordinal_type subcellOrdinal) const {
582
583 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
584 ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
585
586 (void) scratchStorage; //avoid unused variable warning
587
588 const int numPoints = inputPoints.extent(0);
589
590 switch(operatorType) {
591 case OPERATOR_VALUE:
592 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
593 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
594 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
595 using SerialValue = typename Impl::Basis_HGRAD_QUAD_DEG2_FEM<serendipity>::template Serial<OPERATOR_VALUE>;
596 SerialValue::getValues( output, input);
597 });
598 break;
599 case OPERATOR_GRAD:
600 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
601 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
602 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
603 using SerialGrad = typename Impl::Basis_HGRAD_QUAD_DEG2_FEM<serendipity>::template Serial<OPERATOR_GRAD>;
604 SerialGrad::getValues( output, input);
605 });
606 break;
607 case OPERATOR_CURL:
608 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
609 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
610 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
611 using SerialCurl = typename Impl::Basis_HGRAD_QUAD_DEG2_FEM<serendipity>::template Serial<OPERATOR_CURL>;
612 SerialCurl::getValues( output, input);
613 });
614 break;
615 default: {
616 INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getValues), Operator Type not supported.");
617 }
618 }
619 }
620
621}// namespace Intrepid2
622#endif
Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell.
Definition Intrepid2_HGRAD_QUAD_C2_FEM.hpp:158
Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::Basis_HGRAD_QUAD_DEG2_FEM
Basis_HGRAD_QUAD_DEG2_FEM()
Constructor.
Definition Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp:495
Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getValues
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Definition Intrepid2_HGRAD_QUAD_C2_FEM.hpp:178
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125
Generated by doxygen 1.9.8