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 char 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 
178 namespace Lorene {
179 Mtbl_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 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
Lorene::Mtbl_cf::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453