LORENE
poisson_frontiere.C
1/*
2 * Copyright (c) 2000-2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: poisson_frontiere.C,v 1.6 2016/12/05 16:18:10 j_novak Exp $
27 * $Log: poisson_frontiere.C,v $
28 * Revision 1.6 2016/12/05 16:18:10 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.5 2014/10/13 08:53:29 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.4 2014/10/06 15:16:09 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.3 2004/11/23 12:50:44 f_limousin
38 * Intoduce function poisson_dir_neu(...) to solve a scalar poisson
39 * equation with a mixed boundary condition (Dirichlet + Neumann).
40 *
41 * Revision 1.2 2004/09/08 15:12:16 f_limousin
42 * Delete some assert.
43 *
44 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
45 * LORENE
46 *
47 * Revision 2.4 2000/05/22 16:03:32 phil
48 * ajout du cas dzpuis = 3
49 *
50 * Revision 2.3 2000/05/15 15:46:35 phil
51 * *** empty log message ***
52 *
53 * Revision 2.2 2000/04/27 12:28:56 phil
54 * correction pour le raccord des differents domaines
55 *
56 * Revision 2.1 2000/03/20 13:06:32 phil
57 * *** empty log message ***
58 *
59 * Revision 2.0 2000/03/17 17:24:49 phil
60 * *** empty log message ***
61 *
62 *
63 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_frontiere.C,v 1.6 2016/12/05 16:18:10 j_novak Exp $
64 *
65 */
66
67
68// Header C :
69#include <cstdlib>
70#include <cmath>
71
72// Headers Lorene :
73#include "matrice.h"
74#include "tbl.h"
75#include "mtbl_cf.h"
76#include "map.h"
77#include "base_val.h"
78#include "proto.h"
79#include "type_parite.h"
80#include "utilitaires.h"
81#include "valeur.h"
82
83
84
85 //----------------------------------------------
86 // Version Mtbl_cf
87 //----------------------------------------------
88
89/*
90 *
91 * Solution de l'equation de poisson avec Boundary condition a
92 * l'interieur d'une coquille.
93 *
94 * Entree : mapping : le mapping affine
95 * source : les coefficients de la source qui a ete multipliee par
96 * r^4 ou r^2 dans la ZEC.
97 * La base de decomposition doit etre Ylm
98 * limite : la CL (fonction angulaire) sur une frontiere spherique
99 * type_raccord : 1 pour Dirichlet et 2 pour Neumann
100 * num_front : indique la frontiere ou on impose la CL : 1 pour la
101 * frontiere situee entre le domain 1 et 2.
102 * dzpuis : exposant de r dans le factor multiplicatif dans la ZEC
103 * Sortie : renvoie les coefficients de la solution dans la meme base de
104 * decomposition (a savoir Ylm)
105 *
106 */
107
108
109
110namespace Lorene {
111Mtbl_cf sol_poisson_frontiere(const Map_af& mapping, const Mtbl_cf& source,
112 const Mtbl_cf& limite, int type_raccord, int num_front, int dzpuis, double fact_dir, double fact_neu)
113
114{
115
116 // Verifications d'usage sur les zones
117 int nz = source.get_mg()->get_nzone() ;
118 assert (nz>1) ;
119 assert ((num_front>=0) && (num_front<nz-2)) ;
120 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
121 for (int l=num_front+1 ; l<nz-1 ; l++)
122 assert(source.get_mg()->get_type_r(l) == FIN) ;
123
124 assert (source.get_etat() != ETATNONDEF) ;
125 assert (limite.get_etat() != ETATNONDEF) ;
126
127 assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
128 assert ((type_raccord == 1) || (type_raccord == 2)|| (type_raccord == 3)) ;
129
130 // Bases spectrales
131 const Base_val& base = source.base ;
132
133 // donnees sur la zone
134 int nr, nt, np ;
135 int base_r ;
136 double alpha, beta, echelle ;
137 int l_quant, m_quant;
138
139 //Rangement des valeurs intermediaires
140 Tbl *so ;
141 Tbl *sol_hom ;
142 Tbl *sol_part ;
143 Matrice *operateur ;
144 Matrice *nondege ;
145
146
147 // Rangement des solutions, avant raccordement
148 Mtbl_cf solution_part(source.get_mg(), base) ;
149 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
150 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
151 Mtbl_cf resultat(source.get_mg(), base) ;
152
153 solution_part.set_etat_qcq() ;
154 solution_hom_un.set_etat_qcq() ;
155 solution_hom_deux.set_etat_qcq() ;
156 resultat.set_etat_qcq() ;
157
158 for (int l=0 ; l<nz ; l++) {
159 solution_part.t[l]->set_etat_qcq() ;
160 solution_hom_un.t[l]->set_etat_qcq() ;
161 solution_hom_deux.t[l]->set_etat_qcq() ;
162 resultat.t[l]->set_etat_qcq() ;
163 for (int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
164 for (int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
165 for (int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
166 resultat.set(l, k, j, i) = 0 ;
167 }
168
169 // nbre maximum de point en theta et en phi :
170 int np_max = 0 ;
171 int nt_max = 0 ;
172
173
174 //---------------
175 //-- ZEC -----
176 //---------------
177
178 nr = source.get_mg()->get_nr(nz-1) ;
179 nt = source.get_mg()->get_nt(nz-1) ;
180 np = source.get_mg()->get_np(nz-1) ;
181
182 if (np > np_max) np_max = np ;
183 if (nt > nt_max) nt_max = nt ;
184
185 alpha = mapping.get_alpha()[nz-1] ;
186 beta = mapping.get_beta()[nz-1] ;
187
188 for (int k=0 ; k<np+1 ; k++)
189 for (int j=0 ; j<nt ; j++)
190 if (nullite_plm(j, nt, k, np, base) == 1)
191 {
192
193 // calcul des nombres quantiques :
194 donne_lm(nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
195
196 //Construction operateur
197 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis,
198 base_r)) ;
199 (*operateur) = combinaison(*operateur, l_quant, 0., dzpuis, base_r) ;
200
201 // Operateur inversible
202 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0.,
203 dzpuis, base_r)) ;
204
205 // Calcul de la SH
206 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
207
208 // Calcul de la SP
209 so = new Tbl(nr) ;
210 so->set_etat_qcq() ;
211 for (int i=0 ; i<nr ; i++)
212 so->set(i) = source(nz-1, k, j, i) ;
213 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
214 *so, dzpuis, base_r)) ;
215
216 // Rangement
217 for (int i=0 ; i<nr ; i++) {
218 solution_part.set(nz-1, k, j, i) = (*sol_part)(i) ;
219 solution_hom_un.set(nz-1, k, j, i) = (*sol_hom)(i) ;
220 solution_hom_deux.set(nz-1, k, j, i) = 0. ;
221 }
222
223 delete operateur ;
224 delete nondege ;
225 delete so ;
226 delete sol_hom ;
227 delete sol_part ;
228 }
229
230 //---------------
231 //- COQUILLES ---
232 //---------------
233
234 for (int zone=num_front+1 ; zone<nz-1 ; zone++) {
235 nr = source.get_mg()->get_nr(zone) ;
236 nt = source.get_mg()->get_nt(zone) ;
237 np = source.get_mg()->get_np(zone) ;
238
239 if (np > np_max) np_max = np ;
240 if (nt > nt_max) nt_max = nt ;
241
242 alpha = mapping.get_alpha()[zone] ;
243 beta = mapping.get_beta()[zone] ;
244
245 for (int k=0 ; k<np+1 ; k++)
246 for (int j=0 ; j<nt ; j++)
247 if (nullite_plm(j, nt, k, np, base) == 1)
248 {
249 // calcul des nombres quantiques :
250 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
251
252 // Construction de l'operateur
253 operateur = new Matrice(laplacien_mat
254 (nr, l_quant, beta/alpha, 0, base_r)) ;
255
256 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
257 base_r) ;
258
259 // Operateur inversible
260 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
261 beta/alpha, 0, base_r)) ;
262
263 // Calcul DES DEUX SH
264 sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
265
266 // Calcul de la SP
267 so = new Tbl(nr) ;
268 so->set_etat_qcq() ;
269 for (int i=0 ; i<nr ; i++)
270 so->set(i) = source(zone, k, j, i) ;
271
272 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
273 beta, *so, 0, base_r)) ;
274
275
276 // Rangement
277 for (int i=0 ; i<nr ; i++) {
278 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
279 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
280 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
281 }
282
283
284 delete operateur ;
285 delete nondege ;
286 delete so ;
287 delete sol_hom ;
288 delete sol_part ;
289 }
290 }
291
292 //-------------------------------------------
293 // On est parti pour imposer la boundary
294 //-------------------------------------------
295
296 nr = source.get_mg()->get_nr(num_front+1) ;
297 nt = source.get_mg()->get_nt(num_front+1) ;
298 np = source.get_mg()->get_np(num_front+1) ;
299 double facteur ;
300 double somme ;
301
302 alpha = mapping.get_alpha()[num_front+1] ;
303 beta = mapping.get_beta()[num_front+1] ;
304 echelle = beta/alpha ;
305
306 for (int k=0 ; k<np+1 ; k++)
307 for (int j=0 ; j<nt ; j++)
308 if (nullite_plm(j, nt, k, np, base) == 1)
309 {
310 // calcul des nombres quantiques :
311 donne_lm(nz, num_front+1, j, k, base, m_quant, l_quant, base_r) ;
312
313 switch (type_raccord) {
314 case 1 :
315 // Conditions de raccord type Dirichlet :
316 // Pour la sp :
317 somme = 0 ;
318 for (int i=0 ; i<nr ; i++)
319 if (i%2 == 0)
320 somme += solution_part(num_front+1, k, j, i) ;
321 else
322 somme -= solution_part(num_front+1, k, j, i) ;
323 facteur = (limite(num_front, k, j, 0)-somme)
324 * pow(echelle-1, l_quant+1) ;
325
326 for (int i=0 ; i<nr ; i++)
327 solution_part.set(num_front+1, k, j, i) +=
328 facteur*solution_hom_deux(num_front+1, k, j, i) ;
329
330 // pour la solution homogene :
331 facteur = - pow(echelle-1, 2*l_quant+1) ;
332 for (int i=0 ; i<nr ; i++)
333 solution_hom_un.set(num_front+1, k, j, i) +=
334 facteur*solution_hom_deux(num_front+1, k, j, i) ;
335 break ;
336
337
338 case 2 :
339 //Conditions de raccord de type Neumann :
340 // Pour la sp :
341 somme = 0 ;
342 for (int i=0 ; i<nr ; i++)
343 if (i%2 == 0)
344 somme -= i*i/alpha*solution_part(num_front+1, k, j, i) ;
345 else
346 somme += i*i/alpha*solution_part(num_front+1, k, j, i) ;
347 facteur = (somme-limite(num_front, k, j, 0))
348 * alpha*pow(echelle-1, l_quant+2)/(l_quant+1) ;
349 for (int i=0 ; i<nr ; i++)
350 solution_part.set(num_front+1, k, j, i) +=
351 facteur*solution_hom_deux(num_front+1, k, j, i) ;
352
353 // pour la solution homogene :
354 facteur = pow(echelle-1, 2*l_quant+1)*l_quant/(l_quant+1) ;
355 for (int i=0 ; i<nr ; i++)
356 solution_hom_un.set(num_front+1, k, j, i) +=
357 facteur*solution_hom_deux(num_front+1, k, j, i) ;
358 break ;
359
360 case 3 :
361 // Conditions de raccord type Dirichlet-Neumann :
362 somme = 0 ;
363 for (int i=0 ; i<nr ; i++)
364 if (i%2 == 0)
365 somme += solution_part(num_front+1, k, j, i) *
366 fact_dir - fact_neu *
367 i*i/alpha*solution_part(num_front+1, k, j, i) ;
368 else
369 somme += - solution_part(num_front+1, k, j, i) *
370 fact_dir + fact_neu *
371 i*i/alpha*solution_part(num_front+1, k, j, i) ;
372
373 double somme2 ;
374 somme2 = fact_dir * pow(echelle-1, -l_quant-1) -
375 fact_neu/alpha*pow(echelle-1, -l_quant-2)*(l_quant+1) ;
376
377 facteur = (limite(num_front, k, j, 0)- somme) / somme2 ;
378
379 for (int i=0 ; i<nr ; i++)
380 solution_part.set(num_front+1, k, j, i) +=
381 facteur*solution_hom_deux(num_front+1, k, j, i) ;
382
383 // pour la solution homogene :
384 double somme1 ;
385 somme1 = fact_dir * pow(echelle-1, l_quant) +
386 fact_neu / alpha * pow(echelle-1, l_quant-1) *
387 l_quant ;
388 facteur = - somme1 / somme2 ;
389 for (int i=0 ; i<nr ; i++)
390 solution_hom_un.set(num_front+1, k, j, i) +=
391 facteur*solution_hom_deux(num_front+1, k, j, i) ;
392
393 break ;
394
395 default :
396 cout << "Diantres nous ne devrions pas etre ici ! " << endl ;
397 abort() ;
398 break ;
399 }
400
401 // Securite.
402 for (int i=0 ; i<nr ; i++)
403 solution_hom_deux.set(num_front+1, k, j, i) = 0 ;
404 }
405
406
407 //-------------------------------------------
408 // On est parti pour le raccord des solutions
409 //-------------------------------------------
410
411 // On
412 // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
413 // intervient dans le developpement de cette zone.
414 int * indic = new int [nz] ;
415 int taille ;
416 Tbl *sec_membre ; // termes constants du systeme
417 Matrice *systeme ; //le systeme a resoudre
418
419 // des compteurs pour le remplissage du systeme
420 int ligne ;
421 int colonne ;
422
423 // compteurs pour les diagonales du systeme :
424 int sup ;
425 int inf ;
426 int sup_new, inf_new ;
427
428 // on boucle sur le plus grand nombre possible de Plm intervenant...
429 for (int k=0 ; k<np_max+1 ; k++)
430 for (int j=0 ; j<nt_max ; j++)
431 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
432
433 ligne = 0 ;
434 colonne = 0 ;
435 sup = 0 ;
436 inf = 0 ;
437
438 for (int zone=num_front+1 ; zone<nz ; zone++)
439 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
440 k, source.get_mg()->get_np(zone), base);
441
442 // taille du systeme a resoudre pour ce Plm
443 taille = indic[nz-1]+indic[num_front+1] ;
444 for (int zone=num_front+2 ; zone<nz-1 ; zone++)
445 taille+=2*indic[zone] ;
446
447 // on verifie que la taille est non-nulle.
448 // cas pouvant a priori se produire...
449
450 if (taille != 0) {
451
452 sec_membre = new Tbl(taille) ;
453 sec_membre->set_etat_qcq() ;
454 for (int l=0 ; l<taille ; l++)
455 sec_membre->set(l) = 0 ;
456
457 systeme = new Matrice(taille, taille) ;
458 systeme->set_etat_qcq() ;
459 for (int l=0 ; l<taille ; l++)
460 for (int c=0 ; c<taille ; c++)
461 systeme->set(l, c) = 0 ;
462
463 //Calcul des nombres quantiques
464 //base_r est donne dans le noyau, sa valeur dans les autres
465 //zones etant inutile.
466
467 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
468
469
470 //----------------------------------
471 // COQUILLE LA + INTERNE
472 //------------------------------------
473
474 if (indic[num_front+1] == 1) {
475 nr = source.get_mg()->get_nr(num_front+1) ;
476
477 alpha = mapping.get_alpha()[num_front+1] ;
478 beta = mapping.get_beta()[num_front+1] ;
479 echelle = beta/alpha ;
480
481 // valeur de la solhomogene en 1 :
482 somme = 0 ;
483 for (int i=0 ; i<nr ; i++)
484 somme += solution_hom_un (num_front+1, k, j, i) ;
485 systeme->set(ligne, colonne) = somme ;
486
487 // coefficient du Plm dans la solp
488 for (int i=0 ; i<nr ; i++)
489 sec_membre->set(ligne) -= solution_part(num_front+1, k, j, i) ;
490
491 ligne++ ;
492 // on prend les derivees que si Plm existe
493 //dans la zone suivante
494
495 if (indic[num_front+1] == 1) {
496
497 // derivee de la solhomogene en 1 :
498 somme = 0 ;
499 for (int i=0 ; i<nr ; i++)
500 somme += i*i/alpha*
501 solution_hom_un(num_front+1, k, j, i) ;
502 systeme->set(ligne, colonne) = somme ;
503
504 // coefficient de la derivee du Plm dans la solp
505 for (int i=0 ; i<nr ; i++)
506 sec_membre->set(ligne) -= i*i/alpha
507 *solution_part(num_front+1, k, j, i) ;
508
509 // on a au moins un diag inferieure dans ce cas ...
510 inf = 1 ;
511 }
512 colonne++ ;
513 }
514
515 //-----------------------------
516 // COQUILLES "normales"
517 //------------------------------
518
519 for (int zone=num_front+2 ; zone<nz-1 ; zone++)
520 if (indic[zone] == 1) {
521
522 nr = source.get_mg()->get_nr(zone) ;
523 alpha = mapping.get_alpha()[zone] ;
524 echelle = mapping.get_beta()[zone]/alpha ;
525
526 //Frontiere avec la zone precedente :
527 if (indic[zone-1] == 1) ligne -- ;
528
529 //valeur de (x+echelle)^l en -1 :
530 systeme->set(ligne, colonne) =
531 -pow(echelle-1, double(l_quant)) ;
532
533 // valeur de 1/(x+echelle) ^(l+1) en -1 :
534 systeme->set(ligne, colonne+1) =
535 -1/pow(echelle-1, double(l_quant+1)) ;
536
537 // la solution particuliere :
538 for (int i=0 ; i<nr ; i++)
539 if (i%2 == 0)
540 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
541 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
542
543 // les diagonales :
544 sup_new = colonne+1-ligne ;
545 if (sup_new > sup) sup = sup_new ;
546
547 ligne++ ;
548
549 // on prend les derivees si Plm existe dans la zone
550 // precedente :
551
552 if (indic[zone-1] == 1) {
553 // derivee de (x+echelle)^l en -1 :
554 systeme->set(ligne, colonne) =
555 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ;
556 // derivee de 1/(x+echelle)^(l+1) en -1 :
557 systeme->set(ligne, colonne+1) =
558 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ;
559
560
561
562 // la solution particuliere :
563 for (int i=0 ; i<nr ; i++)
564 if (i%2 == 0)
565 sec_membre->set(ligne) -=
566 i*i/alpha*solution_part(zone, k, j, i) ;
567 else
568 sec_membre->set(ligne) +=
569 i*i/alpha*solution_part(zone, k, j, i) ;
570
571 // les diagonales :
572 sup_new = colonne+1-ligne ;
573 if (sup_new > sup) sup = sup_new ;
574
575 ligne++ ;
576 }
577
578 // Frontiere avec la zone suivante :
579 //valeur de (x+echelle)^l en 1 :
580 systeme->set(ligne, colonne) =
581 pow(echelle+1, double(l_quant)) ;
582
583 // valeur de 1/(x+echelle)^(l+1) en 1 :
584 systeme->set(ligne, colonne+1) =
585 1/pow(echelle+1, double(l_quant+1)) ;
586
587 // la solution particuliere :
588 for (int i=0 ; i<nr ; i++)
589 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
590
591 // les diagonales inf :
592 inf_new = ligne-colonne ;
593 if (inf_new > inf) inf = inf_new ;
594
595 ligne ++ ;
596
597 // Utilisation des derivees ssi Plm existe dans la
598 //zone suivante :
599 if (indic[zone+1] == 1) {
600
601 //derivee de (x+echelle)^l en 1 :
602 systeme->set(ligne, colonne) =
603 l_quant*pow(echelle+1, double(l_quant-1))/alpha ;
604
605 //derivee de 1/(echelle+x)^(l+1) en 1 :
606 systeme->set(ligne, colonne+1) =
607 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ;
608
609 // La solution particuliere :
610 for (int i=0 ; i<nr ; i++)
611 sec_membre->set(ligne) -=
612 i*i/alpha*solution_part(zone, k, j, i) ;
613
614 // les diagonales inf :
615 inf_new = ligne-colonne ;
616 if (inf_new > inf) inf = inf_new ;
617
618 }
619 colonne += 2 ;
620 }
621
622
623 //--------------------------------
624 // ZEC
625 //---------------------------------
626
627 if (indic[nz-1] == 1) {
628
629 nr = source.get_mg()->get_nr(nz-1) ;
630
631
632 alpha = mapping.get_alpha()[nz-1] ;
633
634 if (indic[nz-2] == 1) ligne -- ;
635
636 //valeur de (x-1)^(l+1) en -1 :
637 systeme->set(ligne, colonne) = -pow(-2, double(l_quant+1)) ;
638 //solution particuliere :
639 for (int i=0 ; i<nr ; i++)
640 if (i%2 == 0)
641 sec_membre->set(ligne) += solution_part(nz-1, k, j, i) ;
642 else sec_membre->set(ligne) -= solution_part(nz-1, k, j, i) ;
643
644 //on prend les derivees ssi Plm existe dans la zone precedente :
645 if (indic[nz-2] == 1) {
646
647 //derivee de (x-1)^(l+1) en -1 :
648 systeme->set(ligne+1, colonne) =
649 alpha*(l_quant+1)*pow(-2., double(l_quant+2)) ;
650
651 // Solution particuliere :
652 for (int i=0 ; i<nr ; i++)
653 if (i%2 == 0)
654 sec_membre->set(ligne+1) -=
655 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
656 else
657 sec_membre->set(ligne+1) +=
658 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
659
660 //les diags :
661 if (sup == 0) sup = 1 ;
662 }
663 }
664
665 //-------------------------
666 // resolution du systeme
667 //--------------------------
668
669 systeme->set_band(sup, inf) ;
670 systeme->set_lu() ;
671
672 Tbl facteurs(systeme->inverse(*sec_membre)) ;
673 int conte = 0 ;
674
675
676 // rangement dans la coquille interne :
677
678 if (indic[num_front+1] == 1) {
679 nr = source.get_mg()->get_nr(num_front+1) ;
680 for (int i=0 ; i<nr ; i++)
681 resultat.set(num_front+1, k, j, i) =
682 solution_part(num_front+1, k, j, i)
683 +facteurs(conte)*solution_hom_un(num_front+1, k, j, i) ;
684 conte++ ;
685 }
686
687 // rangement dans les coquilles :
688 for (int zone=num_front+2 ; zone<nz-1 ; zone++)
689 if (indic[zone] == 1) {
690 nr = source.get_mg()->get_nr(zone) ;
691 for (int i=0 ; i<nr ; i++)
692 resultat.set(zone, k, j, i) =
693 solution_part(zone, k, j, i)
694 +facteurs(conte)*solution_hom_un(zone, k, j, i)
695 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
696 conte+=2 ;
697 }
698
699 //rangement dans la ZEC :
700 if (indic[nz-1] == 1) {
701 nr = source.get_mg()->get_nr(nz-1) ;
702 for (int i=0 ; i<nr ; i++)
703 resultat.set(nz-1, k, j, i) =
704 solution_part(nz-1, k, j, i)
705 +facteurs(conte)*solution_hom_un(nz-1, k, j, i) ;
706 }
707
708 delete sec_membre ;
709 delete systeme ;
710 }
711
712 }
713
714 delete [] indic ;
715
716 // Les trucs les plus internes sont mis a zero ...
717 for (int l=0 ; l<num_front+1 ; l++)
718 resultat.t[l]->set_etat_zero() ;
719 return resultat;
720}
721}
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::Matrice::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition matrice.C:184
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
Lorene
Lorene prototypes.
Definition app_hor.h:67