48 #ifndef __INTREPID2_POINTTOOLS_DEF_HPP__
49 #define __INTREPID2_POINTTOOLS_DEF_HPP__
51 #if defined(_MSC_VER) || defined(_WIN32) && defined(__ICL)
53 #ifndef _USE_MATH_DEFINES
54 #define _USE_MATH_DEFINES
69 const ordinal_type order,
70 const ordinal_type offset ) {
71 #ifdef HAVE_INTREPID2_DEBUG
72 INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
73 std::invalid_argument ,
74 ">>> ERROR (PointTools::getLatticeSize): order and offset must be positive values." );
76 ordinal_type r_val = 0;
77 switch (cellType.getBaseKey()) {
78 case shards::Tetrahedron<>::key: {
79 const auto effectiveOrder = order - 4 * offset;
80 r_val = (effectiveOrder < 0 ? 0 :(effectiveOrder+1)*(effectiveOrder+2)*(effectiveOrder+3)/6);
83 case shards::Triangle<>::key: {
84 const auto effectiveOrder = order - 3 * offset;
85 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1)*(effectiveOrder+2)/2);
88 case shards::Line<>::key: {
89 const auto effectiveOrder = order - 2 * offset;
90 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1));
93 case shards::Quadrilateral<>::key: {
94 const auto effectiveOrder = order - 2 * offset;
95 r_val = std::pow(effectiveOrder < 0 ? 0 : (effectiveOrder+1),2);
98 case shards::Hexahedron<>::key: {
99 const auto effectiveOrder = order - 2 * offset;
100 r_val = std::pow(effectiveOrder < 0 ? 0 : (effectiveOrder+1),3);
103 case shards::Pyramid<>::key: {
104 const auto effectiveOrder = order - 2 * offset;
105 r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1)*(effectiveOrder+2)*(2*effectiveOrder+3)/6);
109 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
110 ">>> ERROR (Intrepid2::PointTools::getLatticeSize): the specified cell type is not supported." );
116 template<
typename pointValueType,
class ...pointProperties>
119 getLattice( Kokkos::DynRankView<pointValueType,pointProperties...> points,
120 const shards::CellTopology cell,
121 const ordinal_type order,
122 const ordinal_type offset,
123 const EPointType pointType ) {
124 #ifdef HAVE_INTREPID2_DEBUG
125 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
126 std::invalid_argument ,
127 ">>> ERROR (PointTools::getLattice): points rank must be 2." );
128 INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
129 std::invalid_argument ,
130 ">>> ERROR (PointTools::getLattice): order and offset must be positive values." );
132 const size_type latticeSize =
getLatticeSize( cell, order, offset );
133 const size_type spaceDim = cell.getDimension();
135 INTREPID2_TEST_FOR_EXCEPTION( points.extent(0) != latticeSize ||
136 points.extent(1) != spaceDim,
137 std::invalid_argument ,
138 ">>> ERROR (PointTools::getLattice): dimension does not match to lattice size." );
141 switch (cell.getBaseKey()) {
143 case shards::Pyramid<>::key:
getLatticePyramid ( points, order, offset, pointType );
break;
144 case shards::Triangle<>::key:
getLatticeTriangle ( points, order, offset, pointType );
break;
145 case shards::Line<>::key:
getLatticeLine ( points, order, offset, pointType );
break;
146 case shards::Quadrilateral<>::key: {
147 auto hostPoints = Kokkos::create_mirror_view(points);
148 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
149 const ordinal_type numPoints =
getLatticeSize( line, order, offset );
150 auto linePoints = getMatchingViewWithLabel(hostPoints,
"linePoints", numPoints, 1);
153 for (ordinal_type j=0; j<numPoints; ++j) {
154 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
155 hostPoints(idx,0) = linePoints(i,0);
156 hostPoints(idx,1) = linePoints(j,0);
159 Kokkos::deep_copy(points,hostPoints);
162 case shards::Hexahedron<>::key: {
163 auto hostPoints = Kokkos::create_mirror_view(points);
164 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
165 const ordinal_type numPoints =
getLatticeSize( line, order, offset );
166 auto linePoints = getMatchingViewWithLabel(hostPoints,
"linePoints", numPoints, 1);
169 for (ordinal_type k=0; k<numPoints; ++k) {
170 for (ordinal_type j=0; j<numPoints; ++j) {
171 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
172 hostPoints(idx,0) = linePoints(i,0);
173 hostPoints(idx,1) = linePoints(j,0);
174 hostPoints(idx,2) = linePoints(k,0);
178 Kokkos::deep_copy(points,hostPoints);
182 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
183 ">>> ERROR (Intrepid2::PointTools::getLattice): the specified cell type is not supported." );
188 template<
typename pointValueType,
class ...pointProperties>
192 const ordinal_type order ) {
193 #ifdef HAVE_INTREPID2_DEBUG
194 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
195 std::invalid_argument ,
196 ">>> ERROR (PointTools::getGaussPoints): points rank must be 1." );
197 INTREPID2_TEST_FOR_EXCEPTION( order < 0,
198 std::invalid_argument ,
199 ">>> ERROR (PointTools::getGaussPoints): order must be positive value." );
201 const ordinal_type np = order + 1;
202 const double alpha = 0.0, beta = 0.0;
205 Kokkos::View<pointValueType*,Kokkos::HostSpace>
206 zHost(
"PointTools::getGaussPoints::z", np),
207 wHost(
"PointTools::getGaussPoints::w", np);
214 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
215 auto pts = Kokkos::subview( points, range_type(0,np), 0 );
217 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data(), np);
218 Kokkos::deep_copy(pts, z);
225 template<
typename pointValueType,
class ...pointProperties>
229 const ordinal_type order,
230 const ordinal_type offset,
231 const EPointType pointType ) {
236 INTREPID2_TEST_FOR_EXCEPTION(
true ,
237 std::invalid_argument ,
238 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
243 template<
typename pointValueType,
class ...pointProperties>
247 const ordinal_type order,
248 const ordinal_type offset,
249 const EPointType pointType ) {
254 INTREPID2_TEST_FOR_EXCEPTION(
true ,
255 std::invalid_argument ,
256 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
261 template<
typename pointValueType,
class ...pointProperties>
265 const ordinal_type order,
266 const ordinal_type offset,
267 const EPointType pointType ) {
272 INTREPID2_TEST_FOR_EXCEPTION(
true ,
273 std::invalid_argument ,
274 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
279 template<
typename pointValueType,
class ...pointProperties>
283 const ordinal_type order,
284 const ordinal_type offset,
285 const EPointType pointType )
290 INTREPID2_TEST_FOR_EXCEPTION(
true ,
291 std::invalid_argument ,
292 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
299 template<
typename pointValueType,
class ...pointProperties>
303 const ordinal_type order,
304 const ordinal_type offset ) {
305 auto pointsHost = Kokkos::create_mirror_view(points);
308 pointsHost(0,0) = 0.0;
310 const pointValueType h = 2.0 / order;
311 const ordinal_type ibeg = offset, iend = order-offset+1;
312 for (ordinal_type i=ibeg;i<iend;++i)
313 pointsHost(i-ibeg, 0) = -1.0 + h * (pointValueType) i;
316 Kokkos::deep_copy(points, pointsHost);
319 template<
typename pointValueType,
class ...pointProperties>
323 const ordinal_type order,
324 const ordinal_type offset ) {
327 const ordinal_type np = order + 1;
328 const double alpha = 0.0, beta = 0.0;
329 const ordinal_type s = np - 2*offset;
333 Kokkos::View<pointValueType*,Kokkos::HostSpace>
334 zHost(
"PointTools::getGaussPoints::z", np),
335 wHost(
"PointTools::getGaussPoints::w", np);
342 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
343 auto pts = Kokkos::subview( points, range_type(0, s), 0 );
346 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data() + offset, np-offset);
348 const auto common_range = range_type(0, std::min(pts.extent(0), z.extent(0)));
349 Kokkos::deep_copy(Kokkos::subview(pts, common_range),
350 Kokkos::subview(z, common_range));
354 template<
typename pointValueType,
class ...pointProperties>
358 const ordinal_type order,
359 const ordinal_type offset ) {
360 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
361 std::invalid_argument ,
362 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTriangle): order must be positive" );
364 auto pointsHost = Kokkos::create_mirror_view(points);
366 const pointValueType h = 1.0 / order;
367 ordinal_type cur = 0;
369 for (ordinal_type i=offset;i<=order-offset;i++) {
370 for (ordinal_type j=offset;j<=order-i-offset;j++) {
371 pointsHost(cur,0) = j * h ;
372 pointsHost(cur,1) = i * h;
376 Kokkos::deep_copy(points, pointsHost);
379 template<
typename pointValueType,
class ...pointProperties>
383 const ordinal_type order ,
384 const ordinal_type offset ) {
385 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
386 std::invalid_argument ,
387 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
389 auto pointsHost = Kokkos::create_mirror_view(points);
391 const pointValueType h = 1.0 / order;
392 ordinal_type cur = 0;
394 for (ordinal_type i=offset;i<=order-offset;i++) {
395 for (ordinal_type j=offset;j<=order-i-offset;j++) {
396 for (ordinal_type k=offset;k<=order-i-j-offset;k++) {
397 pointsHost(cur,0) = k * h;
398 pointsHost(cur,1) = j * h;
399 pointsHost(cur,2) = i * h;
404 Kokkos::deep_copy(points, pointsHost);
407 template<
typename pointValueType,
class ...pointProperties>
411 const ordinal_type order ,
412 const ordinal_type offset ) {
413 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
414 std::invalid_argument ,
415 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticePyramid): order must be positive" );
417 auto pointsHost = Kokkos::create_mirror_view(points);
419 const pointValueType h = 1.0 / order;
420 ordinal_type cur = 0;
422 for (ordinal_type i=offset;i<=order-offset;i++) {
423 pointValueType z = i * h;
424 for (ordinal_type j=offset;j<=order-i-offset;j++) {
425 for (ordinal_type k=offset;k<=order-i-offset;k++) {
426 pointsHost(cur,0) = (2 * k * h - 1.) * (1-z);
427 pointsHost(cur,1) = (2 * j * h - 1.) * (1-z);
428 pointsHost(cur,2) = z;
433 Kokkos::deep_copy(points, pointsHost);
436 template<
typename pointValueType,
class ...pointProperties>
439 warpFactor( Kokkos::DynRankView<pointValueType,pointProperties...> warp,
440 const ordinal_type order,
441 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes ,
442 const Kokkos::DynRankView<pointValueType,pointProperties...> xout
445 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.extent(0) != xout.extent(0) ) ,
446 std::invalid_argument ,
447 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
449 Kokkos::deep_copy(warp, pointValueType(0.0));
451 ordinal_type xout_dim0 = xout.extent(0);
452 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> d(
"d", xout_dim0 );
454 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> xeq_(
"xeq", order + 1 ,1);
456 const auto xeq = Kokkos::subview(xeq_, Kokkos::ALL(),0);
458 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.extent(0) != xnodes.extent(0) ) ,
459 std::invalid_argument ,
460 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
463 for (ordinal_type i=0;i<=order;i++) {
465 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
467 for (ordinal_type j=1;j<order;j++) {
469 for (ordinal_type k=0;k<xout_dim0;k++) {
470 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
477 for (ordinal_type k=0;k<xout_dim0;k++) {
478 d(k) = -d(k) / (xeq(i) - xeq(0));
483 for (ordinal_type k=0;k<xout_dim0;k++) {
484 d(k) = d(k) / (xeq(i) - xeq(order));
488 for (ordinal_type k=0;k<xout_dim0;k++) {
495 template<
typename pointValueType,
class ...pointProperties>
499 const ordinal_type order ,
500 const ordinal_type offset)
504 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX", order + 1 , 1 );
511 pointValueType alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
512 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
514 pointValueType alpha;
516 if (order >= 1 && order < 16) {
517 alpha = alpopt[order-1];
523 const ordinal_type p = order;
524 ordinal_type N = (p+1)*(p+2)/2;
527 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1(
"L1", N );
528 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2(
"L2", N );
529 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3(
"L3", N );
530 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> X(
"X", N);
531 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> Y(
"Y", N);
534 for (ordinal_type n=1;n<=p+1;n++) {
535 for (ordinal_type m=1;m<=p+2-n;m++) {
536 L1(sk) = (n-1) / (pointValueType)p;
537 L3(sk) = (m-1) / (pointValueType)p;
538 L2(sk) = 1.0 - L1(sk) - L3(sk);
543 for (ordinal_type n=0;n<N;n++) {
544 X(n) = -L2(n) + L3(n);
545 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
549 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend1(
"blend1", N);
550 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend2(
"blend2", N);
551 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend3(
"blend3", N);
553 for (ordinal_type n=0;n<N;n++) {
554 blend1(n) = 4.0 * L2(n) * L3(n);
555 blend2(n) = 4.0 * L1(n) * L3(n);
556 blend3(n) = 4.0 * L1(n) * L2(n);
560 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2", N);
561 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3", N);
562 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1", N);
564 for (ordinal_type k=0;k<N;k++) {
565 L3mL2(k) = L3(k)-L2(k);
566 L1mL3(k) = L1(k)-L3(k);
567 L2mL1(k) = L2(k)-L1(k);
570 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1", N);
571 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2", N);
572 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3", N);
574 warpFactor( warpfactor1, order , gaussX , L3mL2 );
575 warpFactor( warpfactor2, order , gaussX , L1mL3 );
576 warpFactor( warpfactor3, order , gaussX , L2mL1 );
578 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp1(
"warp1", N);
579 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp2(
"warp2", N);
580 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp3(
"warp3", N);
582 for (ordinal_type k=0;k<N;k++) {
583 warp1(k) = blend1(k) * warpfactor1(k) *
584 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
585 warp2(k) = blend2(k) * warpfactor2(k) *
586 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
587 warp3(k) = blend3(k) * warpfactor3(k) *
588 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
591 for (ordinal_type k=0;k<N;k++) {
592 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
593 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
596 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warXY(
"warXY", 1, N,2);
598 for (ordinal_type k=0;k<N;k++) {
599 warXY(0, k,0) = X(k);
600 warXY(0, k,1) = Y(k);
605 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warburtonVerts(
"warburtonVerts", 1,3,2);
606 warburtonVerts(0,0,0) = -1.0;
607 warburtonVerts(0,0,1) = -1.0/std::sqrt(3.0);
608 warburtonVerts(0,1,0) = 1.0;
609 warburtonVerts(0,1,1) = -1.0/std::sqrt(3.0);
610 warburtonVerts(0,2,0) = 0.0;
611 warburtonVerts(0,2,1) = 2.0/std::sqrt(3.0);
613 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> refPts(
"refPts", 1, N,2);
618 shards::getCellTopologyData< shards::Triangle<3> >()
622 auto pointsHost = Kokkos::create_mirror_view(points);
624 ordinal_type noffcur = 0;
625 ordinal_type offcur = 0;
626 for (ordinal_type i=0;i<=order;i++) {
627 for (ordinal_type j=0;j<=order-i;j++) {
628 if ( (i >= offset) && (i <= order-offset) &&
629 (j >= offset) && (j <= order-i-offset) ) {
630 pointsHost(offcur,0) = refPts(0, noffcur,0);
631 pointsHost(offcur,1) = refPts(0, noffcur,1);
638 Kokkos::deep_copy(points, pointsHost);
643 template<
typename pointValueType,
class ...pointProperties>
647 const ordinal_type order ,
648 const pointValueType pval ,
649 const Kokkos::DynRankView<pointValueType,pointProperties...> ,
650 const Kokkos::DynRankView<pointValueType,pointProperties...> L2,
651 const Kokkos::DynRankView<pointValueType,pointProperties...> L3,
652 const Kokkos::DynRankView<pointValueType,pointProperties...> L4
659 template<
typename pointValueType,
class ...pointProperties>
662 evalshift( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
663 const ordinal_type order ,
664 const pointValueType pval ,
665 const Kokkos::DynRankView<pointValueType,pointProperties...> L1 ,
666 const Kokkos::DynRankView<pointValueType,pointProperties...> L2 ,
667 const Kokkos::DynRankView<pointValueType,pointProperties...> L3
671 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX",order+1,1);
674 const ordinal_type N = L1.extent(0);
677 for (ordinal_type k=0;k<=order;k++) {
682 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend1(
"blend1",N);
683 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend2(
"blend2",N);
684 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend3(
"blend3",N);
686 for (ordinal_type i=0;i<N;i++) {
687 blend1(i) = L2(i) * L3(i);
688 blend2(i) = L1(i) * L3(i);
689 blend3(i) = L1(i) * L2(i);
693 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1",N);
694 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2",N);
695 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3",N);
698 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2",N);
699 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3",N);
700 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1",N);
702 for (ordinal_type k=0;k<N;k++) {
703 L3mL2(k) = L3(k)-L2(k);
704 L1mL3(k) = L1(k)-L3(k);
705 L2mL1(k) = L2(k)-L1(k);
708 evalwarp( warpfactor1 , order , gaussX , L3mL2 );
709 evalwarp( warpfactor2 , order , gaussX , L1mL3 );
710 evalwarp( warpfactor3 , order , gaussX , L2mL1 );
712 for (ordinal_type k=0;k<N;k++) {
713 warpfactor1(k) *= 4.0;
714 warpfactor2(k) *= 4.0;
715 warpfactor3(k) *= 4.0;
718 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp1(
"warp1",N);
719 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp2(
"warp2",N);
720 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp3(
"warp3",N);
722 for (ordinal_type k=0;k<N;k++) {
723 warp1(k) = blend1(k) * warpfactor1(k) *
724 ( 1.0 + pval * pval * L1(k) * L1(k) );
725 warp2(k) = blend2(k) * warpfactor2(k) *
726 ( 1.0 + pval * pval * L2(k) * L2(k) );
727 warp3(k) = blend3(k) * warpfactor3(k) *
728 ( 1.0 + pval * pval * L3(k) * L3(k) );
731 for (ordinal_type k=0;k<N;k++) {
732 dxy(k,0) = 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos( 4.0*M_PI/3.0 ) * warp3(k);
733 dxy(k,1) = 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4.0*M_PI/3.0 ) * warp3(k);
738 template<
typename pointValueType,
class ...pointProperties>
741 evalwarp(Kokkos::DynRankView<pointValueType,pointProperties...> warp ,
742 const ordinal_type order ,
743 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes,
744 const Kokkos::DynRankView<pointValueType,pointProperties...> xout )
746 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> xeq(
"xeq",order+1);
748 ordinal_type xout_dim0 = xout.extent(0);
749 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> d(
"d",xout_dim0);
753 for (ordinal_type i=0;i<=order;i++) {
754 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
757 for (ordinal_type i=0;i<=order;i++) {
758 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
759 for (ordinal_type j=1;j<order;j++) {
761 for (ordinal_type k=0;k<xout_dim0;k++) {
762 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
767 for (ordinal_type k=0;k<xout_dim0;k++) {
768 d(k) = -d(k)/(xeq(i)-xeq(0));
772 for (ordinal_type k=0;k<xout_dim0;k++) {
773 d(k) = d(k)/(xeq(i)-xeq(order));
777 for (ordinal_type k=0;k<xout_dim0;k++) {
784 template<
typename pointValueType,
class ...pointProperties>
788 const ordinal_type order ,
789 const ordinal_type offset )
791 pointValueType alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
792 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
793 pointValueType alpha;
796 alpha = alphastore[std::max(order-1,0)];
802 const ordinal_type N = (order+1)*(order+2)*(order+3)/6;
803 pointValueType tol = 1.e-10;
805 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> shift(
"shift",N,3);
806 Kokkos::deep_copy(shift,0.0);
809 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> equipoints(
"equipoints", N,3);
811 for (ordinal_type n=0;n<=order;n++) {
812 for (ordinal_type m=0;m<=order-n;m++) {
813 for (ordinal_type q=0;q<=order-n-m;q++) {
814 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
815 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
816 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
824 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1(
"L1",N);
825 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2(
"L2",N);
826 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3(
"L3",N);
827 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L4(
"L4",N);
828 for (ordinal_type i=0;i<N;i++) {
829 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
830 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
831 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
832 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
837 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warVerts_(
"warVerts",1,4,3);
838 auto warVerts = Kokkos::subview(warVerts_,0,Kokkos::ALL(),Kokkos::ALL());
839 warVerts(0,0) = -1.0;
840 warVerts(0,1) = -1.0/sqrt(3.0);
841 warVerts(0,2) = -1.0/sqrt(6.0);
843 warVerts(1,1) = -1.0/sqrt(3.0);
844 warVerts(1,2) = -1.0/sqrt(6.0);
846 warVerts(2,1) = 2.0 / sqrt(3.0);
847 warVerts(2,2) = -1.0/sqrt(6.0);
850 warVerts(3,2) = 3.0 / sqrt(6.0);
854 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t1(
"t1",4,3);
855 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t2(
"t2",4,3);
856 for (ordinal_type i=0;i<3;i++) {
857 t1(0,i) = warVerts(1,i) - warVerts(0,i);
858 t1(1,i) = warVerts(1,i) - warVerts(0,i);
859 t1(2,i) = warVerts(2,i) - warVerts(1,i);
860 t1(3,i) = warVerts(2,i) - warVerts(0,i);
861 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
862 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
863 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
864 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
868 for (ordinal_type n=0;n<4;n++) {
870 pointValueType normt1n = 0.0;
871 pointValueType normt2n = 0.0;
872 for (ordinal_type i=0;i<3;i++) {
873 normt1n += (t1(n,i) * t1(n,i));
874 normt2n += (t2(n,i) * t2(n,i));
876 normt1n = sqrt(normt1n);
877 normt2n = sqrt(normt2n);
879 for (ordinal_type i=0;i<3;i++) {
886 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> XYZ(
"XYZ",N,3);
887 for (ordinal_type i=0;i<N;i++) {
888 for (ordinal_type j=0;j<3;j++) {
889 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
893 for (ordinal_type face=1;face<=4;face++) {
894 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> La, Lb, Lc, Ld;
895 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp(
"warp",N,2);
896 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend(
"blend",N);
897 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> denom(
"denom",N);
900 La = L1; Lb = L2; Lc = L3; Ld = L4;
break;
902 La = L2; Lb = L1; Lc = L3; Ld = L4;
break;
904 La = L3; Lb = L1; Lc = L4; Ld = L2;
break;
906 La = L4; Lb = L1; Lc = L3; Ld = L2;
break;
912 for (ordinal_type k=0;k<N;k++) {
913 blend(k) = Lb(k) * Lc(k) * Ld(k);
916 for (ordinal_type k=0;k<N;k++) {
917 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
920 for (ordinal_type k=0;k<N;k++) {
921 if (denom(k) > tol) {
922 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
928 for (ordinal_type k=0;k<N;k++) {
929 for (ordinal_type j=0;j<3;j++) {
930 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
931 + blend(k) * warp(k,1) * t2(face-1,j);
935 for (ordinal_type k=0;k<N;k++) {
936 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
937 for (ordinal_type j=0;j<3;j++) {
938 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
945 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> updatedPoints(
"updatedPoints",1,N,3);
946 for (ordinal_type k=0;k<N;k++) {
947 for (ordinal_type j=0;j<3;j++) {
948 updatedPoints(0,k,j) = XYZ(k,j) + shift(k,j);
955 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> refPts(
"refPts",1,N,3);
958 shards::getCellTopologyData<shards::Tetrahedron<4> >()
961 auto pointsHost = Kokkos::create_mirror_view(points);
963 ordinal_type noffcur = 0;
964 ordinal_type offcur = 0;
965 for (ordinal_type i=0;i<=order;i++) {
966 for (ordinal_type j=0;j<=order-i;j++) {
967 for (ordinal_type k=0;k<=order-i-j;k++) {
968 if ( (i >= offset) && (i <= order-offset) &&
969 (j >= offset) && (j <= order-i-offset) &&
970 (k >= offset) && (k <= order-i-j-offset) ) {
971 pointsHost(offcur,0) = refPts(0,noffcur,0);
972 pointsHost(offcur,1) = refPts(0,noffcur,1);
973 pointsHost(offcur,2) = refPts(0,noffcur,2);
981 Kokkos::deep_copy(points, pointsHost);
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.