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