LORENE
FFTW3/cfrchebpii.C
1/*
2 * Copyright (c) 1999-2002 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 cfrchebpii_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfrchebpii.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $" ;
24
25
26/*
27 * Transformation de Tchebyshev (cas rare) sur le troisieme indice (indice
28 * correspondant a r) d'un tableau 3-D decrivant une fonction quelconque.
29 * Utilise la bibliotheque fftw.
30 *
31 *
32 * Entree:
33 * -------
34 * int* deg : tableau du nombre effectif de degres de liberte dans chacune
35 * des 3 dimensions: le nombre de points de collocation
36 * en r est nr = deg[2] et doit etre de la forme
37 * nr = 2*p + 1
38 * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
39 * dimensions.
40 * On doit avoir dimf[2] >= deg[2] = nr.
41 * NB: pour dimf[0] = 1 (un seul point en phi), la transformation
42 * est bien effectuee.
43 * pour dimf[0] > 1 (plus d'un point en phi), la
44 * transformation n'est effectuee que pour les indices (en phi)
45 * j != 1 et j != dimf[0]-1 (cf. commentaires sur borne_phi).
46 *
47 * double* ff : tableau des valeurs de la fonction aux nr points de
48 * de collocation
49 *
50 * x_i = sin( pi/2 i/(nr-1) ) 0 <= i <= nr-1
51 *
52 * Les valeurs de la fonction doivent etre stokees dans le
53 * tableau ff comme suit
54 * f( x_i ) = ff[ dimf[1]*dimf[2] * j + dimf[2] * k + i ]
55 * ou j et k sont les indices correspondant a phi et theta
56 * respectivement.
57 * L'espace memoire correspondant a ce pointeur doit etre
58 * dimf[0]*dimf[1]*dimf[2] et doit etre alloue avant l'appel a
59 * la routine.
60 *
61 * int* dimc : tableau du nombre d'elements de cf dans chacune des trois
62 * dimensions.
63 * On doit avoir dimc[2] >= deg[2] = nr.
64 *
65 * Sortie:
66 * -------
67 * double* cf : tableau des nr coefficients c_i de la fonction definis
68 * comme suit (a theta et phi fixes)
69 *
70 * Si l impair f(x) = som_{i=0}^{nr-1} c_i T_{2i}(x) ,
71 * Si l pair f(x) = som_{i=0}^{nr-1} c_i T_{2i+1}(x)
72 *
73 * ou T_{i}(x) designe le polynome de Tchebyshev de degre i.
74 * Les coefficients c_i (0 <= i <= nr-1) sont stokes dans
75 * le tableau cf comme suit
76 * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
77 * ou j et k sont les indices correspondant a phi et theta
78 * respectivement.
79 * L'espace memoire correspondant a ce pointeur doit etre
80 * dimc[0]*dimc[1]*dimc[2] et doit avoir ete alloue avant
81 * l'appel a la routine.
82 *
83 * NB: Si le pointeur ff est egal a cf, la routine ne travaille que sur un
84 * seul tableau, qui constitue une entree/sortie.
85 */
86
87/*
88 * $Id: cfrchebpii.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
89 * $Log: cfrchebpii.C,v $
90 * Revision 1.3 2014/10/13 08:53:18 j_novak
91 * Lorene classes and functions now belong to the namespace Lorene.
92 *
93 * Revision 1.2 2014/10/06 15:18:47 j_novak
94 * Modified #include directives to use c++ syntax.
95 *
96 * Revision 1.1 2004/12/21 17:06:02 j_novak
97 * Added all files for using fftw3.
98 *
99 * Revision 1.1 2004/11/23 15:13:50 m_forot
100 * Added the bases for the cases without any equatorial symmetry
101 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
102 *
103 *
104 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfrchebpii.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
105 *
106 */
107
108
109// headers du C
110#include <cstdlib>
111#include <fftw3.h>
112
113//Lorene prototypes
114#include "tbl.h"
115
116// Prototypage des sous-routines utilisees:
117namespace Lorene {
118fftw_plan prepare_fft(int, Tbl*&) ;
119double* cheb_ini(const int) ;
120double* chebimp_ini(const int ) ;
121
122//*****************************************************************************
123
124void cfrchebpii(const int* deg, const int* dimf, double* ff, const int* dimc,
125 double* cf)
126
127{
128
129int i, j, k ;
130
131// Dimensions des tableaux ff et cf :
132 int n1f = dimf[0] ;
133 int n2f = dimf[1] ;
134 int n3f = dimf[2] ;
135 int n1c = dimc[0] ;
136 int n2c = dimc[1] ;
137 int n3c = dimc[2] ;
138
139// Nombres de degres de liberte en r :
140 int nr = deg[2] ;
141
142// Tests de dimension:
143 if (nr > n3f) {
144 cout << "cfrchebpii: nr > n3f : nr = " << nr << " , n3f = "
145 << n3f << endl ;
146 abort () ;
147 exit(-1) ;
148 }
149 if (nr > n3c) {
150 cout << "cfrchebpii: nr > n3c : nr = " << nr << " , n3c = "
151 << n3c << endl ;
152 abort () ;
153 exit(-1) ;
154 }
155 if (n1f > n1c) {
156 cout << "cfrchebpii: n1f > n1c : n1f = " << n1f << " , n1c = "
157 << n1c << endl ;
158 abort () ;
159 exit(-1) ;
160 }
161 if (n2f > n2c) {
162 cout << "cfrchebpii: n2f > n2c : n2f = " << n2f << " , n2c = "
163 << n2c << endl ;
164 abort () ;
165 exit(-1) ;
166 }
167
168// Nombre de points pour la FFT:
169 int nm1 = nr - 1;
170 int nm1s2 = nm1 / 2;
171
172// Recherche des tables pour la FFT:
173 Tbl* pg = 0x0 ;
174 fftw_plan p = prepare_fft(nm1, pg) ;
175 Tbl& g = *pg ;
176
177// Recherche de la table des sin(psi) :
178 double* sinp = cheb_ini(nr);
179
180// Recherche de la table des points de collocations x_k :
181 double* x = chebimp_ini(nr);
182
183// boucle sur phi et theta
184
185 int n2n3f = n2f * n3f ;
186 int n2n3c = n2c * n3c ;
187
188/*
189 * Borne de la boucle sur phi:
190 * si n1f = 1, on effectue la boucle une fois seulement.
191 * si n1f > 1, on va jusqu'a j = n1f-2 en sautant j = 1 (les coefficients
192 * j=n1f-1 et j=0 ne sont pas consideres car nuls).
193 */
194 int borne_phi = ( n1f > 1 ) ? n1f-1 : 1 ;
195
196 for (j=0; j< borne_phi; j++) {
197
198 if (j==1) continue ; // on ne traite pas le terme en sin(0 phi)
199
200 /************** Cas l impair ***************/
201
202 for (k=1; k<n2f; k+=2) {
203
204 int i0 = n2n3f * j + n3f * k ; // indice de depart
205 double* ff0 = ff + i0 ; // tableau des donnees a transformer
206
207 i0 = n2n3c * j + n3c * k ; // indice de depart
208 double* cf0 = cf + i0 ; // tableau resultat
209
210
211// Valeur en theta=0 de la partie antisymetrique de F, F_ :
212 double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
213
214// Fonction G(theta) = F+(theta) + F_(theta) sin(theta)
215//---------------------------------------------
216 for ( i = 1; i < nm1s2 ; i++ ) {
217// ... indice (dans le tableau g) du pt symetrique de theta par rapport a pi/2:
218 int isym = nm1 - i ;
219// ... indice (dans le tableau ff0) du point x correspondant a theta
220 int ix = nm1 - i ;
221// ... indice (dans le tableau ff0) du point x correspondant a sym(theta)
222 int ixsym = nm1 - isym ;
223
224// ... F+(psi)
225 double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
226// ... F_(psi) sin(psi)
227 double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
228 g.set(i) = fp + fms ;
229 g.set(isym) = fp - fms ;
230 }
231//... cas particuliers:
232 g.set(0) = 0.5 * ( cf0[nm1] + cf0[0] );
233 g.set(nm1s2) = cf0[nm1s2];
234
235// Developpement de G en series de Fourier par une FFT
236//----------------------------------------------------
237
238 fftw_execute(p) ;
239
240// Coefficients pairs du developmt. de Tchebyshev de h
241//----------------------------------------------------
242// Ces coefficients sont egaux aux coefficients en cosinus du developpement
243// de G en series de Fourier (le facteur 2/nm1 vient de la normalisation
244// de fftw; si fftw donnait reellement les coef. en cosinus, il faudrait le
245// remplacer par un +1.) :
246
247 double fac = 2./double(nm1) ;
248 cf0[0] = g(0) / double(nm1) ;
249 for (i=2; i<nm1; i += 2 ) cf0[i] = fac * g(i/2) ;
250 cf0[nm1] = g(nm1s2) / double(nm1) ;
251
252// Coefficients impairs du developmt. de Tchebyshev de h
253//------------------------------------------------------
254// 1. Coef. c'_k (recurrence amorcee a partir de zero)
255// Le 4/nm1 en facteur de g(i) est du a la normalisation de fftw
256// (si fftw donnait reellement les coef. en sinus, il faudrait le
257// remplacer par un -2.)
258 fac *= 2 ;
259 cf0[1] = 0 ;
260 double som = 0;
261 for ( i = 3; i < nr; i += 2 ) {
262 cf0[i] = cf0[i-2] + fac * g(nm1 - i/2) ;
263 som += cf0[i] ;
264 }
265
266// 2. Calcul de c_1 :
267 double c1 = ( fmoins0 - som ) / nm1s2 ;
268
269// 3. Coef. c_k avec k impair:
270 cf0[1] = c1 ;
271 for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
272
273
274 } // fin de la boucle sur theta
275
276 /************** Cas l pair ***************/
277
278 for (k=0; k<n2f; k+=2) {
279 int i0 = n2n3f * j + n3f * k ; // indice de depart
280 double* ff0 = ff + i0 ; // tableau des donnees a transformer
281
282 i0 = n2n3c * j + n3c * k ; // indice de depart
283 double* cf0 = cf + i0 ; // tableau resultat
284
285// Multiplication de la fonction par x (pour la rendre paire)
286// NB: dans les commentaires qui suivent, on note h(x) = x f(x).
287// (Les valeurs de h dans l'espace des configurations sont stokees dans le
288// tableau cf0).
289 cf0[0] = 0 ;
290 for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
291
292// Valeur en psi=0 de la partie antisymetrique de F, F_ :
293 double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
294
295// Fonction G(theta) = F+(theta) + F_(theta) sin(theta)
296//---------------------------------------------
297 for ( i = 1; i < nm1s2 ; i++ ) {
298// ... indice (dans le tableau g) du pt symetrique de theta par rapport a pi/2:
299 int isym = nm1 - i ;
300// ... indice (dans le tableau cf0) du point x correspondant a theta
301 int ix = nm1 - i ;
302// ... indice (dans le tableau cf0) du point x correspondant a sym(theta)
303 int ixsym = nm1 - isym ;
304
305// ... F+(psi)
306 double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
307// ... F_(psi) sin(psi)
308 double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
309 g.set(i) = fp + fms ;
310 g.set(isym) = fp - fms ;
311 }
312//... cas particuliers:
313 g.set(0) = 0.5 * ( ff0[nm1] + ff0[0] );
314 g.set(nm1s2) = ff0[nm1s2];
315
316// Developpement de G en series de Fourier par une FFT
317//----------------------------------------------------
318
319 fftw_execute(p) ;
320
321// Coefficients pairs du developmt. de Tchebyshev de f
322//----------------------------------------------------
323// Ces coefficients sont egaux aux coefficients en cosinus du developpement
324// de G en series de Fourier (le facteur 2/nm1 vient de la normalisation
325// de fftw) :
326
327 double fac = 2./double(nm1) ;
328 cf0[0] = g(0) / double(nm1) ;
329 for (i=2; i<nm1; i += 2) cf0[i] = fac*g(i/2) ;
330 cf0[nm1] = g(nm1s2) / double(nm1) ;
331
332// Coefficients impairs du developmt. de Tchebyshev de f
333//------------------------------------------------------
334// 1. Coef. c'_k (recurrence amorcee a partir de zero)
335// Le 4/nm1 en facteur de g(i) est du a la normalisation de fftw
336// (si fftw donnait reellement les coef. en sinus, il faudrait le
337// remplacer par un -2.)
338 fac *= 2. ;
339 cf0[1] = 0. ;
340 double som = 0.;
341 for ( i = 3; i < nr; i += 2 ) {
342 cf0[i] = cf0[i-2] + fac * g(nm1-i/2) ;
343 som += cf0[i] ;
344 }
345
346// 2. Calcul de c_1 :
347 double c1 = ( fmoins0 - som ) / nm1s2 ;
348
349// 3. Coef. c_k avec k impair:
350 cf0[1] = c1 ;
351 for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
352
353// Coefficients de f en fonction de ceux de h
354//-------------------------------------------
355
356 cf0[0] = 2* cf0[0] ;
357 for (i=1; i<nm1; i++) {
358 cf0[i] = 2 * cf0[i] - cf0[i-1] ;
359 }
360 cf0[nm1] = 0 ;
361
362
363 } // fin de la boucle sur theta
364 } // fin de la boucle sur phi
365
366}
367}
368
Lorene prototypes.
Definition app_hor.h:64