119getLattice( 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 );
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 );
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." );
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);
662evalshift( 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);
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);