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." );
111 template<
typename pointValueType,
class ...pointProperties>
114 getLattice( Kokkos::DynRankView<pointValueType,pointProperties...> points,
115 const shards::CellTopology cell,
116 const ordinal_type order,
117 const ordinal_type offset,
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." );
182 template<
typename pointValueType,
class ...pointProperties>
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);
219 template<
typename pointValueType,
class ...pointProperties>
223 const ordinal_type order,
224 const ordinal_type offset,
230 INTREPID2_TEST_FOR_EXCEPTION(
true ,
231 std::invalid_argument ,
232 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
237 template<
typename pointValueType,
class ...pointProperties>
241 const ordinal_type order,
242 const ordinal_type offset,
248 INTREPID2_TEST_FOR_EXCEPTION(
true ,
249 std::invalid_argument ,
250 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
255 template<
typename pointValueType,
class ...pointProperties>
259 const ordinal_type order,
260 const ordinal_type offset,
266 INTREPID2_TEST_FOR_EXCEPTION(
true ,
267 std::invalid_argument ,
268 ">>> ERROR (PointTools::getLattice): invalid EPointType." );
275 template<
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);
295 template<
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 Kokkos::deep_copy(pts, z);
328 template<
typename pointValueType,
class ...pointProperties>
332 const ordinal_type order,
333 const ordinal_type offset ) {
334 TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
335 std::invalid_argument ,
336 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTriangle): order must be positive" );
338 auto pointsHost = Kokkos::create_mirror_view(points);
340 const pointValueType h = 1.0 / order;
341 ordinal_type cur = 0;
343 for (ordinal_type i=offset;i<=order-offset;i++) {
344 for (ordinal_type j=offset;j<=order-i-offset;j++) {
345 pointsHost(cur,0) = j * h ;
346 pointsHost(cur,1) = i * h;
350 Kokkos::deep_copy(points, pointsHost);
353 template<
typename pointValueType,
class ...pointProperties>
357 const ordinal_type order ,
358 const ordinal_type offset ) {
359 TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
360 std::invalid_argument ,
361 ">>> ERROR (Intrepid2::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
363 auto pointsHost = Kokkos::create_mirror_view(points);
365 const pointValueType h = 1.0 / order;
366 ordinal_type cur = 0;
368 for (ordinal_type i=offset;i<=order-offset;i++) {
369 for (ordinal_type j=offset;j<=order-i-offset;j++) {
370 for (ordinal_type k=offset;k<=order-i-j-offset;k++) {
371 pointsHost(cur,0) = k * h;
372 pointsHost(cur,1) = j * h;
373 pointsHost(cur,2) = i * h;
378 Kokkos::deep_copy(points, pointsHost);
382 template<
typename pointValueType,
class ...pointProperties>
385 warpFactor( Kokkos::DynRankView<pointValueType,pointProperties...> warp,
386 const ordinal_type order,
387 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes ,
388 const Kokkos::DynRankView<pointValueType,pointProperties...> xout
391 TEUCHOS_TEST_FOR_EXCEPTION( ( warp.extent(0) != xout.extent(0) ) ,
392 std::invalid_argument ,
393 ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
395 Kokkos::deep_copy(warp, pointValueType(0.0));
397 ordinal_type xout_dim0 = xout.extent(0);
398 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> d(
"d", xout_dim0 );
400 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> xeq_(
"xeq", order + 1 ,1);
402 const auto xeq = Kokkos::subview(xeq_, Kokkos::ALL(),0);
404 TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.extent(0) != xnodes.extent(0) ) ,
405 std::invalid_argument ,
406 ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
409 for (ordinal_type i=0;i<=order;i++) {
411 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
413 for (ordinal_type j=1;j<order;j++) {
415 for (ordinal_type k=0;k<xout_dim0;k++) {
416 d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
423 for (ordinal_type k=0;k<xout_dim0;k++) {
424 d(k) = -d(k) / (xeq(i) - xeq(0));
429 for (ordinal_type k=0;k<xout_dim0;k++) {
430 d(k) = d(k) / (xeq(i) - xeq(order));
434 for (ordinal_type k=0;k<xout_dim0;k++) {
441 template<
typename pointValueType,
class ...pointProperties>
445 const ordinal_type order ,
446 const ordinal_type offset)
450 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX", order + 1 , 1 );
457 pointValueType alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
458 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
460 pointValueType alpha;
462 if (order >= 1 && order < 16) {
463 alpha = alpopt[order-1];
469 const ordinal_type p = order;
470 ordinal_type N = (p+1)*(p+2)/2;
473 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1(
"L1", N );
474 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2(
"L2", N );
475 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3(
"L3", N );
476 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> X(
"X", N);
477 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> Y(
"Y", N);
480 for (ordinal_type n=1;n<=p+1;n++) {
481 for (ordinal_type m=1;m<=p+2-n;m++) {
482 L1(sk) = (n-1) / (pointValueType)p;
483 L3(sk) = (m-1) / (pointValueType)p;
484 L2(sk) = 1.0 - L1(sk) - L3(sk);
489 for (ordinal_type n=0;n<N;n++) {
490 X(n) = -L2(n) + L3(n);
491 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
495 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend1(
"blend1", N);
496 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend2(
"blend2", N);
497 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> blend3(
"blend3", N);
499 for (ordinal_type n=0;n<N;n++) {
500 blend1(n) = 4.0 * L2(n) * L3(n);
501 blend2(n) = 4.0 * L1(n) * L3(n);
502 blend3(n) = 4.0 * L1(n) * L2(n);
506 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2", N);
507 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3", N);
508 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1", N);
510 for (ordinal_type k=0;k<N;k++) {
511 L3mL2(k) = L3(k)-L2(k);
512 L1mL3(k) = L1(k)-L3(k);
513 L2mL1(k) = L2(k)-L1(k);
516 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1", N);
517 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2", N);
518 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3", N);
520 warpFactor( warpfactor1, order , gaussX , L3mL2 );
521 warpFactor( warpfactor2, order , gaussX , L1mL3 );
522 warpFactor( warpfactor3, order , gaussX , L2mL1 );
524 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp1(
"warp1", N);
525 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp2(
"warp2", N);
526 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warp3(
"warp3", N);
528 for (ordinal_type k=0;k<N;k++) {
529 warp1(k) = blend1(k) * warpfactor1(k) *
530 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
531 warp2(k) = blend2(k) * warpfactor2(k) *
532 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
533 warp3(k) = blend3(k) * warpfactor3(k) *
534 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
537 for (ordinal_type k=0;k<N;k++) {
538 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
539 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
542 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warXY(
"warXY", 1, N,2);
544 for (ordinal_type k=0;k<N;k++) {
545 warXY(0, k,0) = X(k);
546 warXY(0, k,1) = Y(k);
551 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> warburtonVerts(
"warburtonVerts", 1,3,2);
552 warburtonVerts(0,0,0) = -1.0;
553 warburtonVerts(0,0,1) = -1.0/std::sqrt(3.0);
554 warburtonVerts(0,1,0) = 1.0;
555 warburtonVerts(0,1,1) = -1.0/std::sqrt(3.0);
556 warburtonVerts(0,2,0) = 0.0;
557 warburtonVerts(0,2,1) = 2.0/std::sqrt(3.0);
559 Kokkos::DynRankView<pointValueType, Kokkos::DefaultHostExecutionSpace> refPts(
"refPts", 1, N,2);
564 shards::getCellTopologyData< shards::Triangle<3> >()
568 auto pointsHost = Kokkos::create_mirror_view(points);
570 ordinal_type noffcur = 0;
571 ordinal_type offcur = 0;
572 for (ordinal_type i=0;i<=order;i++) {
573 for (ordinal_type j=0;j<=order-i;j++) {
574 if ( (i >= offset) && (i <= order-offset) &&
575 (j >= offset) && (j <= order-i-offset) ) {
576 pointsHost(offcur,0) = refPts(0, noffcur,0);
577 pointsHost(offcur,1) = refPts(0, noffcur,1);
584 Kokkos::deep_copy(points, pointsHost);
589 template<
typename pointValueType,
class ...pointProperties>
593 const ordinal_type order ,
594 const pointValueType pval ,
595 const Kokkos::DynRankView<pointValueType,pointProperties...> ,
596 const Kokkos::DynRankView<pointValueType,pointProperties...> L2,
597 const Kokkos::DynRankView<pointValueType,pointProperties...> L3,
598 const Kokkos::DynRankView<pointValueType,pointProperties...> L4
605 template<
typename pointValueType,
class ...pointProperties>
608 evalshift( Kokkos::DynRankView<pointValueType,pointProperties...> dxy,
609 const ordinal_type order ,
610 const pointValueType pval ,
611 const Kokkos::DynRankView<pointValueType,pointProperties...> L1 ,
612 const Kokkos::DynRankView<pointValueType,pointProperties...> L2 ,
613 const Kokkos::DynRankView<pointValueType,pointProperties...> L3
617 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> gaussX(
"gaussX",order+1,1);
620 const ordinal_type N = L1.extent(0);
623 for (ordinal_type k=0;k<=order;k++) {
628 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend1(
"blend1",N);
629 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend2(
"blend2",N);
630 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend3(
"blend3",N);
632 for (ordinal_type i=0;i<N;i++) {
633 blend1(i) = L2(i) * L3(i);
634 blend2(i) = L1(i) * L3(i);
635 blend3(i) = L1(i) * L2(i);
639 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor1(
"warpfactor1",N);
640 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor2(
"warpfactor2",N);
641 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warpfactor3(
"warpfactor3",N);
644 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3mL2(
"L3mL2",N);
645 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1mL3(
"L1mL3",N);
646 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2mL1(
"L2mL1",N);
648 for (ordinal_type k=0;k<N;k++) {
649 L3mL2(k) = L3(k)-L2(k);
650 L1mL3(k) = L1(k)-L3(k);
651 L2mL1(k) = L2(k)-L1(k);
654 evalwarp( warpfactor1 , order , gaussX , L3mL2 );
655 evalwarp( warpfactor2 , order , gaussX , L1mL3 );
656 evalwarp( warpfactor3 , order , gaussX , L2mL1 );
658 for (ordinal_type k=0;k<N;k++) {
659 warpfactor1(k) *= 4.0;
660 warpfactor2(k) *= 4.0;
661 warpfactor3(k) *= 4.0;
664 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp1(
"warp1",N);
665 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp2(
"warp2",N);
666 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp3(
"warp3",N);
668 for (ordinal_type k=0;k<N;k++) {
669 warp1(k) = blend1(k) * warpfactor1(k) *
670 ( 1.0 + pval * pval * L1(k) * L1(k) );
671 warp2(k) = blend2(k) * warpfactor2(k) *
672 ( 1.0 + pval * pval * L2(k) * L2(k) );
673 warp3(k) = blend3(k) * warpfactor3(k) *
674 ( 1.0 + pval * pval * L3(k) * L3(k) );
677 for (ordinal_type k=0;k<N;k++) {
678 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);
679 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);
684 template<
typename pointValueType,
class ...pointProperties>
687 evalwarp(Kokkos::DynRankView<pointValueType,pointProperties...> warp ,
688 const ordinal_type order ,
689 const Kokkos::DynRankView<pointValueType,pointProperties...> xnodes,
690 const Kokkos::DynRankView<pointValueType,pointProperties...> xout )
692 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> xeq(
"xeq",order+1);
694 ordinal_type xout_dim0 = xout.extent(0);
695 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> d(
"d",xout_dim0);
699 for (ordinal_type i=0;i<=order;i++) {
700 xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
703 for (ordinal_type i=0;i<=order;i++) {
704 Kokkos::deep_copy(d, xnodes(i,0) - xeq(i));
705 for (ordinal_type j=1;j<order;j++) {
707 for (ordinal_type k=0;k<xout_dim0;k++) {
708 d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
713 for (ordinal_type k=0;k<xout_dim0;k++) {
714 d(k) = -d(k)/(xeq(i)-xeq(0));
718 for (ordinal_type k=0;k<xout_dim0;k++) {
719 d(k) = d(k)/(xeq(i)-xeq(order));
723 for (ordinal_type k=0;k<xout_dim0;k++) {
730 template<
typename pointValueType,
class ...pointProperties>
734 const ordinal_type order ,
735 const ordinal_type offset )
737 pointValueType alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
738 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
739 pointValueType alpha;
742 alpha = alphastore[std::max(order-1,0)];
748 const ordinal_type N = (order+1)*(order+2)*(order+3)/6;
749 pointValueType tol = 1.e-10;
751 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> shift(
"shift",N,3);
752 Kokkos::deep_copy(shift,0.0);
755 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> equipoints(
"equipoints", N,3);
757 for (ordinal_type n=0;n<=order;n++) {
758 for (ordinal_type m=0;m<=order-n;m++) {
759 for (ordinal_type q=0;q<=order-n-m;q++) {
760 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
761 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
762 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
770 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L1(
"L1",N);
771 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L2(
"L2",N);
772 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L3(
"L3",N);
773 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> L4(
"L4",N);
774 for (ordinal_type i=0;i<N;i++) {
775 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
776 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
777 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
778 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
783 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warVerts_(
"warVerts",1,4,3);
784 auto warVerts = Kokkos::subview(warVerts_,0,Kokkos::ALL(),Kokkos::ALL());
785 warVerts(0,0) = -1.0;
786 warVerts(0,1) = -1.0/sqrt(3.0);
787 warVerts(0,2) = -1.0/sqrt(6.0);
789 warVerts(1,1) = -1.0/sqrt(3.0);
790 warVerts(1,2) = -1.0/sqrt(6.0);
792 warVerts(2,1) = 2.0 / sqrt(3.0);
793 warVerts(2,2) = -1.0/sqrt(6.0);
796 warVerts(3,2) = 3.0 / sqrt(6.0);
800 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t1(
"t1",4,3);
801 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> t2(
"t2",4,3);
802 for (ordinal_type i=0;i<3;i++) {
803 t1(0,i) = warVerts(1,i) - warVerts(0,i);
804 t1(1,i) = warVerts(1,i) - warVerts(0,i);
805 t1(2,i) = warVerts(2,i) - warVerts(1,i);
806 t1(3,i) = warVerts(2,i) - warVerts(0,i);
807 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
808 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
809 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
810 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
814 for (ordinal_type n=0;n<4;n++) {
816 pointValueType normt1n = 0.0;
817 pointValueType normt2n = 0.0;
818 for (ordinal_type i=0;i<3;i++) {
819 normt1n += (t1(n,i) * t1(n,i));
820 normt2n += (t2(n,i) * t2(n,i));
822 normt1n = sqrt(normt1n);
823 normt2n = sqrt(normt2n);
825 for (ordinal_type i=0;i<3;i++) {
832 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> XYZ(
"XYZ",N,3);
833 for (ordinal_type i=0;i<N;i++) {
834 for (ordinal_type j=0;j<3;j++) {
835 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
839 for (ordinal_type face=1;face<=4;face++) {
840 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> La, Lb, Lc, Ld;
841 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> warp(
"warp",N,2);
842 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> blend(
"blend",N);
843 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> denom(
"denom",N);
846 La = L1; Lb = L2; Lc = L3; Ld = L4;
break;
848 La = L2; Lb = L1; Lc = L3; Ld = L4;
break;
850 La = L3; Lb = L1; Lc = L4; Ld = L2;
break;
852 La = L4; Lb = L1; Lc = L3; Ld = L2;
break;
858 for (ordinal_type k=0;k<N;k++) {
859 blend(k) = Lb(k) * Lc(k) * Ld(k);
862 for (ordinal_type k=0;k<N;k++) {
863 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
866 for (ordinal_type k=0;k<N;k++) {
867 if (denom(k) > tol) {
868 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
874 for (ordinal_type k=0;k<N;k++) {
875 for (ordinal_type j=0;j<3;j++) {
876 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
877 + blend(k) * warp(k,1) * t2(face-1,j);
881 for (ordinal_type k=0;k<N;k++) {
882 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
883 for (ordinal_type j=0;j<3;j++) {
884 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
891 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> updatedPoints(
"updatedPoints",1,N,3);
892 for (ordinal_type k=0;k<N;k++) {
893 for (ordinal_type j=0;j<3;j++) {
894 updatedPoints(0,k,j) = XYZ(k,j) + shift(k,j);
901 Kokkos::DynRankView<pointValueType,Kokkos::DefaultHostExecutionSpace> refPts(
"refPts",1,N,3);
904 shards::getCellTopologyData<shards::Tetrahedron<4> >()
907 auto pointsHost = Kokkos::create_mirror_view(points);
909 ordinal_type noffcur = 0;
910 ordinal_type offcur = 0;
911 for (ordinal_type i=0;i<=order;i++) {
912 for (ordinal_type j=0;j<=order-i;j++) {
913 for (ordinal_type k=0;k<=order-i-j;k++) {
914 if ( (i >= offset) && (i <= order-offset) &&
915 (j >= offset) && (j <= order-i-offset) &&
916 (k >= offset) && (k <= order-i-j-offset) ) {
917 pointsHost(offcur,0) = refPts(0,noffcur,0);
918 pointsHost(offcur,1) = refPts(0,noffcur,1);
919 pointsHost(offcur,2) = refPts(0,noffcur,2);
927 Kokkos::deep_copy(points, pointsHost);
Gauss-Jacobi/Gauss-Radau-Jacobi/Gauss-Lobatto zeros and weights.
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.
EPointType
Enumeration of types of point distributions in Intrepid.