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