LORENE
op_d2sdx2.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 * Copyright (c) 1999-2001 Philippe Grandclement
5 *
6 * This file is part of LORENE.
7 *
8 * LORENE is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * LORENE is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with LORENE; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 *
22 */
23
24
25
26
27/*
28 * Ensemble des routines de base de derivation seconde par rapport a r
29 * (Utilisation interne)
30 *
31 * void _d2sdx2_XXXX(Tbl * t, int & b)
32 * t pointeur sur le Tbl d'entree-sortie
33 * b base spectrale
34 *
35 */
36
37/*
38 * $Id: op_d2sdx2.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
39 * $Log: op_d2sdx2.C,v $
40 * Revision 1.8 2016/12/05 16:18:07 j_novak
41 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42 *
43 * Revision 1.7 2015/03/05 08:49:32 j_novak
44 * Implemented operators with Legendre bases.
45 *
46 * Revision 1.6 2014/10/13 08:53:24 j_novak
47 * Lorene classes and functions now belong to the namespace Lorene.
48 *
49 * Revision 1.5 2013/06/14 15:54:06 j_novak
50 * Inclusion of Legendre bases.
51 *
52 * Revision 1.4 2008/08/27 08:50:10 jl_cornou
53 * Added Jacobi(0,2) polynomials
54 *
55 * Revision 1.3 2007/12/11 15:28:18 jl_cornou
56 * Jacobi(0,2) polynomials partially implemented
57 *
58 * Revision 1.2 2004/11/23 15:16:01 m_forot
59 *
60 * Added the bases for the cases without any equatorial symmetry
61 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
62 *
63 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
64 * LORENE
65 *
66 * Revision 2.7 2000/02/24 16:40:50 eric
67 * Initialisation a zero du tableau xo avant le calcul.
68 *
69 * Revision 2.6 1999/11/15 16:37:53 eric
70 * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
71 *
72 * Revision 2.5 1999/10/11 15:33:31 phil
73 * *** empty log message ***
74 *
75 * Revision 2.4 1999/10/11 14:25:47 phil
76 * realloc -> delete + new
77 *
78 * Revision 2.3 1999/09/13 11:30:23 phil
79 * gestion du cas k==1
80 *
81 * Revision 2.2 1999/03/01 15:06:31 eric
82 * *** empty log message ***
83 *
84 * Revision 2.1 1999/02/23 10:39:13 hyc
85 * *** empty log message ***
86 *
87 *
88 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdx2.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
89 *
90 */
91
92// Fichier includes
93#include "tbl.h"
94
95namespace Lorene {
96void _d2sdx2_1d_r_jaco02(int, double*, double*) ;
97
98// Prototypage
99void c_est_pas_fait(char * ) ;
100
101// Routine pour les cas non prevus
102//--------------------------------
103void _d2sdx2_pas_prevu(Tbl* , int & b) {
104 cout << "d2sdx2 pas prevu..." << endl ;
105 cout << " base: " << b << endl ;
106 abort () ;
107}
108
109// cas R_CHEBU - dzpuis = 0
110//-------------------------
111void _d2sdx2_r_chebu_0(Tbl *tb, int & )
112{
113
114 // Peut-etre rien a faire ?
115 if (tb->get_etat() == ETATZERO) {
116 return ;
117 }
118
119 // Protection
120 assert(tb->get_etat() == ETATQCQ) ;
121
122 // Pour le confort
123 int nr = (tb->dim).dim[0] ; // Nombre
124 int nt = (tb->dim).dim[1] ; // de points
125 int np = (tb->dim).dim[2] ; // physiques REELS
126 np = np - 2 ; // Nombre de points physiques
127
128 // Variables statiques
129 static double* cx1 = 0x0 ;
130 static double* cx2 = 0x0 ;
131 static double* cx3 = 0x0 ;
132 static int nr_pre = 0 ;
133
134 // Test sur np pour initialisation eventuelle
135 if (nr > nr_pre) {
136 nr_pre = nr ;
137 if (cx1 != 0x0) delete [] cx1 ;
138 if (cx2 != 0x0) delete [] cx2 ;
139 if (cx3 != 0x0) delete [] cx3 ;
140 cx1 = new double [nr] ;
141 cx2 = new double [nr] ;
142 cx3 = new double [nr] ;
143 for (int i=0 ; i<nr ; i++) {
144 cx1[i] = (i+2)*(i+2)*(i+2) ;
145 cx2[i] = (i+2) ;
146 cx3[i] = i*i ;
147 }
148 }
149
150 // pt. sur le tableau de double resultat
151 double* xo = new double[(tb->dim).taille] ;
152
153 // Initialisation a zero :
154 for (int i=0; i<(tb->dim).taille; i++) {
155 xo[i] = 0 ;
156 }
157
158 // On y va...
159 double* xi = tb->t ;
160 double* xci = xi ; // Pointeurs
161 double* xco = xo ; // courants
162
163 for (int k=0 ; k<np+1 ; k++)
164 if (k == 1) {
165 xci += nr*nt ;
166 xco += nr*nt ;
167 }
168 else {
169 for (int j=0 ; j<nt ; j++) {
170
171 double som1, som2 ;
172
173 xco[nr-1] = 0 ;
174 som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
175 som2 = (nr-1) * xci[nr-1] ;
176 xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
177 for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
178 som1 += cx1[i] * xci[i+2] ;
179 som2 += cx2[i] * xci[i+2] ;
180 xco[i] = som1 - cx3[i] * som2 ;
181 } // Fin de la premiere boucle sur r
182 xco[nr-2] = 0 ;
183 som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
184 som2 = (nr-2) * xci[nr-2] ;
185 xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
186 for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
187 som1 += cx1[i] * xci[i+2] ;
188 som2 += cx2[i] * xci[i+2] ;
189 xco[i] = som1 - cx3[i] * som2 ;
190 } // Fin de la deuxieme boucle sur r
191 xco[0] *= .5 ;
192
193 xci += nr ;
194 xco += nr ;
195 } // Fin de la boucle sur theta
196 } // Fin de la boucle sur phi
197
198 // On remet les choses la ou il faut
199 delete [] tb->t ;
200 tb->t = xo ;
201
202 // base de developpement
203 // inchangee
204}
205
206// cas R_CHEB
207//-----------
208void _d2sdx2_r_cheb(Tbl *tb, int & )
209{
210
211 // Peut-etre rien a faire ?
212 if (tb->get_etat() == ETATZERO) {
213 return ;
214 }
215
216 // Protection
217 assert(tb->get_etat() == ETATQCQ) ;
218
219 // Pour le confort
220 int nr = (tb->dim).dim[0] ; // Nombre
221 int nt = (tb->dim).dim[1] ; // de points
222 int np = (tb->dim).dim[2] ; // physiques REELS
223 np = np - 2 ; // Nombre de points physiques
224
225 // Variables statiques
226 static double* cx1 = 0x0 ;
227 static double* cx2 = 0x0 ;
228 static double* cx3 = 0x0 ;
229 static int nr_pre = 0 ;
230
231 // Test sur np pour initialisation eventuelle
232 if (nr > nr_pre) {
233 nr_pre = nr ;
234 if (cx1 != 0x0) delete [] cx1 ;
235 if (cx2 != 0x0) delete [] cx2 ;
236 if (cx3 != 0x0) delete [] cx3 ;
237 cx1 = new double [nr] ;
238 cx2 = new double [nr] ;
239 cx3 = new double [nr] ;
240
241 for (int i=0 ; i<nr ; i++) {
242 cx1[i] = (i+2)*(i+2)*(i+2) ;
243 cx2[i] = (i+2) ;
244 cx3[i] = i*i ;
245 }
246 }
247
248 // pt. sur le tableau de double resultat
249 double* xo = new double[(tb->dim).taille] ;
250
251 // Initialisation a zero :
252 for (int i=0; i<(tb->dim).taille; i++) {
253 xo[i] = 0 ;
254 }
255
256 // On y va...
257 double* xi = tb->t ;
258 double* xci = xi ; // Pointeurs
259 double* xco = xo ; // courants
260
261 for (int k=0 ; k<np+1 ; k++)
262 if (k == 1) {
263 xci += nr*nt ;
264 xco += nr*nt ;
265 }
266 else {
267 for (int j=0 ; j<nt ; j++) {
268
269 double som1, som2 ;
270
271 xco[nr-1] = 0 ;
272 som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
273 som2 = (nr-1) * xci[nr-1] ;
274 xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
275 for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
276 som1 += cx1[i] * xci[i+2] ;
277 som2 += cx2[i] * xci[i+2] ;
278 xco[i] = som1 - cx3[i] * som2 ;
279 } // Fin de la premiere boucle sur r
280 xco[nr-2] = 0 ;
281 som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
282 som2 = (nr-2) * xci[nr-2] ;
283 xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
284 for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
285 som1 += cx1[i] * xci[i+2] ;
286 som2 += cx2[i] * xci[i+2] ;
287 xco[i] = som1 - cx3[i] * som2 ;
288 } // Fin de la deuxieme boucle sur r
289 xco[0] *= .5 ;
290
291 xci += nr ;
292 xco += nr ;
293 } // Fin de la boucle sur theta
294 } // Fin de la boucle sur phi
295
296 // On remet les choses la ou il faut
297 delete [] tb->t ;
298 tb->t = xo ;
299
300 // base de developpement
301 // inchangee
302}
303
304// cas R_CHEBP
305//------------
306void _d2sdx2_r_chebp(Tbl *tb, int & )
307{
308
309 // Peut-etre rien a faire ?
310 if (tb->get_etat() == ETATZERO) {
311 return ;
312 }
313
314 // Protection
315 assert(tb->get_etat() == ETATQCQ) ;
316
317 // Pour le confort
318 int nr = (tb->dim).dim[0] ; // Nombre
319 int nt = (tb->dim).dim[1] ; // de points
320 int np = (tb->dim).dim[2] ; // physiques REELS
321 np = np - 2 ; // Nombre de points physiques
322
323 // Variables statiques
324 static double* cx1 = 0x0 ;
325 static double* cx2 = 0x0 ;
326 static double* cx3 = 0x0 ;
327 static int nr_pre = 0 ;
328
329 // Test sur np pour initialisation eventuelle
330 if (nr > nr_pre) {
331 nr_pre = nr ;
332 if (cx1 != 0x0) delete [] cx1 ;
333 if (cx2 != 0x0) delete [] cx2 ;
334 if (cx3 != 0x0) delete [] cx3 ;
335 cx1 = new double [nr] ;
336 cx2 = new double [nr] ;
337 cx3 = new double [nr] ;
338 for (int i=0 ; i<nr ; i++) {
339 cx1[i] = 8*(i+1)*(i+1)*(i+1) ;
340 cx2[i] = 2*(i+1) ;
341 cx3[i] = 4*i*i ;
342 }
343 }
344 // pt. sur le tableau de double resultat
345 double* xo = new double[(tb->dim).taille] ;
346
347 // Initialisation a zero :
348 for (int i=0; i<(tb->dim).taille; i++) {
349 xo[i] = 0 ;
350 }
351
352 // On y va...
353 double* xi = tb->t ;
354 double* xci = xi ; // Pointeurs
355 double* xco = xo ; // courants
356
357 for (int k=0 ; k<np+1 ; k++)
358 if (k == 1) {
359 xci += nr*nt ;
360 xco += nr*nt ;
361 }
362 else {
363 for (int j=0 ; j<nt ; j++) {
364
365 double som1, som2 ;
366
367 xco[nr-1] = 0 ;
368 som1 = 8*(nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
369 som2 = 2*(nr-1) * xci[nr-1] ;
370 xco[nr-2] = som1 - 4*(nr-2)*(nr-2)*som2 ;
371 for (int i = nr-3 ; i >= 0 ; i-- ) {
372 som1 += cx1[i] * xci[i+1] ;
373 som2 += cx2[i] * xci[i+1] ;
374 xco[i] = som1 - cx3[i] * som2 ;
375 } // Fin de la boucle sur r
376 xco[0] *= .5 ;
377
378 xci += nr ;
379 xco += nr ;
380 } // Fin de la boucle sur theta
381 } // Fin de la boucle sur phi
382
383 // On remet les choses la ou il faut
384 delete [] tb->t ;
385 tb->t = xo ;
386
387 // base de developpement
388 // inchangee
389}
390
391// cas R_CHEBI
392//------------
393void _d2sdx2_r_chebi(Tbl *tb, int & )
394{
395
396 // Peut-etre rien a faire ?
397 if (tb->get_etat() == ETATZERO) {
398 return ;
399 }
400
401 // Protection
402 assert(tb->get_etat() == ETATQCQ) ;
403
404 // Pour le confort
405 int nr = (tb->dim).dim[0] ; // Nombre
406 int nt = (tb->dim).dim[1] ; // de points
407 int np = (tb->dim).dim[2] ; // physiques REELS
408 np = np - 2 ; // Nombre de points physiques
409
410 // Variables statiques
411 static double* cx1 = 0x0 ;
412 static double* cx2 = 0x0 ;
413 static double* cx3 = 0x0 ;
414 static int nr_pre = 0 ;
415
416 // Test sur np pour initialisation eventuelle
417 if (nr > nr_pre) {
418 nr_pre = nr ;
419 if (cx1 != 0x0) delete [] cx1 ;
420 if (cx2 != 0x0) delete [] cx2 ;
421 if (cx3 != 0x0) delete [] cx3 ;
422 cx1 = new double [nr] ;
423 cx2 = new double [nr] ;
424 cx3 = new double [nr] ;
425
426 for (int i=0 ; i<nr ; i++) {
427 cx1[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
428 cx2[i] = (2*i+3) ;
429 cx3[i] = (2*i+1)*(2*i+1) ;
430 }
431 }
432
433 // pt. sur le tableau de double resultat
434 double* xo = new double[(tb->dim).taille] ;
435
436 // Initialisation a zero :
437 for (int i=0; i<(tb->dim).taille; i++) {
438 xo[i] = 0 ;
439 }
440
441 // On y va...
442 double* xi = tb->t ;
443 double* xci = xi ; // Pointeurs
444 double* xco = xo ; // courants
445
446 for (int k=0 ; k<np+1 ; k++)
447 if (k == 1) {
448 xci += nr*nt ;
449 xco += nr*nt ;
450 }
451 else {
452 for (int j=0 ; j<nt ; j++) {
453
454 double som1, som2 ;
455
456 xco[nr-1] = 0 ;
457 som1 = (2*nr-1)*(2*nr-1)*(2*nr-1) * xci[nr-1] ;
458 som2 = (2*nr-1) * xci[nr-1] ;
459 xco[nr-2] = som1 - (2*nr-3)*(2*nr-3)*som2 ;
460 for (int i = nr-3 ; i >= 0 ; i-- ) {
461 som1 += cx1[i] * xci[i+1] ;
462 som2 += cx2[i] * xci[i+1] ;
463 xco[i] = som1 - cx3[i] * som2 ;
464 } // Fin de la boucle sur r
465
466 xci += nr ;
467 xco += nr ;
468 } // Fin de la boucle sur theta
469 } // Fin de la boucle sur phi
470
471 // On remet les choses la ou il faut
472 delete [] tb->t ;
473 tb->t = xo ;
474
475 // base de developpement
476 // inchangee
477}
478
479// cas R_CHEBPIM_P
480//----------------
481void _d2sdx2_r_chebpim_p(Tbl *tb, int & )
482{
483
484 // Peut-etre rien a faire ?
485 if (tb->get_etat() == ETATZERO) {
486 return ;
487 }
488
489 // Protection
490 assert(tb->get_etat() == ETATQCQ) ;
491
492 // Pour le confort
493 int nr = (tb->dim).dim[0] ; // Nombre
494 int nt = (tb->dim).dim[1] ; // de points
495 int np = (tb->dim).dim[2] ; // physiques REELS
496 np = np - 2 ; // Nombre de points physiques
497
498 // Variables statiques
499 static double* cx1p = 0x0 ;
500 static double* cx2p = 0x0 ;
501 static double* cx3p = 0x0 ;
502 static double* cx1i = 0x0 ;
503 static double* cx2i = 0x0 ;
504 static double* cx3i = 0x0 ;
505 static int nr_pre = 0 ;
506
507 // Test sur np pour initialisation eventuelle
508 if (nr > nr_pre) {
509 nr_pre = nr ;
510 if (cx1p != 0x0) delete [] cx1p ;
511 if (cx2p != 0x0) delete [] cx2p ;
512 if (cx3p != 0x0) delete [] cx3p ;
513 if (cx1i != 0x0) delete [] cx1i ;
514 if (cx2i != 0x0) delete [] cx2i ;
515 if (cx3i != 0x0) delete [] cx3i ;
516 cx1p = new double[nr] ;
517 cx2p = new double[nr] ;
518 cx3p = new double[nr] ;
519 cx1i = new double[nr] ;
520 cx2i = new double[nr] ;
521 cx3i = new double[nr] ;
522 for (int i=0 ; i<nr ; i++) {
523 cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
524 cx2p[i] = 2*(i+1) ;
525 cx3p[i] = 4*i*i ;
526
527 cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
528 cx2i[i] = (2*i+3) ;
529 cx3i[i] = (2*i+1)*(2*i+1) ;
530 }
531 }
532
533 double* cx1t[2] ;
534 double* cx2t[2] ;
535 double* cx3t[2] ;
536 cx1t[0] = cx1p ; cx1t[1] = cx1i ;
537 cx2t[0] = cx2p ; cx2t[1] = cx2i ;
538 cx3t[0] = cx3p ; cx3t[1] = cx3i ;
539
540 // pt. sur le tableau de double resultat
541 double* xo = new double[(tb->dim).taille] ;
542
543 // Initialisation a zero :
544 for (int i=0; i<(tb->dim).taille; i++) {
545 xo[i] = 0 ;
546 }
547
548 // On y va...
549 double* xi = tb->t ;
550 double* xci ; // Pointeurs
551 double* xco ; // courants
552
553 // Position depart
554 xci = xi ;
555 xco = xo ;
556
557 double *cx1, *cx2, *cx3 ;
558
559 // k = 0
560 cx1 = cx1t[0] ;
561 cx2 = cx2t[0] ;
562 cx3 = cx3t[0] ;
563 for (int j=0 ; j<nt ; j++) {
564
565 double som1 = 0 ;
566 double som2 = 0 ;
567
568 xco[nr-1] = 0 ;
569 for (int i = nr-2 ; i >= 0 ; i-- ) {
570 som1 += cx1[i] * xci[i+1] ;
571 som2 += cx2[i] * xci[i+1] ;
572 xco[i] = som1 - cx3[i] * som2 ;
573 } // Fin de la boucle sur r
574 xco[0] *= .5 ; // normalisation
575 xci += nr ; // au
576 xco += nr ; // suivant
577 } // Fin de la boucle sur theta
578
579 // k = 1
580 xci += nr*nt ;
581 xco += nr*nt ;
582
583 // k >= 2
584 for (int k=2 ; k<np+1 ; k++) {
585 int m = (k/2) % 2 ;
586 cx1 = cx1t[m] ;
587 cx2 = cx2t[m] ;
588 cx3 = cx3t[m] ;
589 for (int j=0 ; j<nt ; j++) {
590
591 double som1 = 0 ;
592 double som2 = 0 ;
593
594 xco[nr-1] = 0 ;
595 for (int i = nr-2 ; i >= 0 ; i-- ) {
596 som1 += cx1[i] * xci[i+1] ;
597 som2 += cx2[i] * xci[i+1] ;
598 xco[i] = som1 - cx3[i] * som2 ;
599 } // Fin de la boucle sur r
600 if (m == 0) xco[0] *= .5 ; // normalisation eventuelle
601 xci += nr ; // au
602 xco += nr ; // suivant
603 } // Fin de la boucle sur theta
604 } // Fin de la boucle sur phi
605
606 // On remet les choses la ou il faut
607 delete [] tb->t ;
608 tb->t = xo ;
609
610 // base de developpement
611 // inchangee
612}
613
614// cas R_CHEBPIM_I
615//----------------
616void _d2sdx2_r_chebpim_i(Tbl *tb, int & )
617{
618
619 // Peut-etre rien a faire ?
620 if (tb->get_etat() == ETATZERO) {
621 return ;
622 }
623
624 // Protection
625 assert(tb->get_etat() == ETATQCQ) ;
626
627 // Pour le confort
628 int nr = (tb->dim).dim[0] ; // Nombre
629 int nt = (tb->dim).dim[1] ; // de points
630 int np = (tb->dim).dim[2] ; // physiques REELS
631 np = np - 2 ; // Nombre de points physiques
632
633 // Variables statiques
634 static double* cx1p = 0x0 ;
635 static double* cx2p = 0x0 ;
636 static double* cx3p = 0x0 ;
637 static double* cx1i = 0x0 ;
638 static double* cx2i = 0x0 ;
639 static double* cx3i = 0x0 ;
640 static int nr_pre = 0 ;
641
642 // Test sur np pour initialisation eventuelle
643 if (nr > nr_pre) {
644 nr_pre = nr ;
645 if (cx1p != 0x0) delete [] cx1p ;
646 if (cx2p != 0x0) delete [] cx2p ;
647 if (cx3p != 0x0) delete [] cx3p ;
648 if (cx1i != 0x0) delete [] cx1i ;
649 if (cx2i != 0x0) delete [] cx2i ;
650 if (cx3i != 0x0) delete [] cx3i ;
651 cx1p = new double[nr] ;
652 cx2p = new double[nr] ;
653 cx3p = new double[nr] ;
654 cx1i = new double[nr] ;
655 cx2i = new double[nr] ;
656 cx3i = new double[nr] ;
657 for (int i=0 ; i<nr ; i++) {
658 cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
659 cx2p[i] = 2*(i+1) ;
660 cx3p[i] = 4*i*i ;
661
662 cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
663 cx2i[i] = (2*i+3) ;
664 cx3i[i] = (2*i+1)*(2*i+1) ;
665 }
666 }
667
668 double* cx1t[2] ;
669 double* cx2t[2] ;
670 double* cx3t[2] ;
671 cx1t[1] = cx1p ; cx1t[0] = cx1i ;
672 cx2t[1] = cx2p ; cx2t[0] = cx2i ;
673 cx3t[1] = cx3p ; cx3t[0] = cx3i ;
674
675 // pt. sur le tableau de double resultat
676 double* xo = new double[(tb->dim).taille] ;
677
678 // Initialisation a zero :
679 for (int i=0; i<(tb->dim).taille; i++) {
680 xo[i] = 0 ;
681 }
682
683 // On y va...
684 double* xi = tb->t ;
685 double* xci ; // Pointeurs
686 double* xco ; // courants
687
688 // Position depart
689 xci = xi ;
690 xco = xo ;
691
692 double *cx1, *cx2, *cx3 ;
693
694 // k = 0
695 cx1 = cx1t[0] ;
696 cx2 = cx2t[0] ;
697 cx3 = cx3t[0] ;
698 for (int j=0 ; j<nt ; j++) {
699
700 double som1 = 0 ;
701 double som2 = 0 ;
702
703 xco[nr-1] = 0 ;
704 for (int i = nr-2 ; i >= 0 ; i-- ) {
705 som1 += cx1[i] * xci[i+1] ;
706 som2 += cx2[i] * xci[i+1] ;
707 xco[i] = som1 - cx3[i] * som2 ;
708 } // Fin de la boucle sur r
709 xci += nr ;
710 xco += nr ;
711 } // Fin de la boucle sur theta
712
713 // k = 1
714 xci += nr*nt ;
715 xco += nr*nt ;
716
717 // k >= 2
718 for (int k=2 ; k<np+1 ; k++) {
719 int m = (k/2) % 2 ;
720 cx1 = cx1t[m] ;
721 cx2 = cx2t[m] ;
722 cx3 = cx3t[m] ;
723 for (int j=0 ; j<nt ; j++) {
724
725 double som1 = 0 ;
726 double som2 = 0 ;
727
728 xco[nr-1] = 0 ;
729 for (int i = nr-2 ; i >= 0 ; i-- ) {
730 som1 += cx1[i] * xci[i+1] ;
731 som2 += cx2[i] * xci[i+1] ;
732 xco[i] = som1 - cx3[i] * som2 ;
733 } // Fin de la boucle sur r
734 if (m == 1) xco[0] *= .5 ;
735 xci += nr ;
736 xco += nr ;
737 } // Fin de la boucle sur theta
738 } // Fin de la boucle sur phi
739
740 // On remet les choses la ou il faut
741 delete [] tb->t ;
742 tb->t = xo ;
743
744 // base de developpement
745 // inchangee
746}
747
748// cas R_CHEBPI_P
749//----------------
750void _d2sdx2_r_chebpi_p(Tbl *tb, int & )
751{
752
753 // Peut-etre rien a faire ?
754 if (tb->get_etat() == ETATZERO) {
755 return ;
756 }
757
758 // Protection
759 assert(tb->get_etat() == ETATQCQ) ;
760
761 // Pour le confort
762 int nr = (tb->dim).dim[0] ; // Nombre
763 int nt = (tb->dim).dim[1] ; // de points
764 int np = (tb->dim).dim[2] ; // physiques REELS
765 np = np - 2 ; // Nombre de points physiques
766
767 // Variables statiques
768 static double* cx1p = 0x0 ;
769 static double* cx2p = 0x0 ;
770 static double* cx3p = 0x0 ;
771 static double* cx1i = 0x0 ;
772 static double* cx2i = 0x0 ;
773 static double* cx3i = 0x0 ;
774 static int nr_pre = 0 ;
775
776 // Test sur np pour initialisation eventuelle
777 if (nr > nr_pre) {
778 nr_pre = nr ;
779 if (cx1p != 0x0) delete [] cx1p ;
780 if (cx2p != 0x0) delete [] cx2p ;
781 if (cx3p != 0x0) delete [] cx3p ;
782 if (cx1i != 0x0) delete [] cx1i ;
783 if (cx2i != 0x0) delete [] cx2i ;
784 if (cx3i != 0x0) delete [] cx3i ;
785 cx1p = new double[nr] ;
786 cx2p = new double[nr] ;
787 cx3p = new double[nr] ;
788 cx1i = new double[nr] ;
789 cx2i = new double[nr] ;
790 cx3i = new double[nr] ;
791 for (int i=0 ; i<nr ; i++) {
792 cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
793 cx2p[i] = 2*(i+1) ;
794 cx3p[i] = 4*i*i ;
795
796 cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
797 cx2i[i] = (2*i+3) ;
798 cx3i[i] = (2*i+1)*(2*i+1) ;
799 }
800 }
801
802 double* cx1t[2] ;
803 double* cx2t[2] ;
804 double* cx3t[2] ;
805 cx1t[0] = cx1p ; cx1t[1] = cx1i ;
806 cx2t[0] = cx2p ; cx2t[1] = cx2i ;
807 cx3t[0] = cx3p ; cx3t[1] = cx3i ;
808
809 // pt. sur le tableau de double resultat
810 double* xo = new double[(tb->dim).taille] ;
811
812 // Initialisation a zero :
813 for (int i=0; i<(tb->dim).taille; i++) {
814 xo[i] = 0 ;
815 }
816
817 // On y va...
818 double* xi = tb->t ;
819 double* xci ; // Pointeurs
820 double* xco ; // courants
821
822 // Position depart
823 xci = xi ;
824 xco = xo ;
825
826 double *cx1, *cx2, *cx3 ;
827
828 // k = 0
829 for (int j=0 ; j<nt ; j++) {
830 int l = j % 2 ;
831 cx1 = cx1t[l] ;
832 cx2 = cx2t[l] ;
833 cx3 = cx3t[l] ;
834 double som1 = 0 ;
835 double som2 = 0 ;
836
837 xco[nr-1] = 0 ;
838 for (int i = nr-2 ; i >= 0 ; i-- ) {
839 som1 += cx1[i] * xci[i+1] ;
840 som2 += cx2[i] * xci[i+1] ;
841 xco[i] = som1 - cx3[i] * som2 ;
842 } // Fin de la boucle sur r
843 if (l == 0) xco[0] *= .5 ; // normalisation
844 xci += nr ; // au
845 xco += nr ; // suivant
846 } // Fin de la boucle sur theta
847
848 // k = 1
849 xci += nr*nt ;
850 xco += nr*nt ;
851
852 // k >= 2
853 for (int k=2 ; k<np+1 ; k++) {
854 for (int j=0 ; j<nt ; j++) {
855 int l = j % 2 ;
856 cx1 = cx1t[l] ;
857 cx2 = cx2t[l] ;
858 cx3 = cx3t[l] ;
859 double som1 = 0 ;
860 double som2 = 0 ;
861
862 xco[nr-1] = 0 ;
863 for (int i = nr-2 ; i >= 0 ; i-- ) {
864 som1 += cx1[i] * xci[i+1] ;
865 som2 += cx2[i] * xci[i+1] ;
866 xco[i] = som1 - cx3[i] * som2 ;
867 } // Fin de la boucle sur r
868 if (l == 0) xco[0] *= .5 ; // normalisation eventuelle
869 xci += nr ; // au
870 xco += nr ; // suivant
871 } // Fin de la boucle sur theta
872 } // Fin de la boucle sur phi
873
874 // On remet les choses la ou il faut
875 delete [] tb->t ;
876 tb->t = xo ;
877
878 // base de developpement
879 // inchangee
880}
881
882// cas R_CHEBPI_I
883//----------------
884void _d2sdx2_r_chebpi_i(Tbl *tb, int & )
885{
886
887 // Peut-etre rien a faire ?
888 if (tb->get_etat() == ETATZERO) {
889 return ;
890 }
891
892 // Protection
893 assert(tb->get_etat() == ETATQCQ) ;
894
895 // Pour le confort
896 int nr = (tb->dim).dim[0] ; // Nombre
897 int nt = (tb->dim).dim[1] ; // de points
898 int np = (tb->dim).dim[2] ; // physiques REELS
899 np = np - 2 ; // Nombre de points physiques
900
901 // Variables statiques
902 static double* cx1p = 0x0 ;
903 static double* cx2p = 0x0 ;
904 static double* cx3p = 0x0 ;
905 static double* cx1i = 0x0 ;
906 static double* cx2i = 0x0 ;
907 static double* cx3i = 0x0 ;
908 static int nr_pre = 0 ;
909
910 // Test sur np pour initialisation eventuelle
911 if (nr > nr_pre) {
912 nr_pre = nr ;
913 if (cx1p != 0x0) delete [] cx1p ;
914 if (cx2p != 0x0) delete [] cx2p ;
915 if (cx3p != 0x0) delete [] cx3p ;
916 if (cx1i != 0x0) delete [] cx1i ;
917 if (cx2i != 0x0) delete [] cx2i ;
918 if (cx3i != 0x0) delete [] cx3i ;
919 cx1p = new double[nr] ;
920 cx2p = new double[nr] ;
921 cx3p = new double[nr] ;
922 cx1i = new double[nr] ;
923 cx2i = new double[nr] ;
924 cx3i = new double[nr] ;
925 for (int i=0 ; i<nr ; i++) {
926 cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
927 cx2p[i] = 2*(i+1) ;
928 cx3p[i] = 4*i*i ;
929
930 cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
931 cx2i[i] = (2*i+3) ;
932 cx3i[i] = (2*i+1)*(2*i+1) ;
933 }
934 }
935
936 double* cx1t[2] ;
937 double* cx2t[2] ;
938 double* cx3t[2] ;
939 cx1t[1] = cx1p ; cx1t[0] = cx1i ;
940 cx2t[1] = cx2p ; cx2t[0] = cx2i ;
941 cx3t[1] = cx3p ; cx3t[0] = cx3i ;
942
943 // pt. sur le tableau de double resultat
944 double* xo = new double[(tb->dim).taille] ;
945
946 // Initialisation a zero :
947 for (int i=0; i<(tb->dim).taille; i++) {
948 xo[i] = 0 ;
949 }
950
951 // On y va...
952 double* xi = tb->t ;
953 double* xci ; // Pointeurs
954 double* xco ; // courants
955
956 // Position depart
957 xci = xi ;
958 xco = xo ;
959
960 double *cx1, *cx2, *cx3 ;
961
962 // k = 0
963 for (int j=0 ; j<nt ; j++) {
964 int l = j % 2 ;
965 cx1 = cx1t[l] ;
966 cx2 = cx2t[l] ;
967 cx3 = cx3t[l] ;
968 double som1 = 0 ;
969 double som2 = 0 ;
970
971 xco[nr-1] = 0 ;
972 for (int i = nr-2 ; i >= 0 ; i-- ) {
973 som1 += cx1[i] * xci[i+1] ;
974 som2 += cx2[i] * xci[i+1] ;
975 xco[i] = som1 - cx3[i] * som2 ;
976 } // Fin de la boucle sur r
977 if (l == 1) xco[0] *= .5 ;
978 xci += nr ;
979 xco += nr ;
980 } // Fin de la boucle sur theta
981
982 // k = 1
983 xci += nr*nt ;
984 xco += nr*nt ;
985
986 // k >= 2
987 for (int k=2 ; k<np+1 ; k++) {
988 for (int j=0 ; j<nt ; j++) {
989 int l = j % 2 ;
990 cx1 = cx1t[l] ;
991 cx2 = cx2t[l] ;
992 cx3 = cx3t[l] ;
993 double som1 = 0 ;
994 double som2 = 0 ;
995
996 xco[nr-1] = 0 ;
997 for (int i = nr-2 ; i >= 0 ; i-- ) {
998 som1 += cx1[i] * xci[i+1] ;
999 som2 += cx2[i] * xci[i+1] ;
1000 xco[i] = som1 - cx3[i] * som2 ;
1001 } // Fin de la boucle sur r
1002 if (l == 1) xco[0] *= .5 ;
1003 xci += nr ;
1004 xco += nr ;
1005 } // Fin de la boucle sur theta
1006 } // Fin de la boucle sur phi
1007
1008 // On remet les choses la ou il faut
1009 delete [] tb->t ;
1010 tb->t = xo ;
1011
1012 // base de developpement
1013 // inchangee
1014}
1015
1016
1017// cas R_LEG
1018//-----------
1019void _d2sdx2_r_leg(Tbl *tb, int & )
1020{
1021
1022 // Peut-etre rien a faire ?
1023 if (tb->get_etat() == ETATZERO) {
1024 return ;
1025 }
1026
1027 // Protection
1028 assert(tb->get_etat() == ETATQCQ) ;
1029
1030 // Pour le confort
1031 int nr = (tb->dim).dim[0] ; // Nombre
1032 int nt = (tb->dim).dim[1] ; // de points
1033 int np = (tb->dim).dim[2] ; // physiques REELS
1034 np = np - 2 ; // Nombre de points physiques
1035
1036 // Variables statiques
1037 static double* cx1 = 0x0 ;
1038 static double* cx2 = 0x0 ;
1039 static double* cx3 = 0x0 ;
1040 static int nr_pre = 0 ;
1041
1042 // Test sur np pour initialisation eventuelle
1043 if (nr > nr_pre) {
1044 nr_pre = nr ;
1045 if (cx1 != 0x0) delete [] cx1 ;
1046 if (cx2 != 0x0) delete [] cx2 ;
1047 if (cx3 != 0x0) delete [] cx3 ;
1048 cx1 = new double [nr] ;
1049 cx2 = new double [nr] ;
1050 cx3 = new double [nr] ;
1051
1052 for (int i=0 ; i<nr ; i++) {
1053 cx1[i] = (i+2)*(i+3) ;
1054 cx2[i] = i*(i+1) ;
1055 cx3[i] = double(i) + 0.5 ;
1056 }
1057 }
1058
1059 // pt. sur le tableau de double resultat
1060 double* xo = new double[(tb->dim).taille] ;
1061
1062 // Initialisation a zero :
1063 for (int i=0; i<(tb->dim).taille; i++) {
1064 xo[i] = 0 ;
1065 }
1066
1067 // On y va...
1068 double* xi = tb->t ;
1069 double* xci = xi ; // Pointeurs
1070 double* xco = xo ; // courants
1071
1072 for (int k=0 ; k<np+1 ; k++)
1073 if (k == 1) {
1074 xci += nr*nt ;
1075 xco += nr*nt ;
1076 }
1077 else {
1078 for (int j=0 ; j<nt ; j++) {
1079
1080 double som1, som2 ;
1081
1082 xco[nr-1] = 0 ;
1083 som1 = (nr-1)*nr * xci[nr-1] ;
1084 som2 = xci[nr-1] ;
1085 if (nr > 2)
1086 xco[nr-3] = (double(nr) -2.5) * (som1 - (nr-3)*(nr-2)*som2) ;
1087 for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
1088 som1 += cx1[i] * xci[i+2] ;
1089 som2 += xci[i+2] ;
1090 xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1091 } // Fin de la premiere boucle sur r
1092 if (nr > 1) xco[nr-2] = 0 ;
1093 if (nr > 3) {
1094 som1 = (nr-2)*(nr-1)* xci[nr-2] ;
1095 som2 = xci[nr-2] ;
1096 xco[nr-4] = (double(nr) - 3.5) * (som1 - (nr-4)*(nr-3)*som2) ;
1097 }
1098 for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
1099 som1 += cx1[i] * xci[i+2] ;
1100 som2 += xci[i+2] ;
1101 xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1102 } // Fin de la deuxieme boucle sur r
1103
1104 xci += nr ;
1105 xco += nr ;
1106 } // Fin de la boucle sur theta
1107 } // Fin de la boucle sur phi
1108
1109 // On remet les choses la ou il faut
1110 delete [] tb->t ;
1111 tb->t = xo ;
1112
1113 // base de developpement
1114 // inchangee
1115}
1116
1117// cas R_LEGP
1118//------------
1119void _d2sdx2_r_legp(Tbl *tb, int & )
1120{
1121
1122 // Peut-etre rien a faire ?
1123 if (tb->get_etat() == ETATZERO) {
1124 return ;
1125 }
1126
1127 // Protection
1128 assert(tb->get_etat() == ETATQCQ) ;
1129
1130 // Pour le confort
1131 int nr = (tb->dim).dim[0] ; // Nombre
1132 int nt = (tb->dim).dim[1] ; // de points
1133 int np = (tb->dim).dim[2] ; // physiques REELS
1134 np = np - 2 ; // Nombre de points physiques
1135
1136 // Variables statiques
1137 static double* cx1 = 0x0 ;
1138 static double* cx2 = 0x0 ;
1139 static double* cx3 = 0x0 ;
1140 static int nr_pre = 0 ;
1141
1142 // Test sur np pour initialisation eventuelle
1143 if (nr > nr_pre) {
1144 nr_pre = nr ;
1145 if (cx1 != 0x0) delete [] cx1 ;
1146 if (cx2 != 0x0) delete [] cx2 ;
1147 if (cx3 != 0x0) delete [] cx3 ;
1148 cx1 = new double [nr] ;
1149 cx2 = new double [nr] ;
1150 cx3 = new double [nr] ;
1151 for (int i=0 ; i<nr ; i++) {
1152 cx1[i] = (2*i+2)*(2*i+3) ;
1153 cx2[i] = 2*i*(2*i+1) ;
1154 cx3[i] = double(2*i) + 0.5 ;
1155 }
1156 }
1157 // pt. sur le tableau de double resultat
1158 double* xo = new double[(tb->dim).taille] ;
1159
1160 // Initialisation a zero :
1161 for (int i=0; i<(tb->dim).taille; i++) {
1162 xo[i] = 0 ;
1163 }
1164
1165 // On y va...
1166 double* xi = tb->t ;
1167 double* xci = xi ; // Pointeurs
1168 double* xco = xo ; // courants
1169
1170 for (int k=0 ; k<np+1 ; k++)
1171 if (k == 1) {
1172 xci += nr*nt ;
1173 xco += nr*nt ;
1174 }
1175 else {
1176 for (int j=0 ; j<nt ; j++) {
1177
1178 double som1, som2 ;
1179
1180 xco[nr-1] = 0 ;
1181 som1 = (2*nr-2)*(2*nr-1)* xci[nr-1] ;
1182 som2 = xci[nr-1] ;
1183 if (nr > 1)
1184 xco[nr-2] = (double(2*nr) - 1.5)*(som1 - 2*(nr-2)*(2*nr-1)*som2) ;
1185 for (int i = nr-3 ; i >= 0 ; i-- ) {
1186 som1 += cx1[i] * xci[i+1] ;
1187 som2 += xci[i+1] ;
1188 xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1189 } // Fin de la boucle sur r
1190
1191 xci += nr ;
1192 xco += nr ;
1193 } // Fin de la boucle sur theta
1194 } // Fin de la boucle sur phi
1195
1196 // On remet les choses la ou il faut
1197 delete [] tb->t ;
1198 tb->t = xo ;
1199
1200 // base de developpement
1201 // inchangee
1202}
1203
1204// cas R_LEGI
1205//------------
1206void _d2sdx2_r_legi(Tbl *tb, int & )
1207{
1208
1209 // Peut-etre rien a faire ?
1210 if (tb->get_etat() == ETATZERO) {
1211 return ;
1212 }
1213
1214 // Protection
1215 assert(tb->get_etat() == ETATQCQ) ;
1216
1217 // Pour le confort
1218 int nr = (tb->dim).dim[0] ; // Nombre
1219 int nt = (tb->dim).dim[1] ; // de points
1220 int np = (tb->dim).dim[2] ; // physiques REELS
1221 np = np - 2 ; // Nombre de points physiques
1222
1223 // Variables statiques
1224 static double* cx1 = 0x0 ;
1225 static double* cx2 = 0x0 ;
1226 static double* cx3 = 0x0 ;
1227 static int nr_pre = 0 ;
1228
1229 // Test sur np pour initialisation eventuelle
1230 if (nr > nr_pre) {
1231 nr_pre = nr ;
1232 if (cx1 != 0x0) delete [] cx1 ;
1233 if (cx2 != 0x0) delete [] cx2 ;
1234 if (cx3 != 0x0) delete [] cx3 ;
1235 cx1 = new double [nr] ;
1236 cx2 = new double [nr] ;
1237 cx3 = new double [nr] ;
1238
1239 for (int i=0 ; i<nr ; i++) {
1240 cx1[i] = (2*i+3)*(2*i+4) ;
1241 cx2[i] = (2*i+1)*(2*i+2) ;
1242 cx3[i] = double(2*i) + 1.5 ;
1243 }
1244 }
1245
1246 // pt. sur le tableau de double resultat
1247 double* xo = new double[(tb->dim).taille] ;
1248
1249 // Initialisation a zero :
1250 for (int i=0; i<(tb->dim).taille; i++) {
1251 xo[i] = 0 ;
1252 }
1253
1254 // On y va...
1255 double* xi = tb->t ;
1256 double* xci = xi ; // Pointeurs
1257 double* xco = xo ; // courants
1258
1259 for (int k=0 ; k<np+1 ; k++)
1260 if (k == 1) {
1261 xci += nr*nt ;
1262 xco += nr*nt ;
1263 }
1264 else {
1265 for (int j=0 ; j<nt ; j++) {
1266
1267 double som1, som2 ;
1268
1269 xco[nr-1] = 0 ;
1270 som1 = (2*nr-1)*(2*nr) * xci[nr-1] ;
1271 som2 = xci[nr-1] ;
1272 if (nr > 1)
1273 xco[nr-2] = (double(nr) - 1.5)*(som1 - (2*nr-3)*(2*nr-2)*som2) ;
1274 for (int i = nr-3 ; i >= 0 ; i-- ) {
1275 som1 += cx1[i] * xci[i+1] ;
1276 som2 += xci[i+1] ;
1277 xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1278 } // Fin de la boucle sur r
1279
1280 xci += nr ;
1281 xco += nr ;
1282 } // Fin de la boucle sur theta
1283 } // Fin de la boucle sur phi
1284
1285 // On remet les choses la ou il faut
1286 delete [] tb->t ;
1287 tb->t = xo ;
1288
1289 // base de developpement
1290 // inchangee
1291}
1292
1293
1294
1295// cas R_JACO02
1296//-----------
1297void _d2sdx2_r_jaco02(Tbl *tb, int & )
1298{
1299
1300 // Peut-etre rien a faire ?
1301 if (tb->get_etat() == ETATZERO) {
1302 return ;
1303 }
1304
1305 // Protection
1306 assert(tb->get_etat() == ETATQCQ) ;
1307
1308 // Pour le confort
1309 int nr = (tb->dim).dim[0] ; // Nombre
1310 int nt = (tb->dim).dim[1] ; // de points
1311 int np = (tb->dim).dim[2] ; // physiques REELS
1312 np = np - 2 ; // Nombre de points physiques
1313
1314 // pt. sur le tableau de double resultat
1315 double* xo = new double[(tb->dim).taille] ;
1316
1317 // Initialisation a zero :
1318 for (int i=0; i<(tb->dim).taille; i++) {
1319 xo[i] = 0 ;
1320 }
1321
1322 // On y va...
1323 double* xi = tb->t ;
1324 double* xci = xi ; // Pointeurs
1325 double* xco = xo ; // courants
1326
1327 for (int k=0 ; k<np+1 ; k++)
1328 if (k == 1) {
1329 xci += nr*nt ;
1330 xco += nr*nt ;
1331 }
1332 else {
1333 for (int j=0 ; j<nt ; j++) {
1334
1335 double* tb1 = new double[nr] ;
1336 for (int m =0 ; m<nr ; m++) {
1337 tb1[m]=xci[m];
1338 }
1339 double* res = new double[nr] ;
1340 _d2sdx2_1d_r_jaco02(nr,tb1,res) ;
1341 for (int i = 0 ; i<nr ; i++ ) {
1342 xco[i] = res[i] ;
1343 }
1344 delete [] res ;
1345 delete [] tb1 ;
1346
1347 xci += nr ;
1348 xco += nr ;
1349 } // Fin de la boucle sur theta
1350 } // Fin de la boucle sur phi
1351
1352 // On remet les choses la ou il faut
1353 delete [] tb->t ;
1354 tb->t = xo ;
1355
1356 // base de developpement
1357 // inchangee
1358}
1359
1360
1361}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::c_est_pas_fait
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Definition c_est_pas_fait.C:76
Lorene
Lorene prototypes.
Definition app_hor.h:67