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);
104 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
105 ">>> ERROR (Intrepid2::PointTools::getLatticeSize): the specified cell type is not supported." );
111template<
typename pointValueType,
class ...pointProperties>
114getLattice( Kokkos::DynRankView<pointValueType,pointProperties...> points,
115 const shards::CellTopology cell,
116 const ordinal_type order,
117 const ordinal_type offset,
118 const EPointType pointType ) {
119#ifdef HAVE_INTREPID2_DEBUG
120 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
121 std::invalid_argument ,
122 ">>> ERROR (PointTools::getLattice): points rank must be 2." );
123 INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
124 std::invalid_argument ,
125 ">>> ERROR (PointTools::getLattice): order and offset must be positive values." );
127 const size_type latticeSize =
getLatticeSize( cell, order, offset );
128 const size_type spaceDim = cell.getDimension();
130 INTREPID2_TEST_FOR_EXCEPTION( points.extent(0) != latticeSize ||
131 points.extent(1) != spaceDim,
132 std::invalid_argument ,
133 ">>> ERROR (PointTools::getLattice): dimension does not match to lattice size." );
136 switch (cell.getBaseKey()) {
138 case shards::Triangle<>::key:
getLatticeTriangle ( points, order, offset, pointType );
break;
139 case shards::Line<>::key:
getLatticeLine ( points, order, offset, pointType );
break;
140 case shards::Quadrilateral<>::key: {
141 auto hostPoints = Kokkos::create_mirror_view(points);
142 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
143 const ordinal_type numPoints =
getLatticeSize( line, order, offset );
147 for (ordinal_type j=0; j<numPoints; ++j) {
148 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
149 hostPoints(idx,0) = linePoints(i,0);
150 hostPoints(idx,1) = linePoints(j,0);
153 Kokkos::deep_copy(points,hostPoints);
156 case shards::Hexahedron<>::key: {
157 auto hostPoints = Kokkos::create_mirror_view(points);
158 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2> >());
159 const ordinal_type numPoints =
getLatticeSize( line, order, offset );
163 for (ordinal_type k=0; k<numPoints; ++k) {
164 for (ordinal_type j=0; j<numPoints; ++j) {
165 for (ordinal_type i=0; i<numPoints; ++i, ++idx) {
166 hostPoints(idx,0) = linePoints(i,0);
167 hostPoints(idx,1) = linePoints(j,0);
168 hostPoints(idx,2) = linePoints(k,0);
172 Kokkos::deep_copy(points,hostPoints);
176 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
177 ">>> ERROR (Intrepid2::PointTools::getLattice): the specified cell type is not supported." );
182template<
typename pointValueType,
class ...pointProperties>
185getGaussPoints( Kokkos::DynRankView<pointValueType,pointProperties...> points,
186 const ordinal_type order ) {
187#ifdef HAVE_INTREPID2_DEBUG
188 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
189 std::invalid_argument ,
190 ">>> ERROR (PointTools::getGaussPoints): points rank must be 1." );
191 INTREPID2_TEST_FOR_EXCEPTION( order < 0,
192 std::invalid_argument ,
193 ">>> ERROR (PointTools::getGaussPoints): order must be positive value." );
195 const ordinal_type np = order + 1;
196 const double alpha = 0.0, beta = 0.0;
199 Kokkos::View<pointValueType*,Kokkos::HostSpace>
200 zHost(
"PointTools::getGaussPoints::z", np),
201 wHost(
"PointTools::getGaussPoints::w", np);
208 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
209 auto pts = Kokkos::subview( points, range_type(0,np), 0 );
211 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data(), np);
212 Kokkos::deep_copy(pts, z);
219template<
typename pointValueType,
class ...pointProperties>
222getLatticeLine( Kokkos::DynRankView<pointValueType,pointProperties...> points,
223 const ordinal_type order,
224 const ordinal_type offset,
225 const EPointType pointType ) {
230 INTREPID2_TEST_FOR_EXCEPTION(
true ,
231 std::invalid_argument ,
232 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
237template<
typename pointValueType,
class ...pointProperties>
241 const ordinal_type order,
242 const ordinal_type offset,
243 const EPointType pointType ) {
248 INTREPID2_TEST_FOR_EXCEPTION(
true ,
249 std::invalid_argument ,
250 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
255template<
typename pointValueType,
class ...pointProperties>
259 const ordinal_type order,
260 const ordinal_type offset,
261 const EPointType pointType ) {
266 INTREPID2_TEST_FOR_EXCEPTION(
true ,
267 std::invalid_argument ,
268 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
275template<
typename pointValueType,
class ...pointProperties>
279 const ordinal_type order,
280 const ordinal_type offset ) {
281 auto pointsHost = Kokkos::create_mirror_view(points);
284 pointsHost(0,0) = 0.0;
286 const pointValueType h = 2.0 / order;
287 const ordinal_type ibeg = offset, iend = order-offset+1;
288 for (ordinal_type i=ibeg;i<iend;++i)
289 pointsHost(i-ibeg, 0) = -1.0 + h * (pointValueType) i;
292 Kokkos::deep_copy(points, pointsHost);
295template<
typename pointValueType,
class ...pointProperties>
299 const ordinal_type order,
300 const ordinal_type offset ) {
303 const ordinal_type np = order + 1;
304 const double alpha = 0.0, beta = 0.0;
305 const ordinal_type s = np - 2*offset;
309 Kokkos::View<pointValueType*,Kokkos::HostSpace>
310 zHost(
"PointTools::getGaussPoints::z", np),
311 wHost(
"PointTools::getGaussPoints::w", np);
318 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
319 auto pts = Kokkos::subview( points, range_type(0, s), 0 );
322 auto z = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data() + offset, np-offset);
324 const auto common_range = range_type(0, std::min(pts.extent(0), z.extent(0)));
325 Kokkos::deep_copy(Kokkos::subview(pts, common_range),
326 Kokkos::subview(z, common_range));
330template<
typename pointValueType,
class ...pointProperties>
334 const ordinal_type order,
335 const ordinal_type offset ) {
336 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
337 std::invalid_argument ,
338 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTriangle): order must be positive" );
340 auto pointsHost = Kokkos::create_mirror_view(points);
342 const pointValueType h = 1.0 / order;
343 ordinal_type cur = 0;
345 for (ordinal_type i=offset;i<=order-offset;i++) {
346 for (ordinal_type j=offset;j<=order-i-offset;j++) {
347 pointsHost(cur,0) = j * h ;
348 pointsHost(cur,1) = i * h;
352 Kokkos::deep_copy(points, pointsHost);
355template<
typename pointValueType,
class ...pointProperties>
359 const ordinal_type order ,
360 const ordinal_type offset ) {
361 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
362 std::invalid_argument ,
363 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
365 auto pointsHost = Kokkos::create_mirror_view(points);
367 const pointValueType h = 1.0 / order;
368 ordinal_type cur = 0;
370 for (ordinal_type i=offset;i<=order-offset;i++) {
371 for (ordinal_type j=offset;j<=order-i-offset;j++) {
372 for (ordinal_type k=offset;k<=order-i-j-offset;k++) {
373 pointsHost(cur,0) = k * h;
374 pointsHost(cur,1) = j * h;
375 pointsHost(cur,2) = i * h;
380 Kokkos::deep_copy(points, pointsHost);
384template<
typename pointValueType,
class ...pointProperties>
387warpFactor( Kokkos::DynRankView<pointValueType,pointProperties...> warp,
388 const ordinal_type order,
389 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes ,
390 const Kokkos::DynRankView<pointValueType,pointProperties...> xout
393 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.extent(0) != xout.extent(0) ) ,
394 std::invalid_argument ,
395 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
397 Kokkos::deep_copy(warp, pointValueType(0.0));
399 ordinal_type xout_dim0 = xout.extent(0);
400 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> d(
"d", xout_dim0 );
402 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> xeq_(
"xeq", order + 1 ,1);
404 const auto xeq = Kokkos::subview(xeq_, Kokkos::ALL(),0);
406 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.extent(0) != xnodes.extent(0) ) ,
407 std::invalid_argument ,
408 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
411 for (ordinal_type i=0;i<=order;i++) {
413 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
415 for (ordinal_type j=1;j<order;j++) {
417 for (ordinal_type k=0;k<xout_dim0;k++) {
418 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
425 for (ordinal_type k=0;k<xout_dim0;k++) {
426 d(k) = -d(k) / (xeq(i) - xeq(0));
431 for (ordinal_type k=0;k<xout_dim0;k++) {
432 d(k) = d(k) / (xeq(i) - xeq(order));
436 for (ordinal_type k=0;k<xout_dim0;k++) {
443template<
typename pointValueType,
class ...pointProperties>
447 const ordinal_type order ,
448 const ordinal_type offset)
452 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX", order + 1 , 1 );
459 pointValueType alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
460 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
462 pointValueType alpha;
464 if (order >= 1 && order < 16) {
465 alpha = alpopt[order-1];
471 const ordinal_type p = order;
472 ordinal_type N = (p+1)*(p+2)/2;
475 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1(
"L1", N );
476 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2(
"L2", N );
477 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3(
"L3", N );
478 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> X(
"X", N);
479 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> Y(
"Y", N);
482 for (ordinal_type n=1;n<=p+1;n++) {
483 for (ordinal_type m=1;m<=p+2-n;m++) {
484 L1(sk) = (n-1) / (pointValueType)p;
485 L3(sk) = (m-1) / (pointValueType)p;
486 L2(sk) = 1.0 - L1(sk) - L3(sk);
491 for (ordinal_type n=0;n<N;n++) {
492 X(n) = -L2(n) + L3(n);
493 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
497 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend1(
"blend1", N);
498 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend2(
"blend2", N);
499 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend3(
"blend3", N);
501 for (ordinal_type n=0;n<N;n++) {
502 blend1(n) = 4.0 * L2(n) * L3(n);
503 blend2(n) = 4.0 * L1(n) * L3(n);
504 blend3(n) = 4.0 * L1(n) * L2(n);
508 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2", N);
509 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3", N);
510 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1", N);
512 for (ordinal_type k=0;k<N;k++) {
513 L3mL2(k) = L3(k)-L2(k);
514 L1mL3(k) = L1(k)-L3(k);
515 L2mL1(k) = L2(k)-L1(k);
518 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1", N);
519 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2", N);
520 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3", N);
522 warpFactor( warpfactor1, order , gaussX , L3mL2 );
523 warpFactor( warpfactor2, order , gaussX , L1mL3 );
524 warpFactor( warpfactor3, order , gaussX , L2mL1 );
526 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp1(
"warp1", N);
527 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp2(
"warp2", N);
528 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp3(
"warp3", N);
530 for (ordinal_type k=0;k<N;k++) {
531 warp1(k) = blend1(k) * warpfactor1(k) *
532 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
533 warp2(k) = blend2(k) * warpfactor2(k) *
534 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
535 warp3(k) = blend3(k) * warpfactor3(k) *
536 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
539 for (ordinal_type k=0;k<N;k++) {
540 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
541 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
544 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warXY(
"warXY", 1, N,2);
546 for (ordinal_type k=0;k<N;k++) {
547 warXY(0, k,0) = X(k);
548 warXY(0, k,1) = Y(k);
553 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warburtonVerts(
"warburtonVerts", 1,3,2);
554 warburtonVerts(0,0,0) = -1.0;
555 warburtonVerts(0,0,1) = -1.0/std::sqrt(3.0);
556 warburtonVerts(0,1,0) = 1.0;
557 warburtonVerts(0,1,1) = -1.0/std::sqrt(3.0);
558 warburtonVerts(0,2,0) = 0.0;
559 warburtonVerts(0,2,1) = 2.0/std::sqrt(3.0);
561 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> refPts(
"refPts", 1, N,2);
566 shards::getCellTopologyData< shards::Triangle<3> >()
570 auto pointsHost = Kokkos::create_mirror_view(points);
572 ordinal_type noffcur = 0;
573 ordinal_type offcur = 0;
574 for (ordinal_type i=0;i<=order;i++) {
575 for (ordinal_type j=0;j<=order-i;j++) {
576 if ( (i >= offset) && (i <= order-offset) &&
577 (j >= offset) && (j <= order-i-offset) ) {
578 pointsHost(offcur,0) = refPts(0, noffcur,0);
579 pointsHost(offcur,1) = refPts(0, noffcur,1);
586 Kokkos::deep_copy(points, pointsHost);
591template<
typename pointValueType,
class ...pointProperties>
594warpShiftFace3D( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
595 const ordinal_type order ,
596 const pointValueType pval ,
597 const Kokkos::DynRankView<pointValueType,pointProperties...> ,
598 const Kokkos::DynRankView<pointValueType,pointProperties...> L2,
599 const Kokkos::DynRankView<pointValueType,pointProperties...> L3,
600 const Kokkos::DynRankView<pointValueType,pointProperties...> L4
607template<
typename pointValueType,
class ...pointProperties>
610evalshift( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
611 const ordinal_type order ,
612 const pointValueType pval ,
613 const Kokkos::DynRankView<pointValueType,pointProperties...> L1 ,
614 const Kokkos::DynRankView<pointValueType,pointProperties...> L2 ,
615 const Kokkos::DynRankView<pointValueType,pointProperties...> L3
619 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX",order+1,1);
622 const ordinal_type N = L1.extent(0);
625 for (ordinal_type k=0;k<=order;k++) {
630 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend1(
"blend1",N);
631 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend2(
"blend2",N);
632 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend3(
"blend3",N);
634 for (ordinal_type i=0;i<N;i++) {
635 blend1(i) = L2(i) * L3(i);
636 blend2(i) = L1(i) * L3(i);
637 blend3(i) = L1(i) * L2(i);
641 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1",N);
642 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2",N);
643 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3",N);
646 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2",N);
647 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3",N);
648 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1",N);
650 for (ordinal_type k=0;k<N;k++) {
651 L3mL2(k) = L3(k)-L2(k);
652 L1mL3(k) = L1(k)-L3(k);
653 L2mL1(k) = L2(k)-L1(k);
656 evalwarp( warpfactor1 , order , gaussX , L3mL2 );
657 evalwarp( warpfactor2 , order , gaussX , L1mL3 );
658 evalwarp( warpfactor3 , order , gaussX , L2mL1 );
660 for (ordinal_type k=0;k<N;k++) {
661 warpfactor1(k) *= 4.0;
662 warpfactor2(k) *= 4.0;
663 warpfactor3(k) *= 4.0;
666 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp1(
"warp1",N);
667 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp2(
"warp2",N);
668 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp3(
"warp3",N);
670 for (ordinal_type k=0;k<N;k++) {
671 warp1(k) = blend1(k) * warpfactor1(k) *
672 ( 1.0 + pval * pval * L1(k) * L1(k) );
673 warp2(k) = blend2(k) * warpfactor2(k) *
674 ( 1.0 + pval * pval * L2(k) * L2(k) );
675 warp3(k) = blend3(k) * warpfactor3(k) *
676 ( 1.0 + pval * pval * L3(k) * L3(k) );
679 for (ordinal_type k=0;k<N;k++) {
680 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);
681 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);
686template<
typename pointValueType,
class ...pointProperties>
689evalwarp(Kokkos::DynRankView<pointValueType,pointProperties...> warp ,
690 const ordinal_type order ,
691 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes,
692 const Kokkos::DynRankView<pointValueType,pointProperties...> xout )
694 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> xeq(
"xeq",order+1);
696 ordinal_type xout_dim0 = xout.extent(0);
697 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> d(
"d",xout_dim0);
701 for (ordinal_type i=0;i<=order;i++) {
702 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
705 for (ordinal_type i=0;i<=order;i++) {
706 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
707 for (ordinal_type j=1;j<order;j++) {
709 for (ordinal_type k=0;k<xout_dim0;k++) {
710 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
715 for (ordinal_type k=0;k<xout_dim0;k++) {
716 d(k) = -d(k)/(xeq(i)-xeq(0));
720 for (ordinal_type k=0;k<xout_dim0;k++) {
721 d(k) = d(k)/(xeq(i)-xeq(order));
725 for (ordinal_type k=0;k<xout_dim0;k++) {
732template<
typename pointValueType,
class ...pointProperties>
736 const ordinal_type order ,
737 const ordinal_type offset )
739 pointValueType alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
740 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
741 pointValueType alpha;
744 alpha = alphastore[std::max(order-1,0)];
750 const ordinal_type N = (order+1)*(order+2)*(order+3)/6;
751 pointValueType tol = 1.e-10;
753 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> shift(
"shift",N,3);
754 Kokkos::deep_copy(shift,0.0);
757 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> equipoints(
"equipoints", N,3);
759 for (ordinal_type n=0;n<=order;n++) {
760 for (ordinal_type m=0;m<=order-n;m++) {
761 for (ordinal_type q=0;q<=order-n-m;q++) {
762 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
763 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
764 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
772 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1(
"L1",N);
773 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2(
"L2",N);
774 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3(
"L3",N);
775 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L4(
"L4",N);
776 for (ordinal_type i=0;i<N;i++) {
777 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
778 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
779 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
780 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
785 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warVerts_(
"warVerts",1,4,3);
786 auto warVerts = Kokkos::subview(warVerts_,0,Kokkos::ALL(),Kokkos::ALL());
787 warVerts(0,0) = -1.0;
788 warVerts(0,1) = -1.0/sqrt(3.0);
789 warVerts(0,2) = -1.0/sqrt(6.0);
791 warVerts(1,1) = -1.0/sqrt(3.0);
792 warVerts(1,2) = -1.0/sqrt(6.0);
794 warVerts(2,1) = 2.0 / sqrt(3.0);
795 warVerts(2,2) = -1.0/sqrt(6.0);
798 warVerts(3,2) = 3.0 / sqrt(6.0);
802 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t1(
"t1",4,3);
803 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t2(
"t2",4,3);
804 for (ordinal_type i=0;i<3;i++) {
805 t1(0,i) = warVerts(1,i) - warVerts(0,i);
806 t1(1,i) = warVerts(1,i) - warVerts(0,i);
807 t1(2,i) = warVerts(2,i) - warVerts(1,i);
808 t1(3,i) = warVerts(2,i) - warVerts(0,i);
809 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
810 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
811 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
812 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
816 for (ordinal_type n=0;n<4;n++) {
818 pointValueType normt1n = 0.0;
819 pointValueType normt2n = 0.0;
820 for (ordinal_type i=0;i<3;i++) {
821 normt1n += (t1(n,i) * t1(n,i));
822 normt2n += (t2(n,i) * t2(n,i));
824 normt1n = sqrt(normt1n);
825 normt2n = sqrt(normt2n);
827 for (ordinal_type i=0;i<3;i++) {
834 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> XYZ(
"XYZ",N,3);
835 for (ordinal_type i=0;i<N;i++) {
836 for (ordinal_type j=0;j<3;j++) {
837 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
841 for (ordinal_type face=1;face<=4;face++) {
842 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> La, Lb, Lc, Ld;
843 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp(
"warp",N,2);
844 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend(
"blend",N);
845 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> denom(
"denom",N);
848 La = L1; Lb = L2; Lc = L3; Ld = L4;
break;
850 La = L2; Lb = L1; Lc = L3; Ld = L4;
break;
852 La = L3; Lb = L1; Lc = L4; Ld = L2;
break;
854 La = L4; Lb = L1; Lc = L3; Ld = L2;
break;
860 for (ordinal_type k=0;k<N;k++) {
861 blend(k) = Lb(k) * Lc(k) * Ld(k);
864 for (ordinal_type k=0;k<N;k++) {
865 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
868 for (ordinal_type k=0;k<N;k++) {
869 if (denom(k) > tol) {
870 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
876 for (ordinal_type k=0;k<N;k++) {
877 for (ordinal_type j=0;j<3;j++) {
878 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
879 + blend(k) * warp(k,1) * t2(face-1,j);
883 for (ordinal_type k=0;k<N;k++) {
884 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
885 for (ordinal_type j=0;j<3;j++) {
886 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
893 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> updatedPoints(
"updatedPoints",1,N,3);
894 for (ordinal_type k=0;k<N;k++) {
895 for (ordinal_type j=0;j<3;j++) {
896 updatedPoints(0,k,j) = XYZ(k,j) + shift(k,j);
903 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> refPts(
"refPts",1,N,3);
906 shards::getCellTopologyData<shards::Tetrahedron<4> >()
909 auto pointsHost = Kokkos::create_mirror_view(points);
911 ordinal_type noffcur = 0;
912 ordinal_type offcur = 0;
913 for (ordinal_type i=0;i<=order;i++) {
914 for (ordinal_type j=0;j<=order-i;j++) {
915 for (ordinal_type k=0;k<=order-i-j;k++) {
916 if ( (i >= offset) && (i <= order-offset) &&
917 (j >= offset) && (j <= order-i-offset) &&
918 (k >= offset) && (k <= order-i-j-offset) ) {
919 pointsHost(offcur,0) = refPts(0,noffcur,0);
920 pointsHost(offcur,1) = refPts(0,noffcur,1);
921 pointsHost(offcur,2) = refPts(0,noffcur,2);
929 Kokkos::deep_copy(points, pointsHost);
Kokkos::DynRankView< typename ViewType::value_type, typename DeduceLayout< ViewType >::result_layout, typename ViewType::device_type > getMatchingViewWithLabel(const ViewType &view, const std::string &label, DimArgs... dims)
Creates and returns a view that matches the provided view in Kokkos Layout.
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.