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 
20 namespace 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.
static KOKKOS_INLINE_FUNCTION ordinal_type getLocalSubTet(const pointValueType x, const pointValueType y, const pointValueType z)
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.
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::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.