49 #ifndef __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
50 #define __INTREPID2_HDIV_HEX_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_HEX_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));
78 const auto input_z = Kokkos::subview(input, Kokkos::ALL(), range_type(2,3));
80 const ordinal_type dim_s = get_dimension_scalar(work);
81 auto ptr0 = work.data();
82 auto ptr1 = work.data()+cardLine*npts*dim_s;
83 auto ptr2 = work.data()+2*cardLine*npts*dim_s;
84 auto ptr3 = work.data()+(2*cardLine+cardBubble)*npts*dim_s;
86 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
87 auto vcprop = Kokkos::common_view_alloc_prop(work);
90 case OPERATOR_VALUE: {
91 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
92 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
93 viewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
94 viewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
99 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
100 getValues(outputLine, input_x, workLine, vinvLine);
102 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
103 getValues(outputBubble_A, input_y, workLine, vinvBubble);
105 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
106 getValues(outputBubble_B, input_z, workLine, vinvBubble);
110 const auto output_x = outputLine;
111 const auto output_y = outputBubble_A;
112 const auto output_z = outputBubble_B;
114 for (ordinal_type k=0;k<cardBubble;++k)
115 for (ordinal_type j=0;j<cardBubble;++j)
116 for (ordinal_type i=0;i<cardLine;++i,++idx)
117 for (ordinal_type l=0;l<npts;++l) {
118 output.access(idx,l,0) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
119 output.access(idx,l,1) = 0.0;
120 output.access(idx,l,2) = 0.0;
124 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
125 getValues(outputBubble_A, input_x, workLine, vinvBubble);
127 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
128 getValues(outputLine, input_y, workLine, vinvLine);
134 const auto output_x = outputBubble_A;
135 const auto output_y = outputLine;
136 const auto output_z = outputBubble_B;
138 for (ordinal_type k=0;k<cardBubble;++k)
139 for (ordinal_type j=0;j<cardLine;++j)
140 for (ordinal_type i=0;i<cardBubble;++i,++idx)
141 for (ordinal_type l=0;l<npts;++l) {
142 output.access(idx,l,0) = 0.0;
143 output.access(idx,l,1) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
144 output.access(idx,l,2) = 0.0;
151 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
152 getValues(outputBubble_B, input_y, workLine, vinvBubble);
154 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
155 getValues(outputLine, input_z, workLine, vinvLine);
158 const auto output_x = outputBubble_A;
159 const auto output_y = outputBubble_B;
160 const auto output_z = outputLine;
162 for (ordinal_type k=0;k<cardLine;++k)
163 for (ordinal_type j=0;j<cardBubble;++j)
164 for (ordinal_type i=0;i<cardBubble;++i,++idx)
165 for (ordinal_type l=0;l<npts;++l) {
166 output.access(idx,l,0) = 0.0;
167 output.access(idx,l,1) = 0.0;
168 output.access(idx,l,2) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
174 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
176 viewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
178 viewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
180 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
183 ordinal_type idx = 0;
186 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
187 getValues(outputLine, input_x, workLine, vinvLine, 1);
189 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
190 getValues(outputBubble_A, input_y, workLine, vinvBubble);
192 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
193 getValues(outputBubble_B, input_z, workLine, vinvBubble);
196 const auto output_dx = outputLine;
197 const auto output_y = outputBubble_A;
198 const auto output_z = outputBubble_B;
200 for (ordinal_type k=0;k<cardBubble;++k)
201 for (ordinal_type j=0;j<cardBubble;++j)
202 for (ordinal_type i=0;i<cardLine;++i,++idx)
203 for (ordinal_type l=0;l<npts;++l)
204 output.access(idx,l) = output_dx.access(i,l,0)*output_y.access (j,l) *output_z.access(k,l);
207 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
208 getValues(outputBubble_A, input_x, workLine, vinvBubble);
210 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
211 getValues(outputLine, input_y, workLine, vinvLine, 1);
217 const auto output_x = outputBubble_A;
218 const auto output_dy = outputLine;
219 const auto output_z = outputBubble_B;
221 for (ordinal_type k=0;k<cardBubble;++k)
222 for (ordinal_type j=0;j<cardLine;++j)
223 for (ordinal_type i=0;i<cardBubble;++i,++idx)
224 for (ordinal_type l=0;l<npts;++l)
225 output.access(idx,l) = output_x.access(i,l)*output_dy.access(j,l,0)*output_z.access(k,l);
231 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
232 getValues(outputBubble_B, input_y, workLine, vinvBubble);
234 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
235 getValues(outputLine, input_z, workLine, vinvLine, 1);
238 const auto output_x = outputBubble_A;
239 const auto output_y = outputBubble_B;
240 const auto output_dz = outputLine;
242 for (ordinal_type k=0;k<cardLine;++k)
243 for (ordinal_type j=0;j<cardBubble;++j)
244 for (ordinal_type i=0;i<cardBubble;++i,++idx)
245 for (ordinal_type l=0;l<npts;++l)
246 output.access(idx,l) = output_x.access(i,l)*output_y.access(j,l)*output_dz.access(k,l,0);
251 INTREPID2_TEST_FOR_ABORT(
true,
252 ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::Serial::getValues) operator is not supported" );
257 template<
typename SpT, ordinal_type numPtsPerEval,
258 typename outputValueValueType,
class ...outputValueProperties,
259 typename inputPointValueType,
class ...inputPointProperties,
260 typename vinvValueType,
class ...vinvProperties>
262 Basis_HDIV_HEX_In_FEM::
263 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
264 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
265 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
266 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
267 const EOperator operatorType ) {
268 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
269 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
270 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
271 typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
274 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
275 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
276 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
277 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
279 typedef typename inputPointViewType::value_type inputPointType;
281 const ordinal_type cardinality = outputValues.extent(0);
283 ordinal_type order = 0;
284 ordinal_type cardBubble;
285 ordinal_type cardLine;
287 cardBubble = Intrepid2::getPnCardinality<1>(order);
288 cardLine = Intrepid2::getPnCardinality<1>(++order);
291 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
292 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
294 switch (operatorType) {
295 case OPERATOR_VALUE: {
296 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
297 workViewType work(Kokkos::view_alloc(
"Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
298 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
299 OPERATOR_VALUE,numPtsPerEval> FunctorType;
300 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
304 auto workSize = Serial<OPERATOR_DIV>::getWorkSizePerPoint(order);
305 workViewType work(Kokkos::view_alloc(
"Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
306 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
307 OPERATOR_DIV,numPtsPerEval> FunctorType;
308 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
312 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
313 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Operator type not implemented" );
322 template<
typename SpT,
typename OT,
typename PT>
325 const EPointType pointType ) {
327 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
328 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
329 ">>> ERROR (Basis_HDIV_HEX_In_FEM): pointType must be either equispaced or warpblend.");
339 this->vinvLine_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hcurl::Hex::In::vinvLine", cardLine, cardLine);
340 this->vinvBubble_ = Kokkos::DynRankView<typename scalarViewType::value_type,SpT>(
"Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
342 lineBasis.getVandermondeInverse(this->vinvLine_);
343 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
345 this->basisCardinality_ = 3*cardLine*cardBubble*cardBubble;
346 this->basisDegree_ = order;
347 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
348 this->basisType_ = BASIS_FEM_FIAT;
349 this->basisCoordinates_ = COORDINATES_CARTESIAN;
354 const ordinal_type tagSize = 4;
355 const ordinal_type posScDim = 0;
356 const ordinal_type posScOrd = 1;
357 const ordinal_type posDfOrd = 2;
361 ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
363 const ordinal_type face_yz[2] = {3, 1};
364 const ordinal_type face_xz[2] = {0, 2};
365 const ordinal_type face_xy[2] = {4, 5};
368 ordinal_type idx = 0;
378 intr_ndofs_per_direction = (cardLine-2)*cardBubble*cardBubble,
379 intr_ndofs = 3*intr_ndofs_per_direction;
382 for (ordinal_type k=0;k<cardBubble;++k) {
383 const auto tag_z = bubbleBasis.
getDofTag(k);
384 for (ordinal_type j=0;j<cardBubble;++j) {
385 const auto tag_y = bubbleBasis.
getDofTag(j);
386 for (ordinal_type i=0;i<cardLine;++i,++idx) {
387 const auto tag_x = lineBasis.
getDofTag(i);
389 if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
392 tags[idx][1] = face_yz[tag_x(1)];
393 tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2);
394 tags[idx][3] = tag_y(3)*tag_z(3);
399 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2);
400 tags[idx][3] = intr_ndofs;
407 for (ordinal_type k=0;k<cardBubble;++k) {
408 const auto tag_z = bubbleBasis.
getDofTag(k);
409 for (ordinal_type j=0;j<cardLine;++j) {
410 const auto tag_y = lineBasis.
getDofTag(j);
411 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
412 const auto tag_x = bubbleBasis.
getDofTag(i);
414 if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
417 tags[idx][1] = face_xz[tag_y(1)];
418 tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2);
419 tags[idx][3] = tag_x(3)*tag_z(3);
424 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2);
425 tags[idx][3] = intr_ndofs;
432 for (ordinal_type k=0;k<cardLine;++k) {
433 const auto tag_z = lineBasis.
getDofTag(k);
434 for (ordinal_type j=0;j<cardBubble;++j) {
435 const auto tag_y = bubbleBasis.
getDofTag(j);
436 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
437 const auto tag_x = bubbleBasis.
getDofTag(i);
439 if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
442 tags[idx][1] = face_xy[tag_z(1)];
443 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2);
444 tags[idx][3] = tag_x(3)*tag_y(3);
449 tags[idx][2] = 2*intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2);
450 tags[idx][3] = intr_ndofs;
456 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
457 ">>> ERROR (Basis_HDIV_HEX_In_FEM): " \
458 "counted tag index is not same as cardinality." );
465 this->setOrdinalTagData(this->tagToOrdinal_,
468 this->basisCardinality_,
476 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
477 dofCoordsHost(
"dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
480 Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
481 dofCoeffsHost(
"dofCoeffsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
483 Kokkos::DynRankView<typename scalarViewType::value_type,SpT>
484 dofCoordsLine(
"dofCoordsLine", cardLine, 1),
485 dofCoordsBubble(
"dofCoordsBubble", cardBubble, 1);
488 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
489 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
492 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
493 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
496 ordinal_type idx = 0;
499 for (ordinal_type k=0;k<cardBubble;++k) {
500 for (ordinal_type j=0;j<cardBubble;++j) {
501 for (ordinal_type i=0;i<cardLine;++i,++idx) {
502 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
503 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
504 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
505 dofCoeffsHost(idx,0) = 1.0;
511 for (ordinal_type k=0;k<cardBubble;++k) {
512 for (ordinal_type j=0;j<cardLine;++j) {
513 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
514 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
515 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
516 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
517 dofCoeffsHost(idx,1) = 1.0;
523 for (ordinal_type k=0;k<cardLine;++k) {
524 for (ordinal_type j=0;j<cardBubble;++j) {
525 for (ordinal_type i=0;i<cardBubble;++i,++idx) {
526 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
527 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
528 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
529 dofCoeffsHost(idx,2) = 1.0;
535 this->dofCoords_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoordsHost);
536 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
538 this->dofCoeffs_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoeffsHost);
539 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
ordinal_type getCardinality() const
Returns cardinality of the basis.
Basis_HDIV_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
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.