49 #ifndef __INTREPID2_HDIV_QUAD_IN_FEM_DEF_HPP__
50 #define __INTREPID2_HDIV_QUAD_IN_FEM_DEF_HPP__
57 template<EOperator opType>
58 template<
typename outputViewType,
59 typename inputViewType,
60 typename workViewType,
61 typename vinvViewType>
62 KOKKOS_INLINE_FUNCTION
64 Basis_HDIV_QUAD_In_FEM::Serial<opType>::
65 getValues( outputViewType output,
66 const inputViewType input,
68 const vinvViewType vinvLine,
69 const vinvViewType vinvBubble) {
70 const ordinal_type cardLine = vinvLine.extent(0);
71 const ordinal_type cardBubble = vinvBubble.extent(0);
73 const ordinal_type npts = input.extent(0);
75 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
76 const auto input_x = Kokkos::subview(input, Kokkos::ALL(), range_type(0,1));
77 const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
79 const int dim_s = get_dimension_scalar(work);
80 auto ptr0 = work.data();
81 auto ptr1 = work.data()+cardLine*npts*dim_s;
82 auto ptr2 = work.data()+2*cardLine*npts*dim_s;
85 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
86 auto vcprop = Kokkos::common_view_alloc_prop(work);
89 case OPERATOR_VALUE: {
90 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
91 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
92 viewType outputBubble(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
97 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
98 getValues(outputBubble, input_x, workLine, vinvBubble);
100 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
101 getValues(outputLine, input_y, workLine, vinvLine);
104 const auto output_x = outputBubble;
105 const auto output_y = outputLine;
107 for (ordinal_type j=0;j<cardLine;++j)
108 for (ordinal_type i=0;i<cardBubble;++i,++idx)
109 for (ordinal_type k=0;k<npts;++k) {
110 output.access(idx,k,0) = 0.0;
111 output.access(idx,k,1) = output_x.access(i,k)*output_y.access(j,k);
115 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
116 getValues(outputBubble, input_y, workLine, vinvBubble);
118 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119 getValues(outputLine, input_x, workLine, vinvLine);
122 const auto output_x = outputLine;
123 const auto output_y = outputBubble;
124 for (ordinal_type j=0;j<cardBubble;++j)
125 for (ordinal_type i=0;i<cardLine;++i,++idx)
126 for (ordinal_type k=0;k<npts;++k) {
127 output.access(idx,k,0) = output_x.access(i,k)*output_y.access(j,k);
128 output.access(idx,k,1) = 0.0;
134 ordinal_type idx = 0;
136 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
138 viewType output_x(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
140 viewType output_y(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
142 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
143 getValues(output_x, input_x, workLine, vinvBubble);
145 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
146 getValues(output_y, input_y, workLine, vinvLine, 1);
149 for (ordinal_type j=0;j<cardLine;++j)
150 for (ordinal_type i=0;i<cardBubble;++i,++idx)
151 for (ordinal_type k=0;k<npts;++k)
152 output.access(idx,k) = output_x.access(i,k)*output_y.access(j,k,0);
155 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
157 viewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
159 viewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
161 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
162 getValues(output_y, input_y, workLine, vinvBubble);
164 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
165 getValues(output_x, input_x, workLine, vinvLine, 1);
168 for (ordinal_type j=0;j<cardBubble;++j)
169 for (ordinal_type i=0;i<cardLine;++i,++idx)
170 for (ordinal_type k=0;k<npts;++k)
171 output.access(idx,k) = output_x.access(i,k,0)*output_y.access(j,k);
176 INTREPID2_TEST_FOR_ABORT(
true,
177 ">>> ERROR: (Intrepid2::Basis_HDIV_QUAD_In_FEM::Serial::getValues) operator is not supported" );
182 template<
typename SpT, ordinal_type numPtsPerEval,
183 typename outputValueValueType,
class ...outputValueProperties,
184 typename inputPointValueType,
class ...inputPointProperties,
185 typename vinvValueType,
class ...vinvProperties>
187 Basis_HDIV_QUAD_In_FEM::
188 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
189 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
190 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
191 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
192 const EOperator operatorType ) {
193 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
194 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
195 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
196 typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
199 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
200 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
201 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
202 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
204 typedef typename inputPointViewType::value_type inputPointType;
206 const ordinal_type cardinality = outputValues.extent(0);
208 ordinal_type order = 0;
209 ordinal_type cardBubble;
210 ordinal_type cardLine;
212 cardBubble = Intrepid2::getPnCardinality<1>(order);
213 cardLine = Intrepid2::getPnCardinality<1>(++order);
216 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
217 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
219 switch (operatorType) {
220 case OPERATOR_VALUE: {
221 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
222 workViewType work(Kokkos::view_alloc(
"Basis_HDIV_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
223 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
224 OPERATOR_VALUE,numPtsPerEval> FunctorType;
225 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
229 auto workSize = Serial<OPERATOR_DIV>::getWorkSizePerPoint(order);
230 workViewType work(Kokkos::view_alloc(
"Basis_HDIV_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
231 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
232 OPERATOR_DIV,numPtsPerEval> FunctorType;
233 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
237 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
238 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Operator type not implemented" );
245 template<
typename SpT,
typename OT,
typename PT>
248 const EPointType pointType ) {
250 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
251 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
252 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): pointType must be either equispaced or warpblend.");
262 this->vinvLine_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hdiv::Quad::In::vinvLine", cardLine, cardLine);
263 this->vinvBubble_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hdiv::Quad::In::vinvBubble", cardBubble, cardBubble);
265 lineBasis.getVandermondeInverse(this->vinvLine_);
266 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
268 this->basisCardinality_ = 2*cardLine*cardBubble;
269 this->basisDegree_ = order;
270 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
271 this->basisType_ = BASIS_FEM_FIAT;
272 this->basisCoordinates_ = COORDINATES_CARTESIAN;
277 const ordinal_type tagSize = 4;
278 const ordinal_type posScDim = 0;
279 const ordinal_type posScOrd = 1;
280 const ordinal_type posDfOrd = 2;
285 ordinal_type tags[2*maxCardLine*maxCardBubble][4];
287 const ordinal_type edge_x[2] = {0,2};
288 const ordinal_type edge_y[2] = {3,1};
290 ordinal_type idx = 0;
300 intr_ndofs_per_direction = (cardLine-2)*cardBubble,
301 intr_ndofs = 2*intr_ndofs_per_direction;
304 for (ordinal_type j=0;j<cardLine;++j) {
305 const auto tag_y = lineBasis.
getDofTag(j);
306 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
307 const auto tag_x = bubbleBasis.
getDofTag(i);
309 if (tag_x(0) == 1 && tag_y(0) == 0) {
312 tags[idx][1] = edge_x[tag_y(1)];
313 tags[idx][2] = tag_x(2);
314 tags[idx][3] = tag_x(3);
319 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
320 tags[idx][3] = intr_ndofs;
326 for (ordinal_type j=0;j<cardBubble;++j) {
327 const auto tag_y = bubbleBasis.
getDofTag(j);
328 for (ordinal_type i=0;i<cardLine;++i,++idx) {
329 const auto tag_x = lineBasis.
getDofTag(i);
331 if (tag_x(0) == 0 && tag_y(0) == 1) {
334 tags[idx][1] = edge_y[tag_x(1)];
335 tags[idx][2] = tag_y(2);
336 tags[idx][3] = tag_y(3);
341 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2);
342 tags[idx][3] = intr_ndofs;
346 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
347 ">>> ERROR (Basis_HDIV_QUAD_In_FEM): " \
348 "counted tag index is not same as cardinality." );
355 this->setOrdinalTagData(this->tagToOrdinal_,
358 this->basisCardinality_,
366 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
367 dofCoordsHost(
"dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
370 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
371 dofCoeffsHost(
"dofCoeffsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
373 Kokkos::DynRankView<typename scalarViewType::value_type,SpT>
374 dofCoordsLine(
"dofCoordsLine", cardLine, 1),
375 dofCoordsBubble(
"dofCoordsBubble", cardBubble, 1);
378 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
379 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
382 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
383 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
386 ordinal_type idx = 0;
389 for (ordinal_type j=0;j<cardLine;++j) {
390 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
391 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
392 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
393 dofCoeffsHost(idx,1) = 1.0;
398 for (ordinal_type j=0;j<cardBubble;++j) {
399 for (ordinal_type i=0;i<cardLine;++i,++idx) {
400 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
401 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
402 dofCoeffsHost(idx,0) = 1.0;
407 this->dofCoords_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoordsHost);
408 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
410 this->dofCoeffs_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoeffsHost);
411 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
ordinal_type getCardinality() const
Returns cardinality of the basis.
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
const ordinal_type_array_stride_1d_host getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Basis_HDIV_QUAD_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
virtual void getDofCoords(scalarViewType dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.