LORENE
FFTW3/circhebpimi.C
1/*
2 * Copyright (c) 1999-2001 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char circhebpimi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/circhebpimi.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $" ;
24
25/*
26 * Transformation de Tchebyshev inverse T_{2k+1}/T_{2k} (suivant la parite de
27 * l'indice m en phi) sur le troisieme indice
28 * (indice correspondant a r) d'un tableau 3-D decrivant par exemple
29 * d/dr d'une fonction symetrique
30 * par rapport au plan equatorial z = 0 et sans aucune autre symetrie,
31 * cad que l'on a effectue
32 * 1/ un developpement en polynomes de Tchebyshev impairs pour m pair
33 * 2/ un developpement en polynomes de Tchebyshev pairs pour m impair
34 *
35 *
36 * Utilise la bibliotheque fftw.
37 *
38 * Entree:
39 * -------
40 * int* deg : tableau du nombre effectif de degres de liberte dans chacune
41 * des 3 dimensions: le nombre de points de collocation
42 * en r est nr = deg[2] et doit etre de la forme
43 * nr = 2*p + 1
44 * int* dimc : tableau du nombre d'elements de cf dans chacune des trois
45 * dimensions.
46 * On doit avoir dimc[2] >= deg[2] = nr.
47 *
48 * double* cf : tableau des coefficients c_i de la fonction definis
49 * comme suit (a theta et phi fixes)
50 *
51 * -- pour m pair (i.e. j = 0, 1, 4, 5, 8, 9, ...) :
52 *
53 * f(x) = som_{i=0}^{nr-2} c_i T_{2i+1}(x)
54 *
55 * ou T_{2i+1}(x) designe le polynome de Tchebyshev de
56 * degre 2i+1.
57 *
58 * -- pour m impair (i.e. j = 2, 3, 6, 7, 10, 11, ...) :
59 *
60 * f(x) = som_{i=0}^{nr-1} c_i T_{2i}(x) ,
61 *
62 * ou T_{2i}(x) designe le polynome de Tchebyshev de
63 * degre 2i.
64 *
65 * Les coefficients c_i (0 <= i <= nr-1) doivent etre stokes
66 * dans le tableau cf comme suit
67 * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
68 * ou j et k sont les indices correspondant a phi et theta
69 * respectivement.
70 * L'espace memoire correspondant a ce pointeur doit etre
71 * dimc[0]*dimc[1]*dimc[2] et doit etre alloue avant l'appel a
72 * la routine.
73 *
74 * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
75 * dimensions.
76 * On doit avoir dimf[2] >= deg[2] = nr.
77 *
78 * Sortie:
79 * -------
80 * double* ff : tableau des valeurs de la fonction aux nr points de
81 * de collocation
82 *
83 * x_i = sin( pi/2 i/(nr-1) ) 0 <= i <= nr-1
84 *
85 * Les valeurs de la fonction sont stokees dans le
86 * tableau ff comme suit
87 * f( x_i ) = ff[ dimf[1]*dimf[2] * j + dimf[2] * k + i ]
88 * ou j et k sont les indices correspondant a phi et theta
89 * respectivement.
90 * L'espace memoire correspondant a ce pointeur doit etre
91 * dimf[0]*dimf[1]*dimf[2] et doit avoir ete alloue avant
92 * l'appel a la routine.
93 *
94 * NB: Si le pointeur cf est egal a ff, la routine ne travaille que sur un
95 * seul tableau, qui constitue une entree/sortie.
96 */
97
98/*
99 * $Id: circhebpimi.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $
100 * $Log: circhebpimi.C,v $
101 * Revision 1.4 2014/10/13 08:53:20 j_novak
102 * Lorene classes and functions now belong to the namespace Lorene.
103 *
104 * Revision 1.3 2014/10/06 15:18:49 j_novak
105 * Modified #include directives to use c++ syntax.
106 *
107 * Revision 1.2 2004/12/27 14:27:28 j_novak
108 * Added forgotten "delete [] t1"
109 *
110 * Revision 1.1 2004/12/21 17:06:02 j_novak
111 * Added all files for using fftw3.
112 *
113 * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
114 * Suppressed the directive #include <malloc.h> for malloc is defined
115 * in <stdlib.h>
116 *
117 * Revision 1.3 2002/10/16 14:36:53 j_novak
118 * Reorganization of #include instructions of standard C++, in order to
119 * use experimental version 3 of gcc.
120 *
121 * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
122 * Modification of declaration of Fortran 77 prototypes for
123 * a better portability (in particular on IBM AIX systems):
124 * All Fortran subroutine names are now written F77_* and are
125 * defined in the new file C++/Include/proto_f77.h.
126 *
127 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
128 * LORENE
129 *
130 * Revision 2.0 1999/02/22 15:43:19 hyc
131 * *** empty log message ***
132 *
133 *
134 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/circhebpimi.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $
135 *
136 */
137
138
139
140// headers du C
141#include <cassert>
142#include <cstdlib>
143#include <fftw3.h>
144
145//Lorene prototypes
146#include "tbl.h"
147
148// Prototypage des sous-routines utilisees:
149namespace Lorene {
150fftw_plan back_fft(int, Tbl*&) ;
151double* cheb_ini(const int) ;
152double* chebimp_ini(const int ) ;
153//*****************************************************************************
154
155void circhebpimi(const int* deg, const int* dimc, double* cf,
156 const int* dimf, double* ff)
157
158{
159int i, j, k ;
160
161// Dimensions des tableaux ff et cf :
162 int n1f = dimf[0] ;
163 int n2f = dimf[1] ;
164 int n3f = dimf[2] ;
165 int n1c = dimc[0] ;
166 int n2c = dimc[1] ;
167 int n3c = dimc[2] ;
168
169// Nombres de degres de liberte en r :
170 int nr = deg[2] ;
171
172// Tests de dimension:
173 if (nr > n3c) {
174 cout << "circhebpimi: nr > n3c : nr = " << nr << " , n3c = "
175 << n3c << endl ;
176 abort () ;
177 exit(-1) ;
178 }
179 if (nr > n3f) {
180 cout << "circhebpimi: nr > n3f : nr = " << nr << " , n3f = "
181 << n3f << endl ;
182 abort () ;
183 exit(-1) ;
184 }
185 if (n1c > n1f) {
186 cout << "circhebpimi: n1c > n1f : n1c = " << n1c << " , n1f = "
187 << n1f << endl ;
188 abort () ;
189 exit(-1) ;
190 }
191 if (n2c > n2f) {
192 cout << "circhebpimi: n2c > n2f : n2c = " << n2c << " , n2f = "
193 << n2f << endl ;
194 abort () ;
195 exit(-1) ;
196 }
197
198// Nombre de points pour la FFT:
199 int nm1 = nr - 1;
200 int nm1s2 = nm1 / 2;
201
202// Recherche des tables pour la FFT:
203 Tbl* pg = 0x0 ;
204 fftw_plan p = back_fft(nm1, pg) ;
205 Tbl& g = *pg ;
206 double* t1 = new double[nr] ;
207
208// Recherche de la table des sin(psi) :
209 double* sinp = cheb_ini(nr);
210
211// Recherche de la table des points de collocations x_k :
212 double* x = chebimp_ini(nr);
213
214// boucle sur phi et theta
215
216 int n2n3f = n2f * n3f ;
217 int n2n3c = n2c * n3c ;
218
219//=======================================================================
220// Cas m pair
221//=======================================================================
222
223 j = 0 ;
224
225 while (j<n1f-1) { //le dernier coef en phi n'est pas considere
226 // (car nul)
227
228
229//------------------------------------------------------------------------
230// partie cos(m phi) avec m pair : developpement en T_{2i+1}(x)
231//------------------------------------------------------------------------
232
233 for (k=0; k<n2c; k++) {
234
235 int i0 = n2n3c * j + n3c * k ; // indice de depart
236 double* cf0 = cf + i0 ; // tableau des donnees a transformer
237
238 i0 = n2n3f * j + n3f * k ; // indice de depart
239 double* ff0 = ff + i0 ; // tableau resultat
240
241// Calcul des coefficients du developpement en T_{2i}(x) de la fonction
242// h(x) := x f(x) a partir des coefficients de f (resultat stoke dans le
243// tableau t1 :
244 t1[0] = .5 * cf0[0] ;
245 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
246 t1[nm1] = .5 * cf0[nr-2] ;
247
248/*
249 * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
250 * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
251 */
252
253// Calcul des coefficients de Fourier de la fonction
254// G(psi) = F+(psi) + F_(psi) sin(psi)
255// en fonction des coefficients de Tchebyshev de f:
256
257// Coefficients impairs de G
258//--------------------------
259
260 double c1 = t1[1] ;
261
262 double som = 0;
263 ff0[1] = 0 ;
264 for ( i = 3; i < nr; i += 2 ) {
265 ff0[i] = t1[i] - c1 ;
266 som += ff0[i] ;
267 }
268
269// Valeur en psi=0 de la partie antisymetrique de F, F_ :
270 double fmoins0 = nm1s2 * c1 + som ;
271
272// Coef. impairs de G
273// NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
274// donnait exactement les coef. des sinus, ce facteur serait -0.5.
275 for ( i = 3; i < nr; i += 2 ) {
276 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
277 }
278
279
280// Coefficients pairs de G
281//------------------------
282// Ces coefficients sont egaux aux coefficients pairs du developpement de
283// f.
284// NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
285// donnait exactement les coef. des cosinus, ce facteur serait 1.
286
287 g.set(0) = t1[0] ;
288 for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * t1[2*i] ;
289 g.set(nm1s2) = t1[nm1] ;
290
291// Transformation de Fourier inverse de G
292//---------------------------------------
293
294// FFT inverse
295 fftw_execute(p) ;
296
297// Valeurs de f deduites de celles de G
298//-------------------------------------
299
300 for ( i = 1; i < nm1s2 ; i++ ) {
301// ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
302 int isym = nm1 - i ;
303// ... indice (dans le tableau ff0) du point x correspondant a psi
304 int ix = nm1 - i ;
305// ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
306 int ixsym = nm1 - isym ;
307
308 double fp = .5 * ( g(i) + g(isym) ) ;
309 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
310
311 ff0[ix] = ( fp + fm ) / x[ix];
312 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
313 }
314
315//... cas particuliers:
316 ff0[0] = 0 ;
317 ff0[nm1] = g(0) + fmoins0 ;
318 ff0[nm1s2] = g(nm1s2) / x[nm1s2] ;
319
320 } // fin de la boucle sur theta
321
322//------------------------------------------------------------------------
323// partie sin(m phi) avec m pair : developpement en T_{2i+1}(x)
324//------------------------------------------------------------------------
325
326 j++ ;
327
328 if ( (j != 1) && (j != n1f-1) ) {
329// on effectue le calcul seulement dans les cas ou les coef en phi ne sont
330// pas nuls
331
332 for (k=0; k<n2c; k++) {
333
334 int i0 = n2n3c * j + n3c * k ; // indice de depart
335 double* cf0 = cf + i0 ; // tableau des donnees a transformer
336
337 i0 = n2n3f * j + n3f * k ; // indice de depart
338 double* ff0 = ff + i0 ; // tableau resultat
339
340// Calcul des coefficients du developpement en T_{2i}(x) de la fonction
341// h(x) := x f(x) a partir des coefficients de f (resultat stoke dans le
342// tableau t1 :
343 t1[0] = .5 * cf0[0] ;
344 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
345 t1[nm1] = .5 * cf0[nr-2] ;
346
347/*
348 * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
349 * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
350 */
351
352// Calcul des coefficients de Fourier de la fonction
353// G(psi) = F+(psi) + F_(psi) sin(psi)
354// en fonction des coefficients de Tchebyshev de f:
355
356// Coefficients impairs de G
357//--------------------------
358
359 double c1 = t1[1] ;
360
361 double som = 0;
362 ff0[1] = 0 ;
363 for ( i = 3; i < nr; i += 2 ) {
364 ff0[i] = t1[i] - c1 ;
365 som += ff0[i] ;
366 }
367
368// Valeur en psi=0 de la partie antisymetrique de F, F_ :
369 double fmoins0 = nm1s2 * c1 + som ;
370
371// Coef. impairs de G
372// NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
373// donnait exactement les coef. des sinus, ce facteur serait -0.5.
374 for ( i = 3; i < nr; i += 2 ) {
375 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
376 }
377
378
379// Coefficients pairs de G
380//------------------------
381// Ces coefficients sont egaux aux coefficients pairs du developpement de
382// f.
383// NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
384// donnait exactement les coef. des cosinus, ce facteur serait 1.
385
386 g.set(0) = t1[0] ;
387 for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * t1[2*i] ;
388 g.set(nm1s2) = t1[nm1] ;
389
390// Transformation de Fourier inverse de G
391//---------------------------------------
392
393// FFT inverse
394 fftw_execute(p) ;
395
396// Valeurs de f deduites de celles de G
397//-------------------------------------
398
399 for ( i = 1; i < nm1s2 ; i++ ) {
400// ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
401 int isym = nm1 - i ;
402// ... indice (dans le tableau ff0) du point x correspondant a psi
403 int ix = nm1 - i ;
404// ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
405 int ixsym = nm1 - isym ;
406
407 double fp = .5 * ( g(i) + g(isym) ) ;
408 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
409
410 ff0[ix] = ( fp + fm ) / x[ix];
411 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
412 }
413
414//... cas particuliers:
415 ff0[0] = 0 ;
416 ff0[nm1] = g(0) + fmoins0 ;
417 ff0[nm1s2] = g(nm1s2) / x[nm1s2] ;
418
419 } // fin de la boucle sur theta
420
421 } // fin du cas ou le calcul etait necessaire (i.e. ou les
422 // coef en phi n'etaient pas nuls)
423
424// On passe au cas m pair suivant:
425 j+=3 ;
426
427 } // fin de la boucle sur les cas m pair
428
429 if (n1f<=3) { // cas m=0 seulement (symetrie axiale)
430 delete [] t1 ;
431 return ;
432 }
433
434//=======================================================================
435// Cas m impair
436//=======================================================================
437
438 j = 2 ;
439
440 while (j<n1f-1) {
441
442//--------------------------------------------------------------------
443// partie cos(m phi) avec m impair : developpement en T_{2i}(x)
444//--------------------------------------------------------------------
445
446 for (k=0; k<n2c; k++) {
447
448 int i0 = n2n3c * j + n3c * k ; // indice de depart
449 double* cf0 = cf + i0 ; // tableau des donnees a transformer
450
451 i0 = n2n3f * j + n3f * k ; // indice de depart
452 double* ff0 = ff + i0 ; // tableau resultat
453
454/*
455 * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
456 * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
457 */
458
459// Calcul des coefficients de Fourier de la fonction
460// G(psi) = F+(psi) + F_(psi) sin(psi)
461// en fonction des coefficients de Tchebyshev de f:
462
463// Coefficients impairs de G
464//--------------------------
465
466 double c1 = cf0[1] ;
467
468 double som = 0;
469 ff0[1] = 0 ;
470 for ( i = 3; i < nr; i += 2 ) {
471 ff0[i] = cf0[i] - c1 ;
472 som += ff0[i] ;
473 }
474
475// Valeur en psi=0 de la partie antisymetrique de F, F_ :
476 double fmoins0 = nm1s2 * c1 + som ;
477
478// Coef. impairs de G
479// NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
480// donnait exactement les coef. des sinus, ce facteur serait -0.5.
481 for ( i = 3; i < nr; i += 2 ) {
482 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
483 }
484
485
486// Coefficients pairs de G
487//------------------------
488// Ces coefficients sont egaux aux coefficients pairs du developpement de
489// f.
490// NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
491// donnait exactement les coef. des cosinus, ce facteur serait 1.
492
493 g.set(0) = cf0[0] ;
494 for (i=1; i<nm1s2; i++) g.set(i) = 0.5 * cf0[2*i] ;
495 g.set(nm1s2) = cf0[nm1] ;
496
497// Transformation de Fourier inverse de G
498//---------------------------------------
499
500// FFT inverse
501 fftw_execute(p) ;
502
503// Valeurs de f deduites de celles de G
504//-------------------------------------
505
506 for ( i = 1; i < nm1s2 ; i++ ) {
507// ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
508 int isym = nm1 - i ;
509// ... indice (dans le tableau ff0) du point x correspondant a psi
510 int ix = nm1 - i ;
511// ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
512 int ixsym = nm1 - isym ;
513
514 double fp = .5 * ( g(i) + g(isym) ) ;
515 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
516
517 ff0[ix] = fp + fm ;
518 ff0[ixsym] = fp - fm ;
519 }
520
521//... cas particuliers:
522 ff0[0] = g(0) - fmoins0 ;
523 ff0[nm1] = g(0) + fmoins0 ;
524 ff0[nm1s2] = g(nm1s2) ;
525
526 } // fin de la boucle sur theta
527
528
529//--------------------------------------------------------------------
530// partie sin(m phi) avec m impair : developpement en T_{2i}(x)
531//--------------------------------------------------------------------
532
533 j++ ;
534
535 if ( j != n1f-1 ) {
536// on effectue le calcul seulement dans les cas ou les coef en phi ne sont
537// pas nuls
538
539 for (k=0; k<n2c; k++) {
540
541 int i0 = n2n3c * j + n3c * k ; // indice de depart
542 double* cf0 = cf + i0 ; // tableau des donnees a transformer
543
544 i0 = n2n3f * j + n3f * k ; // indice de depart
545 double* ff0 = ff + i0 ; // tableau resultat
546
547/*
548 * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
549 * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
550 */
551
552// Calcul des coefficients de Fourier de la fonction
553// G(psi) = F+(psi) + F_(psi) sin(psi)
554// en fonction des coefficients de Tchebyshev de f:
555
556// Coefficients impairs de G
557//--------------------------
558
559 double c1 = cf0[1] ;
560
561 double som = 0;
562 ff0[1] = 0 ;
563 for ( i = 3; i < nr; i += 2 ) {
564 ff0[i] = cf0[i] - c1 ;
565 som += ff0[i] ;
566 }
567
568// Valeur en psi=0 de la partie antisymetrique de F, F_ :
569 double fmoins0 = nm1s2 * c1 + som ;
570
571// Coef. impairs de G
572// NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
573// donnait exactement les coef. des sinus, ce facteur serait -0.5.
574 for ( i = 3; i < nr; i += 2 ) {
575 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
576 }
577
578
579// Coefficients pairs de G
580//------------------------
581// Ces coefficients sont egaux aux coefficients pairs du developpement de
582// f.
583// NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
584// donnait exactement les coef. des cosinus, ce facteur serait 1.
585
586 g.set(0) = cf0[0] ;
587 for (i=1; i<nm1s2; i++) g.set(i) = 0.5 * cf0[2*i] ;
588 g.set(nm1s2) = cf0[nm1] ;
589
590// Transformation de Fourier inverse de G
591//---------------------------------------
592
593// FFT inverse
594 fftw_execute(p) ;
595
596// Valeurs de f deduites de celles de G
597//-------------------------------------
598
599 for ( i = 1; i < nm1s2 ; i++ ) {
600// ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
601 int isym = nm1 - i ;
602// ... indice (dans le tableau ff0) du point x correspondant a psi
603 int ix = nm1 - i ;
604// ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
605 int ixsym = nm1 - isym ;
606
607 double fp = .5 * ( g(i) + g(isym) ) ;
608 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
609
610 ff0[ix] = fp + fm ;
611 ff0[ixsym] = fp - fm ;
612 }
613
614//... cas particuliers:
615 ff0[0] = g(0) - fmoins0 ;
616 ff0[nm1] = g(0) + fmoins0 ;
617 ff0[nm1s2] = g(nm1s2) ;
618
619 } // fin de la boucle sur theta
620
621 } // fin du cas ou le calcul etait necessaire (i.e. ou les
622 // coef en phi n'etaient pas nuls)
623
624// On passe au cas m impair suivant:
625 j+=3 ;
626
627 } // fin de la boucle sur les cas m impair
628
629 delete [] t1 ;
630}
631
632}
Lorene prototypes.
Definition app_hor.h:64