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