LORENE
cfrcheb_interp.C
1 /*
2  * Copyright (c) 2024 Jerome Novak
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 
23 
24 
25 /*
26  * Chebyshev transforms (cases R_CHEB, R_CHEBP & R_CHEBI)
27  * using the fftw3 library in the "intermediate" case
28  * where there is no coefficicent transform in theta nor in phi.
29  *
30  * Entree:
31  * -------
32  * int* deg : tableau du nombre effectif de degres de liberte dans chacune
33  * des 3 dimensions: le nombre de points de collocation
34  * en r est nr = deg[2] et doit etre de la forme
35  * nr = 2*p + 1
36  * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
37  * dimensions.
38  * On doit avoir dimf[2] >= deg[2] = nr.
39  *
40  * double* ff : tableau des valeurs de la fonction aux nr points de
41  * de collocation
42  *
43  * x_i = - cos( pi i/(nr-1) ) 0 <= i <= nr-1
44  *
45  * Les valeurs de la fonction doivent etre stokees dans le
46  * tableau ff comme suit
47  * f( x_i ) = ff[ dimf[1]*dimf[2] * j + dimf[2] * k + i ]
48  * ou j et k sont les indices correspondant a phi et theta
49  * respectivement.
50  * L'espace memoire correspondant a ce pointeur doit etre
51  * dimf[0]*dimf[1]*dimf[2] et doit avoir ete alloue avant
52  * l'appel a la routine.
53  *
54  * int* dimc : tableau du nombre d'elements de cf dans chacune des trois
55  * dimensions.
56  * On doit avoir dimc[2] >= deg[2] = nr.
57  *
58  * Sortie:
59  * -------
60  * double* cf : tableau des coefficients c_i de la fonction definis
61  * comme suit (a theta et phi fixes)
62  *
63  * f(x) = som_{i=0}^{nr-1} c_i T_i(x) ,
64  *
65  * ou T_i(x) designe le polynome de Tchebyshev de degre i.
66  * Les coefficients c_i (0 <= i <= nr-1) sont stokes dans
67  * le tableau cf comme suit
68  * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
69  * ou j et k sont les indices correspondant a phi et theta
70  * respectivement.
71  * L'espace memoire correspondant au pointeur cf doit etre
72  * dimc[0]*dimc[1]*dimc[2] et doit avoir ete alloue avant
73  * l'appel a la routine.
74  *
75  * NB: Si le pointeur ff est egal a cf, la routine ne travaille que sur un
76  * seul tableau, qui constitue une entree/sortie.
77  */
78 
79 /*
80  * $Id: cfrcheb_interp.C,v 1.1 2024/09/24 14:41:38 j_novak Exp $
81  *
82  *
83  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfrcheb_interp.C,v 1.1 2024/09/24 14:41:38 j_novak Exp $
84  *
85  */
86 
87 
88 // headers du C
89 #include <cstdlib>
90 #include <cassert>
91 #include <fftw3.h>
92 
93 //Lorene prototypes
94 #include "tbl.h"
95 
96 // Prototypage des sous-routines utilisees:
97 namespace Lorene {
98  fftw_plan prepare_fft(int, Tbl*&) ;
99  double* cheb_ini(const int) ;
100  double* chebimp_ini(const int ) ;
101 
102  //*****************************************************************************
103 
104  void cfrcheb_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
105  double* cf)
106 
107  {
108  int i, j, k ;
109 
110 // Dimensions des tableaux ff et cf :
111  int n1f = dimf[0] ;
112  int n2f = dimf[1] ;
113  int n3f = dimf[2] ;
114  int n1c = dimc[0] ;
115  int n2c = dimc[1] ;
116  int n3c = dimc[2] ;
117 
118 // Nombres de degres de liberte en r :
119  int nr = deg[2] ;
120 
121 // Tests de dimension:
122  if (nr > n3f) {
123  cout << "cfrcheb: nr > n3f : nr = " << nr << " , n3f = "
124  << n3f << endl ;
125  abort () ;
126  exit(-1) ;
127  }
128  if (nr > n3c) {
129  cout << "cfrcheb: nr > n3c : nr = " << nr << " , n3c = "
130  << n3c << endl ;
131  abort () ;
132  exit(-1) ;
133  }
134  if (n1f > n1c) {
135  cout << "cfrcheb: n1f > n1c : n1f = " << n1f << " , n1c = "
136  << n1c << endl ;
137  abort () ;
138  exit(-1) ;
139  }
140  if (n2f > n2c) {
141  cout << "cfrcheb: n2f > n2c : n2f = " << n2f << " , n2c = "
142  << n2c << endl ;
143  abort () ;
144  exit(-1) ;
145  }
146 
147 
148 // Nombre de points pour la FFT:
149  int nm1 = nr - 1;
150  int nm1s2 = nm1 / 2;
151 
152 // Recherche des tables pour la FFT:
153  Tbl* pg = 0x0 ;
154  fftw_plan p = prepare_fft(nm1, pg) ;
155  Tbl& g = *pg ;
156 
157 // Recherche de la table des sin(psi) :
158  double* sinp = cheb_ini(nr);
159 
160 // boucle sur phi et theta
161 
162  int n2n3f = n2f * n3f ;
163  int n2n3c = n2c * n3c ;
164 
165  for (j=0; j< n1f; j++) {
166 
167  for (k=0; k<n2f; k++) {
168 
169  int i0 = n2n3f * j + n3f * k ; // indice de depart
170  double* ff0 = ff + i0 ; // tableau des donnees a transformer
171 
172  i0 = n2n3c * j + n3c * k ; // indice de depart
173  double* cf0 = cf + i0 ; // tableau resultat
174 
175 /*
176  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
177  * reliee a x par x = - cos(psi) et F(psi) = f(x(psi)).
178  */
179 
180 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
181  double fmoins0 = 0.5 * ( ff0[0] - ff0[nm1] );
182 
183 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
184 //---------------------------------------------
185  for ( i = 1; i < nm1s2 ; i++ ) {
186 // ... indice du pt symetrique de psi par rapport a pi/2:
187  int isym = nm1 - i ;
188 // ... F+(psi)
189  double fp = 0.5 * ( ff0[i] + ff0[isym] ) ;
190 // ... F_(psi) sin(psi)
191  double fms = 0.5 * ( ff0[i] - ff0[isym] ) * sinp[i] ;
192  g.set(i) = fp + fms ;
193  g.set(isym) = fp - fms ;
194  }
195 //... cas particuliers:
196  g.set(0) = 0.5 * ( ff0[0] + ff0[nm1] );
197  g.set(nm1s2) = ff0[nm1s2];
198 
199 // Developpement de G en series de Fourier par une FFT
200 //----------------------------------------------------
201 
202  fftw_execute(p) ;
203 
204 // Coefficients pairs du developmt. de Tchebyshev de f
205 //----------------------------------------------------
206 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
207 // de G en series de Fourier (le facteur 2/nm1 vient de la normalisation
208 // de fftw) :
209 
210  double fac = 2./double(nm1) ;
211  cf0[0] = g(0) / double(nm1) ;
212  for (i=2; i<nm1; i += 2) cf0[i] = fac*g(i/2) ;
213  cf0[nm1] = g(nm1s2) / double(nm1) ;
214 
215 // Coefficients impairs du developmt. de Tchebyshev de f
216 //------------------------------------------------------
217 // 1. Coef. c'_k (recurrence amorcee a partir de zero):
218 // NB: Le 4/nm1 en facteur de g(i) est du a la normalisation de fftw
219 // (si fftw donnait reellement les coef. en sinus, il faudrait le
220 // remplacer par un +2.)
221  fac *= -2. ;
222  cf0[1] = 0 ;
223  double som = 0;
224  for ( i = 3; i < nr; i += 2 ) {
225  cf0[i] = cf0[i-2] + fac * g(nm1-i/2) ;
226  som += cf0[i] ;
227  }
228 
229 // 2. Calcul de c_1 :
230  double c1 = - ( fmoins0 + som ) / nm1s2 ;
231 
232 // 3. Coef. c_k avec k impair:
233  cf0[1] = c1 ;
234  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
235 
236  } // fin de la boucle sur theta
237  } // fin de la boucle sur phi
238 
239  } // end of cfrcheb_interp
240 
241 //*****************************************************************************
242 
243  void cfrchebp_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
244  double* cf)
245 
246  {
247  int i, j, k ;
248 
249 // Dimensions des tableaux ff et cf :
250  int n1f = dimf[0] ;
251  int n2f = dimf[1] ;
252  int n3f = dimf[2] ;
253  int n1c = dimc[0] ;
254  int n2c = dimc[1] ;
255  int n3c = dimc[2] ;
256 
257 // Nombres de degres de liberte en r :
258  int nr = deg[2] ;
259 
260 // Tests de dimension:
261  if (nr > n3f) {
262  cout << "cfrchebp: nr > n3f : nr = " << nr << " , n3f = "
263  << n3f << endl ;
264  abort () ;
265  exit(-1) ;
266  }
267  if (nr > n3c) {
268  cout << "cfrchebp: nr > n3c : nr = " << nr << " , n3c = "
269  << n3c << endl ;
270  abort () ;
271  exit(-1) ;
272  }
273  if (n1f > n1c) {
274  cout << "cfrchebp: n1f > n1c : n1f = " << n1f << " , n1c = "
275  << n1c << endl ;
276  abort () ;
277  exit(-1) ;
278  }
279  if (n2f > n2c) {
280  cout << "cfrchebp: n2f > n2c : n2f = " << n2f << " , n2c = "
281  << n2c << endl ;
282  abort () ;
283  exit(-1) ;
284  }
285 
286 // Nombre de points pour la FFT:
287  int nm1 = nr - 1;
288  int nm1s2 = nm1 / 2;
289 
290 // Recherche des tables pour la FFT:
291  Tbl* pg = 0x0 ;
292  fftw_plan p = prepare_fft(nm1, pg) ;
293  Tbl& g = *pg ;
294 
295 // Recherche de la table des sin(psi) :
296  double* sinp = cheb_ini(nr);
297 
298 // boucle sur phi et theta
299 
300  int n2n3f = n2f * n3f ;
301  int n2n3c = n2c * n3c ;
302 
303 
304  for (j=0; j< n1f; j++) {
305 
306  for (k=0; k<n2f; k++) {
307 
308  int i0 = n2n3f * j + n3f * k ; // indice de depart
309  double* ff0 = ff + i0 ; // tableau des donnees a transformer
310 
311  i0 = n2n3c * j + n3c * k ; // indice de depart
312  double* cf0 = cf + i0 ; // tableau resultat
313 
314 /*
315  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
316  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
317  */
318 
319 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
320  double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
321 
322 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
323 //---------------------------------------------
324  for ( i = 1; i < nm1s2 ; i++ ) {
325 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
326  int isym = nm1 - i ;
327 // ... indice (dans le tableau ff0) du point x correspondant a psi
328  int ix = nm1 - i ;
329 // ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
330  int ixsym = nm1 - isym ;
331 
332 // ... F+(psi)
333  double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
334 // ... F_(psi) sin(psi)
335  double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
336  g.set(i) = fp + fms ;
337  g.set(isym) = fp - fms ;
338  }
339 //... cas particuliers:
340  g.set(0) = 0.5 * ( ff0[nm1] + ff0[0] );
341  g.set(nm1s2) = ff0[nm1s2];
342 
343 // Developpement de G en series de Fourier par une FFT
344 //----------------------------------------------------
345 
346  fftw_execute(p) ;
347 
348 // Coefficients pairs du developmt. de Tchebyshev de f
349 //----------------------------------------------------
350 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
351 // de G en series de Fourier (le facteur 2/nm1 vient de la normalisation
352 // de fftw) :
353 
354  double fac = 2./double(nm1) ;
355  cf0[0] = g(0) / double(nm1) ;
356  for (i=2; i<nm1; i += 2) cf0[i] = fac*g(i/2) ;
357  cf0[nm1] = g(nm1s2) / double(nm1) ;
358 
359 // Coefficients impairs du developmt. de Tchebyshev de f
360 //------------------------------------------------------
361 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
362 // Le 4/nm1 en facteur de g(i) est du a la normalisation de fftw
363 // (si fftw donnait reellement les coef. en sinus, il faudrait le
364 // remplacer par un -2.)
365  fac *= 2. ;
366  cf0[1] = 0. ;
367  double som = 0.;
368  for ( i = 3; i < nr; i += 2 ) {
369  cf0[i] = cf0[i-2] + fac * g(nm1-i/2) ;
370  som += cf0[i] ;
371  }
372 
373 // 2. Calcul de c_1 :
374  double c1 = ( fmoins0 - som ) / nm1s2 ;
375 
376 // 3. Coef. c_k avec k impair:
377  cf0[1] = c1 ;
378  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
379 
380  } // fin de la boucle sur theta
381  } // fin de la boucle sur phi
382 
383 
384  } // end of cfrchebp_interp
385 
386  //*****************************************************************************
387 
388  void cfrchebi_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
389  double* cf)
390 
391  {
392  int i, j, k ;
393 
394 // Dimensions des tableaux ff et cf :
395  int n1f = dimf[0] ;
396  int n2f = dimf[1] ;
397  int n3f = dimf[2] ;
398  int n1c = dimc[0] ;
399  int n2c = dimc[1] ;
400  int n3c = dimc[2] ;
401 
402 // Nombres de degres de liberte en theta et r :
403  int nr = deg[2] ;
404 
405 // Tests de dimension:
406  if (nr > n3f) {
407  cout << "cfrchebi: nr > n3f : nr = " << nr << " , n3f = "
408  << n3f << endl ;
409  abort () ;
410  exit(-1) ;
411  }
412  if (nr > n3c) {
413  cout << "cfrchebi: nr > n3c : nr = " << nr << " , n3c = "
414  << n3c << endl ;
415  abort () ;
416  exit(-1) ;
417  }
418  if (n1f > n1c) {
419  cout << "cfrchebi: n1f > n1c : n1f = " << n1f << " , n1c = "
420  << n1c << endl ;
421  abort () ;
422  exit(-1) ;
423  }
424  if (n2f > n2c) {
425  cout << "cfrchebi: n2f > n2c : n2f = " << n2f << " , n2c = "
426  << n2c << endl ;
427  abort () ;
428  exit(-1) ;
429  }
430 
431 // Nombre de points pour la FFT:
432  int nm1 = nr - 1;
433  int nm1s2 = nm1 / 2;
434 
435 // Recherche des tables pour la FFT:
436  Tbl* pg = 0x0 ;
437  fftw_plan p = prepare_fft(nm1, pg) ;
438  Tbl& g = *pg ;
439 
440 // Recherche de la table des sin(psi) :
441  double* sinp = cheb_ini(nr);
442 
443 // Recherche de la table des points de collocations x_k :
444  double* x = chebimp_ini(nr);
445 
446 // boucle sur phi et theta
447 
448  int n2n3f = n2f * n3f ;
449  int n2n3c = n2c * n3c ;
450 
451  for (j=0; j< n1f; j++) {
452 
453  for (k=0; k<n2f; k++) {
454 
455  int i0 = n2n3f * j + n3f * k ; // indice de depart
456  double* ff0 = ff + i0 ; // tableau des donnees a transformer
457 
458  i0 = n2n3c * j + n3c * k ; // indice de depart
459  double* cf0 = cf + i0 ; // tableau resultat
460 
461 // Multiplication de la fonction par x (pour la rendre paire)
462 // NB: dans les commentaires qui suivent, on note h(x) = x f(x).
463 // (Les valeurs de h dans l'espace des configurations sont stokees dans le
464 // tableau cf0).
465  cf0[0] = 0 ;
466  for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
467 
468 /*
469  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
470  * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
471  */
472 
473 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
474  double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
475 
476 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
477 //---------------------------------------------
478  for ( i = 1; i < nm1s2 ; i++ ) {
479 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
480  int isym = nm1 - i ;
481 // ... indice (dans le tableau cf0) du point x correspondant a psi
482  int ix = nm1 - i ;
483 // ... indice (dans le tableau cf0) du point x correspondant a sym(psi)
484  int ixsym = nm1 - isym ;
485 
486 // ... F+(psi)
487  double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
488 // ... F_(psi) sin(psi)
489  double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
490  g.set(i) = fp + fms ;
491  g.set(isym) = fp - fms ;
492  }
493 //... cas particuliers:
494  g.set(0) = 0.5 * ( cf0[nm1] + cf0[0] );
495  g.set(nm1s2) = cf0[nm1s2];
496 
497 // Developpement de G en series de Fourier par une FFT
498 //----------------------------------------------------
499 
500  fftw_execute(p) ;
501 
502 // Coefficients pairs du developmt. de Tchebyshev de h
503 //----------------------------------------------------
504 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
505 // de G en series de Fourier (le facteur 2/nm1 vient de la normalisation
506 // de fftw; si fftw donnait reellement les coef. en cosinus, il faudrait le
507 // remplacer par un +1.) :
508 
509  double fac = 2./double(nm1) ;
510  cf0[0] = g(0) / double(nm1) ;
511  for (i=2; i<nm1; i += 2 ) cf0[i] = fac * g(i/2) ;
512  cf0[nm1] = g(nm1s2) / double(nm1) ;
513 
514 // Coefficients impairs du developmt. de Tchebyshev de h
515 //------------------------------------------------------
516 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
517 // Le 4/nm1 en facteur de g(i) est du a la normalisation de fftw
518 // (si fftw donnait reellement les coef. en sinus, il faudrait le
519 // remplacer par un -2.)
520  fac *= 2 ;
521  cf0[1] = 0 ;
522  double som = 0;
523  for ( i = 3; i < nr; i += 2 ) {
524  cf0[i] = cf0[i-2] + fac * g(nm1 - i/2) ;
525  som += cf0[i] ;
526  }
527 
528 // 2. Calcul de c_1 :
529  double c1 = ( fmoins0 - som ) / nm1s2 ;
530 
531 // 3. Coef. c_k avec k impair:
532  cf0[1] = c1 ;
533  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
534 
535 // Coefficients de f en fonction de ceux de h
536 //-------------------------------------------
537 
538  cf0[0] = 2* cf0[0] ;
539  for (i=1; i<nm1; i++) {
540  cf0[i] = 2 * cf0[i] - cf0[i-1] ;
541  }
542  cf0[nm1] = 0 ;
543 
544 
545  } // fin de la boucle sur theta
546  } // fin de la boucle sur phi
547 
548  } // end of cfrchebi_interp
549 
550 
551 } // end of namespace Lorene
552 
Lorene
Lorene prototypes.
Definition: app_hor.h:67