|
Intrepid
|
00001 namespace Intrepid { 00002 00003 template<class Scalar> 00004 void OrthogonalBases::jrc(const Scalar &alpha , const Scalar &beta , 00005 const int &n , 00006 Scalar &an , Scalar &bn, Scalar &cn ) 00007 { 00008 an = (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta ) 00009 / ( 2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) ); 00010 bn = (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta) 00011 / ( 2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) ); 00012 cn = (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta) 00013 / ( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) ); 00014 00015 return; 00016 } 00017 00018 00019 template<class Scalar, class ScalarArray1, class ScalarArray2> 00020 void OrthogonalBases::tabulateTriangle( const ScalarArray1& z , 00021 const int n , 00022 ScalarArray2 & poly_val ) 00023 { 00024 const int np = z.dimension( 0 ); 00025 00026 // each point needs to be transformed from Pavel's element 00027 // z(i,0) --> (2.0 * z(i,0) - 1.0) 00028 // z(i,1) --> (2.0 * z(i,1) - 1.0) 00029 00030 // set up constant term 00031 int idx_cur = OrthogonalBases::idxtri(0,0); 00032 int idx_curp1,idx_curm1; 00033 00034 // set D^{0,0} = 1.0 00035 for (int i=0;i<np;i++) { 00036 poly_val(idx_cur,i) = 1.0; 00037 } 00038 00039 Teuchos::Array<Scalar> f1(np),f2(np),f3(np); 00040 00041 for (int i=0;i<np;i++) { 00042 f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0)); 00043 f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0)); 00044 f3[i] = f2[i] * f2[i]; 00045 } 00046 00047 // set D^{1,0} = f1 00048 idx_cur = OrthogonalBases::idxtri(1,0); 00049 for (int i=0;i<np;i++) { 00050 poly_val(idx_cur,i) = f1[i]; 00051 } 00052 00053 // recurrence in p 00054 for (int p=1;p<n;p++) { 00055 idx_cur = OrthogonalBases::idxtri(p,0); 00056 idx_curp1 = OrthogonalBases::idxtri(p+1,0); 00057 idx_curm1 = OrthogonalBases::idxtri(p-1,0); 00058 Scalar a = (2.0*p+1.0)/(1.0+p); 00059 Scalar b = p / (p+1.0); 00060 00061 for (int i=0;i<np;i++) { 00062 poly_val(idx_curp1,i) = a * f1[i] * poly_val(idx_cur,i) 00063 - b * f3[i] * poly_val(idx_curm1,i); 00064 } 00065 } 00066 00067 // D^{p,1} 00068 for (int p=0;p<n;p++) { 00069 int idxp0 = OrthogonalBases::idxtri(p,0); 00070 int idxp1 = OrthogonalBases::idxtri(p,1); 00071 for (int i=0;i<np;i++) { 00072 poly_val(idxp1,i) = poly_val(idxp0,i) 00073 *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0)); 00074 } 00075 } 00076 00077 // recurrence in q 00078 for (int p=0;p<n-1;p++) { 00079 for (int q=1;q<n-p;q++) { 00080 int idxpqp1=OrthogonalBases::idxtri(p,q+1); 00081 int idxpq=OrthogonalBases::idxtri(p,q); 00082 int idxpqm1=OrthogonalBases::idxtri(p,q-1); 00083 Scalar a,b,c; 00084 jrc((Scalar)(2*p+1),(Scalar)0,q,a,b,c); 00085 for (int i=0;i<np;i++) { 00086 poly_val(idxpqp1,i) 00087 = (a*(2.0*z(i,1)-1.0)+b)*poly_val(idxpq,i) 00088 - c*poly_val(idxpqm1,i); 00089 } 00090 } 00091 } 00092 00093 return; 00094 } 00095 00096 template<class Scalar, class ScalarArray1, class ScalarArray2> 00097 void OrthogonalBases::tabulateTetrahedron(const ScalarArray1 &z , 00098 const int n , 00099 ScalarArray2 &poly_val ) 00100 { 00101 const int np = z.dimension( 0 ); 00102 int idxcur; 00103 00104 // each point needs to be transformed from Pavel's element 00105 // z(i,0) --> (2.0 * z(i,0) - 1.0) 00106 // z(i,1) --> (2.0 * z(i,1) - 1.0) 00107 // z(i,2) --> (2.0 * z(i,2) - 1.0) 00108 00109 Teuchos::Array<Scalar> f1(np),f2(np),f3(np),f4(np),f5(np); 00110 00111 for (int i=0;i<np;i++) { 00112 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) ); 00113 f2[i] = pow( 0.5 * ( (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) , 2 ); 00114 f3[i] = 0.5 * ( 1.0 + 2.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ); 00115 f4[i] = 0.5 * ( 1.0 - (2.0*z(i,2)-1.0) ); 00116 f5[i] = f4[i] * f4[i]; 00117 } 00118 00119 // constant term 00120 idxcur = idxtet(0,0,0); 00121 for (int i=0;i<np;i++) { 00122 poly_val(idxcur,i) = 1.0; 00123 } 00124 00125 // D^{1,0,0} 00126 idxcur = idxtet(1,0,0); 00127 for (int i=0;i<np;i++) { 00128 poly_val(idxcur,i) = f1[i]; 00129 } 00130 00131 // p recurrence 00132 for (int p=1;p<n;p++) { 00133 Scalar a1 = (2.0 * p + 1.0) / ( p + 1.0); 00134 Scalar a2 = p / ( p + 1.0 ); 00135 int idxp = idxtet(p,0,0); 00136 int idxpp1 = idxtet(p+1,0,0); 00137 int idxpm1 = idxtet(p-1,0,0); 00138 //cout << idxpm1 << " " << idxp << " " << idxpp1 << endl; 00139 for (int i=0;i<np;i++) { 00140 poly_val(idxpp1,i) = a1 * f1[i] * poly_val(idxp,i) - a2 * f2[i] * poly_val(idxpm1,i); 00141 } 00142 } 00143 // q = 1 00144 for (int p=0;p<n;p++) { 00145 int idx0 = idxtet(p,0,0); 00146 int idx1 = idxtet(p,1,0); 00147 for (int i=0;i<np;i++) { 00148 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) ) ); 00149 } 00150 } 00151 00152 // q recurrence 00153 for (int p=0;p<n-1;p++) { 00154 for (int q=1;q<n-p;q++) { 00155 Scalar aq,bq,cq; 00156 jrc((Scalar)(2.0*p+1.0),(Scalar)(0),q,aq,bq,cq); 00157 int idxpqp1 = idxtet(p,q+1,0); 00158 int idxpq = idxtet(p,q,0); 00159 int idxpqm1 = idxtet(p,q-1,0); 00160 for (int i=0;i<np;i++) { 00161 poly_val(idxpqp1,i) = ( aq * f3[i] + bq * f4[i] ) * poly_val(idxpq,i) 00162 - ( cq * f5[i] ) * poly_val(idxpqm1,i); 00163 } 00164 } 00165 } 00166 00167 // r = 1 00168 for (int p=0;p<n;p++) { 00169 for (int q=0;q<n-p;q++) { 00170 int idxpq1 = idxtet(p,q,1); 00171 int idxpq0 = idxtet(p,q,0); 00172 for (int i=0;i<np;i++) { 00173 poly_val(idxpq1,i) = poly_val(idxpq0,i) * ( 1.0 + p + q + ( 2.0 + q + p ) * (2.0*z(i,2)-1.0) ); 00174 } 00175 } 00176 } 00177 00178 // general r recurrence 00179 for (int p=0;p<n-1;p++) { 00180 for (int q=0;q<n-p-1;q++) { 00181 for (int r=1;r<n-p-q;r++) { 00182 Scalar ar,br,cr; 00183 int idxpqrp1 = idxtet(p,q,r+1); 00184 int idxpqr = idxtet(p,q,r); 00185 int idxpqrm1 = idxtet(p,q,r-1); 00186 jrc(2.0*p+2.0*q+2.0,0.0,r,ar,br,cr); 00187 for (int i=0;i<np;i++) { 00188 poly_val(idxpqrp1,i) = (ar * (2.0*z(i,2)-1.0) + br) * poly_val( idxpqr , i ) - cr * poly_val(idxpqrm1,i); 00189 } 00190 } 00191 } 00192 } 00193 00194 return; 00195 00196 } 00197 00198 } // namespace Intrepid;
1.7.4