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