Intrepid2
Intrepid2_HGRAD_PYR_I2_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_PYR_I2_SERENDIPITY_FEM_DEF_HPP__
17#define __INTREPID2_HGRAD_PYR_I2_SERENDIPITY_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22
23 namespace Impl {
24
25 template<EOperator opType>
26 template<typename OutputViewType,
27 typename inputViewType>
28 KOKKOS_INLINE_FUNCTION
29 void
30 Basis_HGRAD_PYR_I2_FEM::Serial<opType>::
31 getValues( OutputViewType output,
32 const inputViewType input ) {
33 const auto eps = epsilon();
34
35 static_assert(std::is_same<
36 typename OutputViewType::value_type,
37 typename inputViewType::value_type>::value,"Input/output view has different value types");
38
39 typedef typename OutputViewType::value_type value_type;
40
41 const value_type x = input(0);
42 const value_type y = input(1);
43 const value_type ztmp = input(2);
44
45 //be sure that the basis functions are defined when z is very close to 1.
46 const value_type z = ( (value_type(1.0) - ztmp) < value_type(eps) ? value_type(1.0 - eps) : ztmp );
47
48 switch (opType) {
49
50 case OPERATOR_VALUE: {
51 const value_type w = 1.0/(1.0 - z);
52
53 output.access(0) = 0.25 * (-x - y - 1.0)*((1.0-x)*(1.0-y) - z + x*y*z*w);
54 output.access(1) = 0.25 * ( x - y - 1.0)*((1.0+x)*(1.0-y) - z - x*y*z*w);
55 output.access(2) = 0.25 * ( x + y - 1.0)*((1.0+x)*(1.0+y) - z + x*y*z*w);
56 output.access(3) = 0.25 * (-x + y - 1.0)*((1.0-x)*(1.0+y) - z - x*y*z*w);
57
58 output.access(4) = z * (2.0*z - 1.0);
59
60 output.access(5) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 - y - z)*w;
61 output.access(6) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 + x - z)*w;
62 output.access(7) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 + y - z)*w;
63 output.access(8) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 - x - z)*w;
64
65 output.access(9) = z*(1.0 - x - z)*(1.0 - y - z)*w;
66 output.access(10) = z*(1.0 + x - z)*(1.0 - y - z)*w;
67 output.access(11) = z*(1.0 + x - z)*(1.0 + y - z)*w;
68 output.access(12) = z*(1.0 - x - z)*(1.0 + y - z)*w;
69
70 break;
71 }
72 case OPERATOR_GRAD:
73 case OPERATOR_D1: {
74 const value_type w = 1.0/(1.0 - z);
75
76 output.access(0, 0) = 0.25*(-1.0-x-y)*(-1.0+y + y*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
77 output.access(0, 1) = 0.25*(-1.0-x-y)*(-1.0+x + x*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
78 output.access(0, 2) = 0.25*(-1.0-x-y)*(-1.0 + x*y*w + x*y*z*w*w);
79
80 output.access(1, 0) = 0.25*(-1.0+x-y)*( 1.0-y - y*z*w) + 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
81 output.access(1, 1) = 0.25*(-1.0+x-y)*(-1.0-x - x*z*w) - 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
82 output.access(1, 2) = 0.25*(-1.0+x-y)*(-1.0 - x*y*w - x*y*z*w*w);
83
84 output.access(2, 0) = 0.25*(-1.0+x+y)*(1.0+y + y*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
85 output.access(2, 1) = 0.25*(-1.0+x+y)*(1.0+x + x*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
86 output.access(2, 2) = 0.25*(-1.0+x+y)*(-1.0 + x*y*w + x*y*z*w*w);
87
88 output.access(3, 0) = 0.25*(-1.0-x+y)*(-1.0-y - y*z*w) - 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
89 output.access(3, 1) = 0.25*(-1.0-x+y)*( 1.0-x - x*z*w) + 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
90 output.access(3, 2) = 0.25*(-1.0-x+y)*(-1.0 - x*y*w - x*y*z*w*w);
91
92 output.access(4, 0) = 0.0;
93 output.access(4, 1) = 0.0;
94 output.access(4, 2) = -1.0 + 4.0*z;
95
96 output.access(5, 0) = -x*w*(1.0-y-z);
97 output.access(5, 1) = -0.5*(1.0-x-z)*(1.0+x-z)*w;
98 output.access(5, 2) = 0.5*y*x*x*w*w + 0.5*y - 1.0+z;
99
100 output.access(6, 0) = 0.5*(1.0-y-z)*(1.0+y-z)*w;
101 output.access(6, 1) = -y*w*(1.0+x-z);
102 output.access(6, 2) = -0.5*x*y*y*w*w - 0.5*x - 1.0+z;
103
104 output.access(7, 0) = -x*w*(1.0+y-z);
105 output.access(7, 1) = 0.5*(1.0-x-z)*(1.0+x-z)*w;
106 output.access(7, 2) = -0.5*y*x*x*w*w - 0.5*y - 1.0+z;
107
108 output.access(8, 0) = -0.5*(1.0-y-z)*(1.0+y-z)*w;
109 output.access(8, 1) = -y*w*(1.0-x-z);
110 output.access(8, 2) = 0.5*x*y*y*w*w + 0.5*x - 1.0+z;
111
112 output.access(9, 0) = -(1.0-y-z)*z*w;
113 output.access(9, 1) = -(1.0-x-z)*z*w;
114 output.access(9, 2) = x*y*w*w + 1.0 - x - y - 2.0*z;
115
116 output.access(10,0) = (1.0-y-z)*z*w;
117 output.access(10,1) = -(1.0+x-z)*z*w;
118 output.access(10,2) = -x*y*w*w + 1.0 + x - y - 2.0*z;
119
120 output.access(11,0) = (1.0+y-z)*z*w;
121 output.access(11,1) = (1.0+x-z)*z*w;
122 output.access(11,2) = x*y*w*w + 1.0 + x + y - 2.0*z;
123
124 output.access(12,0) = -(1.0+y-z)*z*w;
125 output.access(12,1) = (1.0-x-z)*z*w;
126 output.access(12,2) = -x*y*w*w + 1.0 - x + y - 2.0*z;
127
128 break;
129 }
130 case OPERATOR_D2: {
131 const value_type w = 1.0/(1.0 - z);
132
133 output.access(0, 0) = -0.5*(-1.0+y+y*z*w);
134 output.access(0, 1) = -(-0.25 + 0.5*x + 0.5*y + 0.5*z)*w;
135 output.access(0, 2) = 0.25 + (-0.25*y-0.5*x*y-0.25*y*y)*w*w;
136 output.access(0, 3) = -0.5*(-1.0+x+x*z*w);
137 output.access(0, 4) = 0.25 + (-0.25*x*x-0.25*x-0.5*x*y)*w*w;
138 output.access(0, 5) = 0.5*x*y*(-1.0-x-y)*w*w*w;
139
140 output.access(1, 0) = 0.5*(1.0-y-y*z*w);
141 output.access(1, 1) =-(0.25 + 0.5*x - 0.5*y - 0.5*z)*w;
142 output.access(1, 2) = -0.25 + (0.25*y-0.5*x*y+0.25*y*y)*w*w;
143 output.access(1, 3) = -0.5*(-1.0-x-x*z*w);
144 output.access(1, 4) = 0.25 + (-0.25*x*x + 0.25*x + 0.5*x*y)*w*w;
145 output.access(1, 5) = -0.5*x*y*(-1.0+x-y)*w*w*w;
146
147 output.access(2, 0) = 0.5*(1.0+y+y*z*w);
148 output.access(2, 1) =-(-0.25 - 0.5*x - 0.5*y + 0.5*z)*w;
149 output.access(2, 2) = -0.25 + (-0.25*y+0.5*x*y+0.25*y*y)*w*w;
150 output.access(2, 3) = 0.5*(1.0+x+x*z*w);
151 output.access(2, 4) = -0.25 + (0.25*x*x -0.25*x + 0.5*x*y)*w*w;
152 output.access(2, 5) = 0.5*x*y*(-1.0+x+y)*w*w*w;
153
154 output.access(3, 0) = -0.5*(-1.0-y-y*z*w);
155 output.access(3, 1) =-(0.25 - 0.5*x + 0.5*y - 0.5*z)*w;
156 output.access(3, 2) = 0.25 + (0.25*y+0.5*x*y-0.25*y*y)*w*w;
157 output.access(3, 3) = 0.5*(1.0-x-x*z*w);
158 output.access(3, 4) = -0.25 + (0.25*x + 0.25*x*x - 0.5*x*y)*w*w;
159 output.access(3, 5) = -0.5*x*y*(-1.0-x+y)*w*w*w;
160
161 output.access(4, 0) = 0.0;
162 output.access(4, 1) = 0.0;
163 output.access(4, 2) = 0.0;
164 output.access(4, 3) = 0.0;
165 output.access(4, 4) = 0.0;
166 output.access(4, 5) = 4.0;
167
168 output.access(5, 0) = -(1.0-y-z)*w;
169 output.access(5, 1) = x*w;
170 output.access(5, 2) = x*y*w*w;
171 output.access(5, 3) = 0.0;
172 output.access(5, 4) = 0.5*x*x*w*w + 0.5;
173 output.access(5, 5) = x*x*y*w*w*w + 1.0;
174
175 output.access(6, 0) = 0.0;
176 output.access(6, 1) = -y*w;
177 output.access(6, 2) = -0.5*y*y*w*w - 0.5;
178 output.access(6, 3) =-(1.0+x-z)*w;
179 output.access(6, 4) = -x*y*w*w;
180 output.access(6, 5) = -x*y*y*w*w*w + 1.0;
181
182 output.access(7, 0) = -(1.0+y-z)*w;
183 output.access(7, 1) = -x*w;
184 output.access(7, 2) = -x*y*w*w;
185 output.access(7, 3) = 0.0;
186 output.access(7, 4) = -0.5*x*x*w*w - 0.5;
187 output.access(7, 5) = -x*x*y*w*w*w + 1.0;
188
189 output.access(8, 0) = 0.0;
190 output.access(8, 1) = y*w;
191 output.access(8, 2) = 0.5*y*y*w*w + 0.5;
192 output.access(8, 3) = -(1.0-x-z)*w;
193 output.access(8, 4) = x*y*w*w;
194 output.access(8, 5) = x*y*y*w*w*w + 1.0;
195
196 output.access(9, 0) = 0.0;
197 output.access(9, 1) = z*w;
198 output.access(9, 2) = y*w*w - 1.0;
199 output.access(9, 3) = 0.0;
200 output.access(9, 4) = x*w*w - 1.0;
201 output.access(9, 5) = 2.0*x*y*w*w*w - 2.0;
202
203 output.access(10,0) = 0.0;
204 output.access(10,1) = -z*w;
205 output.access(10,2) = -y*w*w + 1.0;
206 output.access(10,3) = 0.0;
207 output.access(10,4) = -x*w*w - 1.0;
208 output.access(10,5) = -2.0*x*y*w*w*w - 2.0;
209
210 output.access(11,0) = 0.0;
211 output.access(11,1) = z*w;
212 output.access(11,2) = y*w*w + 1.0;
213 output.access(11,3) = 0.0;
214 output.access(11,4) = x*w*w + 1.0;
215 output.access(11,5) = 2.0*x*y*w*w*w - 2.0;
216
217 output.access(12,0) = 0.0;
218 output.access(12,1) = -z*w;
219 output.access(12,2) = -y*w*w - 1.0;
220 output.access(12,3) = 0.0;
221 output.access(12,4) = -x*w*w + 1.0;
222 output.access(12,5) = -2.0*x*y*w*w*w - 2.0;
223
224 break;
225 }
226 default: {
227 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
228 opType != OPERATOR_GRAD &&
229 opType != OPERATOR_D1 &&
230 opType != OPERATOR_D2,
231 ">>> ERROR: (Intrepid2::Basis_HGRAD_PYR_I2_FEM::Serial::getValues) operator is not supported");
232 }
233 }
234 }
235
236 template<typename DT,
237 typename outputValueValueType, class ...outputValueProperties,
238 typename inputPointValueType, class ...inputPointProperties>
239 void
240 Basis_HGRAD_PYR_I2_FEM::
241 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
242 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
243 const EOperator operatorType ) {
244 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
245 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
246 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
247
248 // Number of evaluation points = dim 0 of inputPoints
249 const auto loopSize = inputPoints.extent(0);
250 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
251
252 switch (operatorType) {
253
254 case OPERATOR_VALUE: {
255 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
256 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
257 break;
258 }
259 case OPERATOR_GRAD:
260 case OPERATOR_D1: {
261 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
262 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
263 break;
264 }
265 case OPERATOR_CURL: {
266 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
267 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
268 break;
269 }
270 case OPERATOR_DIV: {
271 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
272 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
273 break;
274 }
275 case OPERATOR_D2: {
276 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
277 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
278 break;
279 }
280 case OPERATOR_D3:
281 case OPERATOR_D4:
282 case OPERATOR_D5:
283 case OPERATOR_D6:
284 case OPERATOR_D7:
285 case OPERATOR_D8:
286 case OPERATOR_D9:
287 case OPERATOR_D10: {
288 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
289 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Operator not implemented yet");
290 break;
291 }
292 default: {
293 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
294 ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Invalid operator type");
295 }
296 }
297 }
298 }
299
300 // -------------------------------------------------------------------------------------
301
302 template<typename DT, typename OT, typename PT>
303 Basis_HGRAD_PYR_I2_FEM<DT,OT,PT>::
304 Basis_HGRAD_PYR_I2_FEM() {
305 const ordinal_type spaceDim = 3;
306 this->basisCardinality_ = 13;
307 this->basisDegree_ = 2;
308 this->basisCellTopologyKey_ = shards::Pyramid<5>::key;
309 this->basisType_ = BASIS_FEM_DEFAULT;
310 this->basisCoordinates_ = COORDINATES_CARTESIAN;
311 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
312
313 // initialize tags
314 {
315 // Basis-dependent intializations
316 const ordinal_type tagSize = 4; // size of DoF tag
317 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
318 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
319 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
320
321 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
322 ordinal_type tags[52] = { 0, 0, 0, 1,
323 0, 1, 0, 1,
324 0, 2, 0, 1,
325 0, 3, 0, 1,
326 0, 4, 0, 1,
327 1, 0, 0, 1,
328 1, 1, 0, 1,
329 1, 2, 0, 1,
330 1, 3, 0, 1,
331 1, 4, 0, 1,
332 1, 5, 0, 1,
333 1, 6, 0, 1,
334 1, 7, 0, 1 };
335
336
337 // host tags
338 OrdinalTypeArray1DHost tagView(&tags[0], 52);
339
340 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
341 this->setOrdinalTagData(this->tagToOrdinal_,
342 this->ordinalToTag_,
343 tagView,
344 this->basisCardinality_,
345 tagSize,
346 posScDim,
347 posScOrd,
348 posDfOrd);
349 }
350
351 // dofCoords on host and create its mirror view to device
352 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
353 dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
354
355 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = 0.0;
356 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = 0.0;
357 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
358 dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = 0.0;
359 dofCoords(4,0) = 0.0; dofCoords(4,1) = 0.0; dofCoords(4,2) = 1.0;
360
361 dofCoords(5,0) = 0.0; dofCoords(5,1) = -1.0; dofCoords(5,2) = 0.0;
362 dofCoords(6,0) = 1.0; dofCoords(6,1) = 0.0; dofCoords(6,2) = 0.0;
363 dofCoords(7,0) = 0.0; dofCoords(7,1) = 1.0; dofCoords(7,2) = 0.0;
364 dofCoords(8,0) = -1.0; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.0;
365 dofCoords(9,0) = -0.5; dofCoords(9,1) = -0.5; dofCoords(9,2) = 0.5;
366 dofCoords(10,0) = 0.5; dofCoords(10,1) = -0.5; dofCoords(10,2) = 0.5;
367 dofCoords(11,0) = 0.5; dofCoords(11,1) = 0.5; dofCoords(11,2) = 0.5;
368 dofCoords(12,0) = -0.5; dofCoords(12,1) = 0.5; dofCoords(12,2) = 0.5;
369
370 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
371 Kokkos::deep_copy(this->dofCoords_, dofCoords);
372 }
373
374 template<typename DT, typename OT, typename PT>
375 void
376 Basis_HGRAD_PYR_I2_FEM<DT,OT,PT>::getScratchSpaceSize(
377 ordinal_type& perTeamSpaceSize,
378 ordinal_type& perThreadSpaceSize,
379 const PointViewType inputPoints,
380 const EOperator operatorType) const {
381 perTeamSpaceSize = 0;
382 perThreadSpaceSize = 0;
383 }
384
385 template<typename DT, typename OT, typename PT>
386 KOKKOS_INLINE_FUNCTION
387 void
388 Basis_HGRAD_PYR_I2_FEM<DT,OT,PT>::getValues(
389 OutputViewType outputValues,
390 const PointViewType inputPoints,
391 const EOperator operatorType,
392 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
393 const typename DT::execution_space::scratch_memory_space & scratchStorage,
394 const ordinal_type subcellDim,
395 const ordinal_type subcellOrdinal) const {
396
397 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
398 ">>> ERROR: (Intrepid2::Basis_HGRAD_PYR_I2_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
399
400 (void) scratchStorage; //avoid unused variable warning
401
402 const int numPoints = inputPoints.extent(0);
403
404 switch(operatorType) {
405 case OPERATOR_VALUE:
406 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
407 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
408 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
409 Impl::Basis_HGRAD_PYR_I2_FEM::Serial<OPERATOR_VALUE>::getValues( output, input);
410 });
411 break;
412 case OPERATOR_GRAD:
413 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
414 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
415 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
416 Impl::Basis_HGRAD_PYR_I2_FEM::Serial<OPERATOR_GRAD>::getValues( output, input);
417 });
418 break;
419 default: {}
420 }
421 }
422
423}// namespace Intrepid2
424#endif
Intrepid2::Basis_HGRAD_PYR_I2_FEM
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
Definition Intrepid2_HGRAD_PYR_I2_FEM.hpp:148
Intrepid2::Basis_HGRAD_PYR_I2_FEM::Basis_HGRAD_PYR_I2_FEM
Basis_HGRAD_PYR_I2_FEM()
Constructor.
Definition Intrepid2_HGRAD_PYR_I2_FEMDef.hpp:304
Generated by doxygen 1.9.8