Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/Basis/Intrepid_OrthogonalBasesDef.hpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //
00005 //                           Intrepid Package
00006 //                 Copyright (2007) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00039 //                    Denis Ridzal  (dridzal@sandia.gov), or
00040 //                    Kara Peterson (kjpeter@sandia.gov)
00041 //
00042 // ************************************************************************
00043 // @HEADER
00044 */
00045 
00046 namespace Intrepid {
00047   
00048   template<class Scalar>
00049   void OrthogonalBases::jrc(const Scalar &alpha , const Scalar &beta , 
00050                             const int &n ,
00051                             Scalar &an , Scalar &bn, Scalar &cn )
00052   {
00053     an = (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta ) 
00054       / ( 2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) );
00055     bn = (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta) 
00056       / ( 2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) );
00057     cn = (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta) 
00058       / ( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) );
00059     
00060     return;
00061   }
00062   
00063   
00064   template<class Scalar, class ScalarArray1, class ScalarArray2>
00065   void OrthogonalBases::tabulateTriangle( const ScalarArray1& z ,
00066                                           const int n ,
00067                                           ScalarArray2 & poly_val )
00068   {
00069     const int np = z.dimension( 0 );
00070 
00071     // each point needs to be transformed from Pavel's element
00072     // z(i,0) --> (2.0 * z(i,0) - 1.0)
00073     // z(i,1) --> (2.0 * z(i,1) - 1.0)
00074 
00075     // set up constant term
00076     int idx_cur = OrthogonalBases::idxtri(0,0);
00077     int idx_curp1,idx_curm1;
00078 
00079     // set D^{0,0} = 1.0
00080     for (int i=0;i<np;i++) {
00081       poly_val(idx_cur,i) = 1.0;
00082     }
00083 
00084     Teuchos::Array<Scalar> f1(np),f2(np),f3(np);
00085     
00086     for (int i=0;i<np;i++) {
00087       f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0));
00088       f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0));
00089       f3[i] = f2[i] * f2[i];
00090     }
00091 
00092     // set D^{1,0} = f1
00093     idx_cur = OrthogonalBases::idxtri(1,0);
00094     for (int i=0;i<np;i++) {
00095       poly_val(idx_cur,i) = f1[i];
00096     }
00097 
00098     // recurrence in p
00099     for (int p=1;p<n;p++) {
00100       idx_cur = OrthogonalBases::idxtri(p,0);
00101       idx_curp1 = OrthogonalBases::idxtri(p+1,0);
00102       idx_curm1 = OrthogonalBases::idxtri(p-1,0);
00103       Scalar a = (2.0*p+1.0)/(1.0+p);
00104       Scalar b = p / (p+1.0);
00105 
00106       for (int i=0;i<np;i++) {
00107         poly_val(idx_curp1,i) = a * f1[i] * poly_val(idx_cur,i)
00108           - b * f3[i] * poly_val(idx_curm1,i);
00109       }
00110     }
00111     
00112     // D^{p,1}
00113     for (int p=0;p<n;p++) {
00114       int idxp0 = OrthogonalBases::idxtri(p,0);
00115       int idxp1 = OrthogonalBases::idxtri(p,1);
00116       for (int i=0;i<np;i++) {
00117         poly_val(idxp1,i) = poly_val(idxp0,i)
00118           *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
00119       }
00120     }
00121 
00122     // recurrence in q
00123     for (int p=0;p<n-1;p++) {
00124       for (int q=1;q<n-p;q++) {
00125         int idxpqp1=OrthogonalBases::idxtri(p,q+1);
00126         int idxpq=OrthogonalBases::idxtri(p,q);
00127         int idxpqm1=OrthogonalBases::idxtri(p,q-1);
00128         Scalar a,b,c;
00129         jrc((Scalar)(2*p+1),(Scalar)0,q,a,b,c);
00130         for (int i=0;i<np;i++) {
00131           poly_val(idxpqp1,i)
00132             = (a*(2.0*z(i,1)-1.0)+b)*poly_val(idxpq,i)
00133             - c*poly_val(idxpqm1,i);
00134         }
00135       }
00136     }
00137     
00138     return;
00139   }
00140 
00141   template<class Scalar, class ScalarArray1, class ScalarArray2>
00142   void OrthogonalBases::tabulateTetrahedron(const ScalarArray1 &z , 
00143                                             const int n ,
00144                                             ScalarArray2 &poly_val )
00145   {
00146     const int np = z.dimension( 0 );
00147     int idxcur;
00148 
00149     // each point needs to be transformed from Pavel's element
00150     // z(i,0) --> (2.0 * z(i,0) - 1.0)
00151     // z(i,1) --> (2.0 * z(i,1) - 1.0)
00152     // z(i,2) --> (2.0 * z(i,2) - 1.0)
00153     
00154     Teuchos::Array<Scalar> f1(np),f2(np),f3(np),f4(np),f5(np);
00155 
00156     for (int i=0;i<np;i++) {
00157       f1[i] = 0.5 * ( 2.0 + 2.0*(2.0*z(i,0)-1.0) + (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
00158       f2[i] = pow( 0.5 * ( (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) , 2 );
00159       f3[i] = 0.5 * ( 1.0 + 2.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
00160       f4[i] = 0.5 * ( 1.0 - (2.0*z(i,2)-1.0) );
00161       f5[i] = f4[i] * f4[i];
00162     }
00163     
00164     // constant term
00165     idxcur = idxtet(0,0,0);
00166     for (int i=0;i<np;i++) {
00167       poly_val(idxcur,i) = 1.0;
00168     }
00169 
00170     // D^{1,0,0}
00171     idxcur = idxtet(1,0,0);
00172     for (int i=0;i<np;i++) {
00173       poly_val(idxcur,i) = f1[i];
00174     }
00175 
00176     // p recurrence
00177     for (int p=1;p<n;p++) {
00178       Scalar a1 = (2.0 * p + 1.0) / ( p + 1.0);
00179       Scalar a2 = p / ( p + 1.0 );
00180       int idxp = idxtet(p,0,0);
00181       int idxpp1 = idxtet(p+1,0,0);
00182       int idxpm1 = idxtet(p-1,0,0);
00183       //cout << idxpm1 << " " << idxp << " " << idxpp1 << endl;
00184       for (int i=0;i<np;i++) {
00185         poly_val(idxpp1,i) = a1 * f1[i] * poly_val(idxp,i) - a2 * f2[i] * poly_val(idxpm1,i);
00186       }
00187     }
00188     // q = 1
00189     for (int p=0;p<n;p++) {
00190       int idx0 = idxtet(p,0,0);
00191       int idx1 = idxtet(p,1,0);
00192       for (int i=0;i<np;i++) {
00193         poly_val(idx1,i) = poly_val(idx0,i) * ( p * ( 1.0 + (2.0*z(i,1)-1.0) ) + 0.5 * ( 2.0 + 3.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) );
00194       }
00195     }
00196 
00197     // q recurrence
00198     for (int p=0;p<n-1;p++) {
00199       for (int q=1;q<n-p;q++) {
00200         Scalar aq,bq,cq;
00201         jrc((Scalar)(2.0*p+1.0),(Scalar)(0),q,aq,bq,cq);
00202         int idxpqp1 = idxtet(p,q+1,0);
00203         int idxpq = idxtet(p,q,0);
00204         int idxpqm1 = idxtet(p,q-1,0);
00205         for (int i=0;i<np;i++) {
00206           poly_val(idxpqp1,i) = ( aq * f3[i] + bq * f4[i] ) * poly_val(idxpq,i) 
00207             - ( cq * f5[i] ) * poly_val(idxpqm1,i);
00208         }
00209       }
00210     }
00211     
00212     // r = 1
00213     for (int p=0;p<n;p++) {
00214       for (int q=0;q<n-p;q++) {
00215         int idxpq1 = idxtet(p,q,1);
00216         int idxpq0 = idxtet(p,q,0);
00217         for (int i=0;i<np;i++) {
00218           poly_val(idxpq1,i) = poly_val(idxpq0,i) * ( 1.0 + p + q + ( 2.0 + q + p ) * (2.0*z(i,2)-1.0) );
00219         }
00220       }
00221     }
00222     
00223     // general r recurrence
00224     for (int p=0;p<n-1;p++) {
00225       for (int q=0;q<n-p-1;q++) {
00226         for (int r=1;r<n-p-q;r++) {
00227           Scalar ar,br,cr;
00228           int idxpqrp1 = idxtet(p,q,r+1);
00229           int idxpqr = idxtet(p,q,r);
00230           int idxpqrm1 = idxtet(p,q,r-1);
00231           jrc(2.0*p+2.0*q+2.0,0.0,r,ar,br,cr);
00232           for (int i=0;i<np;i++) {
00233             poly_val(idxpqrp1,i) = (ar * (2.0*z(i,2)-1.0) + br) * poly_val( idxpqr , i ) - cr * poly_val(idxpqrm1,i);
00234           }
00235         }
00236       }
00237     }
00238     
00239     return;
00240     
00241   }
00242 
00243 } // namespace Intrepid;