Intrepid2
Intrepid2_HGRAD_TET_COMP12_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
17#ifndef __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
18#define __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
19
20namespace Intrepid2 {
21
22 // -------------------------------------------------------------------------------------
23 namespace Impl {
24
25 // assume that subtets are disjoint and a single point belong to one subtet
26 // at the interface, it returns the first one that satisfy the condition
27 template<typename pointValueType>
28 KOKKOS_INLINE_FUNCTION
29 ordinal_type
31 const pointValueType y,
32 const pointValueType z ) {
33
34 const pointValueType
35 xyz = x + y + z,
36 xy = x + y,
37 xz = x + z,
38 yz = y + z;
39
40 // cycle through each subdomain and push back if the point lies within
41
42 // subtet #0 E0 := 0.0 <= r + s + t <= 0.5 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
43 if ( (0.0 <= xyz && xyz <= 0.5) &&
44 (0.0 <= x && x <= 0.5) &&
45 (0.0 <= y && y <= 0.5) &&
46 (0.0 <= z && z <= 0.5) )
47 return 0;
48
49 // subtet #1 E1 := 0.5 <= r + s + t <= 1.0 && 0.5 <= r <= 1.0 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
50 if ( (0.5 <= xyz && xyz <= 1.0) &&
51 (0.5 <= x && x <= 1.0) &&
52 (0.0 <= y && y <= 0.5) &&
53 (0.0 <= z && z <= 0.5) )
54 return 1;
55
56 // subtet #2 E2 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.5 <= s <= 1.0 && 0.0 <= t <= 0.5
57 if ( (0.5 <= xyz && xyz <= 1.0) &&
58 (0.0 <= x && x <= 0.5) &&
59 (0.5 <= y && y <= 1.0) &&
60 (0.0 <= z && z <= 0.5) )
61 return 2;
62
63 // subtet #3 E3 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.5 <= t <= 1.0
64 if ( (0.5 <= xyz && xyz <= 1.0) &&
65 (0.0 <= x && x <= 0.5) &&
66 (0.0 <= y && y <= 0.5) &&
67 (0.5 <= z && z <= 1.0) )
68 return 3;
69
70 // subtet #4 E4 := 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= r + t <= 1.0 && 0.0 <= r <= 0.5
71 if ( (0.0 <= yz && yz <= 0.5) &&
72 (0.5 <= xy && xy <= 1.0) &&
73 (0.5 <= xz && xz <= 1.0) &&
74 (0.0 <= x && x <= 0.5) )
75 return 4;
76
77 // subtet #5 E5 := 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.5 <= r + t <= 1.0 && 0.75 <= r + s + t <= 1.0
78 if ( (0.5 <= xy && xy <= 1.0) &&
79 (0.5 <= yz && yz <= 1.0) &&
80 (0.5 <= xz && xz <= 1.0) &&
81 (0.75 <= xyz && xyz <= 1.0) )
82 return 5;
83
84 // subtet #6 E6 := 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= t <= 0.5
85 if ( (0.5 <= yz && yz <= 1.0) &&
86 (0.0 <= xy && xy <= 0.5) &&
87 (0.5 <= xz && xz <= 1.0) &&
88 (0.0 <= z && z <= 0.5) )
89 return 6;
90
91 // subtet #7 E7 := 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= s <= 0.25
92 if ( (0.0 <= yz && yz <= 0.5) &&
93 (0.0 <= xy && xy <= 0.5) &&
94 (0.5 <= xz && xz <= 1.0) &&
95 (0.0 <= y && y <= 0.25) )
96 return 7;
97
98 // subtet #8 E8 := 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.0 <= t <= 0.25
99 if ( (0.0 <= xz && xz <= 0.5) &&
100 (0.0 <= yz && yz <= 0.5) &&
101 (0.5 <= xy && xy <= 1.0) &&
102 (0.0 <= z && z <= 0.25) )
103 return 8;
104
105 // subtet #9 E9 := 0.0 <= r + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.0 <= s <= 0.5
106 if ( (0.0 <= xz && xz <= 0.5) &&
107 (0.5 <= xy && xy <= 1.0) &&
108 (0.5 <= yz && yz <= 1.0) &&
109 (0.0 <= y && y <= 0.5) )
110 return 9;
111
112 // subtet #10 E10 := 0.0 <= r + t <= 0.5 && 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.0 <= r <= 0.25
113 if ( (0.0 <= xz && xz <= 0.5) &&
114 (0.5 <= yz && yz <= 1.0) &&
115 (0.0 <= xy && xy <= 0.5) &&
116 (0.0 <= x && x <= 0.25) )
117 return 10;
118
119 // subtet #11 E11 := 0.5 <= r + s + t <= 0.75 && 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5
120 if ( (0.5 <= xyz && xyz <= 0.75) &&
121 (0.0 <= xz && xz <= 0.5) &&
122 (0.0 <= yz && yz <= 0.5) &&
123 (0.0 <= xy && xy <= 0.5) )
124 return 11;
125
126 return -1;
127 }
128
129 template<EOperator opType>
130 template<typename outputValueViewType,
131 typename inputPointViewType>
132 KOKKOS_INLINE_FUNCTION
133 void
135 getValues( outputValueViewType output,
136 const inputPointViewType input ) {
137 switch (opType) {
138 case OPERATOR_VALUE: {
139 const typename inputPointViewType::value_type r = input(0);
140 const typename inputPointViewType::value_type s = input(1);
141 const typename inputPointViewType::value_type t = input(2);
142
143 // initialize output
144 for (auto i=0;i<10;++i)
145 output.access(i) = 0.0;
146
147 const auto subtet = getLocalSubTet( r, s, t );
148
149 // idependent verification shows that a given point will produce
150 // the same shape functions for each tet that contains it, so
151 // we only need to use the first one returned.
152 if (subtet != -1) {
153 typename inputPointViewType::value_type aux = 0.0;
154 switch (subtet) {
155 case 0:
156 output.access(0) = 1. - 2. * (r + s + t);
157 output.access(4) = 2. * r;
158 output.access(6) = 2. * s;
159 output.access(7) = 2. * t;
160 break;
161 case 1:
162 output.access(1) = 2. * r - 1.;
163 output.access(4) = 2. - 2. * (r + s + t);
164 output.access(5) = 2. * s;
165 output.access(8) = 2. * t;
166 break;
167 case 2:
168 output.access(2) = 2. * s - 1.;
169 output.access(5) = 2. * r;
170 output.access(6) = 2. - 2. * (r + s + t);
171 output.access(9) = 2. * t;
172 break;
173 case 3:
174 output.access(3) = 2. * t - 1.;
175 output.access(7) = 2. - 2. * (r + s + t);
176 output.access(8) = 2. * r;
177 output.access(9) = 2. * s;
178 break;
179 case 4:
180 output.access(4) = 1. - 2. * (s + t);
181 output.access(5) = 2. * (r + s) - 1.;
182 output.access(8) = 2. * (r + t) - 1.;
183 aux = 2. - 4. * r;
184 break;
185 case 5:
186 output.access(5) = 2. * (r + s) - 1.;
187 output.access(8) = 2. * (r + t) - 1.;
188 output.access(9) = 2. * (s + t) - 1.;
189 aux = 4. - 4. * (r + s + t);
190 break;
191 case 6:
192 output.access(7) = 1. - 2. * (r + s);
193 output.access(8) = 2. * (r + t) - 1.;
194 output.access(9) = 2. * (s + t) - 1.;
195 aux = 2. - 4. * t;
196 break;
197 case 7:
198 output.access(4) = 1. - 2. * (s + t);
199 output.access(7) = 1. - 2. * (r + s);
200 output.access(8) = 2. * (r + t) - 1.;
201 aux = 4. * s;
202 break;
203 case 8:
204 output.access(4) = 1. - 2. * (s + t);
205 output.access(5) = 2. * (r + s) - 1.;
206 output.access(6) = 1. - 2. * (r + t);
207 aux = 4. * t;
208 break;
209 case 9:
210 output.access(5) = 2. * (r + s) - 1.;
211 output.access(6) = 1. - 2. * (r + t);
212 output.access(9) = 2. * (s + t) - 1.;
213 aux = 2. - 4. * s;
214 break;
215 case 10:
216 output.access(6) = 1. - 2. * (r + t);
217 output.access(7) = 1. - 2. * (r + s);
218 output.access(9) = 2. * (s + t) - 1.;
219 aux = 4. * r;
220 break;
221 case 11:
222 output.access(4) = 1. - 2. * (s + t);
223 output.access(6) = 1. - 2. * (r + t);
224 output.access(7) = 1. - 2. * (r + s);
225 aux = 4. * (r + s + t) - 2.;
226 break;
227 }
228 for (auto i=4;i<10;++i)
229 output.access(i) += aux/6.0;
230 }
231 break;
232 }
233 case OPERATOR_GRAD: {
234 const typename inputPointViewType::value_type r = input(0);
235 const typename inputPointViewType::value_type s = input(1);
236 const typename inputPointViewType::value_type t = input(2);
237
238 output.access(0,0) = (-17 + 20*r + 20*s + 20*t)/8.;
239 output.access(0,1) = (-17 + 20*r + 20*s + 20*t)/8.;
240 output.access(0,2) = (-17 + 20*r + 20*s + 20*t)/8.;
241 output.access(1,0) = -0.375 + (5*r)/2.;
242 output.access(1,1) = 0.;
243 output.access(1,2) = 0.;
244 output.access(2,0) = 0.;
245 output.access(2,1) = -0.375 + (5*s)/2.;
246 output.access(2,2) = 0.;
247 output.access(3,0) = 0.;
248 output.access(3,1) = 0.;
249 output.access(3,2) = -0.375 + (5*t)/2.;
250 output.access(4,0) = (-35*(-1 + 2*r + s + t))/12.;
251 output.access(4,1) = (-4 - 35*r + 5*s + 10*t)/12.;
252 output.access(4,2) = (-4 - 35*r + 10*s + 5*t)/12.;
253 output.access(5,0) = (-1 + 5*r + 40*s - 5*t)/12.;
254 output.access(5,1) = (-1 + 40*r + 5*s - 5*t)/12.;
255 output.access(5,2) = (-5*(-1 + r + s + 2*t))/12.;
256 output.access(6,0) = (-4 + 5*r - 35*s + 10*t)/12.;
257 output.access(6,1) = (-35*(-1 + r + 2*s + t))/12.;
258 output.access(6,2) = (-4 + 10*r - 35*s + 5*t)/12.;
259 output.access(7,0) = (-4 + 5*r + 10*s - 35*t)/12.;
260 output.access(7,1) = (-4 + 10*r + 5*s - 35*t)/12.;
261 output.access(7,2) = (-35*(-1 + r + s + 2*t))/12.;
262 output.access(8,0) = (-1 + 5*r - 5*s + 40*t)/12.;
263 output.access(8,1) = (-5*(-1 + r + 2*s + t))/12.;
264 output.access(8,2) = (-1 + 40*r - 5*s + 5*t)/12.;
265 output.access(9,0) = (-5*(-1 + 2*r + s + t))/12.;
266 output.access(9,1) = (-1 - 5*r + 5*s + 40*t)/12.;
267 output.access(9,2) = (-1 - 5*r + 40*s + 5*t)/12.;
268 break;
269 }
270 case OPERATOR_MAX: {
271 const ordinal_type jend = output.extent(1);
272 const ordinal_type iend = output.extent(0);
273
274 for (ordinal_type j=0;j<jend;++j)
275 for (auto i=0;i<iend;++i)
276 output.access(i, j) = 0.0;
277 break;
278 }
279 default: {
280 INTREPID2_TEST_FOR_ABORT( true ,
281 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
282 }
283 }
284 }
285
286 template<typename DT,
287 typename outputValueValueType, class ...outputValueProperties,
288 typename inputPointValueType, class ...inputPointProperties>
289 void
290 Basis_HGRAD_TET_COMP12_FEM::
291 getValues( const typename DT::execution_space& space,
292 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
293 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
294 const EOperator operatorType ) {
295 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
296 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
297 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
298
299 // Number of evaluation points = dim 0 of inputPoints
300 const auto loopSize = inputPoints.extent(0);
301 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
302
303 switch (operatorType) {
304
305 case OPERATOR_VALUE: {
306 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
307 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
308 break;
309 }
310 case OPERATOR_GRAD:
311 case OPERATOR_D1: {
312 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
313 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
314 break;
315 }
316 case OPERATOR_CURL: {
317 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
318 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
319 break;
320 }
321 case OPERATOR_DIV: {
322 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
323 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
324 break;
325 }
326 case OPERATOR_D2:
327 case OPERATOR_D3:
328 case OPERATOR_D4:
329 case OPERATOR_D5:
330 case OPERATOR_D6:
331 case OPERATOR_D7:
332 case OPERATOR_D8:
333 case OPERATOR_D9:
334 case OPERATOR_D10: {
335 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
336 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
337 break;
338 }
339 default: {
340 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
341 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
342 }
343 }
344 }
345
346 }
347
348 // -------------------------------------------------------------------------------------
349 template<typename DT, typename OT, typename PT>
352 const ordinal_type spaceDim = 3;
353 this->basisCardinality_ = 10;
354 this->basisDegree_ = 1;
355 this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
356 this->basisType_ = BASIS_FEM_DEFAULT;
357 this->basisCoordinates_ = COORDINATES_CARTESIAN;
358 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
359
360 {
361 // Basis-dependent intializations
362 const ordinal_type tagSize = 4; // size of DoF tag
363 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
364 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
365 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
366
367 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
368 ordinal_type tags[] = { 0, 0, 0, 1,
369 0, 1, 0, 1,
370 0, 2, 0, 1,
371 0, 3, 0, 1,
372 1, 0, 0, 1,
373 1, 1, 0, 1,
374 1, 2, 0, 1,
375 1, 3, 0, 1,
376 1, 4, 0, 1,
377 1, 5, 0, 1 };
378
379 // host view
380 OrdinalTypeArray1DHost tagView(&tags[0],40);
381
382 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
383 this->setOrdinalTagData(this->tagToOrdinal_,
384 this->ordinalToTag_,
385 tagView,
386 this->basisCardinality_,
387 tagSize,
388 posScDim,
389 posScOrd,
390 posDfOrd);
391 }
392
393 // dofCoords on host and create its mirror view to device
394 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
395 dofCoords("dofCoordsHost", this->basisCardinality_, spaceDim);
396
397 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
398 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
399 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
400 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
401 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
402 dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
403 dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
404 dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
405 dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
406 dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
407
408 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
409 Kokkos::deep_copy(this->dofCoords_, dofCoords);
410 }
411
412 template<typename DT, typename OT, typename PT>
413 void
415 ordinal_type& perTeamSpaceSize,
416 ordinal_type& perThreadSpaceSize,
417 const PointViewType inputPoints,
418 const EOperator operatorType) const {
419 perTeamSpaceSize = 0;
420 perThreadSpaceSize = 0;
421 }
422
423 template<typename DT, typename OT, typename PT>
424 KOKKOS_INLINE_FUNCTION
425 void
427 OutputViewType outputValues,
428 const PointViewType inputPoints,
429 const EOperator operatorType,
430 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
431 const typename DT::execution_space::scratch_memory_space & scratchStorage,
432 const ordinal_type subcellDim,
433 const ordinal_type subcellOrdinal) const {
434
435 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
436 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_COMP12_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
437
438 (void) scratchStorage; //avoid unused variable warning
439
440 const int numPoints = inputPoints.extent(0);
441
442 switch(operatorType) {
443 case OPERATOR_VALUE:
444 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
445 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
446 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
448 });
449 break;
450 case OPERATOR_GRAD:
451 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
452 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
453 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
454 Impl::Basis_HGRAD_TET_COMP12_FEM::Serial<OPERATOR_GRAD>::getValues( output, input);
455 });
456 break;
457 default: {}
458 }
459 }
460
461}// namespace Intrepid2
462#endif
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
FEM basis evaluation on a reference Tetrahedron cell.
virtual void getScratchSpaceSize(ordinal_type &perTeamSpaceSize, ordinal_type &perThreadSpaceSize, const PointViewType inputPointsconst, const EOperator operatorType=OPERATOR_VALUE) const override
Return the size of the scratch space, in bytes, needed for using the team-level implementation of get...
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
static KOKKOS_INLINE_FUNCTION ordinal_type getLocalSubTet(const pointValueType x, const pointValueType y, const pointValueType z)
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.