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