LORENE
poisson.C
1/*
2 * Copyright (c) 1999-2001 Philippe Grandclement
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
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 poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson.C,v 1.9 2014/10/13 08:53:29 j_novak Exp $" ;
25
26/*
27 * $Id: poisson.C,v 1.9 2014/10/13 08:53:29 j_novak Exp $
28 * $Log: poisson.C,v $
29 * Revision 1.9 2014/10/13 08:53:29 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.8 2014/10/06 15:16:09 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.7 2013/06/05 15:10:43 j_novak
36 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
37 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
38 *
39 * Revision 1.6 2007/12/13 15:48:46 jl_cornou
40 * *** empty log message ***
41 *
42 * Revision 1.5 2007/12/11 15:28:23 jl_cornou
43 * Jacobi(0,2) polynomials partially implemented
44 *
45 * Revision 1.4 2004/10/05 15:44:21 j_novak
46 * Minor speed enhancements.
47 *
48 * Revision 1.3 2004/02/20 10:55:23 j_novak
49 * The versions dzpuis 5 -> 3 has been improved and polished. Should be
50 * operational now...
51 *
52 * Revision 1.2 2004/02/06 10:53:54 j_novak
53 * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
54 *
55 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
56 * LORENE
57 *
58 * Revision 2.24 2000/07/18 10:23:09 eric
59 * Suppression d'une erreur sans consequence : dans le noyau, remplacement de
60 * beta = mapping.get_alpha()[0] par
61 * beta = mapping.get_beta()[0]
62 * Cette erreur etait sans consequence car beta n'intervient pas dans la
63 * suite des calculs pour le noyau.
64 *
65 * Revision 2.23 2000/05/22 13:43:46 phil
66 * ajout du cas dzpuis =3
67 *
68 * Revision 2.22 2000/02/09 14:40:49 eric
69 * Ajout de l'argument dzpuis a sol_poisson.
70 * (dzpuis n'est plus lu sur le Mtbl_cf).
71 *
72 * Revision 2.21 1999/12/02 14:29:05 eric
73 * Suppression de la base en argument de la version Mtbl_cf.
74 * La version Map se trouve desormais dans le fichier map_af_poisson.C
75 * La version Cmp se trouve desormais dans le fichier cmp_pde.C
76 *
77 * Revision 2.20 1999/11/30 12:56:45 eric
78 * Valeur::base est desormais du type Base_val et non plus Base_val*.
79 *
80 * Revision 2.19 1999/11/24 14:34:25 eric
81 * Accession des membres alpha et beta de Map_af (desormais prives) par
82 * Map_af::get_alpha() et Map_af::get_beta().
83 *
84 * Revision 2.18 1999/10/27 15:47:19 eric
85 * Suppression du membre Cmp::c.
86 *
87 * Revision 2.17 1999/10/14 14:28:07 eric
88 * Methode const.
89 *
90 * Revision 2.16 1999/10/13 15:52:22 eric
91 * Ajout de la base dans l'appel au constructeur de Mtbl_cf.
92 *
93 * Revision 2.15 1999/10/11 16:27:59 phil
94 * suppression du sur echantillonnage
95 *
96 * Revision 2.14 1999/10/11 14:28:50 phil
97 * & -> &&
98 *
99 * Revision 2.13 1999/09/06 16:25:42 phil
100 * ajout de la version Cmp
101 *
102 * Revision 2.12 1999/07/02 15:05:01 phil
103 * *** empty log message ***
104 *
105 * Revision 2.11 1999/06/23 12:35:18 phil
106 * ajout de dzpuis = 2
107 *
108 * Revision 2.10 1999/04/27 13:12:28 phil
109 * *** empty log message ***
110 *
111 * Revision 2.9 1999/04/14 13:56:03 phil
112 * *** empty log message ***
113 *
114 * Revision 2.8 1999/04/14 10:22:28 phil
115 * *** empty log message ***
116 *
117 * Revision 2.7 1999/04/14 10:20:09 phil
118 * *** empty log message ***
119 *
120 * Revision 2.6 1999/04/13 17:18:50 phil
121 * Ajout de la version Valeur
122 *
123 * Revision 2.5 1999/04/12 15:21:37 phil
124 * correction de la construction du systeme
125 *
126 * Revision 2.4 1999/04/12 14:55:16 phil
127 * correction matrices
128 *
129 * Revision 2.3 1999/04/12 14:28:46 phil
130 * *** empty log message ***
131 *
132 * Revision 2.2 1999/04/07 14:56:33 phil
133 * Changement prototypage
134 *
135 * Revision 2.1 1999/04/07 14:40:55 phil
136 * mise au point des includes.
137 *
138 * Revision 2.0 1999/04/07 14:10:55 phil
139 * *** empty log message ***
140 *
141 *
142 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson.C,v 1.9 2014/10/13 08:53:29 j_novak Exp $
143 *
144 */
145
146// Header C :
147#include <cstdlib>
148#include <cmath>
149
150// Headers Lorene :
151#include "matrice.h"
152#include "map.h"
153#include "proto.h"
154#include "type_parite.h"
155
156
157
158 //----------------------------------------------
159 // Version Mtbl_cf
160 //----------------------------------------------
161
162/*
163 *
164 * Solution de l'equation de poisson
165 *
166 * Entree : mapping : le mapping affine
167 * source : les coefficients de la source qui a ete multipliee par
168 * r^4 ou r^2 dans la ZEC.
169 * La base de decomposition doit etre Ylm
170 * dzpuis : exposant de r dans le factor multiplicatif dans la ZEC
171 * Sortie : renvoie les coefficients de la solution dans la meme base de
172 * decomposition (a savoir Ylm)
173 *
174 */
175
176
177
178namespace Lorene {
179Mtbl_cf sol_poisson(const Map_af& mapping, const Mtbl_cf& source, int dzpuis,
180 bool match)
181{
182
183 // Verifications d'usage sur les zones
184 int nz = source.get_mg()->get_nzone() ;
185 assert (nz>1) ;
186 assert ((source.get_mg()->get_type_r(0) == RARE) || (source.get_mg()->get_type_r(0) == FIN)) ;
187 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
188 for (int l=1 ; l<nz-1 ; l++)
189 assert(source.get_mg()->get_type_r(l) == FIN) ;
190
191 assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3) || (dzpuis==5)) ;
192 assert ((!match) || (dzpuis != 5)) ;
193
194 // Bases spectrales
195 const Base_val& base = source.base ;
196
197
198 // donnees sur la zone
199 int nr, nt, np ;
200 int base_r ;
201 double alpha, beta, echelle ;
202 int l_quant, m_quant;
203
204 //Rangement des valeurs intermediaires
205 Tbl *so = 0x0 ;
206 Tbl *sol_hom = 0x0 ;
207 Tbl *sol_part = 0x0 ;
208 Matrice *operateur = 0x0 ;
209 Matrice *nondege = 0x0 ;
210
211
212 // Rangement des solutions, avant raccordement
213 Mtbl_cf solution_part(source.get_mg(), base) ;
214 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
215 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
216 Mtbl_cf resultat(source.get_mg(), base) ;
217
218 solution_part.set_etat_qcq() ;
219 solution_hom_un.set_etat_qcq() ;
220 solution_hom_deux.set_etat_qcq() ;
221 resultat.annule_hard() ;
222 for (int l=0 ; l<nz ; l++) {
223 solution_part.t[l]->set_etat_qcq() ;
224 solution_hom_un.t[l]->set_etat_qcq() ;
225 solution_hom_deux.t[l]->set_etat_qcq() ;
226 }
227 // nbre maximum de point en theta et en phi :
228 int np_max, nt_max ;
229
230 //---------------
231 //-- NOYAU -----
232 //---------------
233
234 nr = source.get_mg()->get_nr(0) ;
235 nt = source.get_mg()->get_nt(0) ;
236 np = source.get_mg()->get_np(0) ;
237
238 nt_max = nt ;
239 np_max = np ;
240
241 alpha = mapping.get_alpha()[0] ;
242 beta = mapping.get_beta()[0] ;
243
244 for (int k=0 ; k<np+1 ; k++)
245 for (int j=0 ; j<nt ; j++)
246 if (nullite_plm(j, nt, k, np, base) == 1)
247 {
248 // calcul des nombres quantiques :
249 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
250 assert( (source.get_mg()->get_type_r(0) == RARE) ||
251 (base_r == R_JACO02) ) ;
252
253 // Construction operateur
254 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
255 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
256
257 //Operateur inversible
258 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
259
260 if (match) {
261 // Calcul de la SH
262 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
263 }
264
265 //Calcul de la SP
266 so = new Tbl(nr) ;
267 so->set_etat_qcq() ;
268 for (int i=0 ; i<nr ; i++)
269 so->set(i) = source(0, k, j, i) ;
270
271 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
272 *so, 0, base_r)) ;
273
274 // Rangement dans les tableaux globaux ;
275
276 for (int i=0 ; i<nr ; i++) {
277 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
278 if (match) {
279 if (base_r == R_JACO02) {
280 solution_hom_un.set(0, k, j, i) = (*sol_hom)(0, i) ;
281 solution_hom_deux.set(0, k, j, i) = (*sol_hom)(1, i) ;
282 }
283 else {
284 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
285 solution_hom_deux.set(0, k, j, i) = 0. ;
286 }
287 }
288 }
289
290
291
292 delete operateur ;
293 delete nondege ;
294 delete so ;
295 if (match) delete sol_hom ;
296 delete sol_part ;
297 }
298
299
300 //---------------
301 //-- ZEC -----
302 //---------------
303
304 nr = source.get_mg()->get_nr(nz-1) ;
305 nt = source.get_mg()->get_nt(nz-1) ;
306 np = source.get_mg()->get_np(nz-1) ;
307
308 if (np > np_max) np_max = np ;
309 if (nt > nt_max) nt_max = nt ;
310
311 alpha = mapping.get_alpha()[nz-1] ;
312 beta = mapping.get_beta()[nz-1] ;
313
314 for (int k=0 ; k<np+1 ; k++)
315 for (int j=0 ; j<nt ; j++)
316 if (nullite_plm(j, nt, k, np, base) == 1)
317 {
318
319 // calcul des nombres quantiques :
320 donne_lm(nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
321
322 //Construction operateur
323 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis,
324 base_r)) ;
325 (*operateur) = combinaison(*operateur, l_quant, 0., dzpuis, base_r) ;
326 // Operateur inversible
327 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0.,
328 dzpuis, base_r)) ;
329 // Calcul de la SH
330 if (match)
331 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
332
333 // Calcul de la SP
334 so = new Tbl(nr) ;
335 so->set_etat_qcq() ;
336 for (int i=0 ; i<nr ; i++)
337 so->set(i) = source(nz-1, k, j, i) ;
338 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
339 *so, dzpuis, base_r)) ;
340
341 // Rangement
342 for (int i=0 ; i<nr ; i++) {
343 solution_part.set(nz-1, k, j, i) = (*sol_part)(i) ;
344 if (match) {
345 solution_hom_un.set(nz-1, k, j, i) = (*sol_hom)(i) ;
346 solution_hom_deux.set(nz-1, k, j, i) = 0. ;
347 }
348 }
349
350 delete operateur ;
351 delete nondege ;
352 delete so ;
353 if (match) delete sol_hom ;
354 delete sol_part ;
355 }
356
357 //---------------
358 //- COQUILLES ---
359 //---------------
360
361 for (int zone=1 ; zone<nz-1 ; zone++) {
362 nr = source.get_mg()->get_nr(zone) ;
363 nt = source.get_mg()->get_nt(zone) ;
364 np = source.get_mg()->get_np(zone) ;
365
366 if (np > np_max) np_max = np ;
367 if (nt > nt_max) nt_max = nt ;
368
369 alpha = mapping.get_alpha()[zone] ;
370 beta = mapping.get_beta()[zone] ;
371
372 for (int k=0 ; k<np+1 ; k++)
373 for (int j=0 ; j<nt ; j++)
374 if (nullite_plm(j, nt, k, np, base) == 1)
375 {
376 // calcul des nombres quantiques :
377 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
378
379 // Construction de l'operateur
380 operateur = new Matrice(laplacien_mat
381 (nr, l_quant, beta/alpha, 0, base_r)) ;
382
383 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
384 base_r) ;
385
386 // Operateur inversible
387 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
388 beta/alpha, 0, base_r)) ;
389 if (match) {
390 // Calcul DES DEUX SH
391 sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
392 }
393
394 // Calcul de la SP
395 so = new Tbl(nr) ;
396 so->set_etat_qcq() ;
397 for (int i=0 ; i<nr ; i++)
398 so->set(i) = source(zone, k, j, i) ;
399
400 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
401 beta, *so, 0, base_r)) ;
402
403
404 // Rangement
405 for (int i=0 ; i<nr ; i++) {
406 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
407 if (match) {
408 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
409 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
410 }
411 }
412
413
414 delete operateur ;
415 delete nondege ;
416 delete so ;
417 if (match) delete sol_hom ;
418 delete sol_part ;
419 }
420 }
421
422
423 if (match) {
424
425 //-------------------------------------------
426 // On est parti pour le raccord des solutions
427 //-------------------------------------------
428 // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
429 // intervient dans le developpement de cette zone.
430 int * indic = new int [nz] ;
431 int taille ;
432 Tbl *sec_membre ; // termes constants du systeme
433 Matrice *systeme ; //le systeme a resoudre
434
435 // des compteurs pour le remplissage du systeme
436 int ligne ;
437 int colonne ;
438
439 // compteurs pour les diagonales du systeme :
440 int sup ;
441 int inf ;
442 int sup_new, inf_new ;
443
444 // on boucle sur le plus grand nombre possible de Plm intervenant...
445 for (int k=0 ; k<np_max+1 ; k++)
446 for (int j=0 ; j<nt_max ; j++)
447 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
448
449 ligne = 0 ;
450 colonne = 0 ;
451 sup = 0 ;
452 inf = 0 ;
453
454 for (int zone=0 ; zone<nz ; zone++)
455 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
456 k, source.get_mg()->get_np(zone), base);
457
458 // taille du systeme a resoudre pour ce Plm
459 taille = indic[nz-1]+indic[0] ;
460 for (int zone=1 ; zone<nz-1 ; zone++)
461 taille+=2*indic[zone] ;
462
463 // on verifie que la taille est non-nulle.
464 // cas pouvant a priori se produire...
465
466 if (taille != 0) {
467
468 sec_membre = new Tbl(taille) ;
469 sec_membre->set_etat_qcq() ;
470 for (int l=0 ; l<taille ; l++)
471 sec_membre->set(l) = 0 ;
472
473 systeme = new Matrice(taille, taille) ;
474 systeme->set_etat_qcq() ;
475 for (int l=0 ; l<taille ; l++)
476 for (int c=0 ; c<taille ; c++)
477 systeme->set(l, c) = 0 ;
478
479 //Calcul des nombres quantiques
480 //base_r est donne dans le noyau, sa valeur dans les autres
481 //zones etant inutile.
482
483 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
484
485
486 //--------------------------
487 // NOYAU
488 //---------------------------
489
490 if (indic[0] == 1) {
491 nr = source.get_mg()->get_nr(0) ;
492
493 alpha = mapping.get_alpha()[0] ;
494 // valeur de x^l en 1 :
495 systeme->set(ligne, colonne) = 1. ;
496 // coefficient du Plm dans la solp
497 for (int i=0 ; i<nr ; i++)
498 sec_membre->set(ligne) -= solution_part(0, k, j, i) ;
499
500 ligne++ ;
501 // on prend les derivees que si Plm existe
502 //dans la zone suivante
503
504
505
506// Grosses modifications en perspective
507
508 if (indic[1] == 1) {
509 // derivee de x^l en 1 :
510 systeme->set(ligne, colonne) = 1./alpha*l_quant ;
511
512 // coefficient de la derivee du Plm dans la solp
513 if (base_r == R_CHEBP)
514 // cas en R_CHEBP
515 for (int i=0 ; i<nr ; i++)
516 sec_membre->set(ligne) -=
517 4*i*i/alpha
518 *solution_part(0, k, j, i) ;
519 else
520 // cas en R_CHEBI
521 for (int i=0 ; i<nr ; i++)
522 sec_membre->set(ligne) -=
523 (2*i+1)*(2*i+1)/alpha
524 *solution_part(0, k, j, i) ;
525
526 // on a au moins un diag inferieure dans ce cas ...
527 inf = 1 ;
528 }
529 colonne++ ;
530 }
531
532 //-----------------------------
533 // COQUILLES
534 //------------------------------
535
536 for (int zone=1 ; zone<nz-1 ; zone++)
537 if (indic[zone] == 1) {
538
539 nr = source.get_mg()->get_nr(zone) ;
540 alpha = mapping.get_alpha()[zone] ;
541 echelle = mapping.get_beta()[zone]/alpha ;
542
543 //Frontiere avec la zone precedente :
544 if (indic[zone-1] == 1) ligne -- ;
545
546 //valeur de (x+echelle)^l en -1 :
547 systeme->set(ligne, colonne) =
548 -pow(echelle-1, double(l_quant)) ;
549
550 // valeur de 1/(x+echelle) ^(l+1) en -1 :
551 systeme->set(ligne, colonne+1) =
552 -1/pow(echelle-1, double(l_quant+1)) ;
553
554 // la solution particuliere :
555 for (int i=0 ; i<nr ; i++)
556 if (i%2 == 0)
557 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
558 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
559
560 // les diagonales :
561 sup_new = colonne+1-ligne ;
562 if (sup_new > sup) sup = sup_new ;
563
564 ligne++ ;
565
566 // on prend les derivees si Plm existe dans la zone
567 // precedente :
568
569 if (indic[zone-1] == 1) {
570 // derivee de (x+echelle)^l en -1 :
571 systeme->set(ligne, colonne) =
572 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ;
573 // derivee de 1/(x+echelle)^(l+1) en -1 :
574 systeme->set(ligne, colonne+1) =
575 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ;
576
577
578
579 // la solution particuliere :
580 for (int i=0 ; i<nr ; i++)
581 if (i%2 == 0)
582 sec_membre->set(ligne) -=
583 i*i/alpha*solution_part(zone, k, j, i) ;
584 else
585 sec_membre->set(ligne) +=
586 i*i/alpha*solution_part(zone, k, j, i) ;
587
588 // les diagonales :
589 sup_new = colonne+1-ligne ;
590 if (sup_new > sup) sup = sup_new ;
591
592 ligne++ ;
593 }
594
595 // Frontiere avec la zone suivante :
596 //valeur de (x+echelle)^l en 1 :
597 systeme->set(ligne, colonne) =
598 pow(echelle+1, double(l_quant)) ;
599
600 // valeur de 1/(x+echelle)^(l+1) en 1 :
601 systeme->set(ligne, colonne+1) =
602 1/pow(echelle+1, double(l_quant+1)) ;
603
604 // la solution particuliere :
605 for (int i=0 ; i<nr ; i++)
606 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
607
608 // les diagonales inf :
609 inf_new = ligne-colonne ;
610 if (inf_new > inf) inf = inf_new ;
611
612 ligne ++ ;
613
614 // Utilisation des derivees ssi Plm existe dans la
615 //zone suivante :
616 if (indic[zone+1] == 1) {
617
618 //derivee de (x+echelle)^l en 1 :
619 systeme->set(ligne, colonne) =
620 l_quant*pow(echelle+1, double(l_quant-1))/alpha ;
621
622 //derivee de 1/(echelle+x)^(l+1) en 1 :
623 systeme->set(ligne, colonne+1) =
624 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ;
625
626 // La solution particuliere :
627 for (int i=0 ; i<nr ; i++)
628 sec_membre->set(ligne) -=
629 i*i/alpha*solution_part(zone, k, j, i) ;
630
631 // les diagonales inf :
632 inf_new = ligne-colonne ;
633 if (inf_new > inf) inf = inf_new ;
634
635 }
636 colonne += 2 ;
637 }
638
639
640 //--------------------------------
641 // ZEC
642 //---------------------------------
643
644 if (indic[nz-1] == 1) {
645
646 nr = source.get_mg()->get_nr(nz-1) ;
647
648
649 alpha = mapping.get_alpha()[nz-1] ;
650
651 if (indic[nz-2] == 1) ligne -- ;
652
653 //valeur de (x-1)^(l+1) en -1 :
654 systeme->set(ligne, colonne) = -pow(-2, double(l_quant+1)) ;
655 //solution particuliere :
656 for (int i=0 ; i<nr ; i++)
657 if (i%2 == 0)
658 sec_membre->set(ligne) += solution_part(nz-1, k, j, i) ;
659 else sec_membre->set(ligne) -= solution_part(nz-1, k, j, i) ;
660
661 //on prend les derivees ssi Plm existe dans la zone precedente :
662 if (indic[nz-2] == 1) {
663
664 //derivee de (x-1)^(l+1) en -1 :
665 systeme->set(ligne+1, colonne) =
666 alpha*(l_quant+1)*pow(-2., double(l_quant+2)) ;
667
668 // Solution particuliere :
669 for (int i=0 ; i<nr ; i++)
670 if (i%2 == 0)
671 sec_membre->set(ligne+1) -=
672 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
673 else
674 sec_membre->set(ligne+1) +=
675 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
676
677 //les diags :
678 if (sup == 0) sup = 1 ;
679 }
680 }
681
682 //-------------------------
683 // resolution du systeme
684 //--------------------------
685
686 systeme->set_band(sup, inf) ;
687 systeme->set_lu() ;
688
689 Tbl facteurs(systeme->inverse(*sec_membre)) ;
690 int conte = 0 ;
691
692
693 // rangement dans le noyau :
694
695 if (indic[0] == 1) {
696 nr = source.get_mg()->get_nr(0) ;
697 for (int i=0 ; i<nr ; i++)
698 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
699 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
700 conte++ ;
701 }
702
703 // rangement dans les coquilles :
704 for (int zone=1 ; zone<nz-1 ; zone++)
705 if (indic[zone] == 1) {
706 nr = source.get_mg()->get_nr(zone) ;
707 for (int i=0 ; i<nr ; i++)
708 resultat.set(zone, k, j, i) =
709 solution_part(zone, k, j, i)
710 +facteurs(conte)*solution_hom_un(zone, k, j, i)
711 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
712 conte+=2 ;
713 }
714
715 //rangement dans la ZEC :
716 if (indic[nz-1] == 1) {
717 nr = source.get_mg()->get_nr(nz-1) ;
718 for (int i=0 ; i<nr ; i++)
719 resultat.set(nz-1, k, j, i) =
720 solution_part(nz-1, k, j, i)
721 +facteurs(conte)*solution_hom_un(nz-1, k, j, i) ;
722 }
723 delete sec_membre ;
724 delete systeme ;
725 }
726
727 }
728
729 delete [] indic ;
730
731 } // End of the case where the matching has to be done
732 else { //Only a particular solution is given in each domain
733
734 for (int zone = 0; zone < nz; zone++)
735 for (int k=0 ; k<source.get_mg()->get_np(zone)+1 ; k++)
736 for (int j=0 ; j<source.get_mg()->get_nt(zone) ; j++)
737 if (nullite_plm(j,source.get_mg()->get_nt(zone) ,
738 k, source.get_mg()->get_np(zone), base) == 1)
739 for (int i=0; i<source.get_mg()->get_nr(zone); i++)
740 resultat.set(zone, k, j, i) = solution_part(zone, k, j, i) ;
741
742 }
743
744 return resultat;
745}
746}
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_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64