00001 /* 00002 * Copyright (c) 1999-2001 Eric Gourgoulhon 00003 * Copyright (c) 1999-2001 Philippe Grandclement 00004 * 00005 * This file is part of LORENE. 00006 * 00007 * LORENE is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU General Public License as published by 00009 * the Free Software Foundation; either version 2 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * LORENE is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU General Public License 00018 * along with LORENE; if not, write to the Free Software 00019 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00020 * 00021 */ 00022 00023 00024 char map_af_lap_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_lap.C,v 1.5 2005/11/24 09:25:06 j_novak Exp $" ; 00025 00026 00027 00028 /* 00029 * $Id: map_af_lap.C,v 1.5 2005/11/24 09:25:06 j_novak Exp $ 00030 * $Log: map_af_lap.C,v $ 00031 * Revision 1.5 2005/11/24 09:25:06 j_novak 00032 * Added the Scalar version for the Laplacian 00033 * 00034 * Revision 1.4 2003/10/15 16:03:37 j_novak 00035 * Added the angular Laplace operator for Scalar. 00036 * 00037 * Revision 1.3 2003/10/03 15:58:48 j_novak 00038 * Cleaning of some headers 00039 * 00040 * Revision 1.2 2002/10/16 14:36:41 j_novak 00041 * Reorganization of #include instructions of standard C++, in order to 00042 * use experimental version 3 of gcc. 00043 * 00044 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon 00045 * LORENE 00046 * 00047 * Revision 2.16 2000/08/16 10:35:41 eric 00048 * Suppression de Mtbl::dzpuis. 00049 * 00050 * Revision 2.15 2000/02/25 09:01:47 eric 00051 * Remplacement de (uu.get_dzpuis() == 0) par uu.check_dzpuis(0). 00052 * 00053 * Revision 2.14 2000/02/08 14:19:53 phil 00054 * correction annulation derniere zone 00055 * 00056 * Revision 2.13 2000/01/27 17:52:16 phil 00057 * corrections diverses et variees 00058 * 00059 * Revision 2.12 2000/01/26 13:10:34 eric 00060 * Reprototypage complet : 00061 * le resultat est desormais suppose alloue a l'exterieur de la routine 00062 * et est passe en argument (Cmp& resu). 00063 * 00064 * Revision 2.11 1999/11/30 12:49:53 eric 00065 * Valeur::base est desormais du type Base_val et non plus Base_val*. 00066 * 00067 * Revision 2.10 1999/11/29 12:57:42 eric 00068 * REORGANISATION COMPLETE: nouveau prototype : Valeur --> Cmp 00069 * Utilisation de la nouvelle arithmetique des Valeur's. 00070 * 00071 * Revision 2.9 1999/10/27 15:44:00 eric 00072 * Suppression du membre Cmp::c. 00073 * 00074 * Revision 2.8 1999/10/14 14:27:35 eric 00075 * Methode const. 00076 * 00077 * Revision 2.7 1999/09/06 16:26:03 phil 00078 * ajout de la version Cmp 00079 * 00080 * Revision 2.6 1999/09/06 14:51:20 phil 00081 * ajout du laplacien 00082 * 00083 * Revision 2.5 1999/05/03 15:22:00 phil 00084 * Correction des bases 00085 * 00086 * Revision 2.4 1999/04/28 10:33:02 phil 00087 * correction du cas zec_mult_r = 4 00088 * 00089 * Revision 2.3 1999/04/27 09:22:25 phil 00090 * *** empty log message *** 00091 * 00092 * Revision 2.2 1999/04/27 09:17:27 phil 00093 * corrections diverses et variees .... 00094 * 00095 * Revision 2.1 1999/04/26 17:24:21 phil 00096 * correction de gestion de base 00097 * 00098 * Revision 2.0 1999/04/26 16:33:55 phil 00099 * *** empty log message *** 00100 * 00101 * 00102 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_lap.C,v 1.5 2005/11/24 09:25:06 j_novak Exp $ 00103 * 00104 */ 00105 00106 00107 // Fichiers include 00108 // ---------------- 00109 #include <stdlib.h> 00110 #include <math.h> 00111 00112 #include "cmp.h" 00113 #include "tensor.h" 00114 00115 //****************************************************************************** 00116 00117 00118 /* 00119 * Calcul du laplacien d'un Scalar ou Cmp f dans le cas ou le mapping 00120 * (coordonnees-grille) --> (coordonnees physiques) est affine. 00121 * Le Laplacien est calcule suivant l'equation: 00122 * 00123 * Lap(f) = d^2 f/dr^2 + 1/r ( 2 df/dr + 1/r Lap_ang(f) ) (1) 00124 * 00125 * avec 00126 * 00127 * Lap_ang(f) := d^2 f/dtheta^2 + 1/tan(theta) df/dtheta 00128 * + 1/sin^2(theta) d^2 f/dphi^2 (2) 00129 * 00130 * Le laplacien angulaire (2) est calcule par passage aux harmoniques 00131 * spheriques, suivant la formule 00132 * 00133 * Lap_ang( f_{lm} Y_l^m ) = - l(l+1) f_{lm} Y_l^m (3) 00134 * 00135 * Dans la zone externe compactifiee (ZEC), la routine calcule soit Lap(f), 00136 * soit r^2 Lap(f), soit r^4 Lap(f) suivant la valeur du drapeau zec_mult_r : 00137 * 00138 * -- pour zec_mult_r = 0, on calcule Lap(f) suivant la formule 00139 * 00140 * Lap(f) = u^2 [ u^2 d^2 f/du^2 + Lap_ang(f) ] , avec u = 1/r (4) 00141 * 00142 * -- pour zec_mult_r = 2, on calcule r^2 Lap(f) suivant la formule 00143 * 00144 * r^2 Lap(f) = u^2 d^2 f/du^2 + Lap_ang(f) (5) 00145 * 00146 * -- pour zec_mult_r = 4, on calcule 4^2 Lap(f) suivant la formule 00147 * 00148 * r^4 Lap(f) = d^2 f /du^2 + 1/u^2 Lap_ang(f) (6) 00149 * 00150 * 00151 * 00152 * Entree: 00153 * ------ 00154 * const Scalar& / Cmp& uu : champ f dont on veut calculer le laplacien 00155 * 00156 * 00157 * int zec_mult_r : drapeau indiquant la quantite calculee dans la ZEC: 00158 * zec_mult_r = 0 : lapff = Lap(f) 00159 * zec_mult_r = 2 : lapff = r^2 Lap(f) 00160 * zec_mult_r = 4 : lapff = r^4 Lap(f) 00161 * Sortie: 00162 * ------ 00163 * Scalar& / Cmp& resu : Lap(f) 00164 * 00165 */ 00166 00167 //********************** 00168 // Scalar version 00169 //********************** 00170 00171 void Map_af::laplacien(const Scalar& uu, int zec_mult_r, Scalar& resu) const { 00172 00173 assert (uu.get_etat() != ETATNONDEF) ; 00174 assert (uu.get_mp().get_mg() == mg) ; 00175 00176 // Particular case of null value: 00177 00178 if ((uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) { 00179 resu.set_etat_zero() ; 00180 return ; 00181 } 00182 00183 assert( uu.get_etat() == ETATQCQ ) ; 00184 assert( uu.check_dzpuis(0) ) ; 00185 resu.set_etat_qcq() ; 00186 00187 int nz = get_mg()->get_nzone() ; 00188 int nzm1 = nz - 1 ; 00189 00190 // On recupere les bases d'entree : 00191 Base_val base_entree = uu.get_spectral_base() ; 00192 00193 // Separation zones internes / ZEC : 00194 Scalar uuext(uu) ; 00195 00196 Valeur ff = uu.get_spectral_va() ; 00197 00198 if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC 00199 uuext.annule(0, nzm1-1) ; 00200 ff.annule(nzm1) ; 00201 } 00202 else { 00203 uuext.set_etat_zero() ; // pas de ZEC 00204 } 00205 00206 00207 //========================================================================= 00208 // Calcul dans les zones internes (noyau + coquilles) 00209 //========================================================================= 00210 00211 //---------------------------- 00212 // Derivations par rapport a x 00213 //---------------------------- 00214 00215 Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout 00216 00217 Valeur dff = ff.dsdx() ; // d/dx partout 00218 00219 //--------------------------- 00220 // Calcul de 1/x Lap_ang(f) ---> resultat dans ff 00221 //--------------------------- 00222 00223 //... Passage en harmoniques spheriques 00224 ff.coef() ; 00225 ff.ylm() ; 00226 00227 //... Multiplication par -l*(l+1) 00228 ff = ff.lapang() ; 00229 00230 //... Division par x: 00231 ff = ff.sx() ; // 1/x ds. le noyau 00232 // Id ds. les coquilles 00233 00234 //---------------------------------------------- 00235 // On repasse dans l'espace des configurations 00236 // pour effectuer les multiplications par 00237 // les derivees du chgmt. de var. 00238 //---------------------------------------------- 00239 00240 d2ff.coef_i() ; 00241 dff.coef_i() ; 00242 ff.coef_i() ; 00243 00244 00245 //----------------------------------------------- 00246 // Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) ) 00247 // Le resultat est mis dans ff 00248 //----------------------------------------------- 00249 00250 Base_val sauve_base = dff.base ; 00251 ff = double(2) * ( dxdr * dff) + xsr * ff ; 00252 ff.base = sauve_base ; 00253 00254 ff = ff.sx() ; // 1/x ds. le noyau 00255 // Id ds. les coquilles 00256 ff.coef_i() ; 00257 00258 //--------------------------------------------- 00259 // Calcul de Lap(f) suivant l'Eq. (1) 00260 // Le resultat est mis dans ff 00261 //----------------------------------------------- 00262 00263 sauve_base = d2ff.base ; 00264 ff = dxdr * dxdr * d2ff + xsr * ff ; 00265 ff.base = sauve_base ; 00266 00267 00268 if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC 00269 00270 //========================================================================= 00271 // Calcul dans la ZEC 00272 //========================================================================= 00273 00274 Valeur& ffext = uuext.set_spectral_va() ; 00275 00276 //---------------------------- 00277 // Derivation par rapport a x 00278 //---------------------------- 00279 00280 d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout 00281 00282 //--------------------------- 00283 // Calcul de Lap_ang(f) ---> resultat dans ffext 00284 //--------------------------- 00285 00286 //... Passage en harmoniques spheriques 00287 ffext.coef() ; 00288 ffext.ylm() ; 00289 00290 //... Multiplication par -l*(l+1) 00291 00292 ffext = ffext.lapang() ; 00293 00294 switch (zec_mult_r) { 00295 00296 case 0 : { // calcul de Lap(f) suivant l'Eq. (4) 00297 00298 d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2 00299 // par (x-1)^2 00300 00301 sauve_base = ffext.base ; 00302 ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ; 00303 ffext.base = sauve_base ; 00304 00305 ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2 00306 ffext.coef_i() ; 00307 sauve_base = ffext.base ; 00308 ffext = ffext / (xsr*xsr) ; 00309 ffext.base = sauve_base ; 00310 break ; 00311 } 00312 00313 case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5) 00314 00315 d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2 00316 // par (x-1)^2 00317 sauve_base = ffext.base ; 00318 ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ; 00319 ffext.base = sauve_base ; 00320 break ; 00321 } 00322 00323 case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6) 00324 00325 ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2 00326 00327 sauve_base = ffext.base ; 00328 ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ; 00329 ffext.base = sauve_base ; 00330 break ; 00331 } 00332 00333 default : { 00334 cout << "Map_af::laplacien : unknown value of zec_mult_r !" 00335 << endl << " zec_mult_r = " << zec_mult_r << endl ; 00336 abort() ; 00337 } 00338 } // fin des differents cas pour zec_mult_r 00339 00340 // Resultat final 00341 00342 ff = ff + ffext ; 00343 00344 } // fin du cas ou il existe une ZEC 00345 00346 // Les bases de sorties sont egales aux bases d'entree: 00347 ff.base = base_entree ; 00348 resu = ff ; 00349 resu.set_dzpuis(zec_mult_r) ; 00350 00351 } 00352 00353 00354 //********************** 00355 // Cmp version 00356 //********************** 00357 00358 00359 void Map_af::laplacien(const Cmp& uu, int zec_mult_r, Cmp& resu) const { 00360 00361 assert (uu.get_etat() != ETATNONDEF) ; 00362 assert (uu.get_mp()->get_mg() == mg) ; 00363 00364 // Particular case of null value: 00365 00366 if (uu.get_etat() == ETATZERO) { 00367 resu.set_etat_zero() ; 00368 return ; 00369 } 00370 00371 assert( uu.get_etat() == ETATQCQ ) ; 00372 assert( uu.check_dzpuis(0) ) ; 00373 resu.set_etat_qcq() ; 00374 00375 int nz = get_mg()->get_nzone() ; 00376 int nzm1 = nz - 1 ; 00377 00378 // On recupere les bases d'entree : 00379 Base_val base_entree = (uu.va).base ; 00380 00381 // Separation zones internes / ZEC : 00382 Cmp uuext(uu) ; 00383 00384 Valeur ff = uu.va ; 00385 00386 if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC 00387 uuext.annule(0, nzm1-1) ; 00388 ff.annule(nzm1) ; 00389 } 00390 else { 00391 uuext.set_etat_zero() ; // pas de ZEC 00392 } 00393 00394 00395 //========================================================================= 00396 // Calcul dans les zones internes (noyau + coquilles) 00397 //========================================================================= 00398 00399 //---------------------------- 00400 // Derivations par rapport a x 00401 //---------------------------- 00402 00403 Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout 00404 00405 Valeur dff = ff.dsdx() ; // d/dx partout 00406 00407 //--------------------------- 00408 // Calcul de 1/x Lap_ang(f) ---> resultat dans ff 00409 //--------------------------- 00410 00411 //... Passage en harmoniques spheriques 00412 ff.coef() ; 00413 ff.ylm() ; 00414 00415 //... Multiplication par -l*(l+1) 00416 ff = ff.lapang() ; 00417 00418 //... Division par x: 00419 ff = ff.sx() ; // 1/x ds. le noyau 00420 // Id ds. les coquilles 00421 00422 //---------------------------------------------- 00423 // On repasse dans l'espace des configurations 00424 // pour effectuer les multiplications par 00425 // les derivees du chgmt. de var. 00426 //---------------------------------------------- 00427 00428 d2ff.coef_i() ; 00429 dff.coef_i() ; 00430 ff.coef_i() ; 00431 00432 00433 //----------------------------------------------- 00434 // Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) ) 00435 // Le resultat est mis dans ff 00436 //----------------------------------------------- 00437 00438 Base_val sauve_base = dff.base ; 00439 ff = double(2) * ( dxdr * dff) + xsr * ff ; 00440 ff.base = sauve_base ; 00441 00442 ff = ff.sx() ; // 1/x ds. le noyau 00443 // Id ds. les coquilles 00444 ff.coef_i() ; 00445 00446 //--------------------------------------------- 00447 // Calcul de Lap(f) suivant l'Eq. (1) 00448 // Le resultat est mis dans ff 00449 //----------------------------------------------- 00450 00451 sauve_base = d2ff.base ; 00452 ff = dxdr * dxdr * d2ff + xsr * ff ; 00453 ff.base = sauve_base ; 00454 00455 00456 if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC 00457 00458 //========================================================================= 00459 // Calcul dans la ZEC 00460 //========================================================================= 00461 00462 Valeur& ffext = uuext.va ; 00463 00464 //---------------------------- 00465 // Derivation par rapport a x 00466 //---------------------------- 00467 00468 d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout 00469 00470 //--------------------------- 00471 // Calcul de Lap_ang(f) ---> resultat dans ffext 00472 //--------------------------- 00473 00474 //... Passage en harmoniques spheriques 00475 ffext.coef() ; 00476 ffext.ylm() ; 00477 00478 //... Multiplication par -l*(l+1) 00479 00480 ffext = ffext.lapang() ; 00481 00482 switch (zec_mult_r) { 00483 00484 case 0 : { // calcul de Lap(f) suivant l'Eq. (4) 00485 00486 d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2 00487 // par (x-1)^2 00488 00489 sauve_base = ffext.base ; 00490 ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ; 00491 ffext.base = sauve_base ; 00492 00493 ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2 00494 ffext.coef_i() ; 00495 sauve_base = ffext.base ; 00496 ffext = ffext / (xsr*xsr) ; 00497 ffext.base = sauve_base ; 00498 break ; 00499 } 00500 00501 case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5) 00502 00503 d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2 00504 // par (x-1)^2 00505 sauve_base = ffext.base ; 00506 ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ; 00507 ffext.base = sauve_base ; 00508 break ; 00509 } 00510 00511 case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6) 00512 00513 ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2 00514 00515 sauve_base = ffext.base ; 00516 ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ; 00517 ffext.base = sauve_base ; 00518 break ; 00519 } 00520 00521 default : { 00522 cout << "Map_af::laplacien : unknown value of zec_mult_r !" 00523 << endl << " zec_mult_r = " << zec_mult_r << endl ; 00524 abort() ; 00525 } 00526 } // fin des differents cas pour zec_mult_r 00527 00528 // Resultat final 00529 00530 ff = ff + ffext ; 00531 00532 } // fin du cas ou il existe une ZEC 00533 00534 // Les bases de sorties sont egales aux bases d'entree: 00535 ff.base = base_entree ; 00536 resu = ff ; 00537 resu.set_dzpuis(zec_mult_r) ; 00538 00539 } 00540 00541 void Map_af::lapang(const Scalar& uu, Scalar& resu) const { 00542 00543 assert (uu.get_etat() != ETATNONDEF) ; 00544 assert (uu.get_mp().get_mg() == mg) ; 00545 00546 // Particular case of null result: 00547 00548 if ( (uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) { 00549 resu.set_etat_zero() ; 00550 return ; 00551 } 00552 00553 assert( uu.get_etat() == ETATQCQ ) ; 00554 resu.set_etat_qcq() ; 00555 00556 Valeur ff = uu.get_spectral_va() ; 00557 00558 //... Passage en harmoniques spheriques 00559 ff.ylm() ; 00560 00561 //... Multiplication par -l*(l+1) 00562 resu = ff.lapang() ; 00563 00564 resu.set_dzpuis(uu.get_dzpuis()) ; 00565 00566 } 00567 00568 00569 00570
1.4.6