Intrepid2
Intrepid2_HGRAD_QUAD_Cn_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_QUAD_CN_FEM_DEF_HPP__
17 #define __INTREPID2_HGRAD_QUAD_CN_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_HGRAD_QUAD_Cn_FEM::Serial<OpType>::
32  getValues( OutputViewType output,
33  const InputViewType input,
34  WorkViewType work,
35  const VinvViewType vinv,
36  const ordinal_type operatorDn ) {
37  ordinal_type opDn = operatorDn;
38 
39  const ordinal_type cardLine = vinv.extent(0);
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 
46  const int dim_s = get_dimension_scalar(input);
47  auto ptr0 = work.data();
48  auto ptr1 = work.data()+cardLine*npts*dim_s;
49  auto ptr2 = work.data()+2*cardLine*npts*dim_s;
50 
51  typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
52  auto vcprop = Kokkos::common_view_alloc_prop(input);
53 
54  switch (OpType) {
55  case OPERATOR_VALUE: {
56  ViewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
57  ViewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
58  ViewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
59 
60  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
61  getValues(output_x, input_x, work_line, vinv);
62 
63  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
64  getValues(output_y, input_y, work_line, vinv);
65 
66  // tensor product
67  ordinal_type idx = 0;
68  for (ordinal_type j=0;j<cardLine;++j) // y
69  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
70  for (ordinal_type k=0;k<npts;++k)
71  output.access(idx,k) = output_x.access(i,k)*output_y.access(j,k);
72  break;
73  }
74  case OPERATOR_CURL: {
75  for (auto l=0;l<2;++l) {
76  ViewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
77 
78  ViewType output_x, output_y;
79 
80  typename WorkViewType::value_type s = 0.0;
81  if (l) {
82  // l = 1
83  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
84  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
85  getValues(output_x, input_x, work_line, vinv, 1);
86 
87  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
88  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
89  getValues(output_y, input_y, work_line, vinv);
90 
91  s = -1.0;
92  } else {
93  // l = 0
94  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
95  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
96  getValues(output_x, input_x, work_line, vinv);
97 
98  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts, 1);
99  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
100  getValues(output_y, input_y, work_line, vinv, 1);
101 
102  s = 1.0;
103  }
104 
105  // tensor product (extra dimension of ouput x and y are ignored)
106  ordinal_type idx = 0;
107  for (ordinal_type j=0;j<cardLine;++j) // y
108  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
109  for (ordinal_type k=0;k<npts;++k)
110  output.access(idx,k,l) = s*output_x.access(i,k,0)*output_y.access(j,k,0);
111  }
112  break;
113  }
114  case OPERATOR_GRAD:
115  case OPERATOR_D1:
116  case OPERATOR_D2:
117  case OPERATOR_D3:
118  case OPERATOR_D4:
119  case OPERATOR_D5:
120  case OPERATOR_D6:
121  case OPERATOR_D7:
122  case OPERATOR_D8:
123  case OPERATOR_D9:
124  case OPERATOR_D10:
125  opDn = getOperatorOrder(OpType);
126  case OPERATOR_Dn: {
127  const auto dkcard = opDn + 1;
128  for (auto l=0;l<dkcard;++l) {
129  ViewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
130 
131  ViewType output_x, output_y;
132 
133  const auto mult_x = opDn - l;
134  const auto mult_y = l;
135 
136  if (mult_x) {
137  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
138  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
139  getValues(output_x, input_x, work_line, vinv, mult_x);
140  } else {
141  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
142  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
143  getValues(output_x, input_x, work_line, vinv);
144  }
145 
146  if (mult_y) {
147  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts, 1);
148  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
149  getValues(output_y, input_y, work_line, vinv, mult_y);
150  } else {
151  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
152  Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
153  getValues(output_y, input_y, work_line, vinv);
154  }
155 
156  // tensor product (extra dimension of ouput x and y are ignored)
157  ordinal_type idx = 0;
158  for (ordinal_type j=0;j<cardLine;++j) // y
159  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
160  for (ordinal_type k=0;k<npts;++k)
161  output.access(idx,k,l) = output_x.access(i,k,0)*output_y.access(j,k,0);
162  }
163  break;
164  }
165  default: {
166  INTREPID2_TEST_FOR_ABORT( true,
167  ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_Cn_FEM::Serial::getValues) operator is not supported" );
168  }
169  }
170  }
171 
172  template<typename DT, ordinal_type numPtsPerEval,
173  typename outputValueValueType, class ...outputValueProperties,
174  typename inputPointValueType, class ...inputPointProperties,
175  typename vinvValueType, class ...vinvProperties>
176  void
177  Basis_HGRAD_QUAD_Cn_FEM::
178  getValues( const typename DT::execution_space& space,
179  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
180  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
181  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
182  const EOperator operatorType ) {
183  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
184  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
185  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
186  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
187 
188  // loopSize corresponds to cardinality
189  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
190  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
191  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
192  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
193 
194  typedef typename inputPointViewType::value_type inputPointType;
195 
196  const ordinal_type cardinality = outputValues.extent(0);
197  const ordinal_type cardLine = std::sqrt(cardinality);
198  const ordinal_type workSize = 3*cardLine;
199 
200  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
201  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
202  workViewType work(Kokkos::view_alloc(space, "Basis_HGRAD_QUAD_Cn_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
203 
204  switch (operatorType) {
205  case OPERATOR_VALUE: {
206  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
207  OPERATOR_VALUE,numPtsPerEval> FunctorType;
208  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
209  break;
210  }
211  case OPERATOR_CURL: {
212  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
213  OPERATOR_CURL,numPtsPerEval> FunctorType;
214  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
215  break;
216  }
217  case OPERATOR_GRAD:
218  case OPERATOR_D1:
219  case OPERATOR_D2:
220  case OPERATOR_D3:
221  case OPERATOR_D4:
222  case OPERATOR_D5:
223  case OPERATOR_D6:
224  case OPERATOR_D7:
225  case OPERATOR_D8:
226  case OPERATOR_D9:
227  case OPERATOR_D10: {
228  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
229  OPERATOR_Dn,numPtsPerEval> FunctorType;
230  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
231  getOperatorOrder(operatorType)) );
232  break;
233  }
234  default: {
235  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
236  ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented" );
237  // break;commented out because exception
238  }
239  }
240  }
241  }
242 
243  // -------------------------------------------------------------------------------------
244  template<typename DT, typename OT, typename PT>
246  Basis_HGRAD_QUAD_Cn_FEM( const ordinal_type order,
247  const EPointType pointType ) {
248  // INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
249  // pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
250  // ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): pointType must be either equispaced or warpblend." );
251 
252  // this should be in host
253  Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
254  const auto cardLine = lineBasis.getCardinality();
255 
256  this->vinv_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hgrad::Quad::Cn::vinv", cardLine, cardLine);
257  lineBasis.getVandermondeInverse(this->vinv_);
258 
259  const ordinal_type spaceDim = 2;
260  this->basisCardinality_ = cardLine*cardLine;
261  this->basisDegree_ = order;
262  this->basisCellTopologyKey_ = shards::Quadrilateral<4>::key;
263  this->basisType_ = BASIS_FEM_LAGRANGIAN;
264  this->basisCoordinates_ = COORDINATES_CARTESIAN;
265  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
266  pointType_ = pointType;
267 
268  // initialize tags
269  {
270  // Basis-dependent initializations
271  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
272  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
273  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
274  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
275 
276  // 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.
277  // 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.)
278  INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
279  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
280  constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
281  ordinal_type tags[maxCardLine*maxCardLine][4];
282 
283  const ordinal_type vert[2][2] = { {0,1}, {3,2} }; //[y][x]
284 
285  const ordinal_type edge_x[2] = {0,2};
286  const ordinal_type edge_y[2] = {3,1};
287  {
288  ordinal_type idx = 0;
289  for (ordinal_type j=0;j<cardLine;++j) { // y
290  const auto tag_y = lineBasis.getDofTag(j);
291  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
292  const auto tag_x = lineBasis.getDofTag(i);
293 
294  if (tag_x(0) == 0 && tag_y(0) == 0) {
295  // vertices
296  tags[idx][0] = 0; // vertex dof
297  tags[idx][1] = vert[tag_y(1)][tag_x(1)]; // vertex id
298  tags[idx][2] = 0; // local dof id
299  tags[idx][3] = 1; // total number of dofs in this vertex
300  } else if (tag_x(0) == 1 && tag_y(0) == 0) {
301  // edge: x edge, y vert
302  tags[idx][0] = 1; // edge dof
303  tags[idx][1] = edge_x[tag_y(1)];
304  tags[idx][2] = tag_x(2); // local dof id
305  tags[idx][3] = tag_x(3); // total number of dofs in this vertex
306  } else if (tag_x(0) == 0 && tag_y(0) == 1) {
307  // edge: x vert, y edge
308  tags[idx][0] = 1; // edge dof
309  tags[idx][1] = edge_y[tag_x(1)];
310  tags[idx][2] = tag_y(2); // local dof id
311  tags[idx][3] = tag_y(3); // total number of dofs in this vertex
312  } else {
313  // interior
314  tags[idx][0] = 2; // interior dof
315  tags[idx][1] = 0;
316  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
317  tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
318  }
319  }
320  }
321  }
322 
323  OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
324 
325  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
326  // tags are constructed on host
327  this->setOrdinalTagData(this->tagToOrdinal_,
328  this->ordinalToTag_,
329  tagView,
330  this->basisCardinality_,
331  tagSize,
332  posScDim,
333  posScOrd,
334  posDfOrd);
335  }
336 
337  // dofCoords on host and create its mirror view to device
338  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
339  dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
340 
341  Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
342  dofCoordsLine("dofCoordsLine", cardLine, 1);
343 
344  lineBasis.getDofCoords(dofCoordsLine);
345  auto dofCoordsLineHost = Kokkos::create_mirror_view(dofCoordsLine);
346  Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
347  {
348  ordinal_type idx = 0;
349  for (ordinal_type j=0;j<cardLine;++j) { // y
350  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
351  dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
352  dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
353  }
354  }
355  }
356 
357  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
358  Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
359  }
360 
361  template<typename DT, typename OT, typename PT>
362  void
364  ordinal_type& perTeamSpaceSize,
365  ordinal_type& perThreadSpaceSize,
366  const PointViewType inputPoints,
367  const EOperator operatorType) const {
368  perTeamSpaceSize = 0;
369  perThreadSpaceSize = 3*this->vinv_.extent(0)*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
370  }
371 
372  template<typename DT, typename OT, typename PT>
373  KOKKOS_INLINE_FUNCTION
374  void
376  OutputViewType outputValues,
377  const PointViewType inputPoints,
378  const EOperator operatorType,
379  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
380  const typename DT::execution_space::scratch_memory_space & scratchStorage,
381  const ordinal_type subcellDim,
382  const ordinal_type subcellOrdinal) const {
383 
384  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
385  ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_Cn_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
386 
387  const int numPoints = inputPoints.extent(0);
388  using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
389  using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
390  ordinal_type sizePerPoint = 3*this->vinv_.extent(0)*get_dimension_scalar(inputPoints);
391  WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
392  using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
393 
394  switch(operatorType) {
395  case OPERATOR_VALUE:
396  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
397  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
398  const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
399  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
400  Impl::Basis_HGRAD_QUAD_Cn_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinv_ );
401  });
402  break;
403  case OPERATOR_GRAD:
404  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
405  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
406  const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
407  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
408  Impl::Basis_HGRAD_QUAD_Cn_FEM::Serial<OPERATOR_GRAD>::getValues( output, input, work, this->vinv_ );
409  });
410  break;
411  case OPERATOR_CURL:
412  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
413  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
414  const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
415  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
416  Impl::Basis_HGRAD_QUAD_Cn_FEM::Serial<OPERATOR_CURL>::getValues( output, input, work, this->vinv_ );
417  });
418  break;
419  default: {
420  INTREPID2_TEST_FOR_ABORT( true,
421  ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): getValues not implemented for this operator");
422  }
423  }
424  }
425 
426 } // namespace Intrepid2
427 
428 #endif
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
ordinal_type getCardinality() const
Returns cardinality of the basis.
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.
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
Basis_HGRAD_QUAD_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
See Intrepid2::Basis_HGRAD_QUAD_Cn_FEM work is a rank 1 view having the same value_type of inputPoint...
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
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...
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.