Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/src/Shared/Intrepid_PointToolsDef.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00048 #ifdef _MSC_VER
00049 #include "winmath.h"
00050 #endif
00051 
00052 
00053 namespace Intrepid {
00054 
00055   template<class Scalar, class ArrayType>
00056   void PointTools::getLattice(ArrayType &pts ,
00057                               const shards::CellTopology& cellType ,
00058                               const int order ,
00059                               const int offset ,
00060                               const EPointType pointType ) 
00061   {
00062     switch (pointType) {
00063     case POINTTYPE_EQUISPACED:
00064       getEquispacedLattice<Scalar,ArrayType>( cellType , pts , order , offset );
00065       break;
00066     case POINTTYPE_WARPBLEND:
00067       getWarpBlendLattice<Scalar,ArrayType>( cellType , pts , order , offset );
00068       break;
00069     default:
00070       TEUCHOS_TEST_FOR_EXCEPTION( true ,
00071                           std::invalid_argument ,
00072                           "PointTools::getLattice: invalid EPointType" );
00073     }
00074     return;
00075   }
00076 
00077   template<class Scalar, class ArrayType>
00078   void PointTools::getGaussPoints( ArrayType &pts ,
00079                                    const int order )
00080   {
00081 
00082     Scalar *z = new Scalar[order+1];
00083     Scalar *w = new Scalar[order+1];
00084 
00085     IntrepidPolylib::zwgj( z , w , order + 1 , 0.0 , 0.0 );
00086     for (int i=0;i<order+1;i++) {
00087       pts(i,0) = z[i];
00088     }
00089 
00090     delete []z;
00091     delete []w;
00092   }
00093 
00094   
00095   template<class Scalar, class ArrayType>
00096   void PointTools::getEquispacedLattice(const shards::CellTopology& cellType ,
00097                                         ArrayType &points ,
00098                                         const int order ,
00099                                         const int offset )
00100 
00101   {
00102     switch (cellType.getKey()) {
00103     case shards::Tetrahedron<4>::key:
00104     case shards::Tetrahedron<8>::key:
00105     case shards::Tetrahedron<10>::key:
00106       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00107                           || points.dimension(1) != 3 ,
00108                           std::invalid_argument ,
00109                           ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
00110       getEquispacedLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
00111       break;
00112     case shards::Triangle<3>::key:
00113     case shards::Triangle<4>::key:
00114     case shards::Triangle<6>::key:
00115       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00116                           || points.dimension(1) != 2 ,
00117                           std::invalid_argument ,
00118                           ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
00119       getEquispacedLatticeTriangle<Scalar,ArrayType>( points , order , offset );
00120       break;
00121     case shards::Line<2>::key:
00122     case shards::Line<3>::key:
00123       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00124                           || points.dimension(1) != 1 ,
00125                           std::invalid_argument ,
00126                           ">>> ERROR(PointTools::getEquispacedLattice): points argument is ill-sized." );
00127       getEquispacedLatticeLine<Scalar,ArrayType>( points , order , offset );
00128       break;
00129     default:
00130       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
00131                           ">>> ERROR (Intrepid::PointTools::getEquispacedLattice): Illegal cell type" );
00132     }
00133     
00134   }
00135 
00136   template<class Scalar, class ArrayType>
00137   void PointTools::getWarpBlendLattice( const shards::CellTopology& cellType ,
00138                                         ArrayType &points ,
00139                                         const int order ,
00140                                         const int offset )
00141 
00142   {
00143     switch (cellType.getKey()) {
00144     case shards::Tetrahedron<4>::key:
00145     case shards::Tetrahedron<8>::key:
00146     case shards::Tetrahedron<10>::key:
00147       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00148                           || points.dimension(1) != 3 ,
00149                           std::invalid_argument ,
00150                           ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
00151       getWarpBlendLatticeTetrahedron<Scalar,ArrayType>( points , order , offset );
00152       break;
00153     case shards::Triangle<3>::key:
00154     case shards::Triangle<4>::key:
00155     case shards::Triangle<6>::key:
00156       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00157                           || points.dimension(1) != 2 ,
00158                           std::invalid_argument ,
00159                           ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
00160       getWarpBlendLatticeTriangle<Scalar,ArrayType>( points , order , offset );
00161       break;
00162     case shards::Line<2>::key:
00163     case shards::Line<3>::key:
00164       TEUCHOS_TEST_FOR_EXCEPTION( ( points.dimension(0) != getLatticeSize( cellType , order , offset ) ) 
00165                           || points.dimension(1) != 1 ,
00166                           std::invalid_argument ,
00167                           ">>> ERROR(PointTools::getWarpBlendLattice): points argument is ill-sized." );
00168       getWarpBlendLatticeLine<Scalar,ArrayType>( points , order , offset );
00169       break;
00170     default:
00171       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
00172                           ">>> ERROR (Intrepid::PointTools::getWarpBlendLattice): Illegal cell type" );
00173     }
00174     
00175   }
00176 
00177   template<class Scalar, class ArrayType>
00178   void PointTools::getEquispacedLatticeLine( ArrayType &points ,
00179                                              const int order ,
00180                                              const int offset )
00181   {
00182     TEUCHOS_TEST_FOR_EXCEPTION( order < 0 ,
00183                         std::invalid_argument ,
00184                         ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" );
00185     if (order == 0) {
00186       points(0,0) = 0.0;
00187     }
00188     else {
00189       const Scalar h = 2.0 / order;
00190       
00191       for (int i=offset;i<=order-offset;i++) {
00192         points(i-offset,0) = -1.0 + h * (Scalar) i;
00193       }
00194     }
00195 
00196     return;
00197   }
00198 
00199   template<class Scalar, class ArrayType>
00200   void PointTools::getEquispacedLatticeTriangle( ArrayType &points ,
00201                                                  const int order ,
00202                                                  const int offset )
00203   {
00204     TEUCHOS_TEST_FOR_EXCEPTION( order <= 0 ,
00205                         std::invalid_argument ,
00206                         ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeLine): order must be positive" );
00207 
00208     const Scalar h = 1.0 / order;
00209     int cur = 0;
00210 
00211     for (int i=offset;i<=order-offset;i++) {
00212       for (int j=offset;j<=order-i-offset;j++) {
00213         points(cur,0) = (Scalar)0.0 + (Scalar) j * h ;
00214         points(cur,1) = (Scalar)0.0 + (Scalar) i * h;
00215         cur++;
00216       }
00217     }
00218 
00219     return;
00220   }
00221 
00222   template<class Scalar, class ArrayType>
00223   void PointTools::getEquispacedLatticeTetrahedron( ArrayType &points ,
00224                                                     const int order ,
00225                                                     const int offset )
00226   {
00227     TEUCHOS_TEST_FOR_EXCEPTION( (order <= 0) ,
00228                         std::invalid_argument ,
00229                         ">>> ERROR (Intrepid::PointTools::getEquispacedLatticeTetrahedron): order must be positive" );
00230 
00231     const Scalar h = 1.0 / order;
00232     int cur = 0;
00233 
00234     for (int i=offset;i<=order-offset;i++) {
00235       for (int j=offset;j<=order-i-offset;j++) {
00236         for (int k=offset;k<=order-i-j-offset;k++) {
00237           points(cur,0) = (Scalar) k * h;
00238           points(cur,1) = (Scalar) j * h;
00239           points(cur,2) = (Scalar) i * h;
00240           cur++;
00241         }
00242       }
00243     }
00244 
00245     return;
00246   }
00247 
00248   template<class Scalar, class ArrayType>
00249   void PointTools::getWarpBlendLatticeLine( ArrayType &points ,
00250                                             const int order ,
00251                                             const int offset )
00252   {
00253     Scalar *z = new Scalar[order+1];
00254     Scalar *w = new Scalar[order+1];
00255     
00256     // order is order of polynomial degree.  The Gauss-Lobatto points are accurate
00257     // up to degree 2*i-1
00258     IntrepidPolylib::zwglj( z , w , order+1 , 0.0 , 0.0 );
00259 
00260     for (int i=offset;i<=order-offset;i++) {
00261       points(i-offset,0) = z[i];
00262     }
00263 
00264     delete []z;
00265     delete []w;
00266 
00267     return;
00268   }
00269 
00270   template<class Scalar, class ArrayType>
00271   void PointTools::warpFactor( const int order , 
00272                               const ArrayType &xnodes ,
00273                               const ArrayType &xout ,
00274                               ArrayType &warp)
00275   {
00276     TEUCHOS_TEST_FOR_EXCEPTION( ( warp.dimension(0) != xout.dimension(0) ) ,
00277                         std::invalid_argument ,
00278                         ">>> ERROR (PointTools::warpFactor): xout and warp must be same size." );
00279 
00280     warp.initialize();
00281 
00282     FieldContainer<Scalar> d( xout.dimension(0) );
00283     d.initialize();
00284 
00285     FieldContainer<Scalar> xeq( order + 1 ,1);
00286     PointTools::getEquispacedLatticeLine<Scalar,ArrayType>( xeq , order , 0 );
00287     xeq.resize( order + 1 );
00288 
00289     TEUCHOS_TEST_FOR_EXCEPTION( ( xeq.dimension(0) != xnodes.dimension(0) ) ,
00290                         std::invalid_argument ,
00291                         ">>> ERROR (PointTools::warpFactor): xeq and xnodes must be same size." );
00292     
00293     for (int i=0;i<=order;i++) {
00294 
00295       for (int k=0;k<d.dimension(0);k++) {
00296         d(k) = xnodes(i) - xeq(i);
00297       }
00298 
00299       for (int j=1;j<order;j++) {
00300         if (i != j) {
00301           for (int k=0;k<d.dimension(0);k++) {
00302             d(k) = d(k) * ( (xout(k)-xeq(j)) / (xeq(i)-xeq(j)) );
00303           }
00304         }
00305       }
00306       
00307       // deflate end roots
00308       if ( i != 0 ) {
00309         for (int k=0;k<d.dimension(0);k++) {
00310           d(k) = -d(k) / (xeq(i) - xeq(0));
00311         }
00312       }
00313 
00314       if (i != order ) {
00315         for (int k=0;k<d.dimension(0);k++) {
00316           d(k) = d(k) / (xeq(i) - xeq(order));
00317         }
00318       }
00319 
00320       for (int k=0;k<d.dimension(0);k++) {
00321         warp(k) += d(k);
00322       }
00323 
00324     }
00325 
00326 
00327     return;
00328   }
00329 
00330   template<class Scalar, class ArrayType>
00331   void PointTools::getWarpBlendLatticeTriangle( ArrayType &points ,
00332                                                 const int order ,
00333                                                 const int offset  )
00334   {
00335     /* get Gauss-Lobatto points */
00336     Intrepid::FieldContainer<Scalar> gaussX( order + 1 , 1 );
00337     
00338     PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
00339     
00340     gaussX.resize(gaussX.dimension(0));
00341 
00342     Scalar alpopt[] = {0.0000,0.0000,1.4152,0.1001,0.2751,0.9800,1.0999,
00343                         1.2832,1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
00344 
00345     Scalar alpha;
00346 
00347     if (order >= 1 && order < 16) {
00348       alpha = alpopt[order-1];
00349     }
00350     else {
00351       alpha = 5.0 / 3.0;
00352     }
00353 
00354     const int p = order; /* switch to Warburton's notation */
00355     int N = (p+1)*(p+2)/2;
00356     
00357     /* equidistributed nodes on equilateral triangle */
00358     Intrepid::FieldContainer<Scalar> L1( N );
00359     Intrepid::FieldContainer<Scalar> L2( N );
00360     Intrepid::FieldContainer<Scalar> L3( N );     
00361     Intrepid::FieldContainer<Scalar> X(N);
00362     Intrepid::FieldContainer<Scalar> Y(N);
00363     
00364     int sk = 0;
00365     for (int n=1;n<=p+1;n++) {
00366       for (int m=1;m<=p+2-n;m++) {
00367         L1(sk) = (n-1) / (Scalar)p;
00368         L3(sk) = (m-1) / (Scalar)p;
00369         L2(sk) = 1.0 - L1(sk) - L3(sk);
00370         sk++;
00371       }
00372     }
00373     
00374     for (int n=0;n<N;n++) {
00375       X(n) = -L2(n) + L3(n);
00376       Y(n) = (-L2(n) - L3(n) + 2*L1(n))/1.7320508075688772;
00377     }
00378     
00379     /* get blending function for each node at each edge */
00380     Intrepid::FieldContainer<Scalar> blend1(N);
00381     Intrepid::FieldContainer<Scalar> blend2(N);
00382     Intrepid::FieldContainer<Scalar> blend3(N);
00383     
00384     for (int n=0;n<N;n++) {
00385       blend1(n) = 4.0 * L2(n) * L3(n);
00386       blend2(n) = 4.0 * L1(n) * L3(n);
00387       blend3(n) = 4.0 * L1(n) * L2(n);
00388     }
00389     
00390     /* get difference of each barycentric coordinate */
00391     Intrepid::FieldContainer<Scalar> L3mL2(N);
00392     Intrepid::FieldContainer<Scalar> L1mL3(N);
00393     Intrepid::FieldContainer<Scalar> L2mL1(N);
00394 
00395     for (int k=0;k<N;k++) {
00396       L3mL2(k) = L3(k)-L2(k);
00397       L1mL3(k) = L1(k)-L3(k);
00398       L2mL1(k) = L2(k)-L1(k);
00399     }
00400 
00401     FieldContainer<Scalar> warpfactor1(N);
00402     FieldContainer<Scalar> warpfactor2(N);
00403     FieldContainer<Scalar> warpfactor3(N);
00404     
00405     warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L3mL2 , warpfactor1 );
00406     warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L1mL3 , warpfactor2 );
00407     warpFactor<Scalar,FieldContainer<Scalar> >( order , gaussX , L2mL1 , warpfactor3 );
00408 
00409     FieldContainer<Scalar> warp1(N);
00410     FieldContainer<Scalar> warp2(N);
00411     FieldContainer<Scalar> warp3(N);
00412 
00413     for (int k=0;k<N;k++) {
00414       warp1(k) = blend1(k) * warpfactor1(k) *
00415         ( 1.0 + alpha * alpha * L1(k) * L1(k) );
00416       warp2(k) = blend2(k) * warpfactor2(k) *
00417         ( 1.0 + alpha * alpha * L2(k) * L2(k) );
00418       warp3(k) = blend3(k) * warpfactor3(k) *
00419         ( 1.0 + alpha * alpha * L3(k) * L3(k) );
00420     }
00421 
00422     for (int k=0;k<N;k++) {
00423       X(k) += 1.0 * warp1(k) + cos( 2.0 * M_PI / 3.0 ) * warp2(k) + cos(4*M_PI/3.0) * warp3(k);
00424       Y(k) += 0.0 * warp1(k) + sin( 2.0 * M_PI / 3.0 ) * warp2(k) + sin( 4*M_PI/3.0) * warp3(k);
00425     }
00426 
00427     FieldContainer<Scalar> warXY(N,2);
00428     
00429     for (int k=0;k<N;k++) {
00430       warXY(k,0) = X(k);
00431       warXY(k,1) = Y(k);
00432     }
00433 
00434 
00435     // finally, convert the warp-blend points to the correct triangle
00436     FieldContainer<Scalar> warburtonVerts(1,3,2);
00437     warburtonVerts(0,0,0) = -1.0;
00438     warburtonVerts(0,0,1) = -1.0/sqrt(3.0);
00439     warburtonVerts(0,1,0) = 1.0;
00440     warburtonVerts(0,1,1) = -1.0/sqrt(3.0);
00441     warburtonVerts(0,2,0) = 0.0;
00442     warburtonVerts(0,2,1) = 2.0/sqrt(3.0);
00443 
00444     FieldContainer<Scalar> refPts(N,2);
00445 
00446     Intrepid::CellTools<Scalar>::mapToReferenceFrame( refPts ,
00447                                                       warXY ,
00448                                                       warburtonVerts ,
00449                                                       shards::getCellTopologyData< shards::Triangle<3> >(),
00450                                                       0 );
00451     
00452     // now write from refPts into points, taking care of offset
00453     int noffcur = 0;  // index into refPts
00454     int offcur = 0;   // index int points
00455     for (int i=0;i<=order;i++) {
00456       for (int j=0;j<=order-i;j++) {
00457         if ( (i >= offset) && (i <= order-offset) &&
00458               (j >= offset) && (j <= order-i-offset) ) {
00459           points(offcur,0) = refPts(noffcur,0);
00460           points(offcur,1) = refPts(noffcur,1);
00461           offcur++;
00462         }
00463         noffcur++;
00464       }
00465     }
00466 
00467     return;
00468   }
00469   
00470 
00471   template<class Scalar, class ArrayType>
00472   void PointTools::warpShiftFace3D( const int order ,
00473                                     const Scalar pval ,
00474                                     const ArrayType &L1,
00475                                     const ArrayType &L2,
00476                                     const ArrayType &L3,
00477                                     const ArrayType &L4,
00478                                     ArrayType &dxy)
00479   {
00480     evalshift<Scalar,ArrayType>(order,pval,L2,L3,L4,dxy);
00481     return;
00482   }
00483 
00484   template<class Scalar, class ArrayType>
00485   void PointTools::evalshift( const int order ,
00486                               const Scalar pval ,
00487                               const ArrayType &L1 ,
00488                               const ArrayType &L2 ,
00489                               const ArrayType &L3 ,
00490                               ArrayType &dxy )
00491   {
00492     // get Gauss-Lobatto-nodes
00493     FieldContainer<Scalar> gaussX(order+1,1);
00494     PointTools::getWarpBlendLatticeLine<Scalar,FieldContainer<Scalar> >( gaussX , order , 0 );
00495     gaussX.resize(order+1);
00496     const int N = L1.dimension(0);
00497     
00498     // Warburton code reverses them
00499     for (int k=0;k<=order;k++) {
00500       gaussX(k) *= -1.0;
00501     }
00502 
00503     // blending function at each node for each edge
00504     FieldContainer<Scalar> blend1(N);
00505     FieldContainer<Scalar> blend2(N);
00506     FieldContainer<Scalar> blend3(N);
00507 
00508     for (int i=0;i<N;i++) {
00509       blend1(i) = L2(i) * L3(i);
00510       blend2(i) = L1(i) * L3(i);
00511       blend3(i) = L1(i) * L2(i);
00512     }
00513 
00514     // amount of warp for each node for each edge
00515     FieldContainer<Scalar> warpfactor1(N);
00516     FieldContainer<Scalar> warpfactor2(N);
00517     FieldContainer<Scalar> warpfactor3(N);
00518 
00519     // differences of barycentric coordinates 
00520     FieldContainer<Scalar> L3mL2(N);
00521     FieldContainer<Scalar> L1mL3(N);
00522     FieldContainer<Scalar> L2mL1(N);
00523     
00524     for (int k=0;k<N;k++) {
00525       L3mL2(k) = L3(k)-L2(k);
00526       L1mL3(k) = L1(k)-L3(k);
00527       L2mL1(k) = L2(k)-L1(k);
00528     }
00529     
00530     evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor1 , order , gaussX , L3mL2 );
00531     evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor2 , order , gaussX , L1mL3 );
00532     evalwarp<Scalar,FieldContainer<Scalar> >( warpfactor3 , order , gaussX , L2mL1 );
00533     
00534     for (int k=0;k<N;k++) {
00535       warpfactor1(k) *= 4.0;
00536       warpfactor2(k) *= 4.0;
00537       warpfactor3(k) *= 4.0;      
00538     }
00539 
00540     FieldContainer<Scalar> warp1(N);
00541     FieldContainer<Scalar> warp2(N);
00542     FieldContainer<Scalar> warp3(N);
00543     
00544     for (int k=0;k<N;k++) {
00545       warp1(k) = blend1(k) * warpfactor1(k) *
00546         ( 1.0 + pval * pval * L1(k) * L1(k) );
00547       warp2(k) = blend2(k) * warpfactor2(k) *
00548         ( 1.0 + pval * pval * L2(k) * L2(k) );
00549       warp3(k) = blend3(k) * warpfactor3(k) *
00550         ( 1.0 + pval * pval * L3(k) * L3(k) );
00551     }
00552 
00553     for (int k=0;k<N;k++) {
00554       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);
00555       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);
00556     }
00557 
00558     return;
00559 
00560   }
00561 
00562   /* one-d edge warping function */
00563   template<class Scalar, class ArrayType>
00564   void PointTools::evalwarp( ArrayType &warp ,
00565                             const int order ,
00566                             const ArrayType &xnodes ,
00567                             const ArrayType &xout )
00568   {
00569     FieldContainer<Scalar> xeq(order+1);
00570     FieldContainer<Scalar> d(xout.dimension(0));
00571 
00572     d.initialize();
00573 
00574     for (int i=0;i<=order;i++) {
00575       xeq(i) = -1.0 + 2.0 * ( order - i ) / order;
00576     }
00577 
00578 
00579 
00580     for (int i=0;i<=order;i++) {
00581       d.initialize( xnodes(i) - xeq(i) );
00582       for (int j=1;j<order;j++) {
00583         if (i!=j) {
00584           for (int k=0;k<d.dimension(0);k++) {
00585             d(k) = d(k) * (xout(k)-xeq(j))/(xeq(i)-xeq(j));
00586           }
00587         }
00588       }
00589       if (i!=0) {
00590         for (int k=0;k<d.dimension(0);k++) {
00591           d(k) = -d(k)/(xeq(i)-xeq(0));
00592         }
00593       }
00594       if (i!=order) {
00595         for (int k=0;k<d.dimension(0);k++) {
00596           d(k) = d(k)/(xeq(i)-xeq(order));
00597         }
00598       }
00599       
00600       for (int k=0;k<d.dimension(0);k++) {
00601         warp(k) += d(k);
00602       } 
00603     }    
00604 
00605     return;
00606   }
00607 
00608 
00609   template<class Scalar, class ArrayType>
00610   void PointTools::getWarpBlendLatticeTetrahedron(ArrayType &points ,
00611                                                   const int order ,
00612                                                   const int offset  )
00613   {
00614     Scalar alphastore[] = { 0,0,0,0.1002, 1.1332,1.5608,1.3413,1.2577,1.1603,
00615                             1.10153,0.6080,0.4523,0.8856,0.8717,0.9655};
00616     Scalar alpha;
00617 
00618     if (order <= 15) {
00619       alpha = alphastore[order-1]; 
00620     }
00621     else {
00622       alpha = 1.0;
00623     }
00624 
00625     const int N = (order+1)*(order+2)*(order+3)/6;
00626     Scalar tol = 1.e-10;
00627 
00628     FieldContainer<Scalar> shift(N,3);
00629     shift.initialize();
00630 
00631     /* create 3d equidistributed nodes on Warburton tet */
00632     FieldContainer<Scalar> equipoints(N,3);
00633     int sk = 0;
00634     for (int n=0;n<=order;n++) {
00635       for (int m=0;m<=order-n;m++) {
00636         for (int q=0;q<=order-n-m;q++) {
00637           equipoints(sk,0) = -1.0 + (q * 2.0 ) / order;
00638           equipoints(sk,1) = -1.0 + (m * 2.0 ) / order;
00639           equipoints(sk,2) = -1.0 + (n * 2.0 ) / order;
00640           sk++;
00641         }
00642       }
00643     }
00644     
00645 
00646     /* construct barycentric coordinates for equispaced lattice */
00647     FieldContainer<Scalar> L1(N);
00648     FieldContainer<Scalar> L2(N);
00649     FieldContainer<Scalar> L3(N);
00650     FieldContainer<Scalar> L4(N);
00651     for (int i=0;i<N;i++) {
00652       L1(i) = (1.0 + equipoints(i,2)) / 2.0;
00653       L2(i) = (1.0 + equipoints(i,1)) / 2.0;
00654       L3(i) = -(1.0 + equipoints(i,0) + equipoints(i,1) + equipoints(i,2)) / 2.0;
00655       L4(i) = (1.0 + equipoints(i,0)) / 2.0;
00656     }
00657     
00658     
00659     /* vertices of Warburton tet */
00660     FieldContainer<Scalar> warVerts(4,3);
00661     warVerts(0,0) = -1.0;
00662     warVerts(0,1) = -1.0/sqrt(3.0);
00663     warVerts(0,2) = -1.0/sqrt(6.0);
00664     warVerts(1,0) = 1.0;
00665     warVerts(1,1) = -1.0/sqrt(3.0);
00666     warVerts(1,2) = -1.0/sqrt(6.0);
00667     warVerts(2,0) = 0.0;
00668     warVerts(2,1) = 2.0 / sqrt(3.0);
00669     warVerts(2,2) = -1.0/sqrt(6.0);
00670     warVerts(3,0) = 0.0;
00671     warVerts(3,1) = 0.0;
00672     warVerts(3,2) = 3.0 / sqrt(6.0);
00673 
00674 
00675     /* tangents to faces */
00676     FieldContainer<Scalar> t1(4,3);
00677     FieldContainer<Scalar> t2(4,3);
00678     for (int i=0;i<3;i++) {
00679       t1(0,i) = warVerts(1,i) - warVerts(0,i);
00680       t1(1,i) = warVerts(1,i) - warVerts(0,i);
00681       t1(2,i) = warVerts(2,i) - warVerts(1,i);
00682       t1(3,i) = warVerts(2,i) - warVerts(0,i);
00683       t2(0,i) = warVerts(2,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
00684       t2(1,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(1,i) );
00685       t2(2,i) = warVerts(3,i) - 0.5 * ( warVerts(1,i) + warVerts(2,i) );
00686       t2(3,i) = warVerts(3,i) - 0.5 * ( warVerts(0,i) + warVerts(2,i) );
00687     }
00688 
00689     /* normalize tangents */
00690     for (int n=0;n<4;n++) {
00691       /* Compute norm of t1(n) and t2(n) */
00692       Scalar normt1n = 0.0;
00693       Scalar normt2n = 0.0;
00694       for (int i=0;i<3;i++) {
00695         normt1n += (t1(n,i) * t1(n,i));
00696         normt2n += (t2(n,i) * t2(n,i));
00697       }
00698       normt1n = sqrt(normt1n);
00699       normt2n = sqrt(normt2n);
00700       /* normalize each tangent now */
00701       for (int i=0;i<3;i++) {
00702         t1(n,i) /= normt1n;
00703         t2(n,i) /= normt2n;
00704       }
00705     }
00706 
00707     /* undeformed coordinates */
00708     FieldContainer<Scalar> XYZ(N,3);
00709     for (int i=0;i<N;i++) {
00710       for (int j=0;j<3;j++) {
00711         XYZ(i,j) = L3(i)*warVerts(0,j) + L4(i)*warVerts(1,j) + L2(i)*warVerts(2,j) + L1(i)*warVerts(3,j);
00712       }
00713     }
00714 
00715     for (int face=1;face<=4;face++) {
00716       FieldContainer<Scalar> La, Lb, Lc, Ld;
00717       FieldContainer<Scalar> warp(N,2);
00718       FieldContainer<Scalar> blend(N);
00719       FieldContainer<Scalar> denom(N);
00720       switch (face) {
00721       case 1:
00722         La = L1; Lb = L2; Lc = L3; Ld = L4; break;
00723       case 2:
00724         La = L2; Lb = L1; Lc = L3; Ld = L4; break;
00725       case 3:
00726         La = L3; Lb = L1; Lc = L4; Ld = L2; break;
00727       case 4:
00728         La = L4; Lb = L1; Lc = L3; Ld = L2; break;
00729       }
00730       
00731       /* get warp tangential to face */
00732       warpShiftFace3D<Scalar,ArrayType>(order,alpha,La,Lb,Lc,Ld,warp);
00733       
00734       for (int k=0;k<N;k++) {
00735         blend(k) = Lb(k) * Lc(k) * Ld(k);
00736       }
00737 
00738       for (int k=0;k<N;k++) {
00739         denom(k) = (Lb(k) + 0.5 * La(k)) * (Lc(k) + 0.5*La(k)) * (Ld(k) + 0.5 * La(k));
00740       }
00741 
00742       for (int k=0;k<N;k++) {
00743         if (denom(k) > tol) {
00744           blend(k) *= ( 1.0 + alpha * alpha * La(k) * La(k) ) / denom(k);
00745         }
00746       }  
00747 
00748 
00749       // compute warp and blend
00750       for (int k=0;k<N;k++) {
00751         for (int j=0;j<3;j++) {
00752           shift(k,j) = shift(k,j) + blend(k) * warp(k,0) * t1(face-1,j)
00753             + blend(k) * warp(k,1) * t2(face-1,j);
00754         }
00755       }
00756 
00757       for (int k=0;k<N;k++) {
00758         if (La(k) < tol && ( Lb(k) < tol || Lc(k) < tol || Ld(k) < tol )) {
00759           for (int j=0;j<3;j++) {
00760             shift(k,j) = warp(k,0) * t1(face-1,j) + warp(k,1) * t2(face-1,j);
00761           }
00762         }
00763       }
00764       
00765     }
00766 
00767     FieldContainer<Scalar> updatedPoints(N,3);
00768     for (int k=0;k<N;k++) {
00769       for (int j=0;j<3;j++) {
00770         updatedPoints(k,j) = XYZ(k,j) + shift(k,j);
00771       }
00772     }
00773 
00774     warVerts.resize( 1 , 4 , 3 );
00775 
00776     // now we convert to Pavel's reference triangle!
00777     FieldContainer<Scalar> refPts(N,3);
00778     CellTools<Scalar>::mapToReferenceFrame( refPts ,updatedPoints ,
00779                                             warVerts ,
00780                                             shards::getCellTopologyData<shards::Tetrahedron<4> >() ,
00781                                             0 );
00782 
00783     // now write from refPts into points, taking offset into account
00784     int noffcur = 0;
00785     int offcur = 0;
00786     for (int i=0;i<=order;i++) {
00787       for (int j=0;j<=order-i;j++) {
00788         for (int k=0;k<=order-i-j;k++) {
00789           if ( (i >= offset) && (i <= order-offset) &&
00790               (j >= offset) && (j <= order-i-offset) &&
00791               (k >= offset) && (k <= order-i-j-offset) ) {
00792             points(offcur,0) = refPts(noffcur,0);
00793             points(offcur,1) = refPts(noffcur,1);
00794             points(offcur,2) = refPts(noffcur,2);
00795             offcur++;
00796           }
00797           noffcur++;
00798         }
00799       }
00800     }
00801                                             
00802 
00803 
00804   }
00805   
00806 
00807 } // face Intrepid