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