Intrepid2
Intrepid2_HDIV_HEX_In_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_HDIV_HEX_IN_FEM_DEF_HPP__
17 #define __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
18 
19 namespace Intrepid2 {
20 
21  // -------------------------------------------------------------------------------------
22  namespace Impl {
23 
24  template<EOperator OpType>
25  template<typename OutputViewType,
26  typename InputViewType,
27  typename WorkViewType,
28  typename VinvViewType>
29  KOKKOS_INLINE_FUNCTION
30  void
31  Basis_HDIV_HEX_In_FEM::Serial<OpType>::
32  getValues( OutputViewType output,
33  const InputViewType input,
34  WorkViewType work,
35  const VinvViewType vinvLine,
36  const VinvViewType vinvBubble) {
37  const ordinal_type cardLine = vinvLine.extent(0);
38  const ordinal_type cardBubble = vinvBubble.extent(0);
39 
40  const ordinal_type npts = input.extent(0);
41 
42  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
43  const auto input_x = Kokkos::subview(input, Kokkos::ALL(), range_type(0,1));
44  const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
45  const auto input_z = Kokkos::subview(input, Kokkos::ALL(), range_type(2,3));
46 
47  const ordinal_type dim_s = get_dimension_scalar(input);
48  auto ptr0 = work.data();
49  auto ptr1 = work.data()+cardLine*npts*dim_s;
50  auto ptr2 = work.data()+2*cardLine*npts*dim_s;
51  auto ptr3 = work.data()+(2*cardLine+cardBubble)*npts*dim_s;
52 
53  typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
54  auto vcprop = Kokkos::common_view_alloc_prop(input);
55 
56  switch (OpType) {
57  case OPERATOR_VALUE: {
58  ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
59  ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
60  ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
61  ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
62 
63  // tensor product
64  ordinal_type idx = 0;
65  {
66  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
67  getValues(outputLine, input_x, workLine, vinvLine);
68 
69  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
70  getValues(outputBubble_A, input_y, workLine, vinvBubble);
71 
72  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
73  getValues(outputBubble_B, input_z, workLine, vinvBubble);
74 
75 
76  // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
77  const auto output_x = outputLine;
78  const auto output_y = outputBubble_A;
79  const auto output_z = outputBubble_B;
80 
81  for (ordinal_type k=0;k<cardBubble;++k) // z
82  for (ordinal_type j=0;j<cardBubble;++j) // y
83  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
84  for (ordinal_type l=0;l<npts;++l) {
85  output.access(idx,l,0) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
86  output.access(idx,l,1) = 0.0;
87  output.access(idx,l,2) = 0.0;
88  }
89  }
90  {
91  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
92  getValues(outputBubble_A, input_x, workLine, vinvBubble);
93 
94  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
95  getValues(outputLine, input_y, workLine, vinvLine);
96 
97  //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
98  // getValues(outputBubble_B, input_z, workLine, vinvBubble);
99 
100  // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
101  const auto output_x = outputBubble_A;
102  const auto output_y = outputLine;
103  const auto output_z = outputBubble_B;
104 
105  for (ordinal_type k=0;k<cardBubble;++k) // z
106  for (ordinal_type j=0;j<cardLine;++j) // y
107  for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
108  for (ordinal_type l=0;l<npts;++l) {
109  output.access(idx,l,0) = 0.0;
110  output.access(idx,l,1) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
111  output.access(idx,l,2) = 0.0;
112  }
113  }
114  {
115  //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
116  // getValues(outputBubble_A, input_x, workLine, vinvBubble);
117 
118  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119  getValues(outputBubble_B, input_y, workLine, vinvBubble);
120 
121  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
122  getValues(outputLine, input_z, workLine, vinvLine);
123 
124  // z component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
125  const auto output_x = outputBubble_A;
126  const auto output_y = outputBubble_B;
127  const auto output_z = outputLine;
128 
129  for (ordinal_type k=0;k<cardLine;++k) // z
130  for (ordinal_type j=0;j<cardBubble;++j) // y
131  for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
132  for (ordinal_type l=0;l<npts;++l) {
133  output.access(idx,l,0) = 0.0;
134  output.access(idx,l,1) = 0.0;
135  output.access(idx,l,2) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
136  }
137  }
138  break;
139  }
140  case OPERATOR_DIV: {
141  ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
142  // A line value
143  ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
144  // B line value
145  ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
146  // Line grad
147  ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
148 
149  // tensor product
150  ordinal_type idx = 0;
151 
152  { // x - component
153  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
154  getValues(outputLine, input_x, workLine, vinvLine, 1);
155 
156  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
157  getValues(outputBubble_A, input_y, workLine, vinvBubble);
158 
159  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
160  getValues(outputBubble_B, input_z, workLine, vinvBubble);
161 
162  // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
163  const auto output_dx = outputLine;
164  const auto output_y = outputBubble_A;
165  const auto output_z = outputBubble_B;
166 
167  for (ordinal_type k=0;k<cardBubble;++k) // z
168  for (ordinal_type j=0;j<cardBubble;++j) // y
169  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
170  for (ordinal_type l=0;l<npts;++l)
171  output.access(idx,l) = output_dx.access(i,l,0)*output_y.access (j,l) *output_z.access(k,l);
172  }
173  { // y - component
174  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
175  getValues(outputBubble_A, input_x, workLine, vinvBubble);
176 
177  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
178  getValues(outputLine, input_y, workLine, vinvLine, 1);
179 
180  // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
181  // getValues(outputBubble_B, input_z, workLine, vinvBubble);
182 
183  //(bubbleBasis(z) lineBasis(y) bubbleBasis(x))
184  const auto output_x = outputBubble_A;
185  const auto output_dy = outputLine;
186  const auto output_z = outputBubble_B;
187 
188  for (ordinal_type k=0;k<cardBubble;++k) // z
189  for (ordinal_type j=0;j<cardLine;++j) // y
190  for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
191  for (ordinal_type l=0;l<npts;++l)
192  output.access(idx,l) = output_x.access(i,l)*output_dy.access(j,l,0)*output_z.access(k,l);
193  }
194  { // z - component
195  // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
196  // getValues(outputBubble_A, input_x, workLine, vinvBubble);
197 
198  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
199  getValues(outputBubble_B, input_y, workLine, vinvBubble);
200 
201  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
202  getValues(outputLine, input_z, workLine, vinvLine, 1);
203 
204  // (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
205  const auto output_x = outputBubble_A;
206  const auto output_y = outputBubble_B;
207  const auto output_dz = outputLine;
208 
209  for (ordinal_type k=0;k<cardLine;++k) // z
210  for (ordinal_type j=0;j<cardBubble;++j) // y
211  for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
212  for (ordinal_type l=0;l<npts;++l)
213  output.access(idx,l) = output_x.access(i,l)*output_y.access(j,l)*output_dz.access(k,l,0);
214  }
215  break;
216  }
217  default: {
218  INTREPID2_TEST_FOR_ABORT( true,
219  ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::Serial::getValues) operator is not supported" );
220  }
221  }
222  }
223 
224  template<typename DT, ordinal_type numPtsPerEval,
225  typename outputValueValueType, class ...outputValueProperties,
226  typename inputPointValueType, class ...inputPointProperties,
227  typename vinvValueType, class ...vinvProperties>
228  void
229  Basis_HDIV_HEX_In_FEM::
230  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
231  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
232  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
233  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
234  const EOperator operatorType ) {
235  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
236  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
237  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
238  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
239 
240  // loopSize corresponds to cardinality
241  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
242  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
243  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
244  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
245 
246  typedef typename inputPointViewType::value_type inputPointType;
247 
248  const ordinal_type cardinality = outputValues.extent(0);
249  //get basis order based on basis cardinality.
250  ordinal_type order = 0;
251  ordinal_type cardBubble; // = std::cbrt(cardinality/3);
252  ordinal_type cardLine; // = cardBubble+1;
253  do {
254  cardBubble = Intrepid2::getPnCardinality<1>(order);
255  cardLine = Intrepid2::getPnCardinality<1>(++order);
256  } while((3*cardBubble*cardBubble*cardLine != cardinality) && (order != Parameters::MaxOrder));
257 
258  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
259  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
260 
261  switch (operatorType) {
262  case OPERATOR_VALUE: {
263  auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
264  workViewType work(Kokkos::view_alloc("Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
265  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
266  OPERATOR_VALUE,numPtsPerEval> FunctorType;
267  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
268  break;
269  }
270  case OPERATOR_DIV: {
271  auto workSize = Serial<OPERATOR_DIV>::getWorkSizePerPoint(order);
272  workViewType work(Kokkos::view_alloc("Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
273  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
274  OPERATOR_DIV,numPtsPerEval> FunctorType;
275  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
276  break;
277  }
278  default: {
279  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
280  ">>> ERROR (Basis_HDIV_HEX_In_FEM): Operator type not implemented" );
281  // break; commented out since exception is thrown
282  }
283  }
284  }
285  }
286 
287  // -------------------------------------------------------------------------------------
288 
289  template<typename DT, typename OT, typename PT>
290  Basis_HDIV_HEX_In_FEM<DT,OT,PT>::
291  Basis_HDIV_HEX_In_FEM( const ordinal_type order,
292  const EPointType pointType ) {
293 
294  INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
295  pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
296  ">>> ERROR (Basis_HDIV_HEX_In_FEM): pointType must be either equispaced or warpblend.");
297 
298  // this should be created in host and vinv should be deep copied into device space
299  Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
300  Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
301 
302  const ordinal_type
303  cardLine = lineBasis.getCardinality(),
304  cardBubble = bubbleBasis.getCardinality();
305 
306  this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvLine", cardLine, cardLine);
307  this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
308 
309  lineBasis.getVandermondeInverse(this->vinvLine_);
310  bubbleBasis.getVandermondeInverse(this->vinvBubble_);
311 
312  const ordinal_type spaceDim = 3;
313  this->basisCardinality_ = 3*cardLine*cardBubble*cardBubble;
314  this->basisDegree_ = order;
315  this->basisCellTopologyKey_ = shards::Hexahedron<8>::key;
316  this->basisType_ = BASIS_FEM_LAGRANGIAN;
317  this->basisCoordinates_ = COORDINATES_CARTESIAN;
318  this->functionSpace_ = FUNCTION_SPACE_HDIV;
319  pointType_ = pointType;
320 
321  // initialize tags
322  {
323  // Basis-dependent initializations
324  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
325  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
326  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
327  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
328 
329  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
330  constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
331  ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
332 
333  const ordinal_type face_yz[2] = {3, 1};
334  const ordinal_type face_xz[2] = {0, 2};
335  const ordinal_type face_xy[2] = {4, 5};
336 
337  {
338  ordinal_type idx = 0;
339 
343 
344  // since there are x/y components in the interior
345  // dof sum should be computed before the information
346  // is assigned to tags
347  const ordinal_type
348  intr_ndofs_per_direction = (cardLine-2)*cardBubble*cardBubble,
349  intr_ndofs = 3*intr_ndofs_per_direction;
350 
351  // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
352  for (ordinal_type k=0;k<cardBubble;++k) { // z
353  const auto tag_z = bubbleBasis.getDofTag(k);
354  for (ordinal_type j=0;j<cardBubble;++j) { // y
355  const auto tag_y = bubbleBasis.getDofTag(j);
356  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
357  const auto tag_x = lineBasis.getDofTag(i);
358 
359  if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
360  // face, x vert, y edge, z edge
361  tags[idx][0] = 2; // face dof
362  tags[idx][1] = face_yz[tag_x(1)]; // face id
363  tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2); // local dof id
364  tags[idx][3] = tag_y(3)*tag_z(3); // total number of dofs in this vertex
365  } else {
366  // interior
367  tags[idx][0] = 3; // interior dof
368  tags[idx][1] = 0;
369  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
370  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
371  }
372  }
373  }
374  }
375 
376  // y component (bubbleBasis(z) lineBasis(y) bubbleBasis(x))
377  for (ordinal_type k=0;k<cardBubble;++k) { // z
378  const auto tag_z = bubbleBasis.getDofTag(k);
379  for (ordinal_type j=0;j<cardLine;++j) { // y
380  const auto tag_y = lineBasis.getDofTag(j);
381  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
382  const auto tag_x = bubbleBasis.getDofTag(i);
383 
384  if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
385  // face, x edge, y vert, z edge
386  tags[idx][0] = 2; // face dof
387  tags[idx][1] = face_xz[tag_y(1)]; // face id
388  tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2); // local dof id
389  tags[idx][3] = tag_x(3)*tag_z(3); // total number of dofs in this vertex
390  } else {
391  // interior
392  tags[idx][0] = 3; // interior dof
393  tags[idx][1] = 0;
394  tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
395  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
396  }
397  }
398  }
399  }
400 
401  // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
402  for (ordinal_type k=0;k<cardLine;++k) { // y
403  const auto tag_z = lineBasis.getDofTag(k);
404  for (ordinal_type j=0;j<cardBubble;++j) { // z
405  const auto tag_y = bubbleBasis.getDofTag(j);
406  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
407  const auto tag_x = bubbleBasis.getDofTag(i);
408 
409  if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
410  // face, x edge, y edge, z vert
411  tags[idx][0] = 2; // face dof
412  tags[idx][1] = face_xy[tag_z(1)]; // face id
413  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
414  tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
415  } else {
416  // interior
417  tags[idx][0] = 3; // interior dof
418  tags[idx][1] = 0;
419  tags[idx][2] = 2*intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
420  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
421  }
422  }
423  }
424  }
425 
426  INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
427  ">>> ERROR (Basis_HDIV_HEX_In_FEM): " \
428  "counted tag index is not same as cardinality." );
429  }
430 
431  OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
432 
433  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
434  // tags are constructed on host
435  this->setOrdinalTagData(this->tagToOrdinal_,
436  this->ordinalToTag_,
437  tagView,
438  this->basisCardinality_,
439  tagSize,
440  posScDim,
441  posScOrd,
442  posDfOrd);
443  }
444 
445  // dofCoords on host and create its mirror view to device
446  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
447  dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
448 
449  // dofCoeffs on host and create its mirror view to device
450  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
451  dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, spaceDim);
452 
453  Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
454  dofCoordsLine("dofCoordsLine", cardLine, 1),
455  dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
456 
457  lineBasis.getDofCoords(dofCoordsLine);
458  auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
459  Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
460 
461  bubbleBasis.getDofCoords(dofCoordsBubble);
462  auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
463  Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
464 
465  {
466  ordinal_type idx = 0;
467 
468  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
469  for (ordinal_type k=0;k<cardBubble;++k) { // z
470  for (ordinal_type j=0;j<cardBubble;++j) { // y
471  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
472  dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
473  dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
474  dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
475  dofCoeffsHost(idx,0) = 1.0;
476  }
477  }
478  }
479 
480  // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
481  for (ordinal_type k=0;k<cardBubble;++k) { // z
482  for (ordinal_type j=0;j<cardLine;++j) { // y
483  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
484  dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
485  dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
486  dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
487  dofCoeffsHost(idx,1) = 1.0;
488  }
489  }
490  }
491 
492  // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
493  for (ordinal_type k=0;k<cardLine;++k) { // z
494  for (ordinal_type j=0;j<cardBubble;++j) { // y
495  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
496  dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
497  dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
498  dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
499  dofCoeffsHost(idx,2) = 1.0;
500  }
501  }
502  }
503  }
504 
505  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
506  Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
507 
508  this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
509  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
510  }
511 
512  template<typename DT, typename OT, typename PT>
513  void
514  Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
515  ordinal_type& perTeamSpaceSize,
516  ordinal_type& perThreadSpaceSize,
517  const PointViewType inputPoints,
518  const EOperator operatorType) const {
519  perTeamSpaceSize = 0;
520  perThreadSpaceSize = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
521  }
522 
523  template<typename DT, typename OT, typename PT>
524  KOKKOS_INLINE_FUNCTION
525  void
526  Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getValues(
527  OutputViewType outputValues,
528  const PointViewType inputPoints,
529  const EOperator operatorType,
530  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
531  const typename DT::execution_space::scratch_memory_space & scratchStorage,
532  const ordinal_type subcellDim,
533  const ordinal_type subcellOrdinal) const {
534 
535  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
536  ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
537 
538  const int numPoints = inputPoints.extent(0);
539  using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
540  using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
541  ordinal_type sizePerPoint = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints);
542  WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
543  using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
544 
545  switch(operatorType) {
546  case OPERATOR_VALUE:
547  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
548  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
549  const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
550  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
551  Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
552  });
553  break;
554  case OPERATOR_DIV:
555  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
556  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
557  const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
558  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
559  Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_DIV>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
560  });
561  break;
562  default: {
563  INTREPID2_TEST_FOR_ABORT( true,
564  ">>> ERROR (Basis_HDIV_HEX_In_FEM): getValues not implemented for this operator");
565  }
566  }
567  }
568 
569 } // namespace Intrepid2
570 
571 #endif
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::getDofTag
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
Definition: Intrepid2_Basis.hpp:941
Intrepid2::Basis_HVOL_LINE_Cn_FEM::getDofCoords
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Definition: Intrepid2_HVOL_LINE_Cn_FEM.hpp:205
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::getCardinality
ordinal_type getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid2_Basis.hpp:805
Intrepid2::Basis_HDIV_HEX_In_FEM::Basis_HDIV_HEX_In_FEM
Basis_HDIV_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HDIV_HEX_In_FEMDef.hpp:291
Intrepid2::Basis_HGRAD_LINE_Cn_FEM::getDofCoords
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Definition: Intrepid2_HGRAD_LINE_Cn_FEM.hpp:214
Intrepid2::Basis_HGRAD_LINE_Cn_FEM
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Definition: Intrepid2_HGRAD_LINE_Cn_FEM.hpp:143
Intrepid2::Basis_HVOL_LINE_Cn_FEM
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,1] reference line cell, using Lagrange polynomials.
Definition: Intrepid2_HVOL_LINE_Cn_FEM.hpp:144
Intrepid2::Basis_HDIV_HEX_In_FEM
Implementation of the default H(div)-compatible FEM basis on Hexahedron cell.
Definition: Intrepid2_HDIV_HEX_In_FEM.hpp:138
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:102
Generated by   doxygen 1.8.5