Intrepid2
Intrepid2_HVOL_HEX_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 
15 #ifndef __INTREPID2_HVOL_HEX_CN_FEMDEF_HPP__
16 #define __INTREPID2_HVOL_HEX_CN_FEMDEF_HPP__
17 
18 namespace Intrepid2 {
19 
20  // -------------------------------------------------------------------------------------
21  namespace Impl {
22 
23  template<EOperator OpType>
24  template<typename OutputViewType,
25  typename InputViewType,
26  typename WorkViewType,
27  typename VinvViewType>
28  KOKKOS_INLINE_FUNCTION
29  void
30  Basis_HVOL_HEX_Cn_FEM::Serial<OpType>::
31  getValues( OutputViewType output,
32  const InputViewType input,
33  WorkViewType work,
34  const VinvViewType vinv,
35  const ordinal_type operatorDn ) {
36  ordinal_type opDn = operatorDn;
37 
38  const ordinal_type cardLine = vinv.extent(0);
39  const ordinal_type npts = input.extent(0);
40 
41  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
42  const auto input_x = Kokkos::subview(input, Kokkos::ALL(), range_type(0,1));
43  const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
44  const auto input_z = Kokkos::subview(input, Kokkos::ALL(), range_type(2,3));
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  auto ptr3 = work.data()+3*cardLine*npts*dim_s;
51 
52  typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
53  auto vcprop = Kokkos::common_view_alloc_prop(input);
54 
55  switch (OpType) {
56  case OPERATOR_VALUE: {
57  ViewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
58  ViewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
59  ViewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
60  ViewType output_z(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts);
61 
62  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
63  getValues(output_x, input_x, work_line, vinv);
64 
65  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
66  getValues(output_y, input_y, work_line, vinv);
67 
68  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
69  getValues(output_z, input_z, work_line, vinv);
70 
71  // tensor product
72  ordinal_type idx = 0;
73  for (ordinal_type k=0;k<cardLine;++k) // z
74  for (ordinal_type j=0;j<cardLine;++j) // y
75  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
76  for (ordinal_type l=0;l<npts;++l)
77  output.access(idx,l) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
78  break;
79  }
80  case OPERATOR_GRAD:
81  case OPERATOR_D1:
82  case OPERATOR_D2:
83  case OPERATOR_D3:
84  case OPERATOR_D4:
85  case OPERATOR_D5:
86  case OPERATOR_D6:
87  case OPERATOR_D7:
88  case OPERATOR_D8:
89  case OPERATOR_D9:
90  case OPERATOR_D10:
91  opDn = getOperatorOrder(OpType);
92  case OPERATOR_Dn: {
93  const ordinal_type dkcard = opDn + 1;
94 
95  ordinal_type d = 0;
96  for (ordinal_type l1=0;l1<dkcard;++l1)
97  for (ordinal_type l0=0;l0<(l1+1);++l0) {
98  const ordinal_type mult_x = (opDn - l1);
99  const ordinal_type mult_y = l1 - l0;
100  const ordinal_type mult_z = l0;
101 
102  //std::cout << " l0, l1 = " << l0 << " " << l1 << std::endl;
103  //std::cout << " x , y , z = " << mult_x << " " << mult_y << " " << mult_z << std::endl;
104 
105  if (mult_x < 0) {
106  // pass
107  } else {
108  ViewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
109  decltype(work_line) output_x, output_y, output_z;
110 
111  if (mult_x) {
112  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
113  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
114  getValues(output_x, input_x, work_line, vinv, mult_x);
115  } else {
116  output_x = ViewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
117  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
118  getValues(output_x, input_x, work_line, vinv);
119  }
120 
121  if (mult_y) {
122  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts, 1);
123  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
124  getValues(output_y, input_y, work_line, vinv, mult_y);
125  } else {
126  output_y = ViewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
127  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
128  getValues(output_y, input_y, work_line, vinv);
129  }
130 
131  if (mult_z) {
132  output_z = ViewType(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts, 1);
133  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
134  getValues(output_z, input_z, work_line, vinv, mult_z);
135  } else {
136  output_z = ViewType(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts);
137  Impl::Basis_HVOL_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
138  getValues(output_z, input_z, work_line, vinv);
139  }
140 
141  // tensor product (extra dimension of ouput x,y and z are ignored)
142  ordinal_type idx = 0;
143  for (ordinal_type k=0;k<cardLine;++k) // z
144  for (ordinal_type j=0;j<cardLine;++j) // y
145  for (ordinal_type i=0;i<cardLine;++i,++idx) // x
146  for (ordinal_type l=0;l<npts;++l)
147  output.access(idx,l,d) = output_x.access(i,l,0)*output_y.access(j,l,0)*output_z.access(k,l,0);
148  ++d;
149  }
150  }
151  break;
152  }
153  default: {
154  INTREPID2_TEST_FOR_ABORT( true ,
155  ">>> ERROR (Basis_HVOL_HEX_Cn_FEM): Operator type not implemented");
156  break;
157  }
158  }
159  }
160 
161  template<typename DT, ordinal_type numPtsPerEval,
162  typename outputValueValueType, class ...outputValueProperties,
163  typename inputPointValueType, class ...inputPointProperties,
164  typename vinvValueType, class ...vinvProperties>
165  void
166  Basis_HVOL_HEX_Cn_FEM::
167  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
168  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
169  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
170  const EOperator operatorType ) {
171  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
172  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
173  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
174  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
175 
176  // loopSize corresponds to cardinality
177  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
178  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
179  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
180  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
181 
182  typedef typename inputPointViewType::value_type inputPointType;
183 
184  const ordinal_type cardinality = outputValues.extent(0);
185  const ordinal_type cardLine = std::cbrt(cardinality);
186  const ordinal_type workSize = 4*cardLine;
187 
188  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
189  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
190  workViewType work(Kokkos::view_alloc("Basis_HVOL_HEX_Cn_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
191 
192  switch (operatorType) {
193  case OPERATOR_VALUE: {
194  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
195  OPERATOR_VALUE,numPtsPerEval> FunctorType;
196  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
197  break;
198  }
199  case OPERATOR_GRAD:
200  case OPERATOR_D1:
201  case OPERATOR_D2:
202  case OPERATOR_D3:
203  case OPERATOR_D4:
204  case OPERATOR_D5:
205  case OPERATOR_D6:
206  case OPERATOR_D7:
207  case OPERATOR_D8:
208  case OPERATOR_D9:
209  case OPERATOR_D10: {
210  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
211  OPERATOR_Dn,numPtsPerEval> FunctorType;
212  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
213  getOperatorOrder(operatorType)) );
214  break;
215  }
216  default: {
217  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
218  ">>> ERROR (Basis_HVOL_HEX_Cn_FEM): Operator type not implemented" );
219  // break; commented out since exception is thrown
220  }
221  }
222  }
223  }
224 
225  // -------------------------------------------------------------------------------------
226  template<typename DT, typename OT, typename PT>
227  Basis_HVOL_HEX_Cn_FEM<DT,OT,PT>::
228  Basis_HVOL_HEX_Cn_FEM( const ordinal_type order,
229  const EPointType pointType ) {
230 
231  // this should be in host
232  Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
233  const auto cardLine = lineBasis.getCardinality();
234 
235  this->pointType_ = pointType;
236  this->vinv_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("HVOL::HEX::Cn::vinv", cardLine, cardLine);
237  lineBasis.getVandermondeInverse(this->vinv_);
238 
239  const ordinal_type spaceDim = 3;
240  this->basisCardinality_ = cardLine*cardLine*cardLine;
241  this->basisDegree_ = order;
242  this->basisCellTopologyKey_ = shards::Hexahedron<8>::key;
243  this->basisType_ = BASIS_FEM_LAGRANGIAN;
244  this->basisCoordinates_ = COORDINATES_CARTESIAN;
245  this->functionSpace_ = FUNCTION_SPACE_HVOL;
246 
247  // initialize tags
248  {
249  // Basis-dependent initializations
250  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
251  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
252  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
253  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
254 
255  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
256  constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
257  ordinal_type tags[maxCardLine*maxCardLine*maxCardLine][4];
258 
259  {
260  ordinal_type idx = 0;
261  for (auto k=0;k<cardLine;++k) { // z
262  const auto tag_z = lineBasis.getDofTag(k);
263  for (ordinal_type j=0;j<cardLine;++j) { // y
264  const auto tag_y = lineBasis.getDofTag(j);
265  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
266  const auto tag_x = lineBasis.getDofTag(i);
267 
268  // interior
269  tags[idx][0] = 3; // interior dof
270  tags[idx][1] = 0;
271  tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
272  tags[idx][3] = tag_x(3)*tag_y(3)*tag_z(3); // total number of dofs in this vertex
273  }
274  }
275  }
276  }
277 
278  OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
279 
280  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
281  // tags are constructed on host
282  this->setOrdinalTagData(this->tagToOrdinal_,
283  this->ordinalToTag_,
284  tagView,
285  this->basisCardinality_,
286  tagSize,
287  posScDim,
288  posScOrd,
289  posDfOrd);
290  }
291 
292  // dofCoords on host and create its mirror view to device
293  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
294  dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
295 
296  Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
297  dofCoordsLine("dofCoordsLine", cardLine, 1);
298 
299  lineBasis.getDofCoords(dofCoordsLine);
300  auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
301  Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
302  {
303  ordinal_type idx = 0;
304  for (auto k=0;k<cardLine;++k) { // z
305  for (ordinal_type j=0;j<cardLine;++j) { // y
306  for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
307  dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
308  dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
309  dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
310  }
311  }
312  }
313 
314  }
315 
316  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
317  Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
318  }
319 
320  template<typename DT, typename OT, typename PT>
321  void
322  Basis_HVOL_HEX_Cn_FEM<DT,OT,PT>::getScratchSpaceSize(
323  ordinal_type& perTeamSpaceSize,
324  ordinal_type& perThreadSpaceSize,
325  const PointViewType inputPoints,
326  const EOperator operatorType) const {
327  perTeamSpaceSize = 0;
328  perThreadSpaceSize = 4*this->vinv_.extent(0)*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
329  }
330 
331  template<typename DT, typename OT, typename PT>
332  KOKKOS_INLINE_FUNCTION
333  void
334  Basis_HVOL_HEX_Cn_FEM<DT,OT,PT>::getValues(
335  OutputViewType outputValues,
336  const PointViewType inputPoints,
337  const EOperator operatorType,
338  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
339  const typename DT::execution_space::scratch_memory_space & scratchStorage,
340  const ordinal_type subcellDim,
341  const ordinal_type subcellOrdinal) const {
342 
343  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
344  ">>> ERROR: (Intrepid2::Basis_HVOL_HEX_Cn_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
345 
346  const int numPoints = inputPoints.extent(0);
347  using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
348  using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
349  auto sizePerPoint = 4*this->vinv_.extent(0)*get_dimension_scalar(inputPoints);
350  WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
351  using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
352  switch(operatorType) {
353  case OPERATOR_VALUE:
354  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
355  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
356  const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
357  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
358  Impl::Basis_HVOL_HEX_Cn_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinv_, this->basisDegree_);
359  });
360  break;
361  default: {
362  INTREPID2_TEST_FOR_ABORT( true,
363  ">>> ERROR (Basis_HVOL_HEX_Cn_FEM): getValues not implemented for this operator");
364  }
365  }
366  }
367 
368 } // namespace Intrepid2
369 
370 #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< DeviceType, outputValueType, pointValueType >::scalarType
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
Definition: Intrepid2_Basis.hpp:157
Intrepid2::Basis_HVOL_HEX_Cn_FEM::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:125
Intrepid2::Basis_HVOL_HEX_Cn_FEM::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition: Intrepid2_Basis.hpp:338
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_HVOL_HEX_Cn_FEM::getScratchSpaceSize
virtual void getScratchSpaceSize(ordinal_type &perTeamSpaceSize, ordinal_type &perThreadSpaceSize, const PointViewType inputPoints, 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...
Definition: Intrepid2_HVOL_HEX_Cn_FEMDef.hpp:322
Intrepid2::Impl::Basis_HVOL_HEX_Cn_FEM::Serial
See Intrepid2::Basis_HVOL_HEX_Cn_FEM.
Definition: Intrepid2_HVOL_HEX_Cn_FEM.hpp:34
Intrepid2::Basis_HVOL_HEX_Cn_FEM::Basis_HVOL_HEX_Cn_FEM
Basis_HVOL_HEX_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HVOL_HEX_Cn_FEMDef.hpp:228
Intrepid2::Basis_HVOL_HEX_Cn_FEM::getValues
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Definition: Intrepid2_HVOL_HEX_Cn_FEM.hpp:159
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:102
Generated by   doxygen 1.8.5