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
31 Basis_HDIV_QUAD_In_FEM::Serial<OpType>::
32 getValues( OutputViewType output,
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;
52 auto vcprop = Kokkos::common_view_alloc_prop(input);
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);
63 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
64 getValues(outputBubble, input_x, workLine, vinvBubble);
66 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
67 getValues(outputLine, input_y, workLine, vinvLine);
70 const auto output_x = outputBubble;
71 const auto output_y = outputLine;
73 for (ordinal_type j=0;j<cardLine;++j)
74 for (ordinal_type i=0;i<cardBubble;++i,++idx)
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);
81 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
82 getValues(outputBubble, input_y, workLine, vinvBubble);
84 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
85 getValues(outputLine, input_x, workLine, vinvLine);
88 const auto output_x = outputLine;
89 const auto output_y = outputBubble;
90 for (ordinal_type j=0;j<cardBubble;++j)
91 for (ordinal_type i=0;i<cardLine;++i,++idx)
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;
100 ordinal_type idx = 0;
102 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
104 ViewType output_x(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
106 ViewType output_y(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
108 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
109 getValues(output_x, input_x, workLine, vinvBubble);
111 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
112 getValues(output_y, input_y, workLine, vinvLine, 1);
115 for (ordinal_type j=0;j<cardLine;++j)
116 for (ordinal_type i=0;i<cardBubble;++i,++idx)
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);
121 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
123 ViewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
125 ViewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
127 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
128 getValues(output_y, input_y, workLine, vinvBubble);
130 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
131 getValues(output_x, input_x, workLine, vinvLine, 1);
134 for (ordinal_type j=0;j<cardBubble;++j)
135 for (ordinal_type i=0;i<cardLine;++i,++idx)
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);
142 INTREPID2_TEST_FOR_ABORT(
true,
143 ">>> ERROR: (Intrepid2::Basis_HDIV_QUAD_In_FEM::Serial::getValues) operator is not supported" );
148 template<
typename DT, ordinal_type numPtsPerEval,
149 typename outputValueValueType,
class ...outputValueProperties,
150 typename inputPointValueType,
class ...inputPointProperties,
151 typename vinvValueType,
class ...vinvProperties>
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;
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);
170 typedef typename inputPointViewType::value_type inputPointType;
172 const ordinal_type cardinality = outputValues.extent(0);
174 ordinal_type order = 0;
175 ordinal_type cardBubble;
176 ordinal_type cardLine;
178 cardBubble = Intrepid2::getPnCardinality<1>(order);
179 cardLine = Intrepid2::getPnCardinality<1>(++order);
182 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
183 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
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) );
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) );
203 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
204 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Operator type not implemented" );
211 template<
typename DT,
typename OT,
typename PT>
214 const EPointType pointType ) {
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.");
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);
231 lineBasis.getVandermondeInverse(this->vinvLine_);
232 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
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;
246 const ordinal_type tagSize = 4;
247 const ordinal_type posScDim = 0;
248 const ordinal_type posScOrd = 1;
249 const ordinal_type posDfOrd = 2;
254 ordinal_type tags[2*maxCardLine*maxCardBubble][4];
256 const ordinal_type edge_x[2] = {0,2};
257 const ordinal_type edge_y[2] = {3,1};
259 ordinal_type idx = 0;
269 intr_ndofs_per_direction = (cardLine-2)*cardBubble,
270 intr_ndofs = 2*intr_ndofs_per_direction;
273 for (ordinal_type j=0;j<cardLine;++j) {
274 const auto tag_y = lineBasis.
getDofTag(j);
275 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
276 const auto tag_x = bubbleBasis.
getDofTag(i);
278 if (tag_x(0) == 1 && tag_y(0) == 0) {
281 tags[idx][1] = edge_x[tag_y(1)];
282 tags[idx][2] = tag_x(2);
283 tags[idx][3] = tag_x(3);
288 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
289 tags[idx][3] = intr_ndofs;
295 for (ordinal_type j=0;j<cardBubble;++j) {
296 const auto tag_y = bubbleBasis.
getDofTag(j);
297 for (ordinal_type i=0;i<cardLine;++i,++idx) {
298 const auto tag_x = lineBasis.
getDofTag(i);
300 if (tag_x(0) == 0 && tag_y(0) == 1) {
303 tags[idx][1] = edge_y[tag_x(1)];
304 tags[idx][2] = tag_y(2);
305 tags[idx][3] = tag_y(3);
310 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2);
311 tags[idx][3] = intr_ndofs;
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." );
320 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
324 this->setOrdinalTagData(this->tagToOrdinal_,
327 this->basisCardinality_,
335 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
336 dofCoordsHost(
"dofCoordsHost", this->basisCardinality_, spaceDim);
339 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
340 dofCoeffsHost(
"dofCoeffsHost", this->basisCardinality_, spaceDim);
342 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
343 dofCoordsLine(
"dofCoordsLine", cardLine, 1),
344 dofCoordsBubble(
"dofCoordsBubble", cardBubble, 1);
347 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
348 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
351 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
352 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
355 ordinal_type idx = 0;
358 for (ordinal_type j=0;j<cardLine;++j) {
359 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
360 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
361 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
362 dofCoeffsHost(idx,1) = 1.0;
367 for (ordinal_type j=0;j<cardBubble;++j) {
368 for (ordinal_type i=0;i<cardLine;++i,++idx) {
369 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
370 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
371 dofCoeffsHost(idx,0) = 1.0;
376 this->dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoordsHost);
377 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
379 this->dofCoeffs_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoeffsHost);
380 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
383 template<
typename DT,
typename OT,
typename PT>
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);
394 template<
typename DT,
typename OT,
typename PT>
395 KOKKOS_INLINE_FUNCTION
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 {
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.");
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>;
416 switch(operatorType) {
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_ );
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_ );
434 INTREPID2_TEST_FOR_ABORT(
true,
435 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): getValues not implemented for this operator");
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.