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