49 #ifndef __INTREPID2_HCURL_QUAD_IN_FEM_DEF_HPP__
50 #define __INTREPID2_HCURL_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_HCURL_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;
84 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
85 auto vcprop = Kokkos::common_view_alloc_prop(work);
88 case OPERATOR_VALUE: {
89 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
90 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
91 viewType outputBubble(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
96 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
97 getValues(outputBubble, input_x, workLine, vinvBubble);
99 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
100 getValues(outputLine, input_y, workLine, vinvLine);
103 const auto output_x = outputBubble;
104 const auto output_y = outputLine;
106 for (ordinal_type j=0;j<cardLine;++j)
107 for (ordinal_type i=0;i<cardBubble;++i,++idx)
108 for (ordinal_type k=0;k<npts;++k) {
109 output.access(idx,k,0) = output_x.access(i,k)*output_y.access(j,k);
110 output.access(idx,k,1) = 0.0;
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) = 0.0;
128 output.access(idx,k,1) = output_x.access(i,k)*output_y.access(j,k);
134 case OPERATOR_CURL: {
135 ordinal_type idx = 0;
137 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
139 viewType output_x(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
141 viewType output_y(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
143 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
144 getValues(output_x, input_x, workLine, vinvBubble);
146 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
147 getValues(output_y, input_y, workLine, vinvLine, 1);
150 for (ordinal_type j=0;j<cardLine;++j)
151 for (ordinal_type i=0;i<cardBubble;++i,++idx)
152 for (ordinal_type k=0;k<npts;++k)
153 output.access(idx,k) = -output_x.access(i,k)*output_y.access(j,k,0);
156 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
158 viewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
160 viewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
162 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
163 getValues(output_y, input_y, workLine, vinvBubble);
165 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
166 getValues(output_x, input_x, workLine, vinvLine, 1);
169 for (ordinal_type j=0;j<cardBubble;++j)
170 for (ordinal_type i=0;i<cardLine;++i,++idx)
171 for (ordinal_type k=0;k<npts;++k)
172 output.access(idx,k) = output_x.access(i,k,0)*output_y.access(j,k);
177 INTREPID2_TEST_FOR_ABORT(
true,
178 ">>> ERROR: (Intrepid2::Basis_HCURL_QUAD_In_FEM::Serial::getValues) operator is not supported" );
183 template<
typename SpT, ordinal_type numPtsPerEval,
184 typename outputValueValueType,
class ...outputValueProperties,
185 typename inputPointValueType,
class ...inputPointProperties,
186 typename vinvValueType,
class ...vinvProperties>
188 Basis_HCURL_QUAD_In_FEM::
189 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
190 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
191 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
192 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
193 const EOperator operatorType ) {
194 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
195 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
196 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
197 typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
200 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
201 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
202 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
203 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
205 typedef typename inputPointViewType::value_type inputPointType;
207 const ordinal_type cardinality = outputValues.extent(0);
209 ordinal_type order = 0;
210 ordinal_type cardBubble;
211 ordinal_type cardLine;
213 cardBubble = Intrepid2::getPnCardinality<1>(order);
214 cardLine = Intrepid2::getPnCardinality<1>(++order);
217 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
218 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
220 switch (operatorType) {
221 case OPERATOR_VALUE: {
222 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
223 workViewType work(Kokkos::view_alloc(
"Basis_HCURL_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
224 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
225 OPERATOR_VALUE,numPtsPerEval> FunctorType;
226 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
229 case OPERATOR_CURL: {
230 auto workSize = Serial<OPERATOR_CURL>::getWorkSizePerPoint(order);
231 workViewType work(Kokkos::view_alloc(
"Basis_HCURL_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
232 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
233 OPERATOR_CURL,numPtsPerEval> FunctorType;
234 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
238 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
239 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Operator type not implemented" );
246 template<
typename SpT,
typename OT,
typename PT>
249 const EPointType pointType ) {
251 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
252 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
253 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): pointType must be either equispaced or warpblend.");
263 this->vinvLine_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hcurl::Quad::In::vinvLine", cardLine, cardLine);
264 this->vinvBubble_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hcurl::Quad::In::vinvBubble", cardBubble, cardBubble);
266 lineBasis.getVandermondeInverse(this->vinvLine_);
267 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
269 this->basisCardinality_ = 2*cardLine*cardBubble;
270 this->basisDegree_ = order;
271 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
272 this->basisType_ = BASIS_FEM_FIAT;
273 this->basisCoordinates_ = COORDINATES_CARTESIAN;
278 const ordinal_type tagSize = 4;
279 const ordinal_type posScDim = 0;
280 const ordinal_type posScOrd = 1;
281 const ordinal_type posDfOrd = 2;
286 ordinal_type tags[2*maxCardLine*maxCardBubble][4];
288 const ordinal_type edge_x[2] = {0,2};
289 const ordinal_type edge_y[2] = {3,1};
291 ordinal_type idx = 0;
301 intr_ndofs_per_direction = (cardLine-2)*cardBubble,
302 intr_ndofs = 2*intr_ndofs_per_direction;
305 for (ordinal_type j=0;j<cardLine;++j) {
306 const auto tag_y = lineBasis.
getDofTag(j);
307 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
308 const auto tag_x = bubbleBasis.
getDofTag(i);
310 if (tag_x(0) == 1 && tag_y(0) == 0) {
313 tags[idx][1] = edge_x[tag_y(1)];
314 tags[idx][2] = tag_x(2);
315 tags[idx][3] = tag_x(3);
320 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
321 tags[idx][3] = intr_ndofs;
327 for (ordinal_type j=0;j<cardBubble;++j) {
328 const auto tag_y = bubbleBasis.
getDofTag(j);
329 for (ordinal_type i=0;i<cardLine;++i,++idx) {
330 const auto tag_x = lineBasis.
getDofTag(i);
332 if (tag_x(0) == 0 && tag_y(0) == 1) {
335 tags[idx][1] = edge_y[tag_x(1)];
336 tags[idx][2] = tag_y(2);
337 tags[idx][3] = tag_y(3);
342 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2);
343 tags[idx][3] = intr_ndofs;
347 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
348 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): " \
349 "counted tag index is not same as cardinality." );
356 this->setOrdinalTagData(this->tagToOrdinal_,
359 this->basisCardinality_,
367 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
368 dofCoordsHost(
"dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
371 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
372 dofCoeffsHost(
"dofCoeffsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
374 Kokkos::DynRankView<typename scalarViewType::value_type,SpT>
375 dofCoordsLine(
"dofCoordsLine", cardLine, 1),
376 dofCoordsBubble(
"dofCoordsBubble", cardBubble, 1);
379 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
380 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
383 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
384 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
387 ordinal_type idx = 0;
390 for (ordinal_type j=0;j<cardLine;++j) {
391 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
392 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
393 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
394 dofCoeffsHost(idx,0) = 1.0;
399 for (ordinal_type j=0;j<cardBubble;++j) {
400 for (ordinal_type i=0;i<cardLine;++i,++idx) {
401 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
402 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
403 dofCoeffsHost(idx,1) = 1.0;
408 this->dofCoords_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoordsHost);
409 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
411 this->dofCoeffs_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoeffsHost);
412 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
Basis_HCURL_QUAD_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
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.
virtual void getDofCoords(scalarViewType dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.