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." );
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);
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);
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);