Intrepid
Intrepid_OrthogonalBasesDef.hpp
1 /*
2 // @HEADER
3 // ************************************************************************
4 //
5 // Intrepid Package
6 // Copyright (2007) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
39 // Denis Ridzal (dridzal@sandia.gov), or
40 // Kara Peterson (kjpeter@sandia.gov)
41 //
42 // ************************************************************************
43 // @HEADER
44 */
45 
46 namespace Intrepid {
47 
48  template<class Scalar>
49  void OrthogonalBases::jrc(const Scalar &alpha , const Scalar &beta ,
50  const int &n ,
51  Scalar &an , Scalar &bn, Scalar &cn )
52  {
53  an = (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta )
54  / ( 2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) );
55  bn = (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta)
56  / ( 2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) );
57  cn = (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta)
58  / ( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) );
59 
60  return;
61  }
62 
63 
64  template<class Scalar, class ScalarArray1, class ScalarArray2>
65  void OrthogonalBases::tabulateTriangle( const ScalarArray1& z ,
66  const int n ,
67  ScalarArray2 & poly_val )
68  {
69  const int np = z.dimension( 0 );
70 
71  // each point needs to be transformed from Pavel's element
72  // z(i,0) --> (2.0 * z(i,0) - 1.0)
73  // z(i,1) --> (2.0 * z(i,1) - 1.0)
74 
75  // set up constant term
76  int idx_cur = OrthogonalBases::idxtri(0,0);
77  int idx_curp1,idx_curm1;
78 
79  // set D^{0,0} = 1.0
80  for (int i=0;i<np;i++) {
81  poly_val(idx_cur,i) = 1.0;
82  }
83 
84  Teuchos::Array<Scalar> f1(np),f2(np),f3(np);
85 
86  for (int i=0;i<np;i++) {
87  f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0));
88  f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0));
89  f3[i] = f2[i] * f2[i];
90  }
91 
92  // set D^{1,0} = f1
93  idx_cur = OrthogonalBases::idxtri(1,0);
94  for (int i=0;i<np;i++) {
95  poly_val(idx_cur,i) = f1[i];
96  }
97 
98  // recurrence in p
99  for (int p=1;p<n;p++) {
100  idx_cur = OrthogonalBases::idxtri(p,0);
101  idx_curp1 = OrthogonalBases::idxtri(p+1,0);
102  idx_curm1 = OrthogonalBases::idxtri(p-1,0);
103  Scalar a = (2.0*p+1.0)/(1.0+p);
104  Scalar b = p / (p+1.0);
105 
106  for (int i=0;i<np;i++) {
107  poly_val(idx_curp1,i) = a * f1[i] * poly_val(idx_cur,i)
108  - b * f3[i] * poly_val(idx_curm1,i);
109  }
110  }
111 
112  // D^{p,1}
113  for (int p=0;p<n;p++) {
114  int idxp0 = OrthogonalBases::idxtri(p,0);
115  int idxp1 = OrthogonalBases::idxtri(p,1);
116  for (int i=0;i<np;i++) {
117  poly_val(idxp1,i) = poly_val(idxp0,i)
118  *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
119  }
120  }
121 
122  // recurrence in q
123  for (int p=0;p<n-1;p++) {
124  for (int q=1;q<n-p;q++) {
125  int idxpqp1=OrthogonalBases::idxtri(p,q+1);
126  int idxpq=OrthogonalBases::idxtri(p,q);
127  int idxpqm1=OrthogonalBases::idxtri(p,q-1);
128  Scalar a,b,c;
129  jrc((Scalar)(2*p+1),(Scalar)0,q,a,b,c);
130  for (int i=0;i<np;i++) {
131  poly_val(idxpqp1,i)
132  = (a*(2.0*z(i,1)-1.0)+b)*poly_val(idxpq,i)
133  - c*poly_val(idxpqm1,i);
134  }
135  }
136  }
137 
138  return;
139  }
140 
141  template<class Scalar, class ScalarArray1, class ScalarArray2>
142  void OrthogonalBases::tabulateTetrahedron(const ScalarArray1 &z ,
143  const int n ,
144  ScalarArray2 &poly_val )
145  {
146  const int np = z.dimension( 0 );
147  int idxcur;
148 
149  // each point needs to be transformed from Pavel's element
150  // z(i,0) --> (2.0 * z(i,0) - 1.0)
151  // z(i,1) --> (2.0 * z(i,1) - 1.0)
152  // z(i,2) --> (2.0 * z(i,2) - 1.0)
153 
154  Teuchos::Array<Scalar> f1(np),f2(np),f3(np),f4(np),f5(np);
155 
156  for (int i=0;i<np;i++) {
157  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) );
158  f2[i] = pow( 0.5 * ( (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) ) , 2 );
159  f3[i] = 0.5 * ( 1.0 + 2.0 * (2.0*z(i,1)-1.0) + (2.0*z(i,2)-1.0) );
160  f4[i] = 0.5 * ( 1.0 - (2.0*z(i,2)-1.0) );
161  f5[i] = f4[i] * f4[i];
162  }
163 
164  // constant term
165  idxcur = idxtet(0,0,0);
166  for (int i=0;i<np;i++) {
167  poly_val(idxcur,i) = 1.0;
168  }
169 
170  // D^{1,0,0}
171  idxcur = idxtet(1,0,0);
172  for (int i=0;i<np;i++) {
173  poly_val(idxcur,i) = f1[i];
174  }
175 
176  // p recurrence
177  for (int p=1;p<n;p++) {
178  Scalar a1 = (2.0 * p + 1.0) / ( p + 1.0);
179  Scalar a2 = p / ( p + 1.0 );
180  int idxp = idxtet(p,0,0);
181  int idxpp1 = idxtet(p+1,0,0);
182  int idxpm1 = idxtet(p-1,0,0);
183  //cout << idxpm1 << " " << idxp << " " << idxpp1 << endl;
184  for (int i=0;i<np;i++) {
185  poly_val(idxpp1,i) = a1 * f1[i] * poly_val(idxp,i) - a2 * f2[i] * poly_val(idxpm1,i);
186  }
187  }
188  // q = 1
189  for (int p=0;p<n;p++) {
190  int idx0 = idxtet(p,0,0);
191  int idx1 = idxtet(p,1,0);
192  for (int i=0;i<np;i++) {
193  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) ) );
194  }
195  }
196 
197  // q recurrence
198  for (int p=0;p<n-1;p++) {
199  for (int q=1;q<n-p;q++) {
200  Scalar aq,bq,cq;
201  jrc((Scalar)(2.0*p+1.0),(Scalar)(0),q,aq,bq,cq);
202  int idxpqp1 = idxtet(p,q+1,0);
203  int idxpq = idxtet(p,q,0);
204  int idxpqm1 = idxtet(p,q-1,0);
205  for (int i=0;i<np;i++) {
206  poly_val(idxpqp1,i) = ( aq * f3[i] + bq * f4[i] ) * poly_val(idxpq,i)
207  - ( cq * f5[i] ) * poly_val(idxpqm1,i);
208  }
209  }
210  }
211 
212  // r = 1
213  for (int p=0;p<n;p++) {
214  for (int q=0;q<n-p;q++) {
215  int idxpq1 = idxtet(p,q,1);
216  int idxpq0 = idxtet(p,q,0);
217  for (int i=0;i<np;i++) {
218  poly_val(idxpq1,i) = poly_val(idxpq0,i) * ( 1.0 + p + q + ( 2.0 + q + p ) * (2.0*z(i,2)-1.0) );
219  }
220  }
221  }
222 
223  // general r recurrence
224  for (int p=0;p<n-1;p++) {
225  for (int q=0;q<n-p-1;q++) {
226  for (int r=1;r<n-p-q;r++) {
227  Scalar ar,br,cr;
228  int idxpqrp1 = idxtet(p,q,r+1);
229  int idxpqr = idxtet(p,q,r);
230  int idxpqrm1 = idxtet(p,q,r-1);
231  jrc(2.0*p+2.0*q+2.0,0.0,r,ar,br,cr);
232  for (int i=0;i<np;i++) {
233  poly_val(idxpqrp1,i) = (ar * (2.0*z(i,2)-1.0) + br) * poly_val( idxpqr , i ) - cr * poly_val(idxpqrm1,i);
234  }
235  }
236  }
237  }
238 
239  return;
240 
241  }
242 
243 } // namespace Intrepid;
244 
245 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
246 #ifdef __GNUC__
247 #warning "The Intrepid package is deprecated"
248 #endif
249 #endif
250 
static int idxtet(int p, int q, int r)
Given indices p,q,r, computes the linear index of the tetrahedral polynomial D^{p,q,r}.
static void tabulateTriangle(const ScalarArray1 &z, const int n, ScalarArray2 &poly_val)
Calculates triangular orthogonal expansions (e.g. Dubiner basis) at a range of input points...
static int idxtri(int p, int q)
Given indices p,q, computes the linear index of the Dubiner polynomial D^{p,q}.
static void tabulateTetrahedron(const ScalarArray1 &z, const int n, ScalarArray2 &poly_val)
Calculates triangular orthogonal expansions (e.g. Dubiner basis) at a range of input points...
static void jrc(const Scalar &alpha, const Scalar &beta, const int &n, Scalar &an, Scalar &bn, Scalar &cn)
computes Jacobi recurrence coefficients of order n with weights a,b so that P^{alpha,beta}_{n+1}(x) = (an x + bn) P^{alpha,beta}_n(x) - cn P^{alpha,beta}_{n-1}(x)