16#ifndef __INTREPID2_HGRAD_HEX_CN_FEMDEF_HPP__
17#define __INTREPID2_HGRAD_HEX_CN_FEMDEF_HPP__
24 template<EOperator opType>
25 template<
typename OutputViewType,
26 typename inputViewType,
27 typename workViewType,
28 typename vinvViewType>
29 KOKKOS_INLINE_FUNCTION
33 const inputViewType input,
35 const vinvViewType vinv,
36 const ordinal_type operatorDn ) {
37 ordinal_type opDn = operatorDn;
39 const ordinal_type cardLine = vinv.extent(0);
40 const ordinal_type npts = input.extent(0);
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));
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;
53 typedef typename Kokkos::DynRankView<typename inputViewType::value_type, typename workViewType::memory_space> viewType;
56 case OPERATOR_VALUE: {
62 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
63 getValues(output_x, input_x, work_line, vinv);
65 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
66 getValues(output_y, input_y, work_line, vinv);
68 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
69 getValues(output_z, input_z, work_line, vinv);
73 for (ordinal_type k=0;k<cardLine;++k)
74 for (ordinal_type j=0;j<cardLine;++j)
75 for (ordinal_type i=0;i<cardLine;++i,++idx)
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);
93 const ordinal_type dkcard = opDn + 1;
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;
109 decltype(work_line) output_x, output_y, output_z;
113 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
114 getValues(output_x, input_x, work_line, vinv, mult_x);
117 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
118 getValues(output_x, input_x, work_line, vinv);
123 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
124 getValues(output_y, input_y, work_line, vinv, mult_y);
127 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
128 getValues(output_y, input_y, work_line, vinv);
133 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
134 getValues(output_z, input_z, work_line, vinv, mult_z);
137 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
138 getValues(output_z, input_z, work_line, vinv);
142 ordinal_type idx = 0;
143 for (ordinal_type k=0;k<cardLine;++k)
144 for (ordinal_type j=0;j<cardLine;++j)
145 for (ordinal_type i=0;i<cardLine;++i,++idx)
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);
154 INTREPID2_TEST_FOR_ABORT(
true ,
155 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented");
161 template<
typename DT, ordinal_type numPtsPerEval,
162 typename outputValueValueType,
class ...outputValueProperties,
163 typename inputPointValueType,
class ...inputPointProperties,
164 typename vinvValueType,
class ...vinvProperties>
166 Basis_HGRAD_HEX_Cn_FEM::
167 getValues(
const typename DT::execution_space& space,
168 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
169 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
170 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
171 const EOperator operatorType ) {
172 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
173 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
174 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
175 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
178 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
179 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
180 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
181 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
183 const ordinal_type cardinality = outputValues.extent(0);
184 const ordinal_type cardLine = std::cbrt(cardinality);
185 const ordinal_type workSize = 4*cardLine;
187 auto work =
createMatchingDynRankView(inputPoints,
"Basis_HGRAD_HEX_Cn_FEM::getValues::work", workSize, inputPoints.extent(0));
189 switch (operatorType) {
190 case OPERATOR_VALUE: {
191 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,
decltype(work),
192 OPERATOR_VALUE,numPtsPerEval> FunctorType;
193 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
196 case OPERATOR_CURL: {
197 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,
decltype(work),
198 OPERATOR_CURL,numPtsPerEval> FunctorType;
199 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
213 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,
decltype(work),
214 OPERATOR_Dn,numPtsPerEval> FunctorType;
215 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
220 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
221 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented" );
229 template<
typename DT,
typename OT,
typename PT>
232 const EPointType pointType ) {
238 this->
vinv_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>(
"Hgrad::HEX::Cn::vinv", cardLine, cardLine);
239 lineBasis.getVandermondeInverse(this->
vinv_);
241 const ordinal_type spaceDim = 3;
253 const ordinal_type tagSize = 4;
254 const ordinal_type posScDim = 0;
255 const ordinal_type posScOrd = 1;
256 const ordinal_type posDfOrd = 2;
260 INTREPID2_TEST_FOR_EXCEPTION( order >
Parameters::MaxOrder, std::invalid_argument,
"polynomial order exceeds the max supported by this class");
264 ordinal_type tags[maxCardLine*maxCardLine*maxCardLine][4];
266 const ordinal_type vert[2][2][2] = { { {0,1}, {3,2} },
269 const ordinal_type edge_x[2][2] = { {0, 4}, {2, 6} };
270 const ordinal_type edge_y[2][2] = { {3, 7}, {1, 5} };
271 const ordinal_type edge_z[2][2] = { {8,11}, {9,10} };
273 const ordinal_type face_yz[2] = {3, 1};
274 const ordinal_type face_xz[2] = {0, 2};
275 const ordinal_type face_xy[2] = {4, 5};
278 ordinal_type idx = 0;
279 for (
auto k=0;k<cardLine;++k) {
280 const auto tag_z = lineBasis.
getDofTag(k);
281 for (ordinal_type j=0;j<cardLine;++j) {
282 const auto tag_y = lineBasis.
getDofTag(j);
283 for (ordinal_type i=0;i<cardLine;++i,++idx) {
284 const auto tag_x = lineBasis.
getDofTag(i);
286 if (tag_x(0) == 0 && tag_y(0) == 0 && tag_z(0) == 0) {
289 tags[idx][1] = vert[tag_z(1)][tag_y(1)][tag_x(1)];
292 }
else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 0) {
295 tags[idx][1] = edge_x[tag_y(1)][tag_z(1)];
296 tags[idx][2] = tag_x(2);
297 tags[idx][3] = tag_x(3);
298 }
else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 0) {
301 tags[idx][1] = edge_y[tag_x(1)][tag_z(1)];
302 tags[idx][2] = tag_y(2);
303 tags[idx][3] = tag_y(3);
304 }
else if (tag_x(0) == 0 && tag_y(0) == 0 && tag_z(0) == 1) {
307 tags[idx][1] = edge_z[tag_x(1)][tag_y(1)];
308 tags[idx][2] = tag_z(2);
309 tags[idx][3] = tag_z(3);
310 }
else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
313 tags[idx][1] = face_yz[tag_x(1)];
314 tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2);
315 tags[idx][3] = tag_y(3)*tag_z(3);
316 }
else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
319 tags[idx][1] = face_xz[tag_y(1)];
320 tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2);
321 tags[idx][3] = tag_x(3)*tag_z(3);
322 }
else if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
325 tags[idx][1] = face_xy[tag_z(1)];
326 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
327 tags[idx][3] = tag_x(3)*tag_y(3);
332 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2);
333 tags[idx][3] = tag_x(3)*tag_y(3)*tag_z(3);
355 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
358 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
359 dofCoordsLine(
"dofCoordsLine", cardLine, 1);
362 auto dofCoordsLineHost = Kokkos::create_mirror_view(dofCoordsLine);
363 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
365 ordinal_type idx = 0;
366 for (
auto k=0;k<cardLine;++k) {
367 for (ordinal_type j=0;j<cardLine;++j) {
368 for (ordinal_type i=0;i<cardLine;++i,++idx) {
369 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
370 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
371 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
377 this->
dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoordsHost);
378 Kokkos::deep_copy(this->
dofCoords_, dofCoordsHost);
381 template<
typename DT,
typename OT,
typename PT>
383 Basis_HGRAD_HEX_Cn_FEM<DT,OT,PT>::getScratchSpaceSize(
384 ordinal_type& perTeamSpaceSize,
385 ordinal_type& perThreadSpaceSize,
387 const EOperator operatorType)
const {
389 perTeamSpaceSize = 0;
390 perThreadSpaceSize = 4*this->vinv_.extent(0)*get_dimension_scalar(inputPoints)*
sizeof(
typename BasisBase::scalarType);
393 template<
typename DT,
typename OT,
typename PT>
394 KOKKOS_INLINE_FUNCTION
396 Basis_HGRAD_HEX_Cn_FEM<DT,OT,PT>::getValues(
397 OutputViewType outputValues,
398 const PointViewType inputPoints,
399 const EOperator operatorType,
400 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
401 const typename DT::execution_space::scratch_memory_space & scratchStorage,
402 const ordinal_type subcellDim,
403 const ordinal_type subcellOrdinal)
const {
405 INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
406 ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_Cn_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
408 const int numPoints = inputPoints.extent(0);
409 using ScalarType =
typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
410 using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
411 ordinal_type sizePerPoint = 4*this->vinv_.extent(0)*get_dimension_scalar(inputPoints);
412 WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
413 using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
415 switch(operatorType) {
417 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinv_ = this->vinv_] (ordinal_type& pt) {
418 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
419 const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
420 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
421 Impl::Basis_HGRAD_HEX_Cn_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, vinv_ );
425 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinv_ = this->vinv_] (ordinal_type& pt) {
426 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
427 const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
428 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
429 Impl::Basis_HGRAD_HEX_Cn_FEM::Serial<OPERATOR_GRAD>::getValues( output, input, work, vinv_ );
433 INTREPID2_TEST_FOR_ABORT(
true,
434 ">>> ERROR (Basis_HGRAD_TET_Cn_FEM): getValues not implemented for this operator");
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_pointer_v< CtorProp > &&!std::is_convertible_v< CtorProp, constchar * >, OutViewType >::type createMatchingUnmanagedView(const InViewType &view, const CtorProp &data, const Dims... dims)
Creates an unmanaged view that matches the value_type of the provided view The type of the output vie...
DeduceDynRankView< InViewType >::type createMatchingDynRankView(const InViewType &view, const CtorProp &prop, const Dims... dims)
Creates and returns a view that matches the value_type of the provided view The output view type is d...
Basis_HGRAD_HEX_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Kokkos::DynRankView< typename ScalarViewType::value_type, DeviceType > vinv_
inverse of Generalized Vandermonde matrix (isotropic order)
EPointType pointType_
type of lattice used for creating the DoF coordinates
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
ECoordinates basisCoordinates_
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
ordinal_type basisDegree_
void setOrdinalTagData(OrdinalTypeView3D &tagToOrdinal, OrdinalTypeView2D &ordinalToTag, const OrdinalTypeView1D tags, const ordinal_type basisCard, const ordinal_type tagSize, const ordinal_type posScDim, const ordinal_type posScOrd, const ordinal_type posDfOrd)
ordinal_type getCardinality() const
Returns cardinality of the basis.
OrdinalTypeArray2DHost ordinalToTag_
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
ordinal_type basisCardinality_
OrdinalTypeArray3DHost tagToOrdinal_
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
unsigned basisCellTopologyKey_
EFunctionSpace functionSpace_
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
See Intrepid2::Basis_HGRAD_HEX_Cn_FEM.
See Intrepid2::Basis_HGRAD_HEX_Cn_FEM.