LORENE: poisson.C Source File

00001 /*
00002  *   Copyright (c) 1999-2001 Philippe Grandclement
00003  *   Copyright (c) 1999-2001 Eric Gourgoulhon
00004  *
00005  *   This file is part of LORENE.
00006  *
00007  *   LORENE is free software; you can redistribute it and/or modify
00008  *   it under the terms of the GNU General Public License as published by
00009  *   the Free Software Foundation; either version 2 of the License, or
00010  *   (at your option) any later version.
00011  *
00012  *   LORENE is distributed in the hope that it will be useful,
00013  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  *   GNU General Public License for more details.
00016  *
00017  *   You should have received a copy of the GNU General Public License
00018  *   along with LORENE; if not, write to the Free Software
00019  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  */
00022 
00023 
00024 char poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson.C,v 1.6 2007/12/13 15:48:46 jl_cornou Exp $" ;
00025 
00026 /*
00027  * $Id: poisson.C,v 1.6 2007/12/13 15:48:46 jl_cornou Exp $
00028  * $Log: poisson.C,v $
00029  * Revision 1.6  2007/12/13 15:48:46  jl_cornou
00030  * *** empty log message ***
00031  *
00032  * Revision 1.5  2007/12/11 15:28:23  jl_cornou
00033  * Jacobi(0,2) polynomials partially implemented
00034  *
00035  * Revision 1.4  2004/10/05 15:44:21  j_novak
00036  * Minor speed enhancements.
00037  *
00038  * Revision 1.3  2004/02/20 10:55:23  j_novak
00039  * The versions dzpuis 5 -> 3 has been improved and polished. Should be
00040  * operational now...
00041  *
00042  * Revision 1.2  2004/02/06 10:53:54  j_novak
00043  * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
00044  *
00045  * Revision 1.1.1.1  2001/11/20 15:19:28  e_gourgoulhon
00046  * LORENE
00047  *
00048  * Revision 2.24  2000/07/18  10:23:09  eric
00049  * Suppression d'une erreur sans consequence : dans le noyau, remplacement de
00050  *   beta = mapping.get_alpha()[0]    par
00051  *   beta = mapping.get_beta()[0]
00052  * Cette erreur etait sans consequence car beta n'intervient pas dans la
00053  * suite des calculs pour le noyau.
00054  *
00055  * Revision 2.23  2000/05/22  13:43:46  phil
00056  * ajout du cas dzpuis =3
00057  *
00058  * Revision 2.22  2000/02/09  14:40:49  eric
00059  * Ajout de l'argument dzpuis a sol_poisson.
00060  * (dzpuis n'est plus lu sur le Mtbl_cf).
00061  *
00062  * Revision 2.21  1999/12/02  14:29:05  eric
00063  * Suppression de la base en argument de la version Mtbl_cf.
00064  * La version Map se trouve desormais dans le fichier map_af_poisson.C
00065  * La version Cmp se trouve desormais dans le fichier cmp_pde.C
00066  *
00067  * Revision 2.20  1999/11/30  12:56:45  eric
00068  * Valeur::base est desormais du type Base_val et non plus Base_val*.
00069  *
00070  * Revision 2.19  1999/11/24  14:34:25  eric
00071  * Accession des membres alpha et beta de Map_af (desormais prives) par
00072  *   Map_af::get_alpha() et Map_af::get_beta().
00073  *
00074  * Revision 2.18  1999/10/27  15:47:19  eric
00075  * Suppression du membre Cmp::c.
00076  *
00077  * Revision 2.17  1999/10/14  14:28:07  eric
00078  * Methode const.
00079  *
00080  * Revision 2.16  1999/10/13  15:52:22  eric
00081  * Ajout de la base dans l'appel au constructeur de Mtbl_cf.
00082  *
00083  * Revision 2.15  1999/10/11  16:27:59  phil
00084  * suppression du sur echantillonnage
00085  *
00086  * Revision 2.14  1999/10/11  14:28:50  phil
00087  * & -> &&
00088  *
00089  * Revision 2.13  1999/09/06  16:25:42  phil
00090  * ajout de la version Cmp
00091  *
00092  * Revision 2.12  1999/07/02  15:05:01  phil
00093  * *** empty log message ***
00094  *
00095  * Revision 2.11  1999/06/23  12:35:18  phil
00096  * ajout de dzpuis = 2
00097  *
00098  * Revision 2.10  1999/04/27  13:12:28  phil
00099  * *** empty log message ***
00100  *
00101  * Revision 2.9  1999/04/14  13:56:03  phil
00102  * *** empty log message ***
00103  *
00104  * Revision 2.8  1999/04/14  10:22:28  phil
00105  * *** empty log message ***
00106  *
00107  * Revision 2.7  1999/04/14  10:20:09  phil
00108  * *** empty log message ***
00109  *
00110  * Revision 2.6  1999/04/13  17:18:50  phil
00111  * Ajout de la version Valeur
00112  *
00113  * Revision 2.5  1999/04/12  15:21:37  phil
00114  * correction de la construction du systeme
00115  *
00116  * Revision 2.4  1999/04/12  14:55:16  phil
00117  * correction matrices
00118  *
00119  * Revision 2.3  1999/04/12  14:28:46  phil
00120  * *** empty log message ***
00121  *
00122  * Revision 2.2  1999/04/07  14:56:33  phil
00123  * Changement prototypage
00124  *
00125  * Revision 2.1  1999/04/07  14:40:55  phil
00126  * mise au point des includes.
00127  *
00128  * Revision 2.0  1999/04/07  14:10:55  phil
00129  * *** empty log message ***
00130  *
00131  *
00132  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson.C,v 1.6 2007/12/13 15:48:46 jl_cornou Exp $
00133  *
00134  */
00135 
00136 // Header C : 
00137 #include <stdlib.h>
00138 #include <math.h>
00139 
00140 // Headers Lorene :
00141 #include "matrice.h"
00142 #include "map.h"
00143 #include "proto.h"
00144 #include "type_parite.h"
00145 
00146 
00147 
00148         //----------------------------------------------
00149        //       Version Mtbl_cf
00150       //----------------------------------------------
00151 
00152 /*
00153  * 
00154  * Solution de l'equation de poisson
00155  * 
00156  * Entree : mapping :   le mapping affine
00157  *      source : les coefficients de la source qui a ete multipliee par
00158  *          r^4 ou r^2 dans la ZEC.
00159  *          La base de decomposition doit etre Ylm
00160  *      dzpuis : exposant de r dans le factor multiplicatif dans la ZEC
00161  * Sortie : renvoie les coefficients de la solution dans la meme base de 
00162  *      decomposition (a savoir Ylm)
00163  *      
00164  */
00165 
00166 
00167 
00168 Mtbl_cf sol_poisson(const Map_af& mapping, const Mtbl_cf& source, int dzpuis,
00169             bool match)
00170 {
00171     
00172     // Verifications d'usage sur les zones
00173     int nz = source.get_mg()->get_nzone() ;
00174     assert (nz>1) ;
00175     assert ((source.get_mg()->get_type_r(0) == RARE) || (source.get_mg()->get_type_r(0) == FINJAC)) ;
00176     assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
00177     for (int l=1 ; l<nz-1 ; l++)
00178     assert(source.get_mg()->get_type_r(l) == FIN) ;
00179 
00180      assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3) || (dzpuis==5)) ;
00181      assert ((!match) || (dzpuis != 5)) ;
00182     
00183     // Bases spectrales
00184     const Base_val& base = source.base ;
00185     
00186     
00187     // donnees sur la zone
00188     int nr, nt, np ;
00189     int base_r ;
00190     double alpha, beta, echelle ;
00191     int l_quant, m_quant;
00192     
00193     //Rangement des valeurs intermediaires 
00194     Tbl *so = 0x0 ;
00195     Tbl *sol_hom = 0x0 ;
00196     Tbl *sol_part = 0x0 ;
00197     Matrice *operateur = 0x0 ;
00198     Matrice *nondege = 0x0 ;
00199     
00200     
00201     // Rangement des solutions, avant raccordement
00202     Mtbl_cf solution_part(source.get_mg(), base) ;
00203     Mtbl_cf solution_hom_un(source.get_mg(), base) ;
00204     Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
00205     Mtbl_cf resultat(source.get_mg(), base) ;
00206     
00207     solution_part.set_etat_qcq() ;
00208     solution_hom_un.set_etat_qcq() ;
00209     solution_hom_deux.set_etat_qcq() ;
00210     resultat.annule_hard() ;
00211     for (int l=0 ; l<nz ; l++) {
00212     solution_part.t[l]->set_etat_qcq() ;
00213     solution_hom_un.t[l]->set_etat_qcq() ;
00214     solution_hom_deux.t[l]->set_etat_qcq() ;
00215     }    
00216     // nbre maximum de point en theta et en phi :
00217     int np_max, nt_max ;
00218     
00219            //---------------
00220           //--  NOYAU -----
00221          //---------------
00222          
00223     nr = source.get_mg()->get_nr(0) ;
00224     nt = source.get_mg()->get_nt(0) ;
00225     np = source.get_mg()->get_np(0) ;
00226     
00227     nt_max = nt ;
00228     np_max = np ;
00229     
00230     alpha = mapping.get_alpha()[0] ;
00231     beta = mapping.get_beta()[0] ;
00232     
00233     for (int k=0 ; k<np+1 ; k++)
00234     for (int j=0 ; j<nt ; j++) 
00235         if (nullite_plm(j, nt, k, np, base) == 1) 
00236         {
00237         // calcul des nombres quantiques :
00238         donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
00239         
00240         // Construction operateur
00241         operateur = new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
00242         (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
00243         
00244         //Operateur inversible
00245         nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
00246         
00247         if (match) {
00248           // Calcul de la SH
00249           sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
00250         }
00251      
00252         //Calcul de la SP
00253         so = new Tbl(nr) ;
00254         so->set_etat_qcq() ;
00255         for (int i=0 ; i<nr ; i++)
00256         so->set(i) = source(0, k, j, i) ;
00257         
00258         sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta, 
00259                     *so, 0, base_r)) ;
00260         
00261         // Rangement dans les tableaux globaux ;
00262         
00263         for (int i=0 ; i<nr ; i++) {
00264         solution_part.set(0, k, j, i) = (*sol_part)(i) ;
00265         if (match) {
00266             if (base_r == R_JACO02) {
00267                 solution_hom_un.set(0, k, j, i) = (*sol_hom)(0, i) ;
00268                     solution_hom_deux.set(0, k, j, i) = (*sol_hom)(1, i) ;
00269             }
00270             else {
00271                 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
00272                 solution_hom_deux.set(0, k, j, i) = 0. ;
00273             }
00274         }
00275         }       
00276         
00277         
00278         
00279         delete operateur ;
00280         delete nondege ;
00281         delete so ;
00282         if (match) delete sol_hom ;
00283         delete sol_part ;
00284     }
00285 
00286 
00287            //---------------
00288           //--  ZEC   -----
00289          //---------------
00290          
00291     nr = source.get_mg()->get_nr(nz-1) ;
00292     nt = source.get_mg()->get_nt(nz-1) ;
00293     np = source.get_mg()->get_np(nz-1) ;
00294     
00295     if (np > np_max) np_max = np ;
00296     if (nt > nt_max) nt_max = nt ;
00297    
00298     alpha = mapping.get_alpha()[nz-1] ;
00299     beta = mapping.get_beta()[nz-1] ;
00300     
00301     for (int k=0 ; k<np+1 ; k++)
00302     for (int j=0 ; j<nt ; j++) 
00303         if (nullite_plm(j, nt, k, np, base) == 1)
00304         {
00305         
00306         // calcul des nombres quantiques :
00307         donne_lm(nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
00308         
00309         //Construction operateur
00310         operateur = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis,
00311                                  base_r)) ;
00312         (*operateur) = combinaison(*operateur, l_quant, 0., dzpuis, base_r) ;
00313         // Operateur inversible
00314         nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0.,
00315                                  dzpuis, base_r)) ;
00316         // Calcul de la SH
00317         if (match) 
00318           sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
00319        
00320         // Calcul de la SP
00321         so = new Tbl(nr) ;
00322         so->set_etat_qcq() ;
00323         for (int i=0 ; i<nr ; i++)
00324         so->set(i) = source(nz-1, k, j, i) ;
00325         sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta, 
00326                 *so, dzpuis, base_r)) ;
00327 
00328         // Rangement
00329         for (int i=0 ; i<nr ; i++) {
00330         solution_part.set(nz-1, k, j, i) = (*sol_part)(i) ;
00331         if (match) {
00332           solution_hom_un.set(nz-1, k, j, i) = (*sol_hom)(i) ;
00333           solution_hom_deux.set(nz-1, k, j, i) = 0. ;
00334         }
00335         }
00336       
00337         delete operateur ;
00338         delete nondege ;
00339         delete so ;
00340         if (match) delete sol_hom ;
00341         delete sol_part ;
00342     }
00343     
00344            //---------------
00345           //- COQUILLES ---
00346          //---------------
00347 
00348     for (int zone=1 ; zone<nz-1 ; zone++) {
00349     nr = source.get_mg()->get_nr(zone) ;
00350     nt = source.get_mg()->get_nt(zone) ;
00351     np = source.get_mg()->get_np(zone) ;
00352     
00353     if (np > np_max) np_max = np ;
00354     if (nt > nt_max) nt_max = nt ;
00355     
00356     alpha = mapping.get_alpha()[zone] ;
00357     beta = mapping.get_beta()[zone] ;
00358 
00359     for (int k=0 ; k<np+1 ; k++)
00360         for (int j=0 ; j<nt ; j++) 
00361         if (nullite_plm(j, nt, k, np, base) == 1)
00362         {
00363         // calcul des nombres quantiques :
00364         donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
00365         
00366         // Construction de l'operateur
00367         operateur = new Matrice(laplacien_mat
00368                     (nr, l_quant, beta/alpha, 0, base_r)) ;
00369         
00370         (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0, 
00371                                      base_r) ;
00372         
00373          // Operateur inversible
00374         nondege = new Matrice(prepa_nondege(*operateur, l_quant, 
00375                             beta/alpha, 0, base_r)) ;       
00376         if (match) {
00377           // Calcul DES DEUX SH
00378           sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
00379         }
00380         
00381         // Calcul de la SP
00382         so = new Tbl(nr) ;
00383         so->set_etat_qcq() ;
00384         for (int i=0 ; i<nr ; i++)
00385             so->set(i) = source(zone, k, j, i) ;
00386         
00387         sol_part = new Tbl (solp(*operateur, *nondege, alpha,
00388                      beta, *so, 0, base_r)) ;
00389         
00390         
00391         // Rangement
00392         for (int i=0 ; i<nr ; i++) {
00393             solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
00394             if (match) {
00395               solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
00396               solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
00397             }
00398         }
00399             
00400         
00401         delete operateur ;
00402         delete nondege ;
00403         delete so ;
00404         if (match) delete sol_hom ;
00405         delete sol_part ;
00406         }
00407     }
00408 
00409 
00410     if (match) {
00411 
00412          //-------------------------------------------
00413         // On est parti pour le raccord des solutions
00414        //-------------------------------------------
00415     // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
00416     // intervient dans le developpement de cette zone.
00417     int * indic = new int [nz] ;
00418     int taille ;
00419     Tbl *sec_membre ; // termes constants du systeme
00420     Matrice *systeme ; //le systeme a resoudre
00421 
00422     // des compteurs pour le remplissage du systeme
00423     int ligne ;
00424     int colonne ;
00425 
00426     // compteurs pour les diagonales du systeme :
00427     int sup ;
00428     int inf ;
00429     int sup_new, inf_new ;
00430 
00431     // on boucle sur le plus grand nombre possible de Plm intervenant...
00432     for (int k=0 ; k<np_max+1 ; k++)
00433     for (int j=0 ; j<nt_max ; j++)
00434         if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
00435         
00436         ligne = 0 ;
00437         colonne = 0 ;
00438         sup = 0 ;
00439         inf = 0 ;
00440         
00441         for (int zone=0 ; zone<nz ; zone++)
00442             indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone), 
00443                      k,  source.get_mg()->get_np(zone), base);
00444         
00445         // taille du systeme a resoudre pour ce Plm 
00446         taille = indic[nz-1]+indic[0] ;
00447         for (int zone=1 ; zone<nz-1 ; zone++)
00448             taille+=2*indic[zone] ;
00449         
00450         // on verifie que la taille est non-nulle.
00451         // cas pouvant a priori se produire...
00452         
00453         if (taille != 0) {
00454             
00455             sec_membre = new Tbl(taille) ;
00456             sec_membre->set_etat_qcq() ;
00457             for (int l=0 ; l<taille ; l++)
00458             sec_membre->set(l) = 0 ;
00459             
00460             systeme = new Matrice(taille, taille) ;
00461             systeme->set_etat_qcq() ;
00462             for (int l=0 ; l<taille ; l++)
00463             for (int c=0 ; c<taille ; c++)
00464                 systeme->set(l, c) = 0 ;
00465             
00466             //Calcul des nombres quantiques
00467             //base_r est donne dans le noyau, sa valeur dans les autres
00468             //zones etant inutile.
00469             
00470             donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
00471             
00472             
00473              //--------------------------
00474             //      NOYAU 
00475            //---------------------------
00476             
00477             if (indic[0] == 1) {
00478             nr = source.get_mg()->get_nr(0) ;
00479             
00480             alpha = mapping.get_alpha()[0] ;
00481             // valeur de x^l en 1 :
00482             systeme->set(ligne, colonne) = 1. ;
00483             // coefficient du Plm dans la solp
00484             for (int i=0 ; i<nr ; i++)
00485                 sec_membre->set(ligne) -= solution_part(0, k, j, i) ;
00486                 
00487             ligne++ ;
00488             // on prend les derivees que si Plm existe 
00489             //dans la zone suivante
00490 
00491 
00492 
00493 // Grosses modifications en perspective
00494             
00495             if (indic[1] == 1) {
00496                 // derivee de x^l en 1 :
00497                 systeme->set(ligne, colonne) = 1./alpha*l_quant ; 
00498             
00499                 // coefficient de la derivee du Plm dans la solp
00500                 if (base_r == R_CHEBP)
00501                 // cas en R_CHEBP 
00502                 for (int i=0 ; i<nr ; i++)
00503                     sec_membre->set(ligne) -= 
00504                         4*i*i/alpha
00505                         *solution_part(0, k, j, i) ;
00506                 else
00507                 // cas en R_CHEBI
00508                 for (int i=0 ; i<nr ; i++)
00509                     sec_membre->set(ligne) -=
00510                     (2*i+1)*(2*i+1)/alpha
00511                         *solution_part(0, k, j, i) ;
00512                             
00513                 // on a au moins un diag inferieure dans ce cas ...
00514                 inf = 1 ;
00515             }
00516             colonne++ ; 
00517             }
00518             
00519             //-----------------------------
00520            //        COQUILLES
00521           //------------------------------
00522           
00523           for (int zone=1 ; zone<nz-1 ; zone++)
00524             if (indic[zone] == 1) {
00525             
00526             nr = source.get_mg()->get_nr(zone) ;
00527             alpha = mapping.get_alpha()[zone] ;
00528             echelle = mapping.get_beta()[zone]/alpha ;
00529             
00530             //Frontiere avec la zone precedente :
00531             if (indic[zone-1] == 1) ligne -- ;
00532             
00533             //valeur de (x+echelle)^l en -1 :
00534             systeme->set(ligne, colonne) = 
00535                 -pow(echelle-1, double(l_quant)) ;
00536             
00537             // valeur de 1/(x+echelle) ^(l+1) en -1 :
00538             systeme->set(ligne, colonne+1) = 
00539                 -1/pow(echelle-1, double(l_quant+1)) ;
00540             
00541             // la solution particuliere :
00542             for (int i=0 ; i<nr ; i++)
00543                 if (i%2 == 0)
00544                 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
00545                 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
00546             
00547             // les diagonales :
00548             sup_new = colonne+1-ligne ;
00549             if (sup_new > sup) sup = sup_new ;
00550             
00551             ligne++ ;
00552             
00553             // on prend les derivees si Plm existe dans la zone 
00554             // precedente :
00555             
00556             if (indic[zone-1] == 1) {
00557             // derivee de (x+echelle)^l en -1 :
00558                 systeme->set(ligne, colonne) = 
00559                 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ;
00560             // derivee de 1/(x+echelle)^(l+1) en -1 :
00561                 systeme->set(ligne, colonne+1) = 
00562                 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ;
00563             
00564             
00565             
00566             // la solution particuliere :
00567                 for (int i=0 ; i<nr ; i++)
00568                 if (i%2 == 0)
00569                     sec_membre->set(ligne) -= 
00570                     i*i/alpha*solution_part(zone, k, j, i) ;
00571                 else
00572                     sec_membre->set(ligne) +=
00573                     i*i/alpha*solution_part(zone, k, j, i) ;
00574             
00575             // les diagonales :
00576             sup_new = colonne+1-ligne ;
00577             if (sup_new > sup) sup = sup_new ;
00578             
00579             ligne++ ;
00580             }
00581             
00582             // Frontiere avec la zone suivante :
00583             //valeur de (x+echelle)^l en 1 :
00584             systeme->set(ligne, colonne) = 
00585                 pow(echelle+1, double(l_quant)) ;
00586             
00587             // valeur de 1/(x+echelle)^(l+1) en 1 :
00588             systeme->set(ligne, colonne+1) =
00589                 1/pow(echelle+1, double(l_quant+1)) ;
00590                 
00591             // la solution particuliere :
00592             for (int i=0 ; i<nr ; i++)
00593                 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
00594             
00595             // les diagonales inf :
00596             inf_new = ligne-colonne ;
00597             if (inf_new > inf) inf = inf_new ;  
00598             
00599             ligne ++ ;
00600                 
00601             // Utilisation des derivees ssi Plm existe dans la
00602             //zone suivante :
00603             if (indic[zone+1] == 1) {
00604                 
00605                 //derivee de (x+echelle)^l en 1 :
00606                 systeme->set(ligne, colonne) = 
00607                 l_quant*pow(echelle+1, double(l_quant-1))/alpha ;
00608                 
00609                 //derivee de 1/(echelle+x)^(l+1) en 1 :
00610                 systeme->set(ligne, colonne+1) =
00611                 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ;
00612                 
00613                 // La solution particuliere :
00614                 for (int i=0 ; i<nr ; i++)
00615                 sec_membre->set(ligne) -= 
00616                     i*i/alpha*solution_part(zone, k, j, i) ;
00617                     
00618             // les diagonales inf :
00619             inf_new = ligne-colonne ;
00620             if (inf_new > inf) inf = inf_new ;  
00621                 
00622                  }
00623              colonne += 2 ;
00624              }
00625             
00626             
00627             //--------------------------------
00628            //             ZEC
00629           //---------------------------------
00630           
00631           if (indic[nz-1] == 1) {
00632               
00633               nr = source.get_mg()->get_nr(nz-1) ;
00634               
00635               
00636               alpha = mapping.get_alpha()[nz-1] ;
00637               
00638               if (indic[nz-2] == 1) ligne -- ;
00639               
00640               //valeur de (x-1)^(l+1) en -1 :
00641               systeme->set(ligne, colonne) = -pow(-2, double(l_quant+1)) ;
00642               //solution particuliere :
00643               for (int i=0 ; i<nr ; i++)
00644             if (i%2 == 0)
00645                 sec_membre->set(ligne) += solution_part(nz-1, k, j, i) ;
00646             else sec_membre->set(ligne) -= solution_part(nz-1, k, j, i) ;
00647             
00648             //on prend les derivees ssi Plm existe dans la zone precedente :
00649             if (indic[nz-2] == 1) {
00650             
00651             //derivee de (x-1)^(l+1) en -1 :
00652             systeme->set(ligne+1, colonne) = 
00653                 alpha*(l_quant+1)*pow(-2., double(l_quant+2)) ;
00654             
00655             // Solution particuliere :
00656             for (int i=0 ; i<nr ; i++)
00657                 if (i%2 == 0)
00658                 sec_membre->set(ligne+1) -=
00659                     -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
00660                 else
00661                 sec_membre->set(ligne+1) +=
00662                     -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
00663             
00664             //les diags :
00665             if (sup == 0) sup = 1 ;
00666             }
00667           }
00668             
00669             //-------------------------
00670            //  resolution du systeme 
00671           //--------------------------
00672             
00673             systeme->set_band(sup, inf) ;
00674             systeme->set_lu() ;
00675 
00676             Tbl facteurs(systeme->inverse(*sec_membre)) ;
00677             int conte = 0 ;
00678             
00679          
00680         // rangement dans le noyau :
00681             
00682             if (indic[0] == 1) {
00683             nr = source.get_mg()->get_nr(0) ;
00684             for (int i=0 ; i<nr ; i++)
00685                 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
00686                 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
00687             conte++ ;
00688             }
00689             
00690             // rangement dans les coquilles :
00691             for (int zone=1 ; zone<nz-1 ; zone++)
00692             if (indic[zone] == 1) {
00693                 nr = source.get_mg()->get_nr(zone) ;
00694                 for (int i=0 ; i<nr ; i++)
00695                 resultat.set(zone, k, j, i) = 
00696                 solution_part(zone, k, j, i)
00697                 +facteurs(conte)*solution_hom_un(zone, k, j, i) 
00698                 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
00699                 conte+=2 ;
00700             }
00701             
00702             //rangement dans la ZEC :
00703             if (indic[nz-1] == 1) {
00704              nr = source.get_mg()->get_nr(nz-1) ;
00705              for (int i=0 ; i<nr ; i++)
00706              resultat.set(nz-1, k, j, i) = 
00707                solution_part(nz-1, k, j, i)
00708                 +facteurs(conte)*solution_hom_un(nz-1, k, j, i) ;
00709             }
00710             delete sec_membre ;
00711             delete systeme ;
00712         }
00713     
00714         }
00715     
00716     delete [] indic ;
00717 
00718     } // End of the case where the matching has to be done
00719     else { //Only a particular solution is given in each domain
00720 
00721       for (int zone = 0; zone < nz; zone++) 
00722     for (int k=0 ; k<source.get_mg()->get_np(zone)+1 ; k++)
00723       for (int j=0 ; j<source.get_mg()->get_nt(zone) ; j++)
00724         if (nullite_plm(j,source.get_mg()->get_nt(zone) , 
00725                 k, source.get_mg()->get_np(zone), base) == 1) 
00726           for (int i=0; i<source.get_mg()->get_nr(zone); i++) 
00727         resultat.set(zone, k, j, i) = solution_part(zone, k, j, i) ;
00728 
00729     }
00730 
00731     return resultat;
00732 }