Intrepid2
Intrepid2_HCURL_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_HCURL_HEX_IN_FEM_DEF_HPP__
17 #define __INTREPID2_HCURL_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_HCURL_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() + 3*cardLine*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 
59  ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
60  ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
61  ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
62  ViewType outputBubble(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
63 
64  // tensor product
65  ordinal_type idx = 0;
66  {
67  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
68  getValues(outputBubble, input_x, workLine, vinvBubble);
69 
70  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
71  getValues(outputLine_A, input_y, workLine, vinvLine);
72 
73  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
74  getValues(outputLine_B, input_z, workLine, vinvLine);
75 
76  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
77  const auto output_x = outputBubble;
78  const auto output_y = outputLine_A;
79  const auto output_z = outputLine_B;
80 
81  for (ordinal_type k=0;k<cardLine;++k) // z
82  for (ordinal_type j=0;j<cardLine;++j) // y
83  for (ordinal_type i=0;i<cardBubble;++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(outputLine_A, input_x, workLine, vinvLine);
93 
94  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
95  getValues(outputBubble, input_y, workLine, vinvBubble);
96 
97  //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
98  // getValues(outputLine_B, input_z, workLine, vinvLine);
99 
100  // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
101  const auto output_x = outputLine_A;
102  const auto output_y = outputBubble;
103  const auto output_z = outputLine_B;
104 
105  for (ordinal_type k=0;k<cardLine;++k) // z
106  for (ordinal_type j=0;j<cardBubble;++j) // y
107  for (ordinal_type i=0;i<cardLine;++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(outputLine_A, input_x, workLine, vinvLine);
117 
118  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119  getValues(outputLine_B, input_y, workLine, vinvLine);
120 
121  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
122  getValues(outputBubble, input_z, workLine, vinvBubble);
123 
124  // z component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
125  const auto output_x = outputLine_A;
126  const auto output_y = outputLine_B;
127  const auto output_z = outputBubble;
128 
129  for (ordinal_type k=0;k<cardBubble;++k) // z
130  for (ordinal_type j=0;j<cardLine;++j) // y
131  for (ordinal_type i=0;i<cardLine;++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_CURL: {
141 
142  auto ptr4 = work.data() + 4*cardLine*npts*dim_s;
143  auto ptr5 = work.data() + 5*cardLine*npts*dim_s;
144 
145  ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
146  ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
147  ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
148  ViewType outputLine_DA(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts, 1);
149  ViewType outputLine_DB(Kokkos::view_wrap(ptr4, vcprop), cardLine, npts, 1);
150  ViewType outputBubble(Kokkos::view_wrap(ptr5, vcprop), cardBubble, npts);
151 
152  // tensor product
153  ordinal_type idx = 0;
154 
155  { // x - component
156  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
157  getValues(outputBubble, input_x, workLine, vinvBubble);
158 
159  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
160  getValues(outputLine_A, input_y, workLine, vinvLine);
161 
162  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
163  getValues(outputLine_DA, input_y, workLine, vinvLine, 1);
164 
165  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
166  getValues(outputLine_B, input_z, workLine, vinvLine);
167 
168  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
169  getValues(outputLine_DB, input_z, workLine, vinvLine, 1);
170 
171  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
172  const auto output_x = outputBubble;
173  const auto output_y = outputLine_A;
174  const auto output_dy = outputLine_DA;
175  const auto output_z = outputLine_B;
176  const auto output_dz = outputLine_DB;
177 
178  for (ordinal_type k=0;k<cardLine;++k) // z
179  for (ordinal_type j=0;j<cardLine;++j) // y
180  for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
181  for (ordinal_type l=0;l<npts;++l) {
182  output.access(idx,l,0) = 0.0;
183  output.access(idx,l,1) = output_x.access(i,l)*output_y.access (j,l) *output_dz.access(k,l,0);
184  output.access(idx,l,2) = -output_x.access(i,l)*output_dy.access(j,l,0)*output_z.access (k,l);
185  }
186  }
187  { // y - component
188  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
189  getValues(outputLine_A, input_x, workLine, vinvLine);
190 
191  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
192  getValues(outputLine_DA, input_x, workLine, vinvLine, 1);
193 
194  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
195  getValues(outputBubble, input_y, workLine, vinvBubble);
196 
197  //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
198  // getValues(outputLine_B, input_z, workLine, vinvLine);
199 
200  //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
201  // getValues(outputLine_DB, input_z, workLine, vinvLine, 1);
202 
203  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
204  const auto output_x = outputLine_A;
205  const auto output_dx = outputLine_DA;
206  const auto output_y = outputBubble;
207  const auto output_z = outputLine_B;
208  const auto output_dz = outputLine_DB;
209 
210  for (ordinal_type k=0;k<cardLine;++k) // z
211  for (ordinal_type j=0;j<cardBubble;++j) // y
212  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
213  for (ordinal_type l=0;l<npts;++l) {
214  output.access(idx,l,0) = -output_x.access (i,l) *output_y.access(j,l)*output_dz.access(k,l,0);
215  output.access(idx,l,1) = 0.0;
216  output.access(idx,l,2) = output_dx(i,l,0)*output_y.access(j,l)*output_z.access (k,l);
217  }
218  }
219  { // z - component
220  // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
221  // getValues(outputLine_A, input_x, workLine, vinvLine);
222 
223  // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
224  // getValues(outputLine_DA, input_x, workLine, vinvLine, 1);
225 
226  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
227  getValues(outputLine_B, input_y, workLine, vinvLine);
228 
229  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
230  getValues(outputLine_DB, input_y, workLine, vinvLine, 1);
231 
232  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
233  getValues(outputBubble, input_z, workLine, vinvBubble);
234 
235  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
236  const auto output_x = outputLine_A;
237  const auto output_dx = outputLine_DA;
238  const auto output_y = outputLine_B;
239  const auto output_dy = outputLine_DB;
240  const auto output_z = outputBubble;
241 
242  for (ordinal_type k=0;k<cardBubble;++k) // z
243  for (ordinal_type j=0;j<cardLine;++j) // y
244  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
245  for (ordinal_type l=0;l<npts;++l) {
246  output.access(idx,l,0) = output_x.access (i,l) *output_dy.access(j,l,0)*output_z.access(k,l);
247  output.access(idx,l,1) = -output_dx(i,l,0)*output_y.access (j,l) *output_z.access(k,l);
248  output.access(idx,l,2) = 0.0;
249  }
250  }
251  break;
252  }
253  default: {
254  INTREPID2_TEST_FOR_ABORT( true,
255  ">>> ERROR: (Intrepid2::Basis_HCURL_HEX_In_FEM::Serial::getValues) operator is not supported" );
256  }
257  }
258  }
259  template<typename DT, ordinal_type numPtsPerEval,
260  typename outputValueValueType, class ...outputValueProperties,
261  typename inputPointValueType, class ...inputPointProperties,
262  typename vinvValueType, class ...vinvProperties>
263  void
264  Basis_HCURL_HEX_In_FEM::
265  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
266  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
267  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
268  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
269  const EOperator operatorType ) {
270  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
271  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
272  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
273  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
274 
275  // loopSize corresponds to cardinality
276  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
277  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
278  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
279  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
280 
281  typedef typename inputPointViewType::value_type inputPointType;
282 
283  const ordinal_type cardinality = outputValues.extent(0);
284  //get basis order based on basis cardinality.
285  ordinal_type order = 0;
286  ordinal_type cardBubble; // = std::cbrt(cardinality/3);
287  ordinal_type cardLine; // = cardBubble+1;
288  do {
289  cardBubble = Intrepid2::getPnCardinality<1>(order);
290  cardLine = Intrepid2::getPnCardinality<1>(++order);
291  } while((3*cardBubble*cardLine*cardLine != cardinality) && (order != Parameters::MaxOrder));
292  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
293  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
294 
295  switch (operatorType) {
296  case OPERATOR_VALUE: {
297  auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
298  workViewType work(Kokkos::view_alloc("Basis_CURL_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
299  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
300  OPERATOR_VALUE,numPtsPerEval> FunctorType;
301  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
302  break;
303  }
304  case OPERATOR_CURL: {
305  auto workSize = Serial<OPERATOR_CURL>::getWorkSizePerPoint(order);
306  workViewType work(Kokkos::view_alloc("Basis_CURL_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
307  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
308  OPERATOR_CURL,numPtsPerEval> FunctorType;
309  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
310  break;
311  }
312  default: {
313  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
314  ">>> ERROR (Basis_HCURL_HEX_In_FEM): Operator type not implemented" );
315  // break;commented out since exception is thrown
316  }
317  }
318  }
319  }
320 
321  // -------------------------------------------------------------------------------------
322 
323  template<typename DT, typename OT, typename PT>
324  Basis_HCURL_HEX_In_FEM<DT,OT,PT>::
325  Basis_HCURL_HEX_In_FEM( const ordinal_type order,
326  const EPointType pointType ) {
327 
328  INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
329  pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
330  ">>> ERROR (Basis_HCURL_HEX_In_FEM): pointType must be either equispaced or warpblend.");
331 
332  // this should be created in host and vinv should be deep copied into device space
333  Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
334  Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
335 
336  const ordinal_type
337  cardLine = lineBasis.getCardinality(),
338  cardBubble = bubbleBasis.getCardinality();
339 
340  this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvLine", cardLine, cardLine);
341  this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
342 
343  lineBasis.getVandermondeInverse(this->vinvLine_);
344  bubbleBasis.getVandermondeInverse(this->vinvBubble_);
345 
346  const ordinal_type spaceDim = 3;
347  this->basisCardinality_ = 3*cardLine*cardLine*cardBubble;
348  this->basisDegree_ = order;
349  this->basisCellTopologyKey_ = shards::Hexahedron<8>::key;
350  this->basisType_ = BASIS_FEM_LAGRANGIAN;
351  this->basisCoordinates_ = COORDINATES_CARTESIAN;
352  this->functionSpace_ = FUNCTION_SPACE_HCURL;
353  pointType_ = pointType;
354 
355  // initialize tags
356  {
357  // Basis-dependent initializations
358  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
359  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
360  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
361  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
362 
363  // Note: the only reason why equispaced can't support higher order than Parameters::MaxOrder appears to be the fact that the tags below get stored into a fixed-length array.
364  // TODO: relax the maximum order requirement by setting up tags in a different container, perhaps directly into an OrdinalTypeArray1DHost (tagView, below). (As of this writing (1/25/22), looks like other nodal bases do this in a similar way -- those should be fixed at the same time; maybe search for Parameters::MaxOrder.)
365  INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
366 
367  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
368  constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
369  ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
370 
371  const ordinal_type edge_x[2][2] = { {0, 4}, {2, 6} };
372  const ordinal_type edge_y[2][2] = { {3, 7}, {1, 5} };
373  const ordinal_type edge_z[2][2] = { {8,11}, {9,10} };
374 
375  const ordinal_type face_yz[2] = {3, 1};
376  const ordinal_type face_xz[2] = {0, 2};
377  const ordinal_type face_xy[2] = {4, 5};
378 
379  {
380  ordinal_type idx = 0;
381 
385 
386  // since there are x/y components in the interior
387  // dof sum should be computed before the information
388  // is assigned to tags
389  const ordinal_type
390  face_ndofs_per_direction = (cardLine-2)*cardBubble,
391  face_ndofs = 2*face_ndofs_per_direction,
392  intr_ndofs_per_direction = (cardLine-2)*(cardLine-2)*cardBubble,
393  intr_ndofs = 3*intr_ndofs_per_direction;
394 
395  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
396  for (ordinal_type k=0;k<cardLine;++k) { // z
397  const auto tag_z = lineBasis.getDofTag(k);
398  for (ordinal_type j=0;j<cardLine;++j) { // y
399  const auto tag_y = lineBasis.getDofTag(j);
400  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
401  const auto tag_x = bubbleBasis.getDofTag(i);
402 
403  if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 0) {
404  // edge: x edge, y vert, z vert
405  tags[idx][0] = 1; // edge dof
406  tags[idx][1] = edge_x[tag_y(1)][tag_z(1)]; // edge id
407  tags[idx][2] = tag_x(2); // local dof id
408  tags[idx][3] = tag_x(3); // total number of dofs in this edge
409  } else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
410  // face, x edge, y vert, z edge
411  tags[idx][0] = 2; // face dof
412  tags[idx][1] = face_xz[tag_y(1)]; // face id
413  tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2); // local dof id
414  tags[idx][3] = face_ndofs; // total number of dofs in this face
415  } else if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
416  // face, x edge, y edge, z vert
417  tags[idx][0] = 2; // face dof
418  tags[idx][1] = face_xy[tag_z(1)]; // face id
419  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
420  tags[idx][3] = face_ndofs; // total number of dofs in this face
421  } else {
422  // interior
423  tags[idx][0] = 3; // interior dof
424  tags[idx][1] = 0;
425  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
426  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
427  }
428  }
429  }
430  }
431 
432  // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
433  for (ordinal_type k=0;k<cardLine;++k) { // z
434  const auto tag_z = lineBasis.getDofTag(k);
435  for (ordinal_type j=0;j<cardBubble;++j) { // y
436  const auto tag_y = bubbleBasis.getDofTag(j);
437  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
438  const auto tag_x = lineBasis.getDofTag(i);
439 
440  if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 0) {
441  // edge: x vert, y edge, z vert
442  tags[idx][0] = 1; // edge dof
443  tags[idx][1] = edge_y[tag_x(1)][tag_z(1)]; // edge id
444  tags[idx][2] = tag_y(2); // local dof id
445  tags[idx][3] = tag_y(3); // total number of dofs in this vertex
446  } else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
447  // face, x vert, y edge, z edge
448  tags[idx][0] = 2; // face dof
449  tags[idx][1] = face_yz[tag_x(1)]; // face id
450  tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2); // local dof id
451  tags[idx][3] = face_ndofs; // total number of dofs in this vertex
452  } else if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
453  // face, x edge, y edge, z vert
454  tags[idx][0] = 2; // face dof
455  tags[idx][1] = face_xy[tag_z(1)]; // face id
456  tags[idx][2] = face_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2); // local dof id
457  tags[idx][3] = face_ndofs; // total number of dofs in this vertex
458  } else {
459  // interior
460  tags[idx][0] = 3; // interior dof
461  tags[idx][1] = 0;
462  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
463  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
464  }
465  }
466  }
467  }
468 
469  // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
470  for (ordinal_type k=0;k<cardBubble;++k) { // y
471  const auto tag_z = bubbleBasis.getDofTag(k);
472  for (ordinal_type j=0;j<cardLine;++j) { // z
473  const auto tag_y = lineBasis.getDofTag(j);
474  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
475  const auto tag_x = lineBasis.getDofTag(i);
476 
477  if (tag_x(0) == 0 && tag_y(0) == 0 && tag_z(0) == 1) {
478  // edge: x vert, y vert, z edge
479  tags[idx][0] = 1; // edge dof
480  tags[idx][1] = edge_z[tag_x(1)][tag_y(1)]; // edge id
481  tags[idx][2] = tag_z(2); // local dof id
482  tags[idx][3] = tag_z(3); // total number of dofs in this vertex
483  } else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
484  // face, x vert, y edge, z edge
485  tags[idx][0] = 2; // face dof
486  tags[idx][1] = face_yz[tag_x(1)]; // face id
487  tags[idx][2] = face_ndofs_per_direction + tag_y(2) + tag_y(3)*tag_z(2); // local dof id
488  tags[idx][3] = face_ndofs; // total number of dofs in this vertex
489  } else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
490  // face, x edge, y vert, z edge
491  tags[idx][0] = 2; // face dof
492  tags[idx][1] = face_xz[tag_y(1)]; // face id
493  tags[idx][2] = face_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_z(2); // local dof id
494  tags[idx][3] = face_ndofs; // total number of dofs in this vertex
495  } else {
496  // interior
497  tags[idx][0] = 3; // interior dof
498  tags[idx][1] = 0;
499  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
500  tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
501  }
502  }
503  }
504  }
505 
506  INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
507  ">>> ERROR (Basis_HCURL_HEX_In_FEM): " \
508  "counted tag index is not same as cardinality." );
509  }
510 
511  OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
512 
513  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
514  // tags are constructed on host
515  this->setOrdinalTagData(this->tagToOrdinal_,
516  this->ordinalToTag_,
517  tagView,
518  this->basisCardinality_,
519  tagSize,
520  posScDim,
521  posScOrd,
522  posDfOrd);
523  }
524 
525  // dofCoords on host and create its mirror view to device
526  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
527  dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
528 
529  // dofCoords on host and create its mirror view to device
530  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
531  dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, spaceDim);
532 
533  Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
534  dofCoordsLine("dofCoordsLine", cardLine, 1),
535  dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
536 
537  lineBasis.getDofCoords(dofCoordsLine);
538  auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
539  Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
540 
541  bubbleBasis.getDofCoords(dofCoordsBubble);
542  auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
543  Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
544 
545  {
546  ordinal_type idx = 0;
547 
548  // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
549  for (ordinal_type k=0;k<cardLine;++k) { // z
550  for (ordinal_type j=0;j<cardLine;++j) { // y
551  for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
552  dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
553  dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
554  dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
555  dofCoeffsHost(idx,0) = 1.0;
556  }
557  }
558  }
559 
560  // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
561  for (ordinal_type k=0;k<cardLine;++k) { // z
562  for (ordinal_type j=0;j<cardBubble;++j) { // y
563  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
564  dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
565  dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
566  dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
567  dofCoeffsHost(idx,1) = 1.0;
568  }
569  }
570  }
571 
572  // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
573  for (ordinal_type k=0;k<cardBubble;++k) { // z
574  for (ordinal_type j=0;j<cardLine;++j) { // y
575  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
576  dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
577  dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
578  dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
579  dofCoeffsHost(idx,2) = 1.0;
580  }
581  }
582  }
583  }
584 
585  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
586  Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
587 
588  this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
589  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
590  }
591 
592  template<typename DT, typename OT, typename PT>
593  void
594  Basis_HCURL_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
595  ordinal_type& perTeamSpaceSize,
596  ordinal_type& perThreadSpaceSize,
597  const PointViewType inputPoints,
598  const EOperator operatorType) const {
599  perTeamSpaceSize = 0;
600  ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ?
601  3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
602  5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
603  perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
604  }
605 
606  template<typename DT, typename OT, typename PT>
607  KOKKOS_INLINE_FUNCTION
608  void
609  Basis_HCURL_HEX_In_FEM<DT,OT,PT>::getValues(
610  OutputViewType outputValues,
611  const PointViewType inputPoints,
612  const EOperator operatorType,
613  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
614  const typename DT::execution_space::scratch_memory_space & scratchStorage,
615  const ordinal_type subcellDim,
616  const ordinal_type subcellOrdinal) const {
617 
618  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
619  ">>> ERROR: (Intrepid2::Basis_HCURL_HEX_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
620 
621  const int numPoints = inputPoints.extent(0);
622  using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
623  using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
624  ordinal_type scalarSizePerPoint = (operatorType == OPERATOR_VALUE) ?
625  3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
626  5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
627  ordinal_type sizePerPoint = scalarSizePerPoint*get_dimension_scalar(inputPoints);
628  WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
629  using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
630 
631  switch(operatorType) {
632  case OPERATOR_VALUE:
633  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
634  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
635  const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
636  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
637  Impl::Basis_HCURL_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
638  });
639  break;
640  case OPERATOR_CURL:
641  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
642  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
643  const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
644  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
645  Impl::Basis_HCURL_HEX_In_FEM::Serial<OPERATOR_CURL>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
646  });
647  break;
648  default: {
649  INTREPID2_TEST_FOR_ABORT( true,
650  ">>> ERROR (Basis_HCURL_HEX_In_FEM): getValues not implemented for this operator");
651  }
652  }
653  }
654 
655 } // namespace Intrepid2
656 
657 #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_HCURL_HEX_In_FEM
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Definition: Intrepid2_HCURL_HEX_In_FEM.hpp:148
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_HCURL_HEX_In_FEM::Basis_HCURL_HEX_In_FEM
Basis_HCURL_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HCURL_HEX_In_FEMDef.hpp:325
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:102
Generated by   doxygen 1.8.5