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