LORENE
op_lapang.C
1/*
2 * Copyright (c) 1999-2001 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * Ensemble des routines de base pour le calcul du laplacien angulaire,
27 * c'est-a-dire de l'operateur
28 *
29 * d^2/dtheta^2 + cos(theta)/sin(theta) d/dtheta + 1/sin(theta) d^2/dphi^2
30 *
31 * (Utilisation interne)
32 *
33 * void _lapang_XXXX(Mtbl_cf * mt, int l)
34 * mt pointeur sur le Mtbl_cf d'entree-sortie
35 * l indice de la zone ou l'on doit effectuer le calcul
36 *
37 */
38
39/*
40 * $Id: op_lapang.C,v 1.10 2016/12/05 16:18:08 j_novak Exp $
41 * $Log: op_lapang.C,v $
42 * Revision 1.10 2016/12/05 16:18:08 j_novak
43 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
44 *
45 * Revision 1.9 2014/10/13 08:53:25 j_novak
46 * Lorene classes and functions now belong to the namespace Lorene.
47 *
48 * Revision 1.8 2014/10/06 15:16:06 j_novak
49 * Modified #include directives to use c++ syntax.
50 *
51 * Revision 1.7 2009/10/23 12:55:04 j_novak
52 * New base T_LEG_MI
53 *
54 * Revision 1.6 2009/10/13 19:45:01 j_novak
55 * New base T_LEG_MP.
56 *
57 * Revision 1.5 2005/05/18 07:47:36 j_novak
58 * Corrected an error for the T_LEG_II base (ll was set to 2j+1 instead of 2j for
59 * sin(phi)).
60 *
61 * Revision 1.4 2004/12/17 13:20:55 m_forot
62 * Add T_LEG base
63 *
64 * Revision 1.3 2003/09/16 12:11:59 j_novak
65 * Added the base T_LEG_II.
66 *
67 * Revision 1.2 2002/10/16 14:36:58 j_novak
68 * Reorganization of #include instructions of standard C++, in order to
69 * use experimental version 3 of gcc.
70 *
71 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
72 * LORENE
73 *
74 * Revision 2.2 2000/11/14 15:09:08 eric
75 * Traitement du cas np=1 dans T_LEG_PI
76 *
77 * Revision 2.1 2000/10/04 14:54:59 eric
78 * Ajout des bases T_LEG_IP et T_LEG_PI.
79 *
80 * Revision 2.0 1999/04/26 16:42:04 phil
81 * *** empty log message ***
82 *
83 *
84 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_lapang.C,v 1.10 2016/12/05 16:18:08 j_novak Exp $
85 *
86 */
87
88// Headers C
89#include <cstdlib>
90
91// Headers Lorene
92#include "mtbl_cf.h"
93#include "grilles.h"
94#include "type_parite.h"
95
96
97 //--------------------------------------
98 // Routine pour les cas non prevus ----
99 //------------------------------------
100
101namespace Lorene {
102void _lapang_pas_prevu(Mtbl_cf* mt, int l) {
103 cout << "Unknwon theta basis in the operator Mtbl_cf::lapang() !" << endl ;
104 cout << " basis : " << hex << (mt->base).b[l] << endl ;
105 abort () ;
106}
107 //---------------
108 // cas T_LEG --
109 //---------------
110
111void _lapang_t_leg(Mtbl_cf* mt, int l)
112{
113
114 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
115
116 // Peut-etre rien a faire ?
117 if (tb->get_etat() == ETATZERO) {
118 return ;
119 }
120
121 int k, j, i ;
122 // Pour le confort
123 int nr = mt->get_mg()->get_nr(l) ; // Nombre
124 int nt = mt->get_mg()->get_nt(l) ; // de points
125 int np = mt->get_mg()->get_np(l) ; // physiques
126
127 int np1 = ( np == 1 ) ? 1 : np+1 ;
128
129 double* tuu = tb->t ;
130
131 // k = 0 :
132
133 for (j=0 ; j<nt ; j++) {
134 int ll = j ;
135 double xl = - ll*(ll+1) ;
136 for (i=0 ; i<nr ; i++) {
137 tuu[i] *= xl ;
138 } // Fin de boucle sur r
139 tuu += nr ;
140 } // Fin de boucle sur theta
141
142 // On saute k = 1 :
143 tuu += nt*nr ;
144
145 // k=2,...
146 for (k=2 ; k<np1 ; k++) {
147 int m = k/2 ;
148 tuu += (m/2)*nr ;
149 for (j=m/2 ; j<nt ; j++) {
150 int ll = j;
151 double xl = - ll*(ll+1) ;
152 for (i=0 ; i<nr ; i++) {
153 tuu[i] *= xl ;
154 } // Fin de boucle sur r
155 tuu += nr ;
156 } // Fin de boucle sur theta
157 } // Fin de boucle sur phi
158
159 // base de developpement inchangee
160}
161
162 //---------------
163 // cas T_LEG_P --
164 //---------------
165
166void _lapang_t_leg_p(Mtbl_cf* mt, int l)
167{
168
169 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
170
171 // Peut-etre rien a faire ?
172 if (tb->get_etat() == ETATZERO) {
173 return ;
174 }
175
176 int k, j, i ;
177 // Pour le confort
178 int nr = mt->get_mg()->get_nr(l) ; // Nombre
179 int nt = mt->get_mg()->get_nt(l) ; // de points
180 int np = mt->get_mg()->get_np(l) ; // physiques
181
182 int np1 = ( np == 1 ) ? 1 : np+1 ;
183
184 double* tuu = tb->t ;
185
186 // k = 0 :
187
188 for (j=0 ; j<nt ; j++) {
189 int ll = 2*j ;
190 double xl = - ll*(ll+1) ;
191 for (i=0 ; i<nr ; i++) {
192 tuu[i] *= xl ;
193 } // Fin de boucle sur r
194 tuu += nr ;
195 } // Fin de boucle sur theta
196
197 // On saute k = 1 :
198 tuu += nt*nr ;
199
200 // k=2,...
201 for (k=2 ; k<np1 ; k++) {
202 int m = k/2 ;
203 tuu += (m/2)*nr ;
204 for (j=m/2 ; j<nt ; j++) {
205 int ll = 2*j + (m%2) ;
206 double xl = - ll*(ll+1) ;
207 for (i=0 ; i<nr ; i++) {
208 tuu[i] *= xl ;
209 } // Fin de boucle sur r
210 tuu += nr ;
211 } // Fin de boucle sur theta
212 } // Fin de boucle sur phi
213
214 // base de developpement inchangee
215}
216
217 //------------------
218 // cas T_LEG_PP --
219 //----------------
220
221void _lapang_t_leg_pp(Mtbl_cf* mt, int l)
222{
223
224 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
225
226 // Peut-etre rien a faire ?
227 if (tb->get_etat() == ETATZERO) {
228 return ;
229 }
230
231 int k, j, i ;
232 // Pour le confort
233 int nr = mt->get_mg()->get_nr(l) ; // Nombre
234 int nt = mt->get_mg()->get_nt(l) ; // de points
235 int np = mt->get_mg()->get_np(l) ; // physiques
236
237 int np1 = ( np == 1 ) ? 1 : np+1 ;
238
239 double* tuu = tb->t ;
240
241 // k = 0 :
242
243 for (j=0 ; j<nt ; j++) {
244 int ll = 2*j ;
245 double xl = - ll*(ll+1) ;
246 for (i=0 ; i<nr ; i++) {
247 tuu[i] *= xl ;
248 } // Fin de boucle sur r
249 tuu += nr ;
250 } // Fin de boucle sur theta
251
252 // On saute k = 1 :
253 tuu += nt*nr ;
254
255 // k=2,...
256 for (k=2 ; k<np1 ; k++) {
257 int m = 2*(k/2);
258 tuu += (m/2)*nr ;
259 for (j=m/2 ; j<nt ; j++) {
260 int ll = 2*j ;
261 double xl = - ll*(ll+1) ;
262 for (i=0 ; i<nr ; i++) {
263 tuu[i] *= xl ;
264 } // Fin de boucle sur r
265 tuu += nr ;
266 } // Fin de boucle sur theta
267 } // Fin de boucle sur phi
268
269 // base de developpement inchangee
270}
271
272 //----------------
273 // cas T_LEG_I --
274 //---------------
275
276void _lapang_t_leg_i(Mtbl_cf* mt, int l)
277{
278
279 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
280
281 // Peut-etre rien a faire ?
282 if (tb->get_etat() == ETATZERO) {
283 return ;
284 }
285
286 int k, j, i ;
287 // Pour le confort
288 int nr = mt->get_mg()->get_nr(l) ; // Nombre
289 int nt = mt->get_mg()->get_nt(l) ; // de points
290 int np = mt->get_mg()->get_np(l) ; // physiques
291
292 int np1 = ( np == 1 ) ? 1 : np+1 ;
293
294 double* tuu = tb->t ;
295
296 // k = 0 :
297
298 for (j=0 ; j<nt-1 ; j++) {
299 int ll = 2*j+1 ;
300 double xl = - ll*(ll+1) ;
301 for (i=0 ; i<nr ; i++) {
302 tuu[i] *= xl ;
303 } // Fin de boucle sur r
304 tuu += nr ;
305 } // Fin de boucle sur theta
306 tuu += nr ; // On saute j=nt-1
307
308 // On saute k = 1 :
309 tuu += nt*nr ;
310
311 // k=2,...
312 for (k=2 ; k<np1 ; k++) {
313 int m = k/2 ;
314 tuu += ((m+1)/2)*nr ;
315 for (j=(m+1)/2 ; j<nt-1 ; j++) {
316 int ll = 2*j + ((m+1)%2) ;
317 double xl = - ll*(ll+1) ;
318 for (i=0 ; i<nr ; i++) {
319 tuu[i] *= xl ;
320 } // Fin de boucle sur r
321 tuu += nr ;
322 } // Fin de boucle sur theta
323 tuu += nr ; // On saute j=nt-1
324 } // Fin de boucle sur phi
325
326 // base de developpement inchangee
327}
328
329 //------------------
330 // cas T_LEG_IP --
331 //----------------
332
333void _lapang_t_leg_ip(Mtbl_cf* mt, int l)
334{
335
336 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
337
338 // Peut-etre rien a faire ?
339 if (tb->get_etat() == ETATZERO) {
340 return ;
341 }
342
343 int k, j, i ;
344 // Pour le confort
345 int nr = mt->get_mg()->get_nr(l) ; // Nombre
346 int nt = mt->get_mg()->get_nt(l) ; // de points
347 int np = mt->get_mg()->get_np(l) ; // physiques
348
349 int np1 = ( np == 1 ) ? 1 : np+1 ;
350
351 double* tuu = tb->t ;
352
353 // k = 0 :
354
355 for (j=0 ; j<nt-1 ; j++) {
356 int ll = 2*j+1 ;
357 double xl = - ll*(ll+1) ;
358 for (i=0 ; i<nr ; i++) {
359 tuu[i] *= xl ;
360 } // Fin de boucle sur r
361 tuu += nr ;
362 } // Fin de boucle sur theta
363 tuu += nr ; // On saute j=nt-1
364
365 // On saute k = 1 :
366 tuu += nt*nr ;
367
368 // k=2,...
369 for (k=2 ; k<np1 ; k++) {
370 int m = 2*(k/2);
371 tuu += (m/2)*nr ;
372 for (j=m/2 ; j<nt-1 ; j++) {
373 int ll = 2*j+1 ;
374 double xl = - ll*(ll+1) ;
375 for (i=0 ; i<nr ; i++) {
376 tuu[i] *= xl ;
377 } // Fin de boucle sur r
378 tuu += nr ;
379 } // Fin de boucle sur theta
380 tuu += nr ; // On saute j=nt-1
381 } // Fin de boucle sur phi
382
383//## Verif
384 assert (tuu == tb->t + (np+1)*nt*nr) ;
385
386 // base de developpement inchangee
387}
388
389 //------------------
390 // cas T_LEG_PI --
391 //----------------
392
393void _lapang_t_leg_pi(Mtbl_cf* mt, int l)
394{
395
396 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
397
398 // Peut-etre rien a faire ?
399 if (tb->get_etat() == ETATZERO) {
400 return ;
401 }
402
403 int k, j, i ;
404 // Pour le confort
405 int nr = mt->get_mg()->get_nr(l) ; // Nombre
406 int nt = mt->get_mg()->get_nt(l) ; // de points
407 int np = mt->get_mg()->get_np(l) ; // physiques
408
409 double* tuu = tb->t ;
410
411 // k = 0 : cos(phi)
412 // -----
413
414 for (j=0 ; j<nt-1 ; j++) {
415 int ll = 2*j+1 ;
416 double xl = - ll*(ll+1) ;
417 for (i=0 ; i<nr ; i++) {
418 tuu[i] *= xl ;
419 } // Fin de boucle sur r
420 tuu += nr ;
421 } // Fin de boucle sur theta
422 tuu += nr ; // On saute j=nt-1
423
424 if (np==1) {
425 return ;
426 }
427
428 // k = 1 : on saute
429 // -----
430 tuu += nt*nr ;
431
432 // k = 2 : sin(phi)
433 // ------
434 for (j=0 ; j<nt-1 ; j++) {
435 int ll = 2*j+1 ;
436 double xl = - ll*(ll+1) ;
437 for (i=0 ; i<nr ; i++) {
438 tuu[i] *= xl ;
439 } // Fin de boucle sur r
440 tuu += nr ;
441 } // Fin de boucle sur theta
442 tuu += nr ; // On saute j=nt-1
443
444 // 3 <= k <= np
445 // ------------
446 for (k=3 ; k<np+1 ; k++) {
447 int m = (k%2 == 0) ? k-1 : k ;
448 tuu += (m-1)/2*nr ;
449 for (j=(m-1)/2 ; j<nt-1 ; j++) {
450 int ll = 2*j+1 ;
451 double xl = - ll*(ll+1) ;
452 for (i=0 ; i<nr ; i++) {
453 tuu[i] *= xl ;
454 } // Fin de boucle sur r
455 tuu += nr ;
456 } // Fin de boucle sur theta
457 tuu += nr ; // On saute j=nt-1
458 } // Fin de boucle sur phi
459
460//## Verif
461 assert (tuu == tb->t + (np+1)*nt*nr) ;
462
463 // base de developpement inchangee
464}
465
466 //------------------
467 // cas T_LEG_II --
468 //----------------
469
470void _lapang_t_leg_ii(Mtbl_cf* mt, int l)
471{
472
473 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
474
475 // Peut-etre rien a faire ?
476 if (tb->get_etat() == ETATZERO) {
477 return ;
478 }
479
480 int k, j, i ;
481 // Pour le confort
482 int nr = mt->get_mg()->get_nr(l) ; // Nombre
483 int nt = mt->get_mg()->get_nt(l) ; // de points
484 int np = mt->get_mg()->get_np(l) ; // physiques
485
486 double* tuu = tb->t ;
487
488 // k = 0 : cos(phi)
489 // -----
490
491 for (j=0 ; j<nt-1 ; j++) {
492 int ll = 2*j ;
493 double xl = - ll*(ll+1) ;
494 for (i=0 ; i<nr ; i++) {
495 tuu[i] *= xl ;
496 } // Fin de boucle sur r
497 tuu += nr ;
498 } // Fin de boucle sur theta
499 tuu += nr ; // On saute j=nt-1
500
501 if (np==1) {
502 return ;
503 }
504
505 // k = 1 : on saute
506 // -----
507 tuu += nt*nr ;
508
509 // k = 2 : sin(phi)
510 // ------
511 for (j=0 ; j<nt-1 ; j++) {
512 int ll = 2*j ;
513 double xl = - ll*(ll+1) ;
514 for (i=0 ; i<nr ; i++) {
515 tuu[i] *= xl ;
516 } // Fin de boucle sur r
517 tuu += nr ;
518 } // Fin de boucle sur theta
519 tuu += nr ; // On saute j=nt-1
520
521 // 3 <= k <= np
522 // ------------
523 for (k=3 ; k<np+1 ; k++) {
524 int m = (k%2 == 0) ? k-1 : k ;
525 tuu += (m+1)/2*nr ;
526 for (j=(m+1)/2 ; j<nt-1 ; j++) {
527 int ll = 2*j ;
528 double xl = - ll*(ll+1) ;
529 for (i=0 ; i<nr ; i++) {
530 tuu[i] *= xl ;
531 } // Fin de boucle sur r
532 tuu += nr ;
533 } // Fin de boucle sur theta
534 tuu += nr ; // On saute j=nt-1
535 } // Fin de boucle sur phi
536
537//## Verif
538 assert (tuu == tb->t + (np+1)*nt*nr) ;
539
540 // base de developpement inchangee
541}
542
543 //----------------
544 // cas T_LEG_MP --
545 //----------------
546
547void _lapang_t_leg_mp(Mtbl_cf* mt, int l)
548{
549
550 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
551
552 // Peut-etre rien a faire ?
553 if (tb->get_etat() == ETATZERO) {
554 return ;
555 }
556
557 int k, j, i ;
558 // Pour le confort
559 int nr = mt->get_mg()->get_nr(l) ; // Nombre
560 int nt = mt->get_mg()->get_nt(l) ; // de points
561 int np = mt->get_mg()->get_np(l) ; // physiques
562
563 int np1 = ( np == 1 ) ? 1 : np+1 ;
564
565 double* tuu = tb->t ;
566
567 // k = 0 :
568
569 for (j=0 ; j<nt ; j++) {
570 int ll = j ;
571 double xl = - ll*(ll+1) ;
572 for (i=0 ; i<nr ; i++) {
573 tuu[i] *= xl ;
574 } // Fin de boucle sur r
575 tuu += nr ;
576 } // Fin de boucle sur theta
577
578 // On saute k = 1 :
579 tuu += nt*nr ;
580
581 // k=2,...
582 for (k=2 ; k<np1 ; k++) {
583 int m = 2*(k/2);
584 tuu += m*nr ;
585 for (j=m ; j<nt ; j++) {
586 int ll = j ;
587 double xl = - ll*(ll+1) ;
588 for (i=0 ; i<nr ; i++) {
589 tuu[i] *= xl ;
590 } // Fin de boucle sur r
591 tuu += nr ;
592 } // Fin de boucle sur theta
593 } // Fin de boucle sur phi
594
595 // base de developpement inchangee
596}
597 //----------------
598 // cas T_LEG_MI --
599 //----------------
600
601void _lapang_t_leg_mi(Mtbl_cf* mt, int l)
602{
603
604 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
605
606 // Peut-etre rien a faire ?
607 if (tb->get_etat() == ETATZERO) {
608 return ;
609 }
610
611 int k, j, i ;
612 // Pour le confort
613 int nr = mt->get_mg()->get_nr(l) ; // Nombre
614 int nt = mt->get_mg()->get_nt(l) ; // de points
615 int np = mt->get_mg()->get_np(l) ; // physiques
616
617 int np1 = ( np == 1 ) ? 1 : np+1 ;
618
619 double* tuu = tb->t ;
620
621 // k = 0 :
622
623 for (j=0 ; j<nt ; j++) {
624 int ll = j ;
625 double xl = - ll*(ll+1) ;
626 for (i=0 ; i<nr ; i++) {
627 tuu[i] *= xl ;
628 } // Fin de boucle sur r
629 tuu += nr ;
630 } // Fin de boucle sur theta
631
632 // On saute k = 1 :
633 tuu += nt*nr ;
634
635 // k=2,...
636 for (k=2 ; k<np1 ; k++) {
637 int m = 2*((k-1)/2) + 1;
638 tuu += m*nr ;
639 for (j=m ; j<nt ; j++) {
640 int ll = j ;
641 double xl = - ll*(ll+1) ;
642 for (i=0 ; i<nr ; i++) {
643 tuu[i] *= xl ;
644 } // Fin de boucle sur r
645 tuu += nr ;
646 } // Fin de boucle sur theta
647 } // Fin de boucle sur phi
648
649 // base de developpement inchangee
650}
651}
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