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
24char 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
173namespace Lorene {
174
175 //**********************
176 // Scalar version
177 //**********************
178
179void 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
367void 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
549void 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}
Bases of the spectral expansions.
Definition base_val.h:322
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:304
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
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
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
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
virtual void lapang(const Scalar &uu, Scalar &lap) const
Computes the angular Laplacian of a scalar field.
Definition map_af_lap.C:549
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1560
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Values and coefficients of a (real-value) function.
Definition valeur.h:287
const Valeur & sx2() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition valeur_sx2.C:114
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
void coef_i() const
Computes the physical value of *this.
void coef() const
Computes the coeffcients of *this.
const Valeur & d2sdx2() const
Returns of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
const Valeur & lapang() const
Returns the angular Laplacian of *this.
void mult2_xm1_zec()
Applies the following operator to *this : \ Id (r sampling = RARE, FIN ) \ (r -sampling = UNSURR )
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64