Intrepid2
Intrepid2_HDIV_HEX_In_FEMDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
49#ifndef __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
50#define __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
51
52namespace Intrepid2 {
53
54 // -------------------------------------------------------------------------------------
55 namespace Impl {
56
57 template<EOperator opType>
58 template<typename OutputViewType,
59 typename inputViewType,
60 typename workViewType,
61 typename vinvViewType>
62 KOKKOS_INLINE_FUNCTION
63 void
64 Basis_HDIV_HEX_In_FEM::Serial<opType>::
65 getValues( OutputViewType output,
66 const inputViewType input,
67 workViewType work,
68 const vinvViewType vinvLine,
69 const vinvViewType vinvBubble) {
70 const ordinal_type cardLine = vinvLine.extent(0);
71 const ordinal_type cardBubble = vinvBubble.extent(0);
72
73 const ordinal_type npts = input.extent(0);
74
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));
79
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;
85
86 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
87 auto vcprop = Kokkos::common_view_alloc_prop(work);
88
89 switch (opType) {
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);
95
96 // tensor product
97 ordinal_type idx = 0;
98 {
99 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
100 getValues(outputLine, input_x, workLine, vinvLine);
101
102 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
103 getValues(outputBubble_A, input_y, workLine, vinvBubble);
104
105 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
106 getValues(outputBubble_B, input_z, workLine, vinvBubble);
107
108
109 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
110 const auto output_x = outputLine;
111 const auto output_y = outputBubble_A;
112 const auto output_z = outputBubble_B;
113
114 for (ordinal_type k=0;k<cardBubble;++k) // z
115 for (ordinal_type j=0;j<cardBubble;++j) // y
116 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
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;
121 }
122 }
123 {
124 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
125 getValues(outputBubble_A, input_x, workLine, vinvBubble);
126
127 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
128 getValues(outputLine, input_y, workLine, vinvLine);
129
130 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
131 // getValues(outputBubble_B, input_z, workLine, vinvBubble);
132
133 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
134 const auto output_x = outputBubble_A;
135 const auto output_y = outputLine;
136 const auto output_z = outputBubble_B;
137
138 for (ordinal_type k=0;k<cardBubble;++k) // z
139 for (ordinal_type j=0;j<cardLine;++j) // y
140 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
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;
145 }
146 }
147 {
148 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
149 // getValues(outputBubble_A, input_x, workLine, vinvBubble);
150
151 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
152 getValues(outputBubble_B, input_y, workLine, vinvBubble);
153
154 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
155 getValues(outputLine, input_z, workLine, vinvLine);
156
157 // z component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
158 const auto output_x = outputBubble_A;
159 const auto output_y = outputBubble_B;
160 const auto output_z = outputLine;
161
162 for (ordinal_type k=0;k<cardLine;++k) // z
163 for (ordinal_type j=0;j<cardBubble;++j) // y
164 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
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);
169 }
170 }
171 break;
172 }
173 case OPERATOR_DIV: {
174 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
175 // A line value
176 viewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
177 // B line value
178 viewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
179 // Line grad
180 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
181
182 // tensor product
183 ordinal_type idx = 0;
184
185 { // x - component
186 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
187 getValues(outputLine, input_x, workLine, vinvLine, 1);
188
189 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
190 getValues(outputBubble_A, input_y, workLine, vinvBubble);
191
192 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
193 getValues(outputBubble_B, input_z, workLine, vinvBubble);
194
195 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
196 const auto output_dx = outputLine;
197 const auto output_y = outputBubble_A;
198 const auto output_z = outputBubble_B;
199
200 for (ordinal_type k=0;k<cardBubble;++k) // z
201 for (ordinal_type j=0;j<cardBubble;++j) // y
202 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
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);
205 }
206 { // y - component
207 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
208 getValues(outputBubble_A, input_x, workLine, vinvBubble);
209
210 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
211 getValues(outputLine, input_y, workLine, vinvLine, 1);
212
213 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
214 // getValues(outputBubble_B, input_z, workLine, vinvBubble);
215
216 //(bubbleBasis(z) lineBasis(y) bubbleBasis(x))
217 const auto output_x = outputBubble_A;
218 const auto output_dy = outputLine;
219 const auto output_z = outputBubble_B;
220
221 for (ordinal_type k=0;k<cardBubble;++k) // z
222 for (ordinal_type j=0;j<cardLine;++j) // y
223 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
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);
226 }
227 { // z - component
228 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
229 // getValues(outputBubble_A, input_x, workLine, vinvBubble);
230
231 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
232 getValues(outputBubble_B, input_y, workLine, vinvBubble);
233
234 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
235 getValues(outputLine, input_z, workLine, vinvLine, 1);
236
237 // (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
238 const auto output_x = outputBubble_A;
239 const auto output_y = outputBubble_B;
240 const auto output_dz = outputLine;
241
242 for (ordinal_type k=0;k<cardLine;++k) // z
243 for (ordinal_type j=0;j<cardBubble;++j) // y
244 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
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);
247 }
248 break;
249 }
250 default: {
251 INTREPID2_TEST_FOR_ABORT( true,
252 ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::Serial::getValues) operator is not supported" );
253 }
254 }
255 }
256
257 template<typename DT, ordinal_type numPtsPerEval,
258 typename outputValueValueType, class ...outputValueProperties,
259 typename inputPointValueType, class ...inputPointProperties,
260 typename vinvValueType, class ...vinvProperties>
261 void
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,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
272
273 // loopSize corresponds to cardinality
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);
278
279 typedef typename inputPointViewType::value_type inputPointType;
280
281 const ordinal_type cardinality = outputValues.extent(0);
282 //get basis order based on basis cardinality.
283 ordinal_type order = 0;
284 ordinal_type cardBubble; // = std::cbrt(cardinality/3);
285 ordinal_type cardLine; // = cardBubble+1;
286 do {
287 cardBubble = Intrepid2::getPnCardinality<1>(order);
288 cardLine = Intrepid2::getPnCardinality<1>(++order);
289 } while((3*cardBubble*cardBubble*cardLine != cardinality) && (order != Parameters::MaxOrder));
290
291 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
292 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
293
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) );
301 break;
302 }
303 case OPERATOR_DIV: {
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) );
309 break;
310 }
311 default: {
312 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
313 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Operator type not implemented" );
314 // break; commented out since exception is thrown
315 }
316 }
317 }
318 }
319
320 // -------------------------------------------------------------------------------------
321
322 template<typename DT, typename OT, typename PT>
324 Basis_HDIV_HEX_In_FEM( const ordinal_type order,
325 const EPointType pointType ) {
326
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.");
330
331 // this should be created in host and vinv should be deep copied into device space
332 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
333 Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
334
335 const ordinal_type
336 cardLine = lineBasis.getCardinality(),
337 cardBubble = bubbleBasis.getCardinality();
338
339 this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvLine", cardLine, cardLine);
340 this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
341
342 lineBasis.getVandermondeInverse(this->vinvLine_);
343 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
344
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_LAGRANGIAN;
349 this->basisCoordinates_ = COORDINATES_CARTESIAN;
350 this->functionSpace_ = FUNCTION_SPACE_HDIV;
351 pointType_ = pointType;
352
353 // initialize tags
354 {
355 // Basis-dependent initializations
356 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
357 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
358 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
359 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
360
361 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
362 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
363 ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
364
365 const ordinal_type face_yz[2] = {3, 1};
366 const ordinal_type face_xz[2] = {0, 2};
367 const ordinal_type face_xy[2] = {4, 5};
368
369 {
370 ordinal_type idx = 0;
371
375
376 // since there are x/y components in the interior
377 // dof sum should be computed before the information
378 // is assigned to tags
379 const ordinal_type
380 intr_ndofs_per_direction = (cardLine-2)*cardBubble*cardBubble,
381 intr_ndofs = 3*intr_ndofs_per_direction;
382
383 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
384 for (ordinal_type k=0;k<cardBubble;++k) { // z
385 const auto tag_z = bubbleBasis.getDofTag(k);
386 for (ordinal_type j=0;j<cardBubble;++j) { // y
387 const auto tag_y = bubbleBasis.getDofTag(j);
388 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
389 const auto tag_x = lineBasis.getDofTag(i);
390
391 if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
392 // face, x vert, y edge, z edge
393 tags[idx][0] = 2; // face dof
394 tags[idx][1] = face_yz[tag_x(1)]; // face id
395 tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2); // local dof id
396 tags[idx][3] = tag_y(3)*tag_z(3); // total number of dofs in this vertex
397 } else {
398 // interior
399 tags[idx][0] = 3; // interior dof
400 tags[idx][1] = 0;
401 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
402 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
403 }
404 }
405 }
406 }
407
408 // y component (bubbleBasis(z) lineBasis(y) bubbleBasis(x))
409 for (ordinal_type k=0;k<cardBubble;++k) { // z
410 const auto tag_z = bubbleBasis.getDofTag(k);
411 for (ordinal_type j=0;j<cardLine;++j) { // y
412 const auto tag_y = lineBasis.getDofTag(j);
413 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
414 const auto tag_x = bubbleBasis.getDofTag(i);
415
416 if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
417 // face, x edge, y vert, z edge
418 tags[idx][0] = 2; // face dof
419 tags[idx][1] = face_xz[tag_y(1)]; // face id
420 tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2); // local dof id
421 tags[idx][3] = tag_x(3)*tag_z(3); // total number of dofs in this vertex
422 } else {
423 // interior
424 tags[idx][0] = 3; // interior dof
425 tags[idx][1] = 0;
426 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); // local dof id
427 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
428 }
429 }
430 }
431 }
432
433 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
434 for (ordinal_type k=0;k<cardLine;++k) { // y
435 const auto tag_z = lineBasis.getDofTag(k);
436 for (ordinal_type j=0;j<cardBubble;++j) { // z
437 const auto tag_y = bubbleBasis.getDofTag(j);
438 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
439 const auto tag_x = bubbleBasis.getDofTag(i);
440
441 if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
442 // face, x edge, y edge, z vert
443 tags[idx][0] = 2; // face dof
444 tags[idx][1] = face_xy[tag_z(1)]; // face id
445 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
446 tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
447 } else {
448 // interior
449 tags[idx][0] = 3; // interior dof
450 tags[idx][1] = 0;
451 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); // local dof id
452 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
453 }
454 }
455 }
456 }
457
458 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
459 ">>> ERROR (Basis_HDIV_HEX_In_FEM): " \
460 "counted tag index is not same as cardinality." );
461 }
462
463 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
464
465 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
466 // tags are constructed on host
467 this->setOrdinalTagData(this->tagToOrdinal_,
468 this->ordinalToTag_,
469 tagView,
470 this->basisCardinality_,
471 tagSize,
472 posScDim,
473 posScOrd,
474 posDfOrd);
475 }
476
477 // dofCoords on host and create its mirror view to device
478 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
479 dofCoordsHost("dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
480
481 // dofCoeffs on host and create its mirror view to device
482 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
483 dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
484
485 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
486 dofCoordsLine("dofCoordsLine", cardLine, 1),
487 dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
488
489 lineBasis.getDofCoords(dofCoordsLine);
490 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
491 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
492
493 bubbleBasis.getDofCoords(dofCoordsBubble);
494 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
495 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
496
497 {
498 ordinal_type idx = 0;
499
500 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
501 for (ordinal_type k=0;k<cardBubble;++k) { // z
502 for (ordinal_type j=0;j<cardBubble;++j) { // y
503 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
504 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
505 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
506 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
507 dofCoeffsHost(idx,0) = 1.0;
508 }
509 }
510 }
511
512 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
513 for (ordinal_type k=0;k<cardBubble;++k) { // z
514 for (ordinal_type j=0;j<cardLine;++j) { // y
515 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
516 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
517 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
518 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
519 dofCoeffsHost(idx,1) = 1.0;
520 }
521 }
522 }
523
524 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
525 for (ordinal_type k=0;k<cardLine;++k) { // z
526 for (ordinal_type j=0;j<cardBubble;++j) { // y
527 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
528 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
529 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
530 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
531 dofCoeffsHost(idx,2) = 1.0;
532 }
533 }
534 }
535 }
536
537 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
538 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
539
540 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
541 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
542 }
543}
544
545#endif
Basis_HDIV_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
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.
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
ordinal_type getCardinality() const
Returns cardinality of the basis.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.