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
46namespace 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;
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 idxtet(int p, int q, int r)
Given indices p,q,r, computes the linear index of the tetrahedral polynomial D^{p,...
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,...
static int idxtri(int p, int q)
Given indices p,q, computes the linear index of the Dubiner polynomial D^{p,q}.