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 
101 namespace Lorene {
102 void _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 
111 void _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 
166 void _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 
221 void _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 
276 void _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 
333 void _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 
393 void _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 
470 void _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 
547 void _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 
601 void _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
Lorene prototypes.
Definition: app_hor.h:67