LORENE
op_mult_x.C
1  /*
2  * Copyright (c) 1999-2001 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  * $Id: op_mult_x.C,v 1.8 2023/05/24 09:53:46 g_servignat Exp $
27  * $Log: op_mult_x.C,v $
28  * Revision 1.8 2023/05/24 09:53:46 g_servignat
29  * Added multiplication by \xi for R_CHEB base
30  *
31  * Revision 1.7 2016/12/05 16:18:08 j_novak
32  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
33  *
34  * Revision 1.6 2015/03/05 08:49:32 j_novak
35  * Implemented operators with Legendre bases.
36  *
37  * Revision 1.5 2014/10/13 08:53:25 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.4 2008/08/27 11:22:25 j_novak
41  * Minor modifications
42  *
43  * Revision 1.3 2008/08/27 08:50:10 jl_cornou
44  * Added Jacobi(0,2) polynomials
45  *
46  * Revision 1.2 2004/11/23 15:16:01 m_forot
47  *
48  * Added the bases for the cases without any equatorial symmetry
49  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
50  *
51  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
52  * LORENE
53  *
54  * Revision 1.3 2000/09/07 12:49:53 phil
55  * *** empty log message ***
56  *
57  * Revision 1.2 2000/02/24 16:42:18 eric
58  * Initialisation a zero du tableau xo avant le calcul.
59  *
60  * Revision 1.1 1999/11/16 13:37:41 novak
61  * Initial revision
62  *
63  *
64  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_x.C,v 1.8 2023/05/24 09:53:46 g_servignat Exp $
65  *
66  */
67 
68  /*
69  * Ensemble des routines de base de multiplication par x
70  * (Utilisation interne)
71  *
72  * void _mult_x_XXXX(Tbl * t, int & b)
73  * t pointeur sur le Tbl d'entree-sortie
74  * b base spectrale
75  *
76  */
77 
78  // Fichier includes
79  #include "tbl.h"
80 
81 
82  //-----------------------------------
83  // Routine pour les cas non prevus --
84  //-----------------------------------
85 
86  namespace Lorene {
87  void _mult_x_pas_prevu(Tbl * tb, int& base) {
88  cout << "mult_x pas prevu..." << endl ;
89  cout << "Tbl: " << tb << " base: " << base << endl ;
90  abort () ;
91  exit(-1) ;
92  }
93 
94  //-------------
95  // Identite ---
96  //-------------
97 
98  void _mult_x_identite(Tbl* , int& ) {
99  return ;
100  }
101 
102  //---------------
103  // cas R_CHEBP --
104  //--------------
105 
106  void _mult_x_r_chebp(Tbl* tb, int& base)
107  {
108  // Peut-etre rien a faire ?
109  if (tb->get_etat() == ETATZERO) {
110  int base_t = base & MSQ_T ;
111  int base_p = base & MSQ_P ;
112  base = base_p | base_t | R_CHEBI ;
113  return ;
114  }
115 
116  // Pour le confort
117  int nr = (tb->dim).dim[0] ; // Nombre
118  int nt = (tb->dim).dim[1] ; // de points
119  int np = (tb->dim).dim[2] ; // physiques REELS
120  np = np - 2 ; // Nombre de points physiques
121 
122  // pt. sur le tableau de double resultat
123  double* xo = new double [tb->get_taille()];
124 
125  // Initialisation a zero :
126  for (int i=0; i<tb->get_taille(); i++) {
127  xo[i] = 0 ;
128  }
129 
130  // On y va...
131  double* xi = tb->t ;
132  double* xci = xi ; // Pointeurs
133  double* xco = xo ; // courants
134 
135  int borne_phi = np + 1 ;
136  if (np == 1) {
137  borne_phi = 1 ;
138  }
139 
140  for (int k=0 ; k< borne_phi ; k++)
141  if (k==1) {
142  xci += nr*nt ;
143  xco += nr*nt ;
144  }
145  else {
146  for (int j=0 ; j<nt ; j++) {
147 
148  xco[0] = xci[0] + 0.5*xci[1] ;
149  for (int i = 1 ; i < nr-1 ; i++ ) {
150  xco[i] = 0.5*(xci[i]+xci[i+1]) ;
151  } // Fin de la boucle sur r
152  xco[nr-1] = 0 ;
153 
154  xci += nr ;
155  xco += nr ;
156  } // Fin de la boucle sur theta
157  } // Fin de la boucle sur phi
158 
159  // On remet les choses la ou il faut
160  delete [] tb->t ;
161  tb->t = xo ;
162 
163  // base de developpement
164  // pair -> impair
165  int base_t = base & MSQ_T ;
166  int base_p = base & MSQ_P ;
167  base = base_p | base_t | R_CHEBI ;
168 
169  }
170 
171  //----------------
172  // cas R_CHEBI ---
173  //----------------
174 
175  void _mult_x_r_chebi(Tbl* tb, int& base)
176  {
177 
178  // Peut-etre rien a faire ?
179  if (tb->get_etat() == ETATZERO) {
180  int base_t = base & MSQ_T ;
181  int base_p = base & MSQ_P ;
182  base = base_p | base_t | R_CHEBP ;
183  return ;
184  }
185 
186  // Pour le confort
187  int nr = (tb->dim).dim[0] ; // Nombre
188  int nt = (tb->dim).dim[1] ; // de points
189  int np = (tb->dim).dim[2] ; // physiques REELS
190  np = np - 2 ; // Nombre de points physiques
191 
192  // pt. sur le tableau de double resultat
193  double* xo = new double [tb->get_taille()];
194 
195  // Initialisation a zero :
196  for (int i=0; i<tb->get_taille(); i++) {
197  xo[i] = 0 ;
198  }
199 
200  // On y va...
201  double* xi = tb->t ;
202  double* xci = xi ; // Pointeurs
203  double* xco = xo ; // courants
204 
205  int borne_phi = np + 1 ;
206  if (np == 1) {
207  borne_phi = 1 ;
208  }
209 
210  for (int k=0 ; k< borne_phi ; k++)
211  if (k == 1) {
212  xci += nr*nt ;
213  xco += nr*nt ;
214  }
215  else {
216  for (int j=0 ; j<nt ; j++) {
217 
218  xco[0] = 0.5*xci[0] ;
219  for (int i = 1 ; i < nr-1 ; i++ ) {
220  xco[i] = (xci[i]+xci[i-1])*0.5 ;
221  } // Fin de la premiere boucle sur r
222  xco[nr-1] = 0.5*xci[nr-2] ;
223 
224  xci += nr ;
225  xco += nr ;
226  } // Fin de la boucle sur theta
227  } // Fin de la boucle sur phi
228 
229  // On remet les choses la ou il faut
230  delete [] tb->t ;
231  tb->t = xo ;
232 
233  // base de developpement
234  // impair -> pair
235  int base_t = base & MSQ_T ;
236  int base_p = base & MSQ_P ;
237  base = base_p | base_t | R_CHEBP ;
238  }
239 
240  //--------------------
241  // cas R_CHEBPIM_P --
242  //------------------
243 
244  void _mult_x_r_chebpim_p(Tbl* tb, int& base)
245  {
246 
247  // Peut-etre rien a faire ?
248  if (tb->get_etat() == ETATZERO) {
249  int base_t = base & MSQ_T ;
250  int base_p = base & MSQ_P ;
251  base = base_p | base_t | R_CHEBPIM_I ;
252  return ;
253  }
254 
255  // Pour le confort
256  int nr = (tb->dim).dim[0] ; // Nombre
257  int nt = (tb->dim).dim[1] ; // de points
258  int np = (tb->dim).dim[2] ; // physiques REELS
259  np = np - 2 ;
260 
261  // pt. sur le tableau de double resultat
262  double* xo = new double [tb->get_taille()];
263 
264  // Initialisation a zero :
265  for (int i=0; i<tb->get_taille(); i++) {
266  xo[i] = 0 ;
267  }
268 
269 
270  // On y va...
271  double* xi = tb->t ;
272  double* xci ; // Pointeurs
273  double* xco ; // courants
274 
275 
276  // Partie paire
277  xci = xi ;
278  xco = xo ;
279 
280  int auxiliaire ;
281 
282  for (int k=0 ; k<np+1 ; k += 4) {
283 
284  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
285  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
286 
287 
288  // On evite le sin (0*phi)
289 
290  if ((k==0) && (kmod == 1)) {
291  xco += nr*nt ;
292  xci += nr*nt ;
293  }
294 
295  else
296  for (int j=0 ; j<nt ; j++) {
297 
298  xco[0] = xci[0] + 0.5*xci[1] ;
299  for (int i = 1 ; i < nr-1 ; i++ ) {
300  xco[i] = 0.5*(xci[i]+xci[i+1]) ;
301  } // Fin de la boucle sur r
302  xco[nr-1] = 0 ;
303 
304  xci += nr ; // au
305  xco += nr ; // suivant
306  } // Fin de la boucle sur theta
307  } // Fin de la boucle sur kmod
308  xci += 2*nr*nt ; // saute
309  xco += 2*nr*nt ; // 2 phis
310  } // Fin de la boucle sur phi
311 
312  // Partie impaire
313  xci = xi + 2*nr*nt ;
314  xco = xo + 2*nr*nt ;
315  for (int k=2 ; k<np+1 ; k += 4) {
316 
317  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
318  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
319  for (int j=0 ; j<nt ; j++) {
320 
321  xco[0] = 0.5*xci[0] ;
322  for (int i = 1 ; i < nr-1 ; i++ ) {
323  xco[i] = (xci[i]+xci[i-1])*0.5 ;
324  } // Fin de la premiere boucle sur r
325  xco[nr-1] = 0.5*xci[nr-2] ;
326 
327  xci += nr ;
328  xco += nr ;
329  } // Fin de la boucle sur theta
330  } // Fin de la boucle sur kmod
331  xci += 2*nr*nt ;
332  xco += 2*nr*nt ;
333  } // Fin de la boucle sur phi
334 
335  // On remet les choses la ou il faut
336  delete [] tb->t ;
337  tb->t = xo ;
338 
339  // base de developpement
340  // (pair,impair) -> (impair,pair)
341  int base_t = base & MSQ_T ;
342  int base_p = base & MSQ_P ;
343  base = base_p | base_t | R_CHEBPIM_I ;
344  }
345 
346  //-------------------
347  // cas R_CHEBPIM_I --
348  //-------------------
349 
350  void _mult_x_r_chebpim_i(Tbl* tb, int& base)
351  {
352 
353  // Peut-etre rien a faire ?
354  if (tb->get_etat() == ETATZERO) {
355  int base_t = base & MSQ_T ;
356  int base_p = base & MSQ_P ;
357  base = base_p | base_t | R_CHEBPIM_P ;
358  return ;
359  }
360 
361  // Pour le confort
362  int nr = (tb->dim).dim[0] ; // Nombre
363  int nt = (tb->dim).dim[1] ; // de points
364  int np = (tb->dim).dim[2] ; // physiques REELS
365  np = np - 2 ;
366 
367  // pt. sur le tableau de double resultat
368  double* xo = new double [tb->get_taille()];
369 
370  // Initialisation a zero :
371  for (int i=0; i<tb->get_taille(); i++) {
372  xo[i] = 0 ;
373  }
374 
375  // On y va...
376  double* xi = tb->t ;
377  double* xci ; // Pointeurs
378  double* xco ; // courants
379 
380  // Partie impaire
381  xci = xi ;
382  xco = xo ;
383 
384 
385  int auxiliaire ;
386 
387  for (int k=0 ; k<np+1 ; k += 4) {
388 
389  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
390  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
391 
392 
393  // On evite le sin (0*phi)
394 
395  if ((k==0) && (kmod == 1)) {
396  xco += nr*nt ;
397  xci += nr*nt ;
398  }
399 
400  else
401  for (int j=0 ; j<nt ; j++) {
402 
403  xco[0] = 0.5*xci[0] ;
404  for (int i = 1 ; i < nr-1 ; i++ ) {
405  xco[i] = (xci[i]+xci[i-1])*0.5 ;
406  } // Fin de la premiere boucle sur r
407  xco[nr-1] = 0.5*xci[nr-2] ;
408 
409  xci += nr ;
410  xco += nr ;
411  } // Fin de la boucle sur theta
412  } // Fin de la boucle sur kmod
413  xci += 2*nr*nt ;
414  xco += 2*nr*nt ;
415  } // Fin de la boucle sur phi
416 
417  // Partie paire
418  xci = xi + 2*nr*nt ;
419  xco = xo + 2*nr*nt ;
420  for (int k=2 ; k<np+1 ; k += 4) {
421 
422  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
423  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
424  for (int j=0 ; j<nt ; j++) {
425 
426  xco[0] = xci[0] + 0.5*xci[1] ;
427  for (int i = 1 ; i < nr-1 ; i++ ) {
428  xco[i] = 0.5*(xci[i]+xci[i+1]) ;
429  } // Fin de la boucle sur r
430  xco[nr-1] = 0 ;
431 
432  xci += nr ;
433  xco += nr ;
434  } // Fin de la boucle sur theta
435  } // Fin de la boucle sur kmod
436  xci += 2*nr*nt ;
437  xco += 2*nr*nt ;
438  } // Fin de la boucle sur phi
439 
440  // On remet les choses la ou il faut
441  delete [] tb->t ;
442  tb->t = xo ;
443 
444  // base de developpement
445  // (impair,pair) -> (pair,impair)
446  int base_t = base & MSQ_T ;
447  int base_p = base & MSQ_P ;
448  base = base_p | base_t | R_CHEBPIM_P ;
449  }
450 
451  //---------------
452  // cas R_CHEBPI_P --
453  //--------------
454 
455  void _mult_x_r_chebpi_p(Tbl* tb, int& base)
456  {
457  // Peut-etre rien a faire ?
458  if (tb->get_etat() == ETATZERO) {
459  int base_t = base & MSQ_T ;
460  int base_p = base & MSQ_P ;
461  base = base_p | base_t | R_CHEBPI_I ;
462  return ;
463  }
464 
465  // Pour le confort
466  int nr = (tb->dim).dim[0] ; // Nombre
467  int nt = (tb->dim).dim[1] ; // de points
468  int np = (tb->dim).dim[2] ; // physiques REELS
469  np = np - 2 ; // Nombre de points physiques
470 
471  // pt. sur le tableau de double resultat
472  double* xo = new double [tb->get_taille()];
473 
474  // Initialisation a zero :
475  for (int i=0; i<tb->get_taille(); i++) {
476  xo[i] = 0 ;
477  }
478 
479  // On y va...
480  double* xi = tb->t ;
481  double* xci = xi ; // Pointeurs
482  double* xco = xo ; // courants
483 
484  int borne_phi = np + 1 ;
485  if (np == 1) {
486  borne_phi = 1 ;
487  }
488 
489  for (int k=0 ; k< borne_phi ; k++)
490  if (k==1) {
491  xci += nr*nt ;
492  xco += nr*nt ;
493  }
494  else {
495  for (int j=0 ; j<nt ; j++) {
496  int l = j%2 ;
497  if(l==0){
498  xco[0] = xci[0] + 0.5*xci[1] ;
499  for (int i = 1 ; i < nr-1 ; i++ ) {
500  xco[i] = 0.5*(xci[i]+xci[i+1]) ;
501  } // Fin de la boucle sur r
502  xco[nr-1] = 0 ;
503  } else {
504  xco[0] = 0.5*xci[0] ;
505  for (int i = 1 ; i < nr-1 ; i++ ) {
506  xco[i] = (xci[i]+xci[i-1])*0.5 ;
507  } // Fin de la premiere boucle sur r
508  xco[nr-1] = 0.5*xci[nr-2] ;
509  }
510  xci += nr ;
511  xco += nr ;
512  } // Fin de la boucle sur theta
513  } // Fin de la boucle sur phi
514 
515  // On remet les choses la ou il faut
516  delete [] tb->t ;
517  tb->t = xo ;
518 
519  // base de developpement
520  // pair -> impair
521  int base_t = base & MSQ_T ;
522  int base_p = base & MSQ_P ;
523  base = base_p | base_t | R_CHEBPI_I ;
524 
525  }
526 
527  //----------------
528  // cas R_CHEBPI_I ---
529  //----------------
530 
531  void _mult_x_r_chebpi_i(Tbl* tb, int& base)
532  {
533 
534  // Peut-etre rien a faire ?
535  if (tb->get_etat() == ETATZERO) {
536  int base_t = base & MSQ_T ;
537  int base_p = base & MSQ_P ;
538  base = base_p | base_t | R_CHEBPI_P ;
539  return ;
540  }
541 
542  // Pour le confort
543  int nr = (tb->dim).dim[0] ; // Nombre
544  int nt = (tb->dim).dim[1] ; // de points
545  int np = (tb->dim).dim[2] ; // physiques REELS
546  np = np - 2 ; // Nombre de points physiques
547 
548  // pt. sur le tableau de double resultat
549  double* xo = new double [tb->get_taille()];
550 
551  // Initialisation a zero :
552  for (int i=0; i<tb->get_taille(); i++) {
553  xo[i] = 0 ;
554  }
555 
556  // On y va...
557  double* xi = tb->t ;
558  double* xci = xi ; // Pointeurs
559  double* xco = xo ; // courants
560 
561  int borne_phi = np + 1 ;
562  if (np == 1) {
563  borne_phi = 1 ;
564  }
565 
566  for (int k=0 ; k< borne_phi ; k++)
567  if (k == 1) {
568  xci += nr*nt ;
569  xco += nr*nt ;
570  }
571  else {
572  for (int j=0 ; j<nt ; j++) {
573  int l = j%2 ;
574  if(l==1){
575  xco[0] = xci[0] + 0.5*xci[1] ;
576  for (int i = 1 ; i < nr-1 ; i++ ) {
577  xco[i] = 0.5*(xci[i]+xci[i+1]) ;
578  } // Fin de la boucle sur r
579  xco[nr-1] = 0 ;
580  } else {
581  xco[0] = 0.5*xci[0] ;
582  for (int i = 1 ; i < nr-1 ; i++ ) {
583  xco[i] = (xci[i]+xci[i-1])*0.5 ;
584  } // Fin de la premiere boucle sur r
585  xco[nr-1] = 0.5*xci[nr-2] ;
586  }
587  xci += nr ;
588  xco += nr ;
589  } // Fin de la boucle sur theta
590  } // Fin de la boucle sur phi
591 
592  // On remet les choses la ou il faut
593  delete [] tb->t ;
594  tb->t = xo ;
595 
596  // base de developpement
597  // impair -> pair
598  int base_t = base & MSQ_T ;
599  int base_p = base & MSQ_P ;
600  base = base_p | base_t | R_CHEBPI_P ;
601  }
602 
603  //---------------
604  // cas R_JACO02 -
605  //---------------
606 
607  void _mult_x_r_jaco02(Tbl* tb, int&)
608  {
609  // Peut-etre rien a faire ?
610  if (tb->get_etat() == ETATZERO) {
611  return ;
612  }
613 
614  // Pour le confort
615  int nr = (tb->dim).dim[0] ; // Nombre
616  int nt = (tb->dim).dim[1] ; // de points
617  int np = (tb->dim).dim[2] ; // physiques REELS
618  np = np - 2 ; // Nombre de points physiques
619 
620  // pt. sur le tableau de double resultat
621  double* xo = new double [tb->get_taille()];
622 
623  // Initialisation a zero :
624  for (int i=0; i<tb->get_taille(); i++) {
625  xo[i] = 0 ;
626  }
627 
628  // On y va...
629  double* xi = tb->t ;
630  double* xci = xi ; // Pointeurs
631  double* xco = xo ; // courants
632 
633  int borne_phi = np + 1 ;
634  if (np == 1) {
635  borne_phi = 1 ;
636  }
637 
638  for (int k=0 ; k< borne_phi ; k++)
639  if (k==1) {
640  xci += nr*nt ;
641  xco += nr*nt ;
642  }
643  else {
644  for (int j=0 ; j<nt ; j++) {
645 
646  xco[0] = 1.5*xci[0] + 0.3*xci[1] ;
647  for (int i = 1 ; i < nr-1 ; i++) {
648  xco[i] = i*(i+2)/double((i+1)*(2*i+1))*xci[i-1] + (i*i+3*i+3)/double((i+1)*(i+2))*xci[i] + (i+1)*(i+3)/double((i+2)*(2*i+5))*xci[i+1] ;
649  }
650  xco[nr-1] = (nr*nr-1)/double((nr)*(2*nr-1))*xci[nr-2] + (1+1/double((nr)*(nr+1)))*xci[nr-1] ;
651 
652  xci += nr ;
653  xco += nr ;
654  } // Fin de la boucle sur theta
655  } // Fin de la boucle sur phi
656 
657  // On remet les choses la ou il faut
658  delete [] tb->t ;
659  tb->t = xo ;
660 
661  // base de developpement
662  // inchangée
663 
664  }
665 
666  //--------------
667  // cas R_LEGP --
668  //--------------
669 
670  void _mult_x_r_legp(Tbl* tb, int& base)
671  {
672  // Peut-etre rien a faire ?
673  if (tb->get_etat() == ETATZERO) {
674  int base_t = base & MSQ_T ;
675  int base_p = base & MSQ_P ;
676  base = base_p | base_t | R_LEGI ;
677  return ;
678  }
679 
680  // Pour le confort
681  int nr = (tb->dim).dim[0] ; // Nombre
682  int nt = (tb->dim).dim[1] ; // de points
683  int np = (tb->dim).dim[2] ; // physiques REELS
684  np = np - 2 ; // Nombre de points physiques
685 
686  // pt. sur le tableau de double resultat
687  double* xo = new double [tb->get_taille()];
688 
689  // Initialisation a zero :
690  for (int i=0; i<tb->get_taille(); i++) {
691  xo[i] = 0 ;
692  }
693 
694  // On y va...
695  double* xi = tb->t ;
696  double* xci = xi ; // Pointeurs
697  double* xco = xo ; // courants
698 
699  int borne_phi = np + 1 ;
700  if (np == 1) {
701  borne_phi = 1 ;
702  }
703 
704  for (int k=0 ; k< borne_phi ; k++)
705  if (k==1) {
706  xci += nr*nt ;
707  xco += nr*nt ;
708  }
709  else {
710  for (int j=0 ; j<nt ; j++) {
711 
712  xco[0] = xci[0] + 0.4*xci[1] ;
713  for (int i = 1 ; i < nr-1 ; i++ ) {
714  xco[i] = double(2*i+1)*xci[i]/double(4*i+1)
715  + double(2*i+2)*xci[i+1]/double(4*i+5) ;
716  } // Fin de la boucle sur r
717  xco[nr-1] = 0 ;
718 
719  xci += nr ;
720  xco += nr ;
721  } // Fin de la boucle sur theta
722  } // Fin de la boucle sur phi
723 
724  // On remet les choses la ou il faut
725  delete [] tb->t ;
726  tb->t = xo ;
727 
728  // base de developpement
729  // pair -> impair
730  int base_t = base & MSQ_T ;
731  int base_p = base & MSQ_P ;
732  base = base_p | base_t | R_LEGI ;
733 
734  }
735 
736  //----------------
737  // cas R_LEGI ---
738  //----------------
739 
740  void _mult_x_r_legi(Tbl* tb, int& base)
741  {
742 
743  // Peut-etre rien a faire ?
744  if (tb->get_etat() == ETATZERO) {
745  int base_t = base & MSQ_T ;
746  int base_p = base & MSQ_P ;
747  base = base_p | base_t | R_LEGP ;
748  return ;
749  }
750 
751  // Pour le confort
752  int nr = (tb->dim).dim[0] ; // Nombre
753  int nt = (tb->dim).dim[1] ; // de points
754  int np = (tb->dim).dim[2] ; // physiques REELS
755  np = np - 2 ; // Nombre de points physiques
756 
757  // pt. sur le tableau de double resultat
758  double* xo = new double [tb->get_taille()];
759 
760  // Initialisation a zero :
761  for (int i=0; i<tb->get_taille(); i++) {
762  xo[i] = 0 ;
763  }
764 
765  // On y va...
766  double* xi = tb->t ;
767  double* xci = xi ; // Pointeurs
768  double* xco = xo ; // courants
769 
770  int borne_phi = np + 1 ;
771  if (np == 1) {
772  borne_phi = 1 ;
773  }
774 
775  for (int k=0 ; k< borne_phi ; k++)
776  if (k == 1) {
777  xci += nr*nt ;
778  xco += nr*nt ;
779  }
780  else {
781  for (int j=0 ; j<nt ; j++) {
782 
783  xco[0] = xci[0]/3. ;
784  for (int i = 1 ; i < nr-1 ; i++ ) {
785  xco[i] = double(2*i+1)*xci[i]/double(4*i+3)
786  + double(2*i)*xci[i-1]/double(4*i-1) ;
787  } // Fin de la premiere boucle sur r
788  xco[nr-1] = double(2*nr-2)*xci[nr-2]/double(4*nr-5) ;
789 
790  xci += nr ;
791  xco += nr ;
792  } // Fin de la boucle sur theta
793  } // Fin de la boucle sur phi
794 
795  // On remet les choses la ou il faut
796  delete [] tb->t ;
797  tb->t = xo ;
798 
799  // base de developpement
800  // impair -> pair
801  int base_t = base & MSQ_T ;
802  int base_p = base & MSQ_P ;
803  base = base_p | base_t | R_LEGP ;
804  }
805 
806  //--------------
807  // cas R_CHEB --
808  //--------------
809 
810 void _mult_x_r_cheb(Tbl* tb, int& base)
811  {
812  // Peut-etre rien a faire ?
813  if (tb->get_etat() == ETATZERO) {
814  int base_t = base & MSQ_T ;
815  int base_p = base & MSQ_P ;
816  base = base_p | base_t | R_CHEB ;
817  return ;
818  }
819 
820  // Pour le confort
821  int nr = (tb->dim).dim[0] ; // Nombre
822  int nt = (tb->dim).dim[1] ; // de points
823  int np = (tb->dim).dim[2] ; // physiques REELS
824  np = np - 2 ; // Nombre de points physiques
825 
826  // pt. sur le tableau de double resultat
827  double* xo = new double [tb->get_taille()];
828 
829  // Initialisation a zero :
830  for (int i=0; i<tb->get_taille(); i++) {
831  xo[i] = 0 ;
832  }
833 
834  // On y va...
835  double* xi = tb->t ;
836  double* xci = xi ; // Pointeurs
837  double* xco = xo ; // courants
838 
839  int borne_phi = np + 1 ;
840  if (np == 1) {
841  borne_phi = 1 ;
842  }
843 
844  for (int k=0 ; k< borne_phi ; k++)
845  if (k==1) {
846  xci += nr*nt ;
847  xco += nr*nt ;
848  }
849  else {
850  for (int j=0 ; j<nt ; j++) {
851 
852  xco[0] = .5*xci[1] ;
853  xco[1] = xci[0] + 0.5*xci[2] ;
854  for (int i = 2 ; i < nr-1 ; i++ ) {
855  xco[i] = 0.5*(xci[i-1]+xci[i+1]) ;
856  } // Fin de la boucle sur r
857  xco[nr-1] = .5*xci[nr-2] ;
858 
859  xci += nr ;
860  xco += nr ;
861  } // Fin de la boucle sur theta
862  } // Fin de la boucle sur phi
863 
864  // On remet les choses la ou il faut
865  delete [] tb->t ;
866  tb->t = xo ;
867 
868  // base de developpement
869  // pair -> impair
870  int base_t = base & MSQ_T ;
871  int base_p = base & MSQ_P ;
872  base = base_p | base_t | R_CHEB ;
873 
874  }
875 
876 
877 
878  }
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
MSQ_T
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166