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 
181 namespace Lorene {
182 Mtbl_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
Lorene prototypes.
Definition: app_hor.h:67
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
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
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:465
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724