Intrepid
Intrepid_PolylibDef.hpp
Go to the documentation of this file.
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
47//
48// File: Intrepid_PolylibDef.hpp
49//
50// For more information, please see: http://www.nektar.info
51//
52// The MIT License
53//
54// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
55// Department of Aeronautics, Imperial College London (UK), and Scientific
56// Computing and Imaging Institute, University of Utah (USA).
57//
58// License for the specific language governing rights and limitations under
59// Permission is hereby granted, free of charge, to any person obtaining a
60// copy of this software and associated documentation files (the "Software"),
61// to deal in the Software without restriction, including without limitation
62// the rights to use, copy, modify, merge, publish, distribute, sublicense,
63// and/or sell copies of the Software, and to permit persons to whom the
64// Software is furnished to do so, subject to the following conditions:
65//
66// The above copyright notice and this permission notice shall be included
67// in all copies or substantial portions of the Software.
68//
69// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
70// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
71// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
72// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
73// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
74// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
75// DEALINGS IN THE SOFTWARE.
76//
77// Description:
78// This file is redistributed with the Intrepid package. It should be used
79// in accordance with the above MIT license, at the request of the authors.
80// This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
81//
82// Origin: Nektar++ library, http://www.nektar.info, downloaded on
83// March 10, 2009.
84//
86
87
95namespace Intrepid {
96
98#define INTREPID_POLYLIB_STOP 50
99
101#define INTREPID_POLYLIB_POLYNOMIAL_DEFLATION 0
102
103#ifdef INTREPID_POLYLIB_POLYNOMIAL_DEFLATION
105#define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
106#else
108#define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
109#endif
110
111
112template <class Scalar>
113void IntrepidPolylib::zwgj (Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
114 int i;
115 Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
116
117 IntrepidPolylib::jacobz (np,z,alpha,beta);
118 IntrepidPolylib::jacobd (np,z,w,np,alpha,beta);
119
120 fac = std::pow(two,apb + one)*IntrepidPolylib::gammaF(alpha + np + one)*IntrepidPolylib::gammaF(beta + np + one);
121 fac /= IntrepidPolylib::gammaF((Scalar)(np + one))*IntrepidPolylib::gammaF(apb + np + one);
122
123 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
124
125 return;
126}
127
128
129template <class Scalar>
130void IntrepidPolylib::zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
131
132 if(np == 1){
133 z[0] = 0.0;
134 w[0] = 2.0;
135 }
136 else{
137 int i;
138 Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
139
140 z[0] = -one;
141 IntrepidPolylib::jacobz (np-1,z+1,alpha,beta+1);
142 IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
143
144 fac = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
145 fac /= IntrepidPolylib::gammaF((Scalar)np)*(beta + np)*IntrepidPolylib::gammaF(apb + np + 1);
146
147 for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
148 w[0] *= (beta + one);
149 }
150
151 return;
152}
153
154
155template <class Scalar>
156void IntrepidPolylib::zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
157
158 if(np == 1){
159 z[0] = 0.0;
160 w[0] = 2.0;
161 }
162 else{
163 int i;
164 Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
165
166 IntrepidPolylib::jacobz (np-1,z,alpha+1,beta);
167 z[np-1] = one;
168 IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
169
170 fac = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
171 fac /= IntrepidPolylib::gammaF((Scalar)np)*(alpha + np)*IntrepidPolylib::gammaF(apb + np + 1);
172
173 for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
174 w[np-1] *= (alpha + one);
175 }
176
177 return;
178}
179
180
181template <class Scalar>
182void IntrepidPolylib::zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
183
184 if( np == 1 ){
185 z[0] = 0.0;
186 w[0] = 2.0;
187 }
188 else{
189 int i;
190 Scalar fac, one = 1.0, apb = alpha + beta, two = 2.0;
191
192 z[0] = -one;
193 z[np-1] = one;
194 IntrepidPolylib::jacobz (np-2,z + 1,alpha + one,beta + one);
195 IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
196
197 fac = std::pow(two,apb + 1)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
198 fac /= (np-1)*IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha + beta + np + one);
199
200 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
201 w[0] *= (beta + one);
202 w[np-1] *= (alpha + one);
203 }
204
205 return;
206}
207
208
209template <class Scalar>
210void IntrepidPolylib::Dgj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
211{
212
213 Scalar one = 1.0, two = 2.0;
214
215 if (np <= 0){
216 D[0] = 0.0;
217 }
218 else{
219 int i,j;
220 Scalar *pd;
221
222 pd = (Scalar *)malloc(np*sizeof(Scalar));
223 IntrepidPolylib::jacobd(np,z,pd,np,alpha,beta);
224
225 for (i = 0; i < np; i++){
226 for (j = 0; j < np; j++){
227
228 if (i != j)
229 //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
230 D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
231 else
232 D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
233 (two*(one - z[j]*z[j]));
234 }
235 }
236 free(pd);
237 }
238 return;
239}
240
241
242template <class Scalar>
243void IntrepidPolylib::Dgrjm(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
244{
245
246 if (np <= 0){
247 D[0] = 0.0;
248 }
249 else{
250 int i, j;
251 Scalar one = 1.0, two = 2.0;
252 Scalar *pd;
253
254 pd = (Scalar *)malloc(np*sizeof(Scalar));
255
256 pd[0] = std::pow(-one,np-1)*IntrepidPolylib::gammaF((Scalar)np+beta+one);
257 pd[0] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(beta+two);
258 IntrepidPolylib::jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
259 for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
260
261 for (i = 0; i < np; i++) {
262 for (j = 0; j < np; j++){
263 if (i != j)
264 //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
265 D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
266 else {
267 if(j == 0)
268 D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
269 (two*(beta + two));
270 else
271 D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
272 (two*(one - z[j]*z[j]));
273 }
274 }
275 }
276 free(pd);
277 }
278 return;
279}
280
281
282template <class Scalar>
283void IntrepidPolylib::Dgrjp(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
284{
285
286 if (np <= 0){
287 D[0] = 0.0;
288 }
289 else{
290 int i, j;
291 Scalar one = 1.0, two = 2.0;
292 Scalar *pd;
293
294 pd = (Scalar *)malloc(np*sizeof(Scalar));
295
296
297 IntrepidPolylib::jacobd(np-1,z,pd,np-1,alpha+1,beta);
298 for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
299 pd[np-1] = -IntrepidPolylib::gammaF((Scalar)np+alpha+one);
300 pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha+two);
301
302 for (i = 0; i < np; i++) {
303 for (j = 0; j < np; j++){
304 if (i != j)
305 //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
306 D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
307 else {
308 if(j == np-1)
309 D[i*np+j] = (np + alpha + beta + one)*(np - one)/
310 (two*(alpha + two));
311 else
312 D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
313 (two*(one - z[j]*z[j]));
314 }
315 }
316 }
317 free(pd);
318 }
319 return;
320}
321
322
323template <class Scalar>
324void IntrepidPolylib::Dglj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
325{
326
327 if (np <= 0){
328 D[0] = 0.0;
329 }
330 else{
331 int i, j;
332 Scalar one = 1.0, two = 2.0;
333 Scalar *pd;
334
335 pd = (Scalar *)malloc(np*sizeof(Scalar));
336
337 pd[0] = two*std::pow(-one,np)*IntrepidPolylib::gammaF((Scalar)np + beta);
338 pd[0] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(beta + two);
339 IntrepidPolylib::jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
340 for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
341 pd[np-1] = -two*IntrepidPolylib::gammaF((Scalar)np + alpha);
342 pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(alpha + two);
343
344 for (i = 0; i < np; i++) {
345 for (j = 0; j < np; j++){
346 if (i != j)
347 //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
348 D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
349 else {
350 if (j == 0)
351 D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
352 else if (j == np-1)
353 D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
354 else
355 D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
356 (two*(one - z[j]*z[j]));
357 }
358 }
359 }
360 free(pd);
361 }
362 return;
363}
364
365
366template <class Scalar>
367Scalar IntrepidPolylib::hgj (const int i, const Scalar z, const Scalar *zgj,
368 const int np, const Scalar alpha, const Scalar beta)
369{
370
371 Scalar zi, dz, p, pd, h;
372
373 zi = *(zgj+i);
374 dz = z - zi;
375 if (std::abs(dz) < INTREPID_TOL) return 1.0;
376
377 IntrepidPolylib::jacobd (1, &zi, &pd , np, alpha, beta);
378 IntrepidPolylib::jacobfd(1, &z , &p, (Scalar*)0 , np, alpha, beta);
379 h = p/(pd*dz);
380
381 return h;
382}
383
384
385template <class Scalar>
386Scalar IntrepidPolylib::hgrjm (const int i, const Scalar z, const Scalar *zgrj,
387 const int np, const Scalar alpha, const Scalar beta)
388{
389
390 Scalar zi, dz, p, pd, h;
391
392 zi = *(zgrj+i);
393 dz = z - zi;
394 if (std::abs(dz) < INTREPID_TOL) return 1.0;
395
396 IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha, beta + 1);
397 // need to use this routine in case zi = -1 or 1
398 IntrepidPolylib::jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
399 h = (1.0 + zi)*pd + p;
400 IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0, np-1, alpha, beta + 1);
401 h = (1.0 + z )*p/(h*dz);
402
403 return h;
404}
405
406
407template <class Scalar>
408Scalar IntrepidPolylib::hgrjp (const int i, const Scalar z, const Scalar *zgrj,
409 const int np, const Scalar alpha, const Scalar beta)
410{
411
412 Scalar zi, dz, p, pd, h;
413
414 zi = *(zgrj+i);
415 dz = z - zi;
416 if (std::abs(dz) < INTREPID_TOL) return 1.0;
417
418 IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha+1, beta );
419 // need to use this routine in case z = -1 or 1
420 IntrepidPolylib::jacobd (1, &zi, &pd, np-1, alpha+1, beta );
421 h = (1.0 - zi)*pd - p;
422 IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0, np-1, alpha+1, beta);
423 h = (1.0 - z )*p/(h*dz);
424
425 return h;
426}
427
428
429template <class Scalar>
430Scalar IntrepidPolylib::hglj (const int i, const Scalar z, const Scalar *zglj,
431 const int np, const Scalar alpha, const Scalar beta)
432{
433 Scalar one = 1., two = 2.;
434 Scalar zi, dz, p, pd, h;
435
436 zi = *(zglj+i);
437 dz = z - zi;
438 if (std::abs(dz) < INTREPID_TOL) return 1.0;
439
440 IntrepidPolylib::jacobfd(1, &zi, &p , (Scalar*)0, np-2, alpha + one, beta + one);
441 // need to use this routine in case z = -1 or 1
442 IntrepidPolylib::jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
443 h = (one - zi*zi)*pd - two*zi*p;
444 IntrepidPolylib::jacobfd(1, &z, &p, (Scalar*)0, np-2, alpha + one, beta + one);
445 h = (one - z*z)*p/(h*dz);
446
447 return h;
448}
449
450
451template <class Scalar>
452void IntrepidPolylib::Imgj(Scalar *im, const Scalar *zgj, const Scalar *zm, const int nz,
453 const int mz, const Scalar alpha, const Scalar beta){
454 Scalar zp;
455 int i, j;
456
457 /* old Polylib code */
458 for (i = 0; i < mz; ++i) {
459 zp = zm[i];
460 for (j = 0; j < nz; ++j) {
461 im[i*nz+j] = IntrepidPolylib::hgj(j, zp, zgj, nz, alpha, beta);
462 }
463 }
464
465 /* original Nektar++ code: is this correct???
466 for (i = 0; i < nz; ++i) {
467 for (j = 0; j < mz; ++j)
468 {
469 zp = zm[j];
470 im [i*mz+j] = IntrepidPolylib::hgj(i, zp, zgj, nz, alpha, beta);
471 }
472 }
473 */
474
475 return;
476}
477
478
479template <class Scalar>
480void IntrepidPolylib::Imgrjm(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz,
481 const int mz, const Scalar alpha, const Scalar beta)
482{
483 Scalar zp;
484 int i, j;
485
486 for (i = 0; i < mz; i++) {
487 zp = zm[i];
488 for (j = 0; j < nz; j++) {
489 im[i*nz+j] = IntrepidPolylib::hgrjm(j, zp, zgrj, nz, alpha, beta);
490 }
491 }
492
493 /* original Nektar++ code: is this correct???
494 for (i = 0; i < nz; i++) {
495 for (j = 0; j < mz; j++)
496 {
497 zp = zm[j];
498 im [i*mz+j] = IntrepidPolylib::hgrjm(i, zp, zgrj, nz, alpha, beta);
499 }
500 }
501 */
502
503 return;
504}
505
506
507template <class Scalar>
508void IntrepidPolylib::Imgrjp(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz,
509 const int mz, const Scalar alpha, const Scalar beta)
510{
511 Scalar zp;
512 int i, j;
513
514 for (i = 0; i < mz; i++) {
515 zp = zm[i];
516 for (j = 0; j < nz; j++) {
517 im [i*nz+j] = IntrepidPolylib::hgrjp(j, zp, zgrj, nz, alpha, beta);
518 }
519 }
520
521 /* original Nektar++ code: is this correct?
522 for (i = 0; i < nz; i++) {
523 for (j = 0; j < mz; j++)
524 {
525 zp = zm[j];
526 im [i*mz+j] = IntrepidPolylib::hgrjp(i, zp, zgrj, nz, alpha, beta);
527 }
528 }
529 */
530
531 return;
532}
533
534
535template <class Scalar>
536void IntrepidPolylib::Imglj(Scalar *im, const Scalar *zglj, const Scalar *zm, const int nz,
537 const int mz, const Scalar alpha, const Scalar beta)
538{
539 Scalar zp;
540 int i, j;
541
542 for (i = 0; i < mz; i++) {
543 zp = zm[i];
544 for (j = 0; j < nz; j++) {
545 im[i*nz+j] = IntrepidPolylib::hglj(j, zp, zglj, nz, alpha, beta);
546 }
547 }
548
549 /* original Nektar++ code: is this correct?
550 for (i = 0; i < nz; i++) {
551 for (j = 0; j < mz; j++)
552 {
553 zp = zm[j];
554 im[i*mz+j] = IntrepidPolylib::hglj(i, zp, zglj, nz, alpha, beta);
555 }
556 }
557 */
558
559 return;
560}
561
562
563template <class Scalar>
564void
566jacobfd (const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd,
567 const int n, const Scalar alpha, const Scalar beta)
568{
569 const Scalar zero = 0.0, one = 1.0, two = 2.0;
570
571 if (! np) {
572 return;
573 }
574
575 if (n == 0) {
576 if (poly_in) {
577 for (int i = 0; i < np; ++i) {
578 poly_in[i] = one;
579 }
580 }
581 if (polyd) {
582 for (int i = 0; i < np; ++i) {
583 polyd[i] = zero;
584 }
585 }
586 }
587 else if (n == 1) {
588 if (poly_in) {
589 for (int i = 0; i < np; ++i) {
590 poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
591 }
592 }
593 if (polyd) {
594 for (int i = 0; i < np; ++i) {
595 polyd[i] = 0.5*(alpha + beta + two);
596 }
597 }
598 }
599 else {
600 Scalar a1,a2,a3,a4;
601 Scalar apb = alpha + beta;
602 Scalar *poly, *polyn1,*polyn2;
603
604 if (poly_in) { // switch for case of no poynomial function return
605 polyn1 = (Scalar *)malloc(2*np*sizeof(Scalar));
606 polyn2 = polyn1+np;
607 poly = poly_in;
608 }
609 else{
610 polyn1 = (Scalar *)malloc(3*np*sizeof(Scalar));
611 polyn2 = polyn1+np;
612 poly = polyn2+np;
613 }
614
615 for (int i = 0; i < np; ++i) {
616 polyn2[i] = one;
617 polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
618 }
619
620 for (int k = 2; k <= n; ++k) {
621 a1 = two*k*(k + apb)*(two*k + apb - two);
622 a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
623 a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
624 a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
625
626 a2 /= a1;
627 a3 /= a1;
628 a4 /= a1;
629
630 for (int i = 0; i < np; ++i) {
631 poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
632 polyn2[i] = polyn1[i];
633 polyn1[i] = poly [i];
634 }
635 }
636
637 if (polyd) {
638 a1 = n*(alpha - beta);
639 a2 = n*(two*n + alpha + beta);
640 a3 = two*(n + alpha)*(n + beta);
641 a4 = (two*n + alpha + beta);
642 a1 /= a4; a2 /= a4; a3 /= a4;
643
644 // note polyn2 points to polyn1 at end of poly iterations
645 for (int i = 0; i < np; ++i) {
646 polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
647 polyd[i] /= (one - z[i]*z[i]);
648 }
649 }
650
651 free(polyn1);
652 }
653}
654
655
656template <class Scalar>
657void IntrepidPolylib::jacobd(const int np, const Scalar *z, Scalar *polyd, const int n,
658 const Scalar alpha, const Scalar beta)
659{
660 int i;
661 Scalar one = 1.0;
662 if(n == 0)
663 for(i = 0; i < np; ++i) polyd[i] = 0.0;
664 else{
665 //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
666 IntrepidPolylib::jacobfd(np,z,polyd,(Scalar*)0,n-1,alpha+one,beta+one);
667 for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (Scalar)n + one);
668 }
669 return;
670}
671
672
673template <class Scalar>
674void IntrepidPolylib::Jacobz(const int n, Scalar *z, const Scalar alpha, const Scalar beta){
675 int i,j,k;
676 Scalar dth = M_PI/(2.0*(Scalar)n);
677 Scalar poly,pder,rlast=0.0;
678 Scalar sum,delr,r;
679 Scalar one = 1.0, two = 2.0;
680
681 if(!n)
682 return;
683
684 for(k = 0; k < n; ++k){
685 r = -std::cos((two*(Scalar)k + one) * dth);
686 if(k) r = 0.5*(r + rlast);
687
688 for(j = 1; j < INTREPID_POLYLIB_STOP; ++j){
689 IntrepidPolylib::jacobfd(1,&r,&poly, &pder, n, alpha, beta);
690
691 for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
692
693 delr = -poly / (pder - sum * poly);
694 r += delr;
695 if( std::abs(delr) < INTREPID_TOL ) break;
696 }
697 z[k] = r;
698 rlast = r;
699 }
700 return;
701}
702
703
704template <class Scalar>
705void IntrepidPolylib::JacZeros(const int n, Scalar *a, const Scalar alpha, const Scalar beta){
706 int i;
707 Scalar apb, apbi,a2b2;
708 Scalar *b;
709
710 if(!n)
711 return;
712
713 b = (Scalar *) malloc(n*sizeof(Scalar));
714
715 // generate normalised terms
716 apb = alpha + beta;
717 apbi = 2.0 + apb;
718
719 b[n-1] = std::pow(2.0,apb+1.0)*IntrepidPolylib::gammaF(alpha+1.0)*IntrepidPolylib::gammaF(beta+1.0)/gammaF(apbi);
720 a[0] = (beta-alpha)/apbi;
721 b[0] = std::sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
722
723 a2b2 = beta*beta-alpha*alpha;
724 for(i = 1; i < n-1; ++i){
725 apbi = 2.0*(i+1) + apb;
726 a[i] = a2b2/((apbi-2.0)*apbi);
727 b[i] = std::sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
728 ((apbi*apbi-1)*apbi*apbi));
729 }
730
731 apbi = 2.0*n + apb;
732 //a[n-1] = a2b2/((apbi-2.0)*apbi); // THIS IS A BUG!!!
733 if (n>1) a[n-1] = a2b2/((apbi-2.0)*apbi);
734
735 // find eigenvalues
736 IntrepidPolylib::TriQL(n, a, b);
737
738 free(b);
739 return;
740}
741
742
743template <class Scalar>
744void IntrepidPolylib::TriQL(const int n, Scalar *d,Scalar *e) {
745 int m,l,iter,i,k;
746 Scalar s,r,p,g,f,dd,c,b;
747
748 for (l=0;l<n;l++) {
749 iter=0;
750 do {
751 for (m=l;m<n-1;m++) {
752 dd=std::abs(d[m])+std::abs(d[m+1]);
753 if (std::abs(e[m])+dd == dd) break;
754 }
755 if (m != l) {
756 if (iter++ == 50){
757 TEUCHOS_TEST_FOR_EXCEPTION((1),
758 std::runtime_error,
759 ">>> ERROR (IntrepidPolylib): Too many iterations in TQLI.");
760 }
761 g=(d[l+1]-d[l])/(2.0*e[l]);
762 r=std::sqrt((g*g)+1.0);
763 //g=d[m]-d[l]+e[l]/(g+sign(r,g));
764 g=d[m]-d[l]+e[l]/(g+((g)<0 ? -std::abs(r) : std::abs(r)));
765 s=c=1.0;
766 p=0.0;
767 for (i=m-1;i>=l;i--) {
768 f=s*e[i];
769 b=c*e[i];
770 if (std::abs(f) >= std::abs(g)) {
771 c=g/f;
772 r=std::sqrt((c*c)+1.0);
773 e[i+1]=f*r;
774 c *= (s=1.0/r);
775 } else {
776 s=f/g;
777 r=std::sqrt((s*s)+1.0);
778 e[i+1]=g*r;
779 s *= (c=1.0/r);
780 }
781 g=d[i+1]-p;
782 r=(d[i]-g)*s+2.0*c*b;
783 p=s*r;
784 d[i+1]=g+p;
785 g=c*r-b;
786 }
787 d[l]=d[l]-p;
788 e[l]=g;
789 e[m]=0.0;
790 }
791 } while (m != l);
792 }
793
794 // order eigenvalues
795 for(i = 0; i < n-1; ++i){
796 k = i;
797 p = d[i];
798 for(l = i+1; l < n; ++l)
799 if (d[l] < p) {
800 k = l;
801 p = d[l];
802 }
803 d[k] = d[i];
804 d[i] = p;
805 }
806}
807
808
809template <class Scalar>
810Scalar IntrepidPolylib::gammaF(const Scalar x){
811 Scalar gamma = 1.0;
812
813 if (x == -0.5) gamma = -2.0*std::sqrt(M_PI);
814 else if (!x) return gamma;
815 else if ((x-(int)x) == 0.5){
816 int n = (int) x;
817 Scalar tmp = x;
818
819 gamma = std::sqrt(M_PI);
820 while(n--){
821 tmp -= 1.0;
822 gamma *= tmp;
823 }
824 }
825 else if ((x-(int)x) == 0.0){
826 int n = (int) x;
827 Scalar tmp = x;
828
829 while(--n){
830 tmp -= 1.0;
831 gamma *= tmp;
832 }
833 }
834 else
835 TEUCHOS_TEST_FOR_EXCEPTION((1),
836 std::invalid_argument,
837 ">>> ERROR (IntrepidPolylib): Argument is not of integer or half order.");
838 return gamma;
839}
840
841} // end of namespace Intrepid
842
#define INTREPID_POLYLIB_STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
static void Jacobz(const int n, Scalar *z, const Scalar alpha, const Scalar beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e. .
static void zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
static void JacZeros(const int n, Scalar *a, const Scalar alpha, const Scalar beta)
Zero determination through the eigenvalues of a tridiagonal matrix from the three term recursion rela...
static void Imgj(Scalar *im, const Scalar *zgj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.
static Scalar hgrjp(const int i, const Scalar z, const Scalar *zgrj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
static void jacobfd(const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
static void zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
static Scalar hgrjm(const int i, const Scalar z, const Scalar *zgrj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
static void Imgrjp(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at ...
static void Dgrjp(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.
static void Imglj(Scalar *im, const Scalar *zglj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.
static void zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
static void Dgj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
static void TriQL(const int n, Scalar *d, Scalar *e)
QL algorithm for symmetric tridiagonal matrix.
static void Dgrjm(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
static void Dglj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
static void Imgrjm(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at...
static Scalar hglj(const int i, const Scalar z, const Scalar *zglj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zglj ...
static void jacobd(const int np, const Scalar *z, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Calculate the derivative of Jacobi polynomials.
static Scalar gammaF(const Scalar x)
Calculate the Gamma function , , for integer values x and halves.
static Scalar hgj(const int i, const Scalar z, const Scalar *zgj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the ar...
static void zwgj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Jacobi zeros and weights.