16#ifndef __INTREPID2_HDIV_QUAD_IN_FEM_DEF_HPP__
17#define __INTREPID2_HDIV_QUAD_IN_FEM_DEF_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 vinvLine,
36 const VinvViewType vinvBubble) {
37 const ordinal_type cardLine = vinvLine.extent(0);
38 const ordinal_type cardBubble = vinvBubble.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));
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;
51 typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
54 case OPERATOR_VALUE: {
62 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
63 getValues(outputBubble, input_x, workLine, vinvBubble);
65 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
66 getValues(outputLine, input_y, workLine, vinvLine);
69 const auto output_x = outputBubble;
70 const auto output_y = outputLine;
72 for (ordinal_type j=0;j<cardLine;++j)
73 for (ordinal_type i=0;i<cardBubble;++i,++idx)
74 for (ordinal_type k=0;k<npts;++k) {
75 output.access(idx,k,0) = 0.0;
76 output.access(idx,k,1) = output_x.access(i,k)*output_y.access(j,k);
80 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
81 getValues(outputBubble, input_y, workLine, vinvBubble);
83 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
84 getValues(outputLine, input_x, workLine, vinvLine);
87 const auto output_x = outputLine;
88 const auto output_y = outputBubble;
89 for (ordinal_type j=0;j<cardBubble;++j)
90 for (ordinal_type i=0;i<cardLine;++i,++idx)
91 for (ordinal_type k=0;k<npts;++k) {
92 output.access(idx,k,0) = output_x.access(i,k)*output_y.access(j,k);
93 output.access(idx,k,1) = 0.0;
107 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
108 getValues(output_x, input_x, workLine, vinvBubble);
110 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
111 getValues(output_y, input_y, workLine, vinvLine, 1);
114 for (ordinal_type j=0;j<cardLine;++j)
115 for (ordinal_type i=0;i<cardBubble;++i,++idx)
116 for (ordinal_type k=0;k<npts;++k)
117 output.access(idx,k) = output_x.access(i,k)*output_y.access(j,k,0);
126 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
127 getValues(output_y, input_y, workLine, vinvBubble);
129 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
130 getValues(output_x, input_x, workLine, vinvLine, 1);
133 for (ordinal_type j=0;j<cardBubble;++j)
134 for (ordinal_type i=0;i<cardLine;++i,++idx)
135 for (ordinal_type k=0;k<npts;++k)
136 output.access(idx,k) = output_x.access(i,k,0)*output_y.access(j,k);
141 INTREPID2_TEST_FOR_ABORT(
true,
142 ">>> ERROR: (Intrepid2::Basis_HDIV_QUAD_In_FEM::Serial::getValues) operator is not supported" );
147 template<
typename DT, ordinal_type numPtsPerEval,
148 typename outputValueValueType,
class ...outputValueProperties,
149 typename inputPointValueType,
class ...inputPointProperties,
150 typename vinvValueType,
class ...vinvProperties>
152 Basis_HDIV_QUAD_In_FEM::
153 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
154 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
155 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
156 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
157 const EOperator operatorType ) {
158 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
159 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
160 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
161 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
164 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
165 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
166 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
167 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
171 const ordinal_type cardinality = outputValues.extent(0);
173 ordinal_type order = 0;
174 ordinal_type cardBubble;
175 ordinal_type cardLine;
177 cardBubble = Intrepid2::getPnCardinality<1>(order);
178 cardLine = Intrepid2::getPnCardinality<1>(++order);
181 switch (operatorType) {
182 case OPERATOR_VALUE: {
183 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
184 auto work =
createMatchingDynRankView(inputPoints,
"Basis_HDIV_QUAD_In_FEM::getValues::work", workSize, inputPoints.extent(0));
185 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,
decltype(work),
186 OPERATOR_VALUE,numPtsPerEval> FunctorType;
187 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
191 auto workSize = Serial<OPERATOR_DIV>::getWorkSizePerPoint(order);
192 auto work =
createMatchingDynRankView(inputPoints,
"Basis_HDIV_QUAD_In_FEM::getValues::work", workSize, inputPoints.extent(0));
193 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,
decltype(work),
194 OPERATOR_DIV,numPtsPerEval> FunctorType;
195 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
199 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
200 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Operator type not implemented" );
207 template<
typename DT,
typename OT,
typename PT>
210 const EPointType pointType ) {
212 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
213 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
214 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): pointType must be either equispaced or warpblend.");
224 this->
vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>(
"Hdiv::Quad::In::vinvLine", cardLine, cardLine);
225 this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>(
"Hdiv::Quad::In::vinvBubble", cardBubble, cardBubble);
227 lineBasis.getVandermondeInverse(this->
vinvLine_);
228 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
230 const ordinal_type spaceDim = 2;
237 pointType_ = pointType;
242 const ordinal_type tagSize = 4;
243 const ordinal_type posScDim = 0;
244 const ordinal_type posScOrd = 1;
245 const ordinal_type posDfOrd = 2;
250 ordinal_type tags[2*maxCardLine*maxCardBubble][4];
252 const ordinal_type edge_x[2] = {0,2};
253 const ordinal_type edge_y[2] = {3,1};
255 ordinal_type idx = 0;
265 intr_ndofs_per_direction = (cardLine-2)*cardBubble,
266 intr_ndofs = 2*intr_ndofs_per_direction;
269 for (ordinal_type j=0;j<cardLine;++j) {
270 const auto tag_y = lineBasis.
getDofTag(j);
271 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
272 const auto tag_x = bubbleBasis.
getDofTag(i);
274 if (tag_x(0) == 1 && tag_y(0) == 0) {
277 tags[idx][1] = edge_x[tag_y(1)];
278 tags[idx][2] = tag_x(2);
279 tags[idx][3] = tag_x(3);
284 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
285 tags[idx][3] = intr_ndofs;
291 for (ordinal_type j=0;j<cardBubble;++j) {
292 const auto tag_y = bubbleBasis.
getDofTag(j);
293 for (ordinal_type i=0;i<cardLine;++i,++idx) {
294 const auto tag_x = lineBasis.
getDofTag(i);
296 if (tag_x(0) == 0 && tag_y(0) == 1) {
299 tags[idx][1] = edge_y[tag_x(1)];
300 tags[idx][2] = tag_y(2);
301 tags[idx][3] = tag_y(3);
306 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2);
307 tags[idx][3] = intr_ndofs;
311 INTREPID2_TEST_FOR_EXCEPTION( idx != this->
basisCardinality_ , std::runtime_error,
312 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): " \
313 "counted tag index is not same as cardinality." );
331 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
335 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
338 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
339 dofCoordsLine(
"dofCoordsLine", cardLine, 1),
340 dofCoordsBubble(
"dofCoordsBubble", cardBubble, 1);
343 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
344 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
347 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
348 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
351 ordinal_type idx = 0;
354 for (ordinal_type j=0;j<cardLine;++j) {
355 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
356 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
357 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
358 dofCoeffsHost(idx,1) = 1.0;
363 for (ordinal_type j=0;j<cardBubble;++j) {
364 for (ordinal_type i=0;i<cardLine;++i,++idx) {
365 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
366 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
367 dofCoeffsHost(idx,0) = 1.0;
372 this->
dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoordsHost);
373 Kokkos::deep_copy(this->
dofCoords_, dofCoordsHost);
375 this->
dofCoeffs_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoeffsHost);
376 Kokkos::deep_copy(this->
dofCoeffs_, dofCoeffsHost);
379 template<
typename DT,
typename OT,
typename PT>
381 Basis_HDIV_QUAD_In_FEM<DT,OT,PT>::getScratchSpaceSize(
382 ordinal_type& perTeamSpaceSize,
383 ordinal_type& perThreadSpaceSize,
385 const EOperator operatorType)
const {
386 perTeamSpaceSize = 0;
387 perThreadSpaceSize = (2*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints)*
sizeof(
typename BasisBase::scalarType);
390 template<
typename DT,
typename OT,
typename PT>
391 KOKKOS_INLINE_FUNCTION
393 Basis_HDIV_QUAD_In_FEM<DT,OT,PT>::getValues(
394 OutputViewType outputValues,
395 const PointViewType inputPoints,
396 const EOperator operatorType,
397 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
398 const typename DT::execution_space::scratch_memory_space & scratchStorage,
399 const ordinal_type subcellDim,
400 const ordinal_type subcellOrdinal)
const {
402 INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
403 ">>> ERROR: (Intrepid2::Basis_HDIV_QUAD_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
405 const int numPoints = inputPoints.extent(0);
406 using ScalarType =
typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
407 using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
408 ordinal_type sizePerPoint = (2*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints);
409 WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
410 using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
412 switch(operatorType) {
414 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
415 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
416 const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
417 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
418 Impl::Basis_HDIV_QUAD_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, vinvLine_, vinvBubble_ );
422 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
423 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
424 const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
425 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
426 Impl::Basis_HDIV_QUAD_In_FEM::Serial<OPERATOR_DIV>::getValues( output, input, work, vinvLine_, vinvBubble_ );
430 INTREPID2_TEST_FOR_ABORT(
true,
431 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): getValues not implemented for this 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_HDIV_QUAD_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Kokkos::DynRankView< typename ScalarViewType::value_type, DeviceType > vinvLine_
inverse of Generalized Vandermonde matrix (isotropic order)
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.
Implementation of the locally HVOL-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_
Kokkos::DynRankView< scalarType, DeviceType > dofCoeffs_
ordinal_type basisCardinality_
OrdinalTypeArray3DHost tagToOrdinal_
unsigned basisCellTopologyKey_
EFunctionSpace functionSpace_
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
See Intrepid2::Basis_HDIV_QUAD_In_FEM.
See Intrepid2::Basis_HDIV_QUAD_In_FEM.