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
810void _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 }
Lorene::Tbl
Basic array class.
Definition tbl.h:161
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_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition type_parite.h:154
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene
Lorene prototypes.
Definition app_hor.h:67