LORENE
map_af_lap.C
1 /*
2  * Copyright (c) 1999-2001 Eric Gourgoulhon
3  * Copyright (c) 1999-2001 Philippe Grandclement
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 
25 
26 
27 
28 /*
29  * $Id: map_af_lap.C,v 1.8 2016/12/05 16:17:57 j_novak Exp $
30  * $Log: map_af_lap.C,v $
31  * Revision 1.8 2016/12/05 16:17:57 j_novak
32  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
33  *
34  * Revision 1.7 2014/10/13 08:53:02 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.6 2014/10/06 15:13:12 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.5 2005/11/24 09:25:06 j_novak
41  * Added the Scalar version for the Laplacian
42  *
43  * Revision 1.4 2003/10/15 16:03:37 j_novak
44  * Added the angular Laplace operator for Scalar.
45  *
46  * Revision 1.3 2003/10/03 15:58:48 j_novak
47  * Cleaning of some headers
48  *
49  * Revision 1.2 2002/10/16 14:36:41 j_novak
50  * Reorganization of #include instructions of standard C++, in order to
51  * use experimental version 3 of gcc.
52  *
53  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
54  * LORENE
55  *
56  * Revision 2.16 2000/08/16 10:35:41 eric
57  * Suppression de Mtbl::dzpuis.
58  *
59  * Revision 2.15 2000/02/25 09:01:47 eric
60  * Remplacement de (uu.get_dzpuis() == 0) par uu.check_dzpuis(0).
61  *
62  * Revision 2.14 2000/02/08 14:19:53 phil
63  * correction annulation derniere zone
64  *
65  * Revision 2.13 2000/01/27 17:52:16 phil
66  * corrections diverses et variees
67  *
68  * Revision 2.12 2000/01/26 13:10:34 eric
69  * Reprototypage complet :
70  * le resultat est desormais suppose alloue a l'exterieur de la routine
71  * et est passe en argument (Cmp& resu).
72  *
73  * Revision 2.11 1999/11/30 12:49:53 eric
74  * Valeur::base est desormais du type Base_val et non plus Base_val*.
75  *
76  * Revision 2.10 1999/11/29 12:57:42 eric
77  * REORGANISATION COMPLETE: nouveau prototype : Valeur --> Cmp
78  * Utilisation de la nouvelle arithmetique des Valeur's.
79  *
80  * Revision 2.9 1999/10/27 15:44:00 eric
81  * Suppression du membre Cmp::c.
82  *
83  * Revision 2.8 1999/10/14 14:27:35 eric
84  * Methode const.
85  *
86  * Revision 2.7 1999/09/06 16:26:03 phil
87  * ajout de la version Cmp
88  *
89  * Revision 2.6 1999/09/06 14:51:20 phil
90  * ajout du laplacien
91  *
92  * Revision 2.5 1999/05/03 15:22:00 phil
93  * Correction des bases
94  *
95  * Revision 2.4 1999/04/28 10:33:02 phil
96  * correction du cas zec_mult_r = 4
97  *
98  * Revision 2.3 1999/04/27 09:22:25 phil
99  * *** empty log message ***
100  *
101  * Revision 2.2 1999/04/27 09:17:27 phil
102  * corrections diverses et variees ....
103  *
104  * Revision 2.1 1999/04/26 17:24:21 phil
105  * correction de gestion de base
106  *
107  * Revision 2.0 1999/04/26 16:33:55 phil
108  * *** empty log message ***
109  *
110  *
111  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_lap.C,v 1.8 2016/12/05 16:17:57 j_novak Exp $
112  *
113  */
114 
115 
116 // Fichiers include
117 // ----------------
118 #include <cstdlib>
119 #include <cmath>
120 
121 #include "cmp.h"
122 #include "tensor.h"
123 
124 //******************************************************************************
125 
126 
127 /*
128  * Calcul du laplacien d'un Scalar ou Cmp f dans le cas ou le mapping
129  * (coordonnees-grille) --> (coordonnees physiques) est affine.
130  * Le Laplacien est calcule suivant l'equation:
131  *
132  * Lap(f) = d^2 f/dr^2 + 1/r ( 2 df/dr + 1/r Lap_ang(f) ) (1)
133  *
134  * avec
135  *
136  * Lap_ang(f) := d^2 f/dtheta^2 + 1/tan(theta) df/dtheta
137  * + 1/sin^2(theta) d^2 f/dphi^2 (2)
138  *
139  * Le laplacien angulaire (2) est calcule par passage aux harmoniques
140  * spheriques, suivant la formule
141  *
142  * Lap_ang( f_{lm} Y_l^m ) = - l(l+1) f_{lm} Y_l^m (3)
143  *
144  * Dans la zone externe compactifiee (ZEC), la routine calcule soit Lap(f),
145  * soit r^2 Lap(f), soit r^4 Lap(f) suivant la valeur du drapeau zec_mult_r :
146  *
147  * -- pour zec_mult_r = 0, on calcule Lap(f) suivant la formule
148  *
149  * Lap(f) = u^2 [ u^2 d^2 f/du^2 + Lap_ang(f) ] , avec u = 1/r (4)
150  *
151  * -- pour zec_mult_r = 2, on calcule r^2 Lap(f) suivant la formule
152  *
153  * r^2 Lap(f) = u^2 d^2 f/du^2 + Lap_ang(f) (5)
154  *
155  * -- pour zec_mult_r = 4, on calcule 4^2 Lap(f) suivant la formule
156  *
157  * r^4 Lap(f) = d^2 f /du^2 + 1/u^2 Lap_ang(f) (6)
158  *
159  *
160  *
161  * Entree:
162  * ------
163  * const Scalar& / Cmp& uu : champ f dont on veut calculer le laplacien
164  *
165  *
166  * int zec_mult_r : drapeau indiquant la quantite calculee dans la ZEC:
167  * zec_mult_r = 0 : lapff = Lap(f)
168  * zec_mult_r = 2 : lapff = r^2 Lap(f)
169  * zec_mult_r = 4 : lapff = r^4 Lap(f)
170  * Sortie:
171  * ------
172  * Scalar& / Cmp& resu : Lap(f)
173  *
174  */
175 
176 namespace Lorene {
177 
178  //**********************
179  // Scalar version
180  //**********************
181 
182 void Map_af::laplacien(const Scalar& uu, int zec_mult_r, Scalar& resu) const {
183 
184  assert (uu.get_etat() != ETATNONDEF) ;
185  assert (uu.get_mp().get_mg() == mg) ;
186 
187  // Particular case of null value:
188 
189  if ((uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) {
190  resu.set_etat_zero() ;
191  return ;
192  }
193 
194  assert( uu.get_etat() == ETATQCQ ) ;
195  assert( uu.check_dzpuis(0) ) ;
196  resu.set_etat_qcq() ;
197 
198  int nz = get_mg()->get_nzone() ;
199  int nzm1 = nz - 1 ;
200 
201  // On recupere les bases d'entree :
202  Base_val base_entree = uu.get_spectral_base() ;
203 
204  // Separation zones internes / ZEC :
205  Scalar uuext(uu) ;
206 
207  Valeur ff = uu.get_spectral_va() ;
208 
209  if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
210  uuext.annule(0, nzm1-1) ;
211  ff.annule(nzm1) ;
212  }
213  else {
214  uuext.set_etat_zero() ; // pas de ZEC
215  }
216 
217 
218  //=========================================================================
219  // Calcul dans les zones internes (noyau + coquilles)
220  //=========================================================================
221 
222  //----------------------------
223  // Derivations par rapport a x
224  //----------------------------
225 
226  Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout
227 
228  Valeur dff = ff.dsdx() ; // d/dx partout
229 
230  //---------------------------
231  // Calcul de 1/x Lap_ang(f) ---> resultat dans ff
232  //---------------------------
233 
234  //... Passage en harmoniques spheriques
235  ff.coef() ;
236  ff.ylm() ;
237 
238  //... Multiplication par -l*(l+1)
239  ff = ff.lapang() ;
240 
241  //... Division par x:
242  ff = ff.sx() ; // 1/x ds. le noyau
243  // Id ds. les coquilles
244 
245  //----------------------------------------------
246  // On repasse dans l'espace des configurations
247  // pour effectuer les multiplications par
248  // les derivees du chgmt. de var.
249  //----------------------------------------------
250 
251  d2ff.coef_i() ;
252  dff.coef_i() ;
253  ff.coef_i() ;
254 
255 
256  //-----------------------------------------------
257  // Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) )
258  // Le resultat est mis dans ff
259  //-----------------------------------------------
260 
261  Base_val sauve_base = dff.base ;
262  ff = double(2) * ( dxdr * dff) + xsr * ff ;
263  ff.base = sauve_base ;
264 
265  ff = ff.sx() ; // 1/x ds. le noyau
266  // Id ds. les coquilles
267  ff.coef_i() ;
268 
269  //---------------------------------------------
270  // Calcul de Lap(f) suivant l'Eq. (1)
271  // Le resultat est mis dans ff
272  //-----------------------------------------------
273 
274  sauve_base = d2ff.base ;
275  ff = dxdr * dxdr * d2ff + xsr * ff ;
276  ff.base = sauve_base ;
277 
278 
279  if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
280 
281  //=========================================================================
282  // Calcul dans la ZEC
283  //=========================================================================
284 
285  Valeur& ffext = uuext.set_spectral_va() ;
286 
287  //----------------------------
288  // Derivation par rapport a x
289  //----------------------------
290 
291  d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout
292 
293  //---------------------------
294  // Calcul de Lap_ang(f) ---> resultat dans ffext
295  //---------------------------
296 
297  //... Passage en harmoniques spheriques
298  ffext.coef() ;
299  ffext.ylm() ;
300 
301  //... Multiplication par -l*(l+1)
302 
303  ffext = ffext.lapang() ;
304 
305  switch (zec_mult_r) {
306 
307  case 0 : { // calcul de Lap(f) suivant l'Eq. (4)
308 
309  d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
310  // par (x-1)^2
311 
312  sauve_base = ffext.base ;
313  ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ;
314  ffext.base = sauve_base ;
315 
316  ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2
317  ffext.coef_i() ;
318  sauve_base = ffext.base ;
319  ffext = ffext / (xsr*xsr) ;
320  ffext.base = sauve_base ;
321  break ;
322  }
323 
324  case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5)
325 
326  d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
327  // par (x-1)^2
328  sauve_base = ffext.base ;
329  ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ;
330  ffext.base = sauve_base ;
331  break ;
332  }
333 
334  case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6)
335 
336  ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2
337 
338  sauve_base = ffext.base ;
339  ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ;
340  ffext.base = sauve_base ;
341  break ;
342  }
343 
344  default : {
345  cout << "Map_af::laplacien : unknown value of zec_mult_r !"
346  << endl << " zec_mult_r = " << zec_mult_r << endl ;
347  abort() ;
348  }
349  } // fin des differents cas pour zec_mult_r
350 
351  // Resultat final
352 
353  ff = ff + ffext ;
354 
355  } // fin du cas ou il existe une ZEC
356 
357  // Les bases de sorties sont egales aux bases d'entree:
358  ff.base = base_entree ;
359  resu = ff ;
360  resu.set_dzpuis(zec_mult_r) ;
361 
362 }
363 
364 
365  //**********************
366  // Cmp version
367  //**********************
368 
369 
370 void Map_af::laplacien(const Cmp& uu, int zec_mult_r, Cmp& resu) const {
371 
372  assert (uu.get_etat() != ETATNONDEF) ;
373  assert (uu.get_mp()->get_mg() == mg) ;
374 
375  // Particular case of null value:
376 
377  if (uu.get_etat() == ETATZERO) {
378  resu.set_etat_zero() ;
379  return ;
380  }
381 
382  assert( uu.get_etat() == ETATQCQ ) ;
383  assert( uu.check_dzpuis(0) ) ;
384  resu.set_etat_qcq() ;
385 
386  int nz = get_mg()->get_nzone() ;
387  int nzm1 = nz - 1 ;
388 
389  // On recupere les bases d'entree :
390  Base_val base_entree = (uu.va).base ;
391 
392  // Separation zones internes / ZEC :
393  Cmp uuext(uu) ;
394 
395  Valeur ff = uu.va ;
396 
397  if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
398  uuext.annule(0, nzm1-1) ;
399  ff.annule(nzm1) ;
400  }
401  else {
402  uuext.set_etat_zero() ; // pas de ZEC
403  }
404 
405 
406  //=========================================================================
407  // Calcul dans les zones internes (noyau + coquilles)
408  //=========================================================================
409 
410  //----------------------------
411  // Derivations par rapport a x
412  //----------------------------
413 
414  Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout
415 
416  Valeur dff = ff.dsdx() ; // d/dx partout
417 
418  //---------------------------
419  // Calcul de 1/x Lap_ang(f) ---> resultat dans ff
420  //---------------------------
421 
422  //... Passage en harmoniques spheriques
423  ff.coef() ;
424  ff.ylm() ;
425 
426  //... Multiplication par -l*(l+1)
427  ff = ff.lapang() ;
428 
429  //... Division par x:
430  ff = ff.sx() ; // 1/x ds. le noyau
431  // Id ds. les coquilles
432 
433  //----------------------------------------------
434  // On repasse dans l'espace des configurations
435  // pour effectuer les multiplications par
436  // les derivees du chgmt. de var.
437  //----------------------------------------------
438 
439  d2ff.coef_i() ;
440  dff.coef_i() ;
441  ff.coef_i() ;
442 
443 
444  //-----------------------------------------------
445  // Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) )
446  // Le resultat est mis dans ff
447  //-----------------------------------------------
448 
449  Base_val sauve_base = dff.base ;
450  ff = double(2) * ( dxdr * dff) + xsr * ff ;
451  ff.base = sauve_base ;
452 
453  ff = ff.sx() ; // 1/x ds. le noyau
454  // Id ds. les coquilles
455  ff.coef_i() ;
456 
457  //---------------------------------------------
458  // Calcul de Lap(f) suivant l'Eq. (1)
459  // Le resultat est mis dans ff
460  //-----------------------------------------------
461 
462  sauve_base = d2ff.base ;
463  ff = dxdr * dxdr * d2ff + xsr * ff ;
464  ff.base = sauve_base ;
465 
466 
467  if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
468 
469  //=========================================================================
470  // Calcul dans la ZEC
471  //=========================================================================
472 
473  Valeur& ffext = uuext.va ;
474 
475  //----------------------------
476  // Derivation par rapport a x
477  //----------------------------
478 
479  d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout
480 
481  //---------------------------
482  // Calcul de Lap_ang(f) ---> resultat dans ffext
483  //---------------------------
484 
485  //... Passage en harmoniques spheriques
486  ffext.coef() ;
487  ffext.ylm() ;
488 
489  //... Multiplication par -l*(l+1)
490 
491  ffext = ffext.lapang() ;
492 
493  switch (zec_mult_r) {
494 
495  case 0 : { // calcul de Lap(f) suivant l'Eq. (4)
496 
497  d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
498  // par (x-1)^2
499 
500  sauve_base = ffext.base ;
501  ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ;
502  ffext.base = sauve_base ;
503 
504  ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2
505  ffext.coef_i() ;
506  sauve_base = ffext.base ;
507  ffext = ffext / (xsr*xsr) ;
508  ffext.base = sauve_base ;
509  break ;
510  }
511 
512  case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5)
513 
514  d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
515  // par (x-1)^2
516  sauve_base = ffext.base ;
517  ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ;
518  ffext.base = sauve_base ;
519  break ;
520  }
521 
522  case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6)
523 
524  ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2
525 
526  sauve_base = ffext.base ;
527  ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ;
528  ffext.base = sauve_base ;
529  break ;
530  }
531 
532  default : {
533  cout << "Map_af::laplacien : unknown value of zec_mult_r !"
534  << endl << " zec_mult_r = " << zec_mult_r << endl ;
535  abort() ;
536  }
537  } // fin des differents cas pour zec_mult_r
538 
539  // Resultat final
540 
541  ff = ff + ffext ;
542 
543  } // fin du cas ou il existe une ZEC
544 
545  // Les bases de sorties sont egales aux bases d'entree:
546  ff.base = base_entree ;
547  resu = ff ;
548  resu.set_dzpuis(zec_mult_r) ;
549 
550 }
551 
552 void Map_af::lapang(const Scalar& uu, Scalar& resu) const {
553 
554  assert (uu.get_etat() != ETATNONDEF) ;
555  assert (uu.get_mp().get_mg() == mg) ;
556 
557  // Particular case of null result:
558 
559  if ( (uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) {
560  resu.set_etat_zero() ;
561  return ;
562  }
563 
564  assert( uu.get_etat() == ETATQCQ ) ;
565  resu.set_etat_qcq() ;
566 
567  Valeur ff = uu.get_spectral_va() ;
568 
569  //... Passage en harmoniques spheriques
570  ff.ylm() ;
571 
572  //... Multiplication par -l*(l+1)
573  resu = ff.lapang() ;
574 
575  resu.set_dzpuis(uu.get_dzpuis()) ;
576 
577 }
578 
579 
580 
581 
582 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:75
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Lorene prototypes.
Definition: app_hor.h:67
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
void coef_i() const
Computes the physical value of *this.
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
void mult2_xm1_zec()
Applies the following operator to *this : \ Id (r sampling = RARE, FIN ) \ (r -sampling = UNSURR ) ...
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
const Valeur & d2sdx2() const
Returns of *this.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:292
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1581
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
const Valeur & sx2() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx2.C:117
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1570
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
virtual void lapang(const Scalar &uu, Scalar &lap) const
Computes the angular Laplacian of a scalar field.
Definition: map_af_lap.C:552
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:694
Bases of the spectral expansions.
Definition: base_val.h:325
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: cmp.C:718
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:657
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
Definition: map_af_lap.C:182
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607