LORENE
poisson_angu.C
1/*
2 * Resolution of the angular Poisson equation.
3 *
4 *
5 */
6
7/*
8 * Copyright (c) 2003-2005 Eric Gourgoulhon & Jerome Novak
9 *
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30
31
32
33/*
34 * $Id: poisson_angu.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
35 * $Log: poisson_angu.C,v $
36 * Revision 1.10 2016/12/05 16:18:09 j_novak
37 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38 *
39 * Revision 1.9 2014/10/13 08:53:29 j_novak
40 * Lorene classes and functions now belong to the namespace Lorene.
41 *
42 * Revision 1.8 2014/10/06 15:16:09 j_novak
43 * Modified #include directives to use c++ syntax.
44 *
45 * Revision 1.7 2009/10/23 12:55:04 j_novak
46 * New base T_LEG_MI
47 *
48 * Revision 1.6 2009/10/13 19:45:01 j_novak
49 * New base T_LEG_MP.
50 *
51 * Revision 1.5 2005/04/08 07:36:20 f_limousin
52 * Add #include <math.h> to avoid error in the compilation with gcc 3.3.1
53 * (problem with fabs).
54 *
55 * Revision 1.4 2005/04/05 08:34:10 e_gourgoulhon
56 * Treatment case l(l+1) = lambda.
57 *
58 * Revision 1.3 2005/04/04 21:33:37 e_gourgoulhon
59 * Solving now for the generalized angular Poisson equation
60 * Lap_ang u + lambda u = source
61 * with new parameter lambda.
62 *
63 * Revision 1.2 2004/12/17 13:35:03 m_forot
64 * Add the case T_LEG
65 *
66 * Revision 1.1 2003/10/15 21:13:28 e_gourgoulhon
67 * First version.
68 *
69 *
70 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_angu.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
71 *
72 */
73
74// Headers C
75#include <cstdlib>
76#include <cmath>
77
78// Headers Lorene
79#include "mtbl_cf.h"
80#include "grilles.h"
81#include "type_parite.h"
82
83
84 //--------------------------------------
85 // Routine pour les cas non prevus ----
86 //------------------------------------
87
88namespace Lorene {
89void _poisangu_pas_prevu(Mtbl_cf* mt, int l, double) {
90 cout << "Unknwon theta basis in the operator Mtbl_cf::poisson_angu() !" << endl ;
91 cout << " basis : " << hex << (mt->base).b[l] << endl ;
92 abort () ;
93}
94
95 //---------------
96 // cas T_LEG --
97 //---------------
98
99void _poisangu_t_leg(Mtbl_cf* mt, int l, double lambda)
100{
101
102 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
103
104 // Peut-etre rien a faire ?
105 if (tb->get_etat() == ETATZERO) {
106 return ;
107 }
108
109 int k, j, i ;
110 // Pour le confort
111 int nr = mt->get_mg()->get_nr(l) ; // Nombre
112 int nt = mt->get_mg()->get_nt(l) ; // de points
113 int np = mt->get_mg()->get_np(l) ; // physiques
114
115 int np1 = ( np == 1 ) ? 1 : np+1 ;
116
117 double* tuu = tb->t ;
118
119 // k = 0 :
120
121 for (j=0 ; j<nt ; j++) {
122 int ll = j ;
123 double xl = - ll*(ll+1) + lambda ;
124
125 if (fabs(xl) < 1.e-14) {
126 for (i=0 ; i<nr ; i++) {
127 tuu[i] = 0 ;
128 }
129 }
130 else {
131 for (i=0 ; i<nr ; i++) {
132 tuu[i] /= xl ;
133 }
134 }
135 tuu += nr ;
136 } // Fin de boucle sur theta
137
138 // On saute k = 1 :
139 tuu += nt*nr ;
140
141 // k=2,...
142 for (k=2 ; k<np1 ; k++) {
143 int m = k/2 ;
144 tuu += (m/2)*nr ;
145 for (j=m/2 ; j<nt ; j++) {
146 int ll = j ;
147 double xl = - ll*(ll+1) + lambda ;
148
149 if (fabs(xl) < 1.e-14) {
150 for (i=0 ; i<nr ; i++) {
151 tuu[i] = 0 ;
152 }
153 }
154 else {
155 for (i=0 ; i<nr ; i++) {
156 tuu[i] /= xl ;
157 }
158 }
159 tuu += nr ;
160 } // Fin de boucle sur theta
161 } // Fin de boucle sur phi
162
163 // base de developpement inchangee
164}
165
166 //---------------
167 // cas T_LEG_P --
168 //---------------
169
170void _poisangu_t_leg_p(Mtbl_cf* mt, int l, double lambda)
171{
172
173 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
174
175 // Peut-etre rien a faire ?
176 if (tb->get_etat() == ETATZERO) {
177 return ;
178 }
179
180 int k, j, i ;
181 // Pour le confort
182 int nr = mt->get_mg()->get_nr(l) ; // Nombre
183 int nt = mt->get_mg()->get_nt(l) ; // de points
184 int np = mt->get_mg()->get_np(l) ; // physiques
185
186 int np1 = ( np == 1 ) ? 1 : np+1 ;
187
188 double* tuu = tb->t ;
189
190 // k = 0 :
191
192 for (j=0 ; j<nt ; j++) {
193 int ll = 2*j ;
194 double xl = - ll*(ll+1) + lambda ;
195
196 if (fabs(xl) < 1.e-14) {
197 for (i=0 ; i<nr ; i++) {
198 tuu[i] = 0 ;
199 }
200 }
201 else {
202 for (i=0 ; i<nr ; i++) {
203 tuu[i] /= xl ;
204 }
205 }
206 tuu += nr ;
207 } // Fin de boucle sur theta
208
209 // On saute k = 1 :
210 tuu += nt*nr ;
211
212 // k=2,...
213 for (k=2 ; k<np1 ; k++) {
214 int m = k/2 ;
215 tuu += (m/2)*nr ;
216 for (j=m/2 ; j<nt ; j++) {
217 int ll = 2*j + (m%2) ;
218 double xl = - ll*(ll+1) + lambda ;
219
220 if (fabs(xl) < 1.e-14) {
221 for (i=0 ; i<nr ; i++) {
222 tuu[i] = 0 ;
223 }
224 }
225 else {
226 for (i=0 ; i<nr ; i++) {
227 tuu[i] /= xl ;
228 }
229 }
230 tuu += nr ;
231 } // Fin de boucle sur theta
232 } // Fin de boucle sur phi
233
234 // base de developpement inchangee
235}
236
237 //------------------
238 // cas T_LEG_PP --
239 //----------------
240
241void _poisangu_t_leg_pp(Mtbl_cf* mt, int l, double lambda)
242{
243
244 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
245
246 // Peut-etre rien a faire ?
247 if (tb->get_etat() == ETATZERO) {
248 return ;
249 }
250
251 int k, j, i ;
252 // Pour le confort
253 int nr = mt->get_mg()->get_nr(l) ; // Nombre
254 int nt = mt->get_mg()->get_nt(l) ; // de points
255 int np = mt->get_mg()->get_np(l) ; // physiques
256
257 int np1 = ( np == 1 ) ? 1 : np+1 ;
258
259 double* tuu = tb->t ;
260
261 // k = 0 :
262
263 for (j=0 ; j<nt ; j++) {
264 int ll = 2*j ;
265 double xl = - ll*(ll+1) + lambda ;
266
267 if (fabs(xl) < 1.e-14) {
268 for (i=0 ; i<nr ; i++) {
269 tuu[i] = 0 ;
270 }
271 }
272 else {
273 for (i=0 ; i<nr ; i++) {
274 tuu[i] /= xl ;
275 }
276 }
277 tuu += nr ;
278 } // Fin de boucle sur theta
279
280 // On saute k = 1 :
281 tuu += nt*nr ;
282
283 // k=2,...
284 for (k=2 ; k<np1 ; k++) {
285 int m = 2*(k/2);
286 tuu += (m/2)*nr ;
287 for (j=m/2 ; j<nt ; j++) {
288 int ll = 2*j ;
289 double xl = - ll*(ll+1) + lambda ;
290
291 if (fabs(xl) < 1.e-14) {
292 for (i=0 ; i<nr ; i++) {
293 tuu[i] = 0 ;
294 }
295 }
296 else {
297 for (i=0 ; i<nr ; i++) {
298 tuu[i] /= xl ;
299 }
300 }
301 tuu += nr ;
302 } // Fin de boucle sur theta
303 } // Fin de boucle sur phi
304
305 // base de developpement inchangee
306}
307
308 //----------------
309 // cas T_LEG_I --
310 //---------------
311
312void _poisangu_t_leg_i(Mtbl_cf* mt, int l, double lambda)
313{
314
315 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
316
317 // Peut-etre rien a faire ?
318 if (tb->get_etat() == ETATZERO) {
319 return ;
320 }
321
322 int k, j, i ;
323 // Pour le confort
324 int nr = mt->get_mg()->get_nr(l) ; // Nombre
325 int nt = mt->get_mg()->get_nt(l) ; // de points
326 int np = mt->get_mg()->get_np(l) ; // physiques
327
328 int np1 = ( np == 1 ) ? 1 : np+1 ;
329
330 double* tuu = tb->t ;
331
332 // k = 0 :
333
334 for (j=0 ; j<nt-1 ; j++) {
335 int ll = 2*j+1 ;
336 double xl = - ll*(ll+1) + lambda ;
337
338 if (fabs(xl) < 1.e-14) {
339 for (i=0 ; i<nr ; i++) {
340 tuu[i] = 0 ;
341 }
342 }
343 else {
344 for (i=0 ; i<nr ; i++) {
345 tuu[i] /= xl ;
346 }
347 }
348 tuu += nr ;
349 } // Fin de boucle sur theta
350 tuu += nr ; // On saute j=nt-1
351
352 // On saute k = 1 :
353 tuu += nt*nr ;
354
355 // k=2,...
356 for (k=2 ; k<np1 ; k++) {
357 int m = k/2 ;
358 tuu += ((m+1)/2)*nr ;
359 for (j=(m+1)/2 ; j<nt-1 ; j++) {
360 int ll = 2*j + ((m+1)%2) ;
361 double xl = - ll*(ll+1) + lambda ;
362
363 if (fabs(xl) < 1.e-14) {
364 for (i=0 ; i<nr ; i++) {
365 tuu[i] = 0 ;
366 }
367 }
368 else {
369 for (i=0 ; i<nr ; i++) {
370 tuu[i] /= xl ;
371 }
372 }
373 tuu += nr ;
374 } // Fin de boucle sur theta
375 tuu += nr ; // On saute j=nt-1
376 } // Fin de boucle sur phi
377
378 // base de developpement inchangee
379}
380
381 //------------------
382 // cas T_LEG_IP --
383 //----------------
384
385void _poisangu_t_leg_ip(Mtbl_cf* mt, int l, double lambda)
386{
387
388 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
389
390 // Peut-etre rien a faire ?
391 if (tb->get_etat() == ETATZERO) {
392 return ;
393 }
394
395 int k, j, i ;
396 // Pour le confort
397 int nr = mt->get_mg()->get_nr(l) ; // Nombre
398 int nt = mt->get_mg()->get_nt(l) ; // de points
399 int np = mt->get_mg()->get_np(l) ; // physiques
400
401 int np1 = ( np == 1 ) ? 1 : np+1 ;
402
403 double* tuu = tb->t ;
404
405 // k = 0 :
406
407 for (j=0 ; j<nt-1 ; j++) {
408 int ll = 2*j+1 ;
409 double xl = - ll*(ll+1) + lambda ;
410
411 if (fabs(xl) < 1.e-14) {
412 for (i=0 ; i<nr ; i++) {
413 tuu[i] = 0 ;
414 }
415 }
416 else {
417 for (i=0 ; i<nr ; i++) {
418 tuu[i] /= xl ;
419 }
420 }
421 tuu += nr ;
422 } // Fin de boucle sur theta
423 tuu += nr ; // On saute j=nt-1
424
425 // On saute k = 1 :
426 tuu += nt*nr ;
427
428 // k=2,...
429 for (k=2 ; k<np1 ; k++) {
430 int m = 2*(k/2);
431 tuu += (m/2)*nr ;
432 for (j=m/2 ; j<nt-1 ; j++) {
433 int ll = 2*j+1 ;
434 double xl = - ll*(ll+1) + lambda ;
435
436 if (fabs(xl) < 1.e-14) {
437 for (i=0 ; i<nr ; i++) {
438 tuu[i] = 0 ;
439 }
440 }
441 else {
442 for (i=0 ; i<nr ; i++) {
443 tuu[i] /= xl ;
444 }
445 }
446 tuu += nr ;
447 } // Fin de boucle sur theta
448 tuu += nr ; // On saute j=nt-1
449 } // Fin de boucle sur phi
450
451//## Verif
452 assert (tuu == tb->t + (np+1)*nt*nr) ;
453
454 // base de developpement inchangee
455}
456
457 //------------------
458 // cas T_LEG_PI --
459 //----------------
460
461void _poisangu_t_leg_pi(Mtbl_cf* mt, int l, double lambda)
462{
463
464 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
465
466 // Peut-etre rien a faire ?
467 if (tb->get_etat() == ETATZERO) {
468 return ;
469 }
470
471 int k, j, i ;
472 // Pour le confort
473 int nr = mt->get_mg()->get_nr(l) ; // Nombre
474 int nt = mt->get_mg()->get_nt(l) ; // de points
475 int np = mt->get_mg()->get_np(l) ; // physiques
476
477 double* tuu = tb->t ;
478
479 // k = 0 : cos(phi)
480 // -----
481
482 for (j=0 ; j<nt-1 ; j++) {
483 int ll = 2*j+1 ;
484 double xl = - ll*(ll+1) + lambda ;
485
486 if (fabs(xl) < 1.e-14) {
487 for (i=0 ; i<nr ; i++) {
488 tuu[i] = 0 ;
489 }
490 }
491 else {
492 for (i=0 ; i<nr ; i++) {
493 tuu[i] /= xl ;
494 }
495 }
496 tuu += nr ;
497 } // Fin de boucle sur theta
498 tuu += nr ; // On saute j=nt-1
499
500 if (np==1) {
501 return ;
502 }
503
504 // k = 1 : on saute
505 // -----
506 tuu += nt*nr ;
507
508 // k = 2 : sin(phi)
509 // ------
510 for (j=0 ; j<nt-1 ; j++) {
511 int ll = 2*j+1 ;
512 double xl = - ll*(ll+1) + lambda ;
513
514 if (fabs(xl) < 1.e-14) {
515 for (i=0 ; i<nr ; i++) {
516 tuu[i] = 0 ;
517 }
518 }
519 else {
520 for (i=0 ; i<nr ; i++) {
521 tuu[i] /= xl ;
522 }
523 }
524 tuu += nr ;
525 } // Fin de boucle sur theta
526 tuu += nr ; // On saute j=nt-1
527
528 // 3 <= k <= np
529 // ------------
530 for (k=3 ; k<np+1 ; k++) {
531 int m = (k%2 == 0) ? k-1 : k ;
532 tuu += (m-1)/2*nr ;
533 for (j=(m-1)/2 ; j<nt-1 ; j++) {
534 int ll = 2*j+1 ;
535 double xl = - ll*(ll+1) + lambda ;
536
537 if (fabs(xl) < 1.e-14) {
538 for (i=0 ; i<nr ; i++) {
539 tuu[i] = 0 ;
540 }
541 }
542 else {
543 for (i=0 ; i<nr ; i++) {
544 tuu[i] /= xl ;
545 }
546 }
547 tuu += nr ;
548 } // Fin de boucle sur theta
549 tuu += nr ; // On saute j=nt-1
550 } // Fin de boucle sur phi
551
552//## Verif
553 assert (tuu == tb->t + (np+1)*nt*nr) ;
554
555 // base de developpement inchangee
556}
557
558 //------------------
559 // cas T_LEG_II --
560 //----------------
561
562void _poisangu_t_leg_ii(Mtbl_cf* mt, int l, double lambda)
563{
564
565 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
566
567 // Peut-etre rien a faire ?
568 if (tb->get_etat() == ETATZERO) {
569 return ;
570 }
571
572 int k, j, i ;
573 // Pour le confort
574 int nr = mt->get_mg()->get_nr(l) ; // Nombre
575 int nt = mt->get_mg()->get_nt(l) ; // de points
576 int np = mt->get_mg()->get_np(l) ; // physiques
577
578 double* tuu = tb->t ;
579
580 // k = 0 : cos(phi)
581 // -----
582
583 for (j=0 ; j<nt-1 ; j++) {
584 int ll = 2*j ;
585 double xl = - ll*(ll+1) + lambda ;
586
587 if (fabs(xl) < 1.e-14) {
588 for (i=0 ; i<nr ; i++) {
589 tuu[i] = 0 ;
590 }
591 }
592 else {
593 for (i=0 ; i<nr ; i++) {
594 tuu[i] /= xl ;
595 }
596 }
597 tuu += nr ;
598 } // Fin de boucle sur theta
599 tuu += nr ; // On saute j=nt-1
600
601 if (np==1) {
602 return ;
603 }
604
605 // k = 1 : on saute
606 // -----
607 tuu += nt*nr ;
608
609 // k = 2 : sin(phi)
610 // ------
611 for (j=0 ; j<nt-1 ; j++) {
612 int ll = 2*j+1 ;
613 double xl = - ll*(ll+1) + lambda ;
614
615 if (fabs(xl) < 1.e-14) {
616 for (i=0 ; i<nr ; i++) {
617 tuu[i] = 0 ;
618 }
619 }
620 else {
621 for (i=0 ; i<nr ; i++) {
622 tuu[i] /= xl ;
623 }
624 }
625 tuu += nr ;
626 } // Fin de boucle sur theta
627 tuu += nr ; // On saute j=nt-1
628
629 // 3 <= k <= np
630 // ------------
631 for (k=3 ; k<np+1 ; k++) {
632 int m = (k%2 == 0) ? k-1 : k ;
633 tuu += (m+1)/2*nr ;
634 for (j=(m+1)/2 ; j<nt-1 ; j++) {
635 int ll = 2*j ;
636 double xl = - ll*(ll+1) + lambda ;
637
638 if (fabs(xl) < 1.e-14) {
639 for (i=0 ; i<nr ; i++) {
640 tuu[i] = 0 ;
641 }
642 }
643 else {
644 for (i=0 ; i<nr ; i++) {
645 tuu[i] /= xl ;
646 }
647 }
648 tuu += nr ;
649 } // Fin de boucle sur theta
650 tuu += nr ; // On saute j=nt-1
651 } // Fin de boucle sur phi
652
653//## Verif
654 assert (tuu == tb->t + (np+1)*nt*nr) ;
655
656 // base de developpement inchangee
657}
658
659 //------------------
660 // cas T_LEG_MP --
661 //----------------
662
663void _poisangu_t_leg_mp(Mtbl_cf* mt, int l, double lambda)
664{
665
666 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
667
668 // Peut-etre rien a faire ?
669 if (tb->get_etat() == ETATZERO) {
670 return ;
671 }
672
673 int k, j, i ;
674 // Pour le confort
675 int nr = mt->get_mg()->get_nr(l) ; // Nombre
676 int nt = mt->get_mg()->get_nt(l) ; // de points
677 int np = mt->get_mg()->get_np(l) ; // physiques
678
679 int np1 = ( np == 1 ) ? 1 : np+1 ;
680
681 double* tuu = tb->t ;
682
683 // k = 0 :
684
685 for (j=0 ; j<nt ; j++) {
686 int ll = j ;
687 double xl = - ll*(ll+1) + lambda ;
688
689 if (fabs(xl) < 1.e-14) {
690 for (i=0 ; i<nr ; i++) {
691 tuu[i] = 0 ;
692 }
693 }
694 else {
695 for (i=0 ; i<nr ; i++) {
696 tuu[i] /= xl ;
697 }
698 }
699 tuu += nr ;
700 } // Fin de boucle sur theta
701
702 // On saute k = 1 :
703 tuu += nt*nr ;
704
705 // k=2,...
706 for (k=2 ; k<np1 ; k++) {
707 int m = 2*(k/2);
708 tuu += m*nr ;
709 for (j=m ; j<nt ; j++) {
710 int ll = j ;
711 double xl = - ll*(ll+1) + lambda ;
712
713 if (fabs(xl) < 1.e-14) {
714 for (i=0 ; i<nr ; i++) {
715 tuu[i] = 0 ;
716 }
717 }
718 else {
719 for (i=0 ; i<nr ; i++) {
720 tuu[i] /= xl ;
721 }
722 }
723 tuu += nr ;
724 } // Fin de boucle sur theta
725 } // Fin de boucle sur phi
726
727 // base de developpement inchangee
728}
729
730 //----------------
731 // cas T_LEG_MI --
732 //----------------
733
734void _poisangu_t_leg_mi(Mtbl_cf* mt, int l, double lambda)
735{
736
737 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
738
739 // Peut-etre rien a faire ?
740 if (tb->get_etat() == ETATZERO) {
741 return ;
742 }
743
744 int k, j, i ;
745 // Pour le confort
746 int nr = mt->get_mg()->get_nr(l) ; // Nombre
747 int nt = mt->get_mg()->get_nt(l) ; // de points
748 int np = mt->get_mg()->get_np(l) ; // physiques
749
750 int np1 = ( np == 1 ) ? 1 : np+1 ;
751
752 double* tuu = tb->t ;
753
754 // k = 0 :
755
756 for (j=0 ; j<nt-1 ; j++) {
757 int ll = j ;
758 double xl = - ll*(ll+1) + lambda ;
759
760 if (fabs(xl) < 1.e-14) {
761 for (i=0 ; i<nr ; i++) {
762 tuu[i] = 0 ;
763 }
764 }
765 else {
766 for (i=0 ; i<nr ; i++) {
767 tuu[i] /= xl ;
768 }
769 }
770 tuu += nr ;
771 } // Fin de boucle sur theta
772 tuu += nr ; // On saute j=nt-1
773
774 // On saute k = 1 :
775 tuu += nt*nr ;
776
777 // k=2,...
778 for (k=2 ; k<np1 ; k++) {
779 int m = 2*((k-1)/2) + 1 ;
780 tuu += m*nr ;
781 for (j=m ; j<nt-1 ; j++) {
782 int ll = j ;
783 double xl = - ll*(ll+1) + lambda ;
784
785 if (fabs(xl) < 1.e-14) {
786 for (i=0 ; i<nr ; i++) {
787 tuu[i] = 0 ;
788 }
789 }
790 else {
791 for (i=0 ; i<nr ; i++) {
792 tuu[i] /= xl ;
793 }
794 }
795 tuu += nr ;
796 } // Fin de boucle sur theta
797 tuu += nr ; // On saute j=nt-1
798 } // Fin de boucle sur phi
799
800 // base de developpement inchangee
801}
802
803}
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene
Lorene prototypes.
Definition app_hor.h:67