99 switch (cellType.getKey()) {
100 case shards::Tetrahedron<4>::key:
101 case shards::Tetrahedron<8>::key:
102 case shards::Tetrahedron<10>::key:
103 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
104 || points.dimension(1) != 3 ,
105 std::invalid_argument ,
106 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
107 getEquispacedLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
109 case shards::Triangle<3>::key:
110 case shards::Triangle<4>::key:
111 case shards::Triangle<6>::key:
112 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
113 || points.dimension(1) != 2 ,
114 std::invalid_argument ,
115 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
116 getEquispacedLatticeTriangle<Scalar,ArrayType>( points , order , offset );
118 case shards::Line<2>::key:
119 case shards::Line<3>::key:
120 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
121 || points.dimension(1) != 1 ,
122 std::invalid_argument ,
123 ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
124 getEquispacedLatticeLine<Scalar,ArrayType>( points , order , offset );
127 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
128 ">>> ERROR (Intrepid::PointTools::getEquispacedLattice): Illegal cell type" );
141 switch (cellType.getKey()) {
142 case shards::Tetrahedron<4>::key:
143 case shards::Tetrahedron<8>::key:
144 case shards::Tetrahedron<10>::key:
145 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
146 || points.dimension(1) != 3 ,
147 std::invalid_argument ,
148 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
150 getWarpBlendLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
152 case shards::Triangle<3>::key:
153 case shards::Triangle<4>::key:
154 case shards::Triangle<6>::key:
155 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
156 || points.dimension(1) != 2 ,
157 std::invalid_argument ,
158 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
160 getWarpBlendLatticeTriangle<Scalar,ArrayType>( points , order , offset );
162 case shards::Line<2>::key:
163 case shards::Line<3>::key:
164 TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) !=
getLatticeSize( cellType , order , offset ) )
165 || points.dimension(1) != 1 ,
166 std::invalid_argument ,
167 ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
169 getWarpBlendLatticeLine<Scalar,ArrayType>( points , order , offset );
172 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument ,
173 ">>> ERROR (Intrepid::PointTools::getWarpBlendLattice): Illegal cell type" );
340 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
342 gaussX.resize(gaussX.dimension(0));
344 Scalar alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
345 1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
349 if (order >= 1 && order < 16) {
350 alpha = alpopt[order-1];
357 int N = (p+1)*(p+2)/2;
367 for (
int n=1;n<=p+1;n++) {
368 for (
int m=1;m<=p+2-n;m++) {
369 L1(sk) = (n-1) / (Scalar)p;
370 L3(sk) = (m-1) / (Scalar)p;
371 L2(sk) = 1.0 - L1(sk) - L3(sk);
376 for (
int n=0;n<N;n++) {
377 X(n) = -L2(n) + L3(n);
378 Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
386 for (
int n=0;n<N;n++) {
387 blend1(n) = 4.0 * L2(n) * L3(n);
388 blend2(n) = 4.0 * L1(n) * L3(n);
389 blend3(n) = 4.0 * L1(n) * L2(n);
397 for (
int k=0;k<N;k++) {
398 L3mL2(k) = L3(k)-L2(k);
399 L1mL3(k) = L1(k)-L3(k);
400 L2mL1(k) = L2(k)-L1(k);
403 FieldContainer<Scalar> warpfactor1(N);
404 FieldContainer<Scalar> warpfactor2(N);
405 FieldContainer<Scalar> warpfactor3(N);
407 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L3mL2 , warpfactor1 );
408 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L1mL3 , warpfactor2 );
409 warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L2mL1 , warpfactor3 );
411 FieldContainer<Scalar> warp1(N);
412 FieldContainer<Scalar> warp2(N);
413 FieldContainer<Scalar> warp3(N);
415 for (
int k=0;k<N;k++) {
416 warp1(k) = blend1(k) * warpfactor1(k) *
417 ( 1.0 + alpha * alpha * L1(k) * L1(k) );
418 warp2(k) = blend2(k) * warpfactor2(k) *
419 ( 1.0 + alpha * alpha * L2(k) * L2(k) );
420 warp3(k) = blend3(k) * warpfactor3(k) *
421 ( 1.0 + alpha * alpha * L3(k) * L3(k) );
424 for (
int k=0;k<N;k++) {
425 X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
426 Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
429 FieldContainer<Scalar> warXY(N,2);
431 for (
int k=0;k<N;k++) {
438 FieldContainer<Scalar> warburtonVerts(1,3,2);
439 warburtonVerts(0,0,0) = -1.0;
440 warburtonVerts(0,0,1) = -1.0/sqrt(3.0);
441 warburtonVerts(0,1,0) = 1.0;
442 warburtonVerts(0,1,1) = -1.0/sqrt(3.0);
443 warburtonVerts(0,2,0) = 0.0;
444 warburtonVerts(0,2,1) = 2.0/sqrt(3.0);
446 FieldContainer<Scalar> refPts(N,2);
451 shards::getCellTopologyData< shards::Triangle<3> >(),
457 for (
int i=0;i<=order;i++) {
458 for (
int j=0;j<=order-i;j++) {
459 if ( (i >= offset) && (i <= order-offset) &&
460 (j >= offset) && (j <= order-i-offset) ) {
461 points(offcur,0) = refPts(noffcur,0);
462 points(offcur,1) = refPts(noffcur,1);
489 const ArrayType &L1 ,
490 const ArrayType &L2 ,
491 const ArrayType &L3 ,
495 FieldContainer<Scalar> gaussX(order+1,1);
496 PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
497 gaussX.resize(order+1);
498 const int N = L1.dimension(0);
501 for (
int k=0;k<=order;k++) {
506 FieldContainer<Scalar> blend1(N);
507 FieldContainer<Scalar> blend2(N);
508 FieldContainer<Scalar> blend3(N);
510 for (
int i=0;i<N;i++) {
511 blend1(i) = L2(i) * L3(i);
512 blend2(i) = L1(i) * L3(i);
513 blend3(i) = L1(i) * L2(i);
517 FieldContainer<Scalar> warpfactor1(N);
518 FieldContainer<Scalar> warpfactor2(N);
519 FieldContainer<Scalar> warpfactor3(N);
522 FieldContainer<Scalar> L3mL2(N);
523 FieldContainer<Scalar> L1mL3(N);
524 FieldContainer<Scalar> L2mL1(N);
526 for (
int k=0;k<N;k++) {
527 L3mL2(k) = L3(k)-L2(k);
528 L1mL3(k) = L1(k)-L3(k);
529 L2mL1(k) = L2(k)-L1(k);
532 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor1 , order , gaussX , L3mL2 );
533 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor2 , order , gaussX , L1mL3 );
534 evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor3 , order , gaussX , L2mL1 );
536 for (
int k=0;k<N;k++) {
537 warpfactor1(k) *= 4.0;
538 warpfactor2(k) *= 4.0;
539 warpfactor3(k) *= 4.0;
542 FieldContainer<Scalar> warp1(N);
543 FieldContainer<Scalar> warp2(N);
544 FieldContainer<Scalar> warp3(N);
546 for (
int k=0;k<N;k++) {
547 warp1(k) = blend1(k) * warpfactor1(k) *
548 ( 1.0 + pval * pval * L1(k) * L1(k) );
549 warp2(k) = blend2(k) * warpfactor2(k) *
550 ( 1.0 + pval * pval * L2(k) * L2(k) );
551 warp3(k) = blend3(k) * warpfactor3(k) *
552 ( 1.0 + pval * pval * L3(k) * L3(k) );
555 for (
int k=0;k<N;k++) {
556 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);
557 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);
616 Scalar alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
617 1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
621 alpha = alphastore[order-1];
627 const int N = (order+1)*(order+2)*(order+3)/6;
630 FieldContainer<Scalar> shift(N,3);
634 FieldContainer<Scalar> equipoints(N,3);
636 for (
int n=0;n<=order;n++) {
637 for (
int m=0;m<=order-n;m++) {
638 for (
int q=0;q<=order-n-m;q++) {
639 equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
640 equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
641 equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
649 FieldContainer<Scalar> L1(N);
650 FieldContainer<Scalar> L2(N);
651 FieldContainer<Scalar> L3(N);
652 FieldContainer<Scalar> L4(N);
653 for (
int i=0;i<N;i++) {
654 L1(i) = (1.0 + equipoints(i,2)) / 2.0;
655 L2(i) = (1.0 + equipoints(i,1)) / 2.0;
656 L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
657 L4(i) = (1.0 + equipoints(i,0)) / 2.0;
662 FieldContainer<Scalar> warVerts(4,3);
663 warVerts(0,0) = -1.0;
664 warVerts(0,1) = -1.0/sqrt(3.0);
665 warVerts(0,2) = -1.0/sqrt(6.0);
667 warVerts(1,1) = -1.0/sqrt(3.0);
668 warVerts(1,2) = -1.0/sqrt(6.0);
670 warVerts(2,1) = 2.0 / sqrt(3.0);
671 warVerts(2,2) = -1.0/sqrt(6.0);
674 warVerts(3,2) = 3.0 / sqrt(6.0);
678 FieldContainer<Scalar> t1(4,3);
679 FieldContainer<Scalar> t2(4,3);
680 for (
int i=0;i<3;i++) {
681 t1(0,i) = warVerts(1,i) - warVerts(0,i);
682 t1(1,i) = warVerts(1,i) - warVerts(0,i);
683 t1(2,i) = warVerts(2,i) - warVerts(1,i);
684 t1(3,i) = warVerts(2,i) - warVerts(0,i);
685 t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
686 t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
687 t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
688 t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
692 for (
int n=0;n<4;n++) {
694 Scalar normt1n = 0.0;
695 Scalar normt2n = 0.0;
696 for (
int i=0;i<3;i++) {
697 normt1n += (t1(n,i) * t1(n,i));
698 normt2n += (t2(n,i) * t2(n,i));
700 normt1n = sqrt(normt1n);
701 normt2n = sqrt(normt2n);
703 for (
int i=0;i<3;i++) {
710 FieldContainer<Scalar> XYZ(N,3);
711 for (
int i=0;i<N;i++) {
712 for (
int j=0;j<3;j++) {
713 XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
717 for (
int face=1;face<=4;face++) {
718 FieldContainer<Scalar> La, Lb, Lc, Ld;
719 FieldContainer<Scalar> warp(N,2);
720 FieldContainer<Scalar> blend(N);
721 FieldContainer<Scalar> denom(N);
724 La = L1; Lb = L2; Lc = L3; Ld = L4;
break;
726 La = L2; Lb = L1; Lc = L3; Ld = L4;
break;
728 La = L3; Lb = L1; Lc = L4; Ld = L2;
break;
730 La = L4; Lb = L1; Lc = L3; Ld = L2;
break;
734 warpShiftFace3D<Scalar,ArrayType>(order,alpha,La,Lb,Lc,Ld,warp);
736 for (
int k=0;k<N;k++) {
737 blend(k) = Lb(k) * Lc(k) * Ld(k);
740 for (
int k=0;k<N;k++) {
741 denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
744 for (
int k=0;k<N;k++) {
745 if (denom(k) > tol) {
746 blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
752 for (
int k=0;k<N;k++) {
753 for (
int j=0;j<3;j++) {
754 shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
755 + blend(k) * warp(k,1) * t2(face-1,j);
759 for (
int k=0;k<N;k++) {
760 if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
761 for (
int j=0;j<3;j++) {
762 shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
769 FieldContainer<Scalar> updatedPoints(N,3);
770 for (
int k=0;k<N;k++) {
771 for (
int j=0;j<3;j++) {
772 updatedPoints(k,j) = XYZ(k,j) + shift(k,j);
776 warVerts.resize( 1 , 4 , 3 );
779 FieldContainer<Scalar> refPts(N,3);
782 shards::getCellTopologyData<shards::Tetrahedron<4> >() ,
788 for (
int i=0;i<=order;i++) {
789 for (
int j=0;j<=order-i;j++) {
790 for (
int k=0;k<=order-i-j;k++) {
791 if ( (i >= offset) && (i <= order-offset) &&
792 (j >= offset) && (j <= order-i-offset) &&
793 (k >= offset) && (k <= order-i-j-offset) ) {
794 points(offcur,0) = refPts(noffcur,0);
795 points(offcur,1) = refPts(noffcur,1);
796 points(offcur,2) = refPts(noffcur,2);