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 char 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 
107 namespace Lorene {
108 Mtbl_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