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 
19 namespace 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.8.5