|
Intrepid
|
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;
1.7.6.1