LORENE
poisson_falloff.C
1/*
2 * Method of Non_class_members for solving Poisson equations
3 * with a falloff condition at the outer boundary
4 *
5 */
6
7/*
8 * Copyright (c) 2004 Joshua A. Faber
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27char poisson_falloff_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_falloff.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $" ;
28
29/*
30 * $Id: poisson_falloff.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $
31 * $Log: poisson_falloff.C,v $
32 * Revision 1.3 2014/10/13 08:53:29 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.2 2014/10/06 15:16:09 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.1 2004/11/30 20:55:03 k_taniguchi
39 * *** empty log message ***
40 *
41 *
42 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_falloff.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $
43 *
44 */
45
46// Header C :
47#include <cstdlib>
48#include <cmath>
49
50// Headers Lorene :
51#include "matrice.h"
52#include "tbl.h"
53#include "mtbl_cf.h"
54#include "map.h"
55#include "base_val.h"
56#include "proto.h"
57#include "type_parite.h"
58
59
60
61 //----------------------------------------------
62 // Version Mtbl_cf
63 //----------------------------------------------
64
65/*
66 *
67 * Solution de l'equation de poisson
68 *
69 * Entree : mapping : le mapping affine
70 * source : les coefficients de la source qui a ete multipliee par
71 * r^4 ou r^2 dans la ZEC.
72 * La base de decomposition doit etre Ylm
73 * k_falloff: exponent in radial dependence of field: phi \propto r^{-k}
74 * Sortie : renvoie les coefficients de la solution dans la meme base de
75 * decomposition (a savoir Ylm)
76 *
77 */
78
79
80
81namespace Lorene {
82Mtbl_cf sol_poisson_falloff(const Map_af& mapping, const Mtbl_cf& source, const int k_falloff)
83{
84
85 // Verifications d'usage sur les zones
86 int nz = source.get_mg()->get_nzone() ;
87 assert (nz>1) ;
88 assert (source.get_mg()->get_type_r(0) == RARE) ;
89 // assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
90 for (int l=1 ; l<nz ; l++)
91 assert(source.get_mg()->get_type_r(l) == FIN) ;
92
93 // assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
94
95
96 // Bases spectrales
97 const Base_val& base = source.base ;
98
99
100 // donnees sur la zone
101 int nr, nt, np ;
102 int base_r ;
103 double alpha, beta, echelle ;
104 int l_quant, m_quant;
105
106 //Rangement des valeurs intermediaires
107 Tbl *so ;
108 Tbl *sol_hom ;
109 Tbl *sol_part ;
110 Matrice *operateur ;
111 Matrice *nondege ;
112
113
114 // Rangement des solutions, avant raccordement
115 Mtbl_cf solution_part(source.get_mg(), base) ;
116 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
117 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
118 Mtbl_cf resultat(source.get_mg(), base) ;
119
120 solution_part.set_etat_qcq() ;
121 solution_hom_un.set_etat_qcq() ;
122 solution_hom_deux.set_etat_qcq() ;
123 resultat.set_etat_qcq() ;
124
125 for (int l=0 ; l<nz ; l++) {
126 solution_part.t[l]->set_etat_qcq() ;
127 solution_hom_un.t[l]->set_etat_qcq() ;
128 solution_hom_deux.t[l]->set_etat_qcq() ;
129 resultat.t[l]->set_etat_qcq() ;
130 for (int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
131 for (int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
132 for (int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
133 resultat.set(l, k, j, i) = 0 ;
134 }
135
136 // nbre maximum de point en theta et en phi :
137 int np_max, nt_max ;
138
139 //---------------
140 //-- NOYAU -----
141 //---------------
142
143 nr = source.get_mg()->get_nr(0) ;
144 nt = source.get_mg()->get_nt(0) ;
145 np = source.get_mg()->get_np(0) ;
146
147 nt_max = nt ;
148 np_max = np ;
149
150 alpha = mapping.get_alpha()[0] ;
151 beta = mapping.get_beta()[0] ;
152
153 for (int k=0 ; k<np+1 ; k++)
154 for (int j=0 ; j<nt ; j++)
155 if (nullite_plm(j, nt, k, np, base) == 1)
156 {
157 // calcul des nombres quantiques :
158 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
159
160 // Construction operateur
161 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
162 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
163
164 //Operateur inversible
165 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
166
167 // Calcul de la SH
168 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
169
170 //Calcul de la SP
171 so = new Tbl(nr) ;
172 so->set_etat_qcq() ;
173 for (int i=0 ; i<nr ; i++)
174 so->set(i) = source(0, k, j, i) ;
175
176 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
177 *so, 0, base_r)) ;
178
179 // Rangement dans les tableaux globaux ;
180
181 for (int i=0 ; i<nr ; i++) {
182 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
183 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
184 solution_hom_deux.set(0, k, j, i) = 0. ;
185 }
186
187
188
189 delete operateur ;
190 delete nondege ;
191 delete so ;
192 delete sol_hom ;
193 delete sol_part ;
194 }
195
196
197
198 //---------------
199 //- COQUILLES ---
200 //---------------
201
202 for (int zone=1 ; zone<nz ; zone++) {
203 nr = source.get_mg()->get_nr(zone) ;
204 nt = source.get_mg()->get_nt(zone) ;
205 np = source.get_mg()->get_np(zone) ;
206
207 if (np > np_max) np_max = np ;
208 if (nt > nt_max) nt_max = nt ;
209
210 alpha = mapping.get_alpha()[zone] ;
211 beta = mapping.get_beta()[zone] ;
212
213 for (int k=0 ; k<np+1 ; k++)
214 for (int j=0 ; j<nt ; j++)
215 if (nullite_plm(j, nt, k, np, base) == 1)
216 {
217 // calcul des nombres quantiques :
218 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
219
220 // Construction de l'operateur
221 operateur = new Matrice(laplacien_mat
222 (nr, l_quant, beta/alpha, 0, base_r)) ;
223
224 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
225 base_r) ;
226
227 // Operateur inversible
228 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
229 beta/alpha, 0, base_r)) ;
230
231 // Calcul DES DEUX SH
232 sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
233
234 // Calcul de la SP
235 so = new Tbl(nr) ;
236 so->set_etat_qcq() ;
237 for (int i=0 ; i<nr ; i++)
238 so->set(i) = source(zone, k, j, i) ;
239
240 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
241 beta, *so, 0, base_r)) ;
242
243
244 // Rangement
245 for (int i=0 ; i<nr ; i++) {
246 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
247 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
248 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
249 }
250
251
252 delete operateur ;
253 delete nondege ;
254 delete so ;
255 delete sol_hom ;
256 delete sol_part ;
257 }
258 }
259
260 //-------------------------------------------
261 // On est parti pour le raccord des solutions
262 //-------------------------------------------
263 // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
264 // intervient dans le developpement de cette zone.
265 int * indic = new int [nz] ;
266 int taille ;
267 Tbl *sec_membre ; // termes constants du systeme
268 Matrice *systeme ; //le systeme a resoudre
269
270 // des compteurs pour le remplissage du systeme
271 int ligne ;
272 int colonne ;
273
274 // compteurs pour les diagonales du systeme :
275 int sup ;
276 int inf ;
277 int sup_new, inf_new ;
278
279 // on boucle sur le plus grand nombre possible de Plm intervenant...
280 for (int k=0 ; k<np_max+1 ; k++)
281 for (int j=0 ; j<nt_max ; j++)
282 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
283
284 ligne = 0 ;
285 colonne = 0 ;
286 sup = 0 ;
287 inf = 0 ;
288
289 for (int zone=0 ; zone<nz ; zone++)
290 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
291 k, source.get_mg()->get_np(zone), base);
292
293 // taille du systeme a resoudre pour ce Plm
294 taille = indic[0] ;
295 for (int zone=1 ; zone<nz ; zone++)
296 taille+=2*indic[zone] ;
297
298 // on verifie que la taille est non-nulle.
299 // cas pouvant a priori se produire...
300
301 if (taille != 0) {
302
303 sec_membre = new Tbl(taille) ;
304 sec_membre->set_etat_qcq() ;
305 for (int l=0 ; l<taille ; l++)
306 sec_membre->set(l) = 0 ;
307
308 systeme = new Matrice(taille, taille) ;
309 systeme->set_etat_qcq() ;
310 for (int l=0 ; l<taille ; l++)
311 for (int c=0 ; c<taille ; c++)
312 systeme->set(l, c) = 0 ;
313
314 //Calcul des nombres quantiques
315 //base_r est donne dans le noyau, sa valeur dans les autres
316 //zones etant inutile.
317
318 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
319
320
321 //--------------------------
322 // NOYAU
323 //---------------------------
324
325 if (indic[0] == 1) {
326 nr = source.get_mg()->get_nr(0) ;
327
328 alpha = mapping.get_alpha()[0] ;
329 // valeur de x^l en 1 :
330 systeme->set(ligne, colonne) = 1. ; /* ligne=0, colonne=0 */
331 // coefficient du Plm dans la solp
332 for (int i=0 ; i<nr ; i++)
333 sec_membre->set(ligne) -= solution_part(0, k, j, i) ; /* ligne=0 */
334
335 ligne++ ; /* ligne=1, colonne=0 */
336 // on prend les derivees que si Plm existe
337 //dans la zone suivante
338
339 if (indic[1] == 1) {
340 // derivee de x^l en 1 :
341 systeme->set(ligne, colonne) = 1./alpha*l_quant ; /* ligne=1, colonne=0 */
342
343 // coefficient de la derivee du Plm dans la solp
344 if (base_r == R_CHEBP)
345 // cas en R_CHEBP
346 for (int i=0 ; i<nr ; i++)
347 sec_membre->set(ligne) -=
348 4*i*i/alpha
349 *solution_part(0, k, j, i) ; /* ligne=1 */
350 else
351 // cas en R_CHEBI
352 for (int i=0 ; i<nr ; i++)
353 sec_membre->set(ligne) -=
354 (2*i+1)*(2*i+1)/alpha
355 *solution_part(0, k, j, i) ; /* ligne=1 */
356
357 // on a au moins un diag inferieure dans ce cas ...
358 inf = 1 ;
359 }
360 colonne++ ; /* ligne=1, colonne=1 */
361 }
362
363 //-----------------------------
364 // COQUILLES
365 //------------------------------
366
367 for (int zone=1 ; zone<nz ; zone++)
368 if (indic[zone] == 1) {
369
370 nr = source.get_mg()->get_nr(zone) ;
371 alpha = mapping.get_alpha()[zone] ;
372 echelle = mapping.get_beta()[zone]/alpha ;
373
374 //Frontiere avec la zone precedente :
375 if (indic[zone-1] == 1) ligne -- ; /* ligne=0, colonne=1 */
376
377 //valeur de (x+echelle)^l en -1 :
378 systeme->set(ligne, colonne) =
379 -pow(echelle-1, double(l_quant)) ; /* ligne=0, colonne=1 */
380
381 // valeur de 1/(x+echelle) ^(l+1) en -1 :
382 systeme->set(ligne, colonne+1) =
383 -1/pow(echelle-1, double(l_quant+1)) ; /* ligne=0, colonne=1, colonne+1=2 */
384
385 // la solution particuliere :
386 for (int i=0 ; i<nr ; i++)
387 if (i%2 == 0)
388 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
389 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=0 */
390
391 // les diagonales :
392 sup_new = colonne+1-ligne ; /* ligne=0, colonne=1, colonne+1-ligne=2 */
393 if (sup_new > sup) sup = sup_new ;
394
395 ligne++ ; /* ligne=1 */
396
397 // on prend les derivees si Plm existe dans la zone
398 // precedente :
399
400 if (indic[zone-1] == 1) {
401 // derivee de (x+echelle)^l en -1 :
402 systeme->set(ligne, colonne) =
403 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ; /* ligne=1, colonne=1 */
404 // derivee de 1/(x+echelle)^(l+1) en -1 :
405 systeme->set(ligne, colonne+1) =
406 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ; /* ligne=1, colonne=1, colonne+1=2 */
407
408
409
410 // la solution particuliere :
411 for (int i=0 ; i<nr ; i++)
412 if (i%2 == 0)
413 sec_membre->set(ligne) -=
414 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
415 else
416 sec_membre->set(ligne) +=
417 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
418
419 // les diagonales :
420 sup_new = colonne+1-ligne ; /* ligne=1, colonne=1, colonne+1-ligne=1 */
421 if (sup_new > sup) sup = sup_new ;
422
423 ligne++ ; /* ligne=2 */
424 }
425
426
427 if(zone < nz-1) {
428
429 // Frontiere avec la zone suivante :
430 //valeur de (x+echelle)^l en 1 :
431 systeme->set(ligne, colonne) =
432 pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
433
434 // valeur de 1/(x+echelle)^(l+1) en 1 :
435 systeme->set(ligne, colonne+1) =
436 1/pow(echelle+1, double(l_quant+1)) ; /* ligne=2, colonne=1, colonne+1=2 */
437
438 // la solution particuliere :
439 for (int i=0 ; i<nr ; i++)
440 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=2 */
441
442 // les diagonales inf :
443 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
444 if (inf_new > inf) inf = inf_new ;
445
446 ligne ++ ; /* ligne=3 */
447
448 // Utilisation des derivees ssi Plm existe dans la
449 //zone suivante :
450 if (indic[zone+1] == 1) {
451
452 //derivee de (x+echelle)^l en 1 :
453 systeme->set(ligne, colonne) =
454 l_quant*pow(echelle+1, double(l_quant-1))/alpha ; /* ligne=3, colonne=1 */
455
456 //derivee de 1/(echelle+x)^(l+1) en 1 :
457 systeme->set(ligne, colonne+1) =
458 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ; /* ligne=3, colonne=1, colonne+1=2 */
459
460 // La solution particuliere :
461 for (int i=0 ; i<nr ; i++)
462 sec_membre->set(ligne) -=
463 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=3 */
464
465 // les diagonales inf :
466 inf_new = ligne-colonne ; /* ligne=3, colonne=1, ligne-colonne=2 */
467 if (inf_new > inf) inf = inf_new ;
468
469 }
470 colonne += 2 ; /* ligne=colonne=3 -> ligne=colonne=1 next block of two */
471 } else {
472
473
474
475 //-------------------------
476 // Falloff outer boundary
477 //---------------------------
478
479 /* ligne=2, colonne=1 */
480
481 // The coefficient for A_n is (k+l)(echelle+1)^l
482 systeme->set(ligne, colonne) =
483 double(k_falloff+l_quant)*pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
484
485 // The coefficient for B_n is (k-(l+1))(echelle+1)^{-(l+1)}
486 systeme->set(ligne, colonne+1) =
487 double(k_falloff-l_quant-1)/pow(echelle+1, double(l_quant+1)) ; /* ligne=2, colonne=1, colonne+1=2 */
488
489 // Here we have -(1+echelle)*F'(x=1)-k*F(x=1)
490 // La solution particuliere :
491 for (int i=0 ; i<nr ; i++)
492 sec_membre->set(ligne) -=
493 (k_falloff+(echelle+1)*i*i)*solution_part(zone, k, j, i) ; /* ligne=2 */
494
495 // les diagonales inf :
496 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
497 if (inf_new > inf) inf = inf_new ;
498
499 }
500 }
501
502
503 //-------------------------
504 // resolution du systeme
505 //--------------------------
506
507 systeme->set_band(sup, inf) ;
508 systeme->set_lu() ;
509
510 Tbl facteurs(systeme->inverse(*sec_membre)) ;
511 int conte = 0 ;
512
513
514 // rangement dans le noyau :
515
516 if (indic[0] == 1) {
517 nr = source.get_mg()->get_nr(0) ;
518 for (int i=0 ; i<nr ; i++)
519 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
520 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
521 conte++ ;
522 }
523
524 // rangement dans les coquilles :
525 for (int zone=1 ; zone<nz ; zone++)
526 if (indic[zone] == 1) {
527 nr = source.get_mg()->get_nr(zone) ;
528 for (int i=0 ; i<nr ; i++)
529 resultat.set(zone, k, j, i) =
530 solution_part(zone, k, j, i)
531 +facteurs(conte)*solution_hom_un(zone, k, j, i)
532 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
533 conte+=2 ;
534 }
535
536 delete sec_membre ;
537 delete systeme ;
538 }
539
540 }
541
542 delete [] indic ;
543
544 return resultat;
545}
546}
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:453
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64