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