op_d2sdx2.C

00001 /*
00002  *   Copyright (c) 1999-2000 Jean-Alain Marck
00003  *   Copyright (c) 1999-2001 Eric Gourgoulhon
00004  *   Copyright (c) 1999-2001 Philippe Grandclement
00005  *
00006  *   This file is part of LORENE.
00007  *
00008  *   LORENE is free software; you can redistribute it and/or modify
00009  *   it under the terms of the GNU General Public License as published by
00010  *   the Free Software Foundation; either version 2 of the License, or
00011  *   (at your option) any later version.
00012  *
00013  *   LORENE is distributed in the hope that it will be useful,
00014  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00015  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016  *   GNU General Public License for more details.
00017  *
00018  *   You should have received a copy of the GNU General Public License
00019  *   along with LORENE; if not, write to the Free Software
00020  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00021  *
00022  */
00023 
00024 
00025 char op_d2sdx2_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdx2.C,v 1.4 2008/08/27 08:50:10 jl_cornou Exp $" ;
00026 
00027 /* 
00028  * Ensemble des routines de base de derivation seconde par rapport a r
00029  * (Utilisation interne)
00030  * 
00031  *  void _d2sdx2_XXXX(Tbl * t, int & b)
00032  *  t   pointeur sur le Tbl d'entree-sortie
00033  *  b   base spectrale
00034  * 
00035  */
00036 
00037 /*
00038  * $Id: op_d2sdx2.C,v 1.4 2008/08/27 08:50:10 jl_cornou Exp $
00039  * $Log: op_d2sdx2.C,v $
00040  * Revision 1.4  2008/08/27 08:50:10  jl_cornou
00041  * Added Jacobi(0,2) polynomials
00042  *
00043  * Revision 1.3  2007/12/11 15:28:18  jl_cornou
00044  * Jacobi(0,2) polynomials partially implemented
00045  *
00046  * Revision 1.2  2004/11/23 15:16:01  m_forot
00047  *
00048  * Added the bases for the cases without any equatorial symmetry
00049  *  (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
00050  *
00051  * Revision 1.1.1.1  2001/11/20 15:19:29  e_gourgoulhon
00052  * LORENE
00053  *
00054  * Revision 2.7  2000/02/24  16:40:50  eric
00055  * Initialisation a zero du tableau xo avant le calcul.
00056  *
00057  * Revision 2.6  1999/11/15  16:37:53  eric
00058  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
00059  *
00060  * Revision 2.5  1999/10/11  15:33:31  phil
00061  * *** empty log message ***
00062  *
00063  * Revision 2.4  1999/10/11  14:25:47  phil
00064  * realloc -> delete + new
00065  *
00066  * Revision 2.3  1999/09/13  11:30:23  phil
00067  * gestion du cas k==1
00068  *
00069  * Revision 2.2  1999/03/01  15:06:31  eric
00070  * *** empty log message ***
00071  *
00072  * Revision 2.1  1999/02/23  10:39:13  hyc
00073  * *** empty log message ***
00074  *
00075  *
00076  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdx2.C,v 1.4 2008/08/27 08:50:10 jl_cornou Exp $
00077  *
00078  */
00079 
00080 // Fichier includes
00081 #include "tbl.h"
00082 
00083 void d2sdx2_1d(int, double**, int);
00084 void _d2sdx2_1d_r_jaco02(int, double*, double*) ;
00085 
00086 // Prototypage
00087 void c_est_pas_fait(char * ) ;
00088 
00089 // Routine pour les cas non prevus
00090 //--------------------------------
00091 void _d2sdx2_pas_prevu(Tbl* , int & b) {
00092     cout << "d2sdx2 pas prevu..." << endl ;
00093     cout << " base: " << b << endl ;
00094     abort () ;
00095 }
00096 
00097 // cas R_CHEBU - dzpuis = 0
00098 //-------------------------
00099 void _d2sdx2_r_chebu_0(Tbl *tb, int & )
00100 {
00101 
00102     // Peut-etre rien a faire ?
00103     if (tb->get_etat() == ETATZERO) {
00104     return ;
00105     }
00106     
00107     // Protection
00108     assert(tb->get_etat() == ETATQCQ) ;
00109 
00110     // Pour le confort
00111     int nr = (tb->dim).dim[0] ;     // Nombre
00112     int nt = (tb->dim).dim[1] ;     //   de points
00113     int np = (tb->dim).dim[2] ;     //      physiques REELS
00114     np = np - 2 ;           // Nombre de points physiques
00115     
00116     // Variables statiques
00117     static double* cx1 = 0x0 ;
00118     static double* cx2 = 0x0 ;
00119     static double* cx3 = 0x0 ;
00120     static int nr_pre = 0 ;
00121 
00122     // Test sur np pour initialisation eventuelle
00123     if (nr > nr_pre) {
00124     nr_pre = nr ;
00125     if (cx1 != 0x0) delete [] cx1 ;
00126     if (cx2 != 0x0) delete [] cx2 ;
00127     if (cx3 != 0x0) delete [] cx3 ;
00128     cx1 = new double [nr] ;
00129     cx2 = new double [nr] ;
00130     cx3 = new double [nr] ;
00131     for (int i=0 ; i<nr ; i++) {
00132         cx1[i] =  (i+2)*(i+2)*(i+2) ;
00133         cx2[i] =  (i+2) ;
00134         cx3[i] =  i*i ;
00135     }
00136     }
00137     
00138     // pt. sur le tableau de double resultat
00139     double* xo = new double[(tb->dim).taille] ;
00140     
00141     // Initialisation a zero :
00142     for (int i=0; i<(tb->dim).taille; i++) {
00143     xo[i] = 0 ; 
00144     }
00145     
00146     // On y va...
00147     double* xi = tb->t ;
00148     double* xci = xi ;  // Pointeurs
00149     double* xco = xo ;  //  courants
00150     
00151     for (int k=0 ; k<np+1 ; k++) 
00152     if (k == 1)  {
00153         xci += nr*nt ;
00154         xco += nr*nt ;
00155         }
00156     else {
00157     for (int j=0 ; j<nt ; j++) {
00158 
00159         double som1, som2 ;
00160        
00161         xco[nr-1] = 0 ;
00162         som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
00163         som2 = (nr-1) * xci[nr-1] ;
00164         xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
00165         for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
00166         som1 += cx1[i] * xci[i+2] ;
00167         som2 += cx2[i] * xci[i+2] ;
00168         xco[i] = som1 - cx3[i] * som2 ;
00169         }   // Fin de la premiere boucle sur r
00170         xco[nr-2] = 0 ;
00171         som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
00172         som2 = (nr-2) * xci[nr-2] ;
00173         xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
00174         for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
00175         som1 += cx1[i] * xci[i+2] ;
00176         som2 += cx2[i] * xci[i+2] ;
00177         xco[i] = som1 - cx3[i] * som2 ;
00178         }   // Fin de la deuxieme boucle sur r
00179         xco[0] *= .5 ;
00180         
00181         xci += nr ;
00182         xco += nr ;
00183     }   // Fin de la boucle sur theta
00184     }   // Fin de la boucle sur phi
00185     
00186     // On remet les choses la ou il faut
00187     delete [] tb->t ;
00188     tb->t = xo ;
00189     
00190     // base de developpement
00191     // inchangee
00192 }
00193 
00194 // cas R_CHEB
00195 //-----------
00196 void _d2sdx2_r_cheb(Tbl *tb, int & )
00197 {
00198 
00199     // Peut-etre rien a faire ?
00200     if (tb->get_etat() == ETATZERO) {
00201     return ;
00202     }
00203     
00204     // Protection
00205     assert(tb->get_etat() == ETATQCQ) ;
00206     
00207     // Pour le confort
00208     int nr = (tb->dim).dim[0] ;     // Nombre
00209     int nt = (tb->dim).dim[1] ;     //   de points
00210     int np = (tb->dim).dim[2] ;     //      physiques REELS
00211     np = np - 2 ;           // Nombre de points physiques
00212     
00213     // Variables statiques
00214     static double* cx1 = 0x0 ;
00215     static double* cx2 = 0x0 ;
00216     static double* cx3 = 0x0 ;
00217     static int nr_pre = 0 ;
00218 
00219     // Test sur np pour initialisation eventuelle
00220     if (nr > nr_pre) {
00221     nr_pre = nr ;
00222     if (cx1 != 0x0) delete [] cx1 ;
00223     if (cx2 != 0x0) delete [] cx2 ;
00224     if (cx3 != 0x0) delete [] cx3 ;
00225     cx1 = new double [nr] ;
00226     cx2 = new double [nr] ;
00227     cx3 = new double [nr] ;  
00228     
00229     for (int i=0 ; i<nr ; i++) {
00230         cx1[i] =  (i+2)*(i+2)*(i+2) ;
00231         cx2[i] =  (i+2) ;
00232         cx3[i] =  i*i ;
00233     }
00234     }
00235     
00236     // pt. sur le tableau de double resultat
00237     double* xo = new double[(tb->dim).taille] ;
00238     
00239     // Initialisation a zero :
00240     for (int i=0; i<(tb->dim).taille; i++) {
00241     xo[i] = 0 ; 
00242     }
00243     
00244     // On y va...
00245     double* xi = tb->t ;
00246     double* xci = xi ;  // Pointeurs
00247     double* xco = xo ;  //  courants
00248     
00249     for (int k=0 ; k<np+1 ; k++) 
00250     if (k == 1)  {
00251         xci += nr*nt ;
00252         xco += nr*nt ;
00253         }
00254     else {
00255     for (int j=0 ; j<nt ; j++) {
00256 
00257         double som1, som2 ;
00258         
00259         xco[nr-1] = 0 ;
00260         som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
00261         som2 = (nr-1) * xci[nr-1] ;
00262         xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
00263         for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
00264         som1 += cx1[i] * xci[i+2] ;
00265         som2 += cx2[i] * xci[i+2] ;
00266         xco[i] = som1 - cx3[i] * som2 ;
00267         }   // Fin de la premiere boucle sur r
00268         xco[nr-2] = 0 ;
00269         som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
00270         som2 = (nr-2) * xci[nr-2] ;
00271         xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
00272         for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
00273         som1 += cx1[i] * xci[i+2] ;
00274         som2 += cx2[i] * xci[i+2] ;
00275         xco[i] = som1 - cx3[i] * som2 ;
00276         }   // Fin de la deuxieme boucle sur r
00277         xco[0] *= .5 ;
00278         
00279         xci += nr ;
00280         xco += nr ;
00281     }   // Fin de la boucle sur theta
00282     }   // Fin de la boucle sur phi
00283     
00284     // On remet les choses la ou il faut
00285     delete [] tb->t ;
00286     tb->t = xo ;
00287     
00288     // base de developpement
00289     // inchangee
00290 }
00291 
00292 // cas R_CHEBP
00293 //------------
00294 void _d2sdx2_r_chebp(Tbl *tb, int & )
00295 {
00296 
00297     // Peut-etre rien a faire ?
00298     if (tb->get_etat() == ETATZERO) {
00299     return ;
00300     }
00301     
00302     // Protection
00303     assert(tb->get_etat() == ETATQCQ) ;
00304     
00305     // Pour le confort
00306     int nr = (tb->dim).dim[0] ;     // Nombre
00307     int nt = (tb->dim).dim[1] ;     //   de points
00308     int np = (tb->dim).dim[2] ;     //      physiques REELS
00309     np = np - 2 ;           // Nombre de points physiques
00310     
00311     // Variables statiques
00312     static double* cx1 = 0x0 ;
00313     static double* cx2 = 0x0 ;
00314     static double* cx3 = 0x0 ;
00315     static int nr_pre = 0 ;
00316 
00317     // Test sur np pour initialisation eventuelle
00318     if (nr > nr_pre) {
00319     nr_pre = nr ;
00320     if (cx1 != 0x0) delete [] cx1 ;
00321     if (cx2 != 0x0) delete [] cx2 ;
00322     if (cx3 != 0x0) delete [] cx3 ;
00323     cx1 = new double [nr] ;
00324     cx2 = new double [nr] ;
00325     cx3 = new double [nr] ;  
00326     for (int i=0 ; i<nr ; i++) {
00327         cx1[i] =  8*(i+1)*(i+1)*(i+1) ;
00328         cx2[i] =  2*(i+1) ;
00329         cx3[i] =  4*i*i ;
00330     }
00331     }
00332     // pt. sur le tableau de double resultat
00333     double* xo = new double[(tb->dim).taille] ;
00334     
00335     // Initialisation a zero :
00336     for (int i=0; i<(tb->dim).taille; i++) {
00337     xo[i] = 0 ; 
00338     }
00339     
00340     // On y va...
00341     double* xi = tb->t ;
00342     double* xci = xi ;  // Pointeurs
00343     double* xco = xo ;  //  courants
00344     
00345     for (int k=0 ; k<np+1 ; k++) 
00346     if (k == 1)  {
00347         xci += nr*nt ;
00348         xco += nr*nt ;
00349         }
00350     else {
00351     for (int j=0 ; j<nt ; j++) {
00352 
00353         double som1, som2 ;
00354         
00355         xco[nr-1] = 0 ;
00356         som1 = 8*(nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
00357         som2 = 2*(nr-1) * xci[nr-1] ;
00358         xco[nr-2] = som1 - 4*(nr-2)*(nr-2)*som2 ;
00359         for (int i = nr-3 ; i >= 0 ; i-- ) {
00360         som1 += cx1[i] * xci[i+1] ;
00361         som2 += cx2[i] * xci[i+1] ;
00362         xco[i] = som1 - cx3[i] * som2 ;
00363         }   // Fin de la boucle sur r
00364         xco[0] *= .5 ;
00365         
00366         xci += nr ;
00367         xco += nr ;
00368     }   // Fin de la boucle sur theta
00369     }   // Fin de la boucle sur phi
00370     
00371     // On remet les choses la ou il faut
00372     delete [] tb->t ;
00373     tb->t = xo ;
00374     
00375     // base de developpement
00376     // inchangee
00377 }
00378 
00379 // cas R_CHEBI
00380 //------------
00381 void _d2sdx2_r_chebi(Tbl *tb, int & )
00382 {
00383 
00384     // Peut-etre rien a faire ?
00385     if (tb->get_etat() == ETATZERO) {
00386     return ;
00387     }
00388     
00389     // Protection
00390     assert(tb->get_etat() == ETATQCQ) ;
00391     
00392     // Pour le confort
00393     int nr = (tb->dim).dim[0] ;     // Nombre
00394     int nt = (tb->dim).dim[1] ;     //   de points
00395     int np = (tb->dim).dim[2] ;     //      physiques REELS
00396     np = np - 2 ;           // Nombre de points physiques
00397     
00398     // Variables statiques
00399     static double* cx1 = 0x0 ;
00400     static double* cx2 = 0x0 ;
00401     static double* cx3 = 0x0 ;
00402     static int nr_pre = 0 ;
00403 
00404     // Test sur np pour initialisation eventuelle
00405     if (nr > nr_pre) {
00406     nr_pre = nr ;
00407     if (cx1 != 0x0) delete [] cx1 ;
00408     if (cx2 != 0x0) delete [] cx2 ;
00409     if (cx3 != 0x0) delete [] cx3 ;
00410     cx1 = new double [nr] ;
00411     cx2 = new double [nr] ;
00412     cx3 = new double [nr] ;   
00413    
00414     for (int i=0 ; i<nr ; i++) {
00415         cx1[i] =  (2*i+3)*(2*i+3)*(2*i+3) ;
00416         cx2[i] =  (2*i+3) ;
00417         cx3[i] =  (2*i+1)*(2*i+1) ;
00418     }
00419     }
00420     
00421     // pt. sur le tableau de double resultat
00422     double* xo = new double[(tb->dim).taille] ;
00423     
00424     // Initialisation a zero :
00425     for (int i=0; i<(tb->dim).taille; i++) {
00426     xo[i] = 0 ; 
00427     }
00428     
00429     // On y va...
00430     double* xi = tb->t ;
00431     double* xci = xi ;  // Pointeurs
00432     double* xco = xo ;  //  courants
00433     
00434     for (int k=0 ; k<np+1 ; k++)  
00435     if (k == 1)  {
00436         xci += nr*nt ;
00437         xco += nr*nt ;
00438         }
00439     else {
00440     for (int j=0 ; j<nt ; j++) {
00441 
00442         double som1, som2 ;
00443         
00444         xco[nr-1] = 0 ;
00445         som1 = (2*nr-1)*(2*nr-1)*(2*nr-1) * xci[nr-1] ;
00446         som2 = (2*nr-1) * xci[nr-1] ;
00447         xco[nr-2] = som1 - (2*nr-3)*(2*nr-3)*som2 ;
00448         for (int i = nr-3 ; i >= 0 ; i-- ) {
00449         som1 += cx1[i] * xci[i+1] ;
00450         som2 += cx2[i] * xci[i+1] ;
00451         xco[i] = som1 - cx3[i] * som2 ;
00452         }   // Fin de la boucle sur r
00453         
00454         xci += nr ;
00455         xco += nr ;
00456     }   // Fin de la boucle sur theta
00457     }   // Fin de la boucle sur phi
00458     
00459     // On remet les choses la ou il faut
00460     delete [] tb->t ;
00461     tb->t = xo ;
00462     
00463     // base de developpement
00464     // inchangee
00465 }
00466 
00467 // cas R_CHEBPIM_P
00468 //----------------
00469 void _d2sdx2_r_chebpim_p(Tbl *tb, int & )
00470 {
00471 
00472     // Peut-etre rien a faire ?
00473     if (tb->get_etat() == ETATZERO) {
00474     return ;
00475     }
00476     
00477     // Protection
00478     assert(tb->get_etat() == ETATQCQ) ;
00479     
00480     // Pour le confort
00481     int nr = (tb->dim).dim[0] ;     // Nombre
00482     int nt = (tb->dim).dim[1] ;     //   de points
00483     int np = (tb->dim).dim[2] ;     //      physiques REELS
00484     np = np - 2 ;           // Nombre de points physiques
00485     
00486     // Variables statiques
00487     static double* cx1p = 0x0 ;
00488     static double* cx2p = 0x0 ;
00489     static double* cx3p = 0x0 ;
00490     static double* cx1i = 0x0 ;
00491     static double* cx2i = 0x0 ;
00492     static double* cx3i = 0x0 ;
00493     static int nr_pre = 0 ;
00494 
00495     // Test sur np pour initialisation eventuelle
00496     if (nr > nr_pre) {
00497     nr_pre = nr ;
00498     if (cx1p != 0x0) delete [] cx1p ;
00499     if (cx2p != 0x0) delete [] cx2p ;
00500     if (cx3p != 0x0) delete [] cx3p ;
00501     if (cx1i != 0x0) delete [] cx1i ;
00502     if (cx2i != 0x0) delete [] cx2i ;
00503     if (cx3i != 0x0) delete [] cx3i ;
00504     cx1p = new double[nr] ;
00505     cx2p = new double[nr] ;
00506     cx3p = new double[nr] ;
00507     cx1i = new double[nr] ;
00508     cx2i = new double[nr] ;
00509     cx3i = new double[nr] ;
00510     for (int i=0 ; i<nr ; i++) {
00511         cx1p[i] =  8*(i+1)*(i+1)*(i+1) ;
00512         cx2p[i] =  2*(i+1) ;
00513         cx3p[i] =  4*i*i ;
00514 
00515         cx1i[i] =  (2*i+3)*(2*i+3)*(2*i+3) ;
00516         cx2i[i] =  (2*i+3) ;
00517         cx3i[i] =  (2*i+1)*(2*i+1) ;
00518     }
00519     }
00520     
00521     double* cx1t[2] ;
00522     double* cx2t[2] ;
00523     double* cx3t[2] ;
00524     cx1t[0] = cx1p ; cx1t[1] = cx1i ;
00525     cx2t[0] = cx2p ; cx2t[1] = cx2i ;
00526     cx3t[0] = cx3p ; cx3t[1] = cx3i ;
00527 
00528     // pt. sur le tableau de double resultat
00529     double* xo = new double[(tb->dim).taille] ;
00530     
00531     // Initialisation a zero :
00532     for (int i=0; i<(tb->dim).taille; i++) {
00533     xo[i] = 0 ; 
00534     }
00535     
00536     // On y va...
00537     double* xi = tb->t ;
00538     double* xci ;   // Pointeurs
00539     double* xco ;   //  courants
00540     
00541     // Position depart
00542     xci = xi ;
00543     xco = xo ;
00544     
00545     double *cx1, *cx2, *cx3 ;
00546 
00547     // k = 0
00548     cx1 = cx1t[0] ;
00549     cx2 = cx2t[0] ;
00550     cx3 = cx3t[0] ;
00551         for (int j=0 ; j<nt ; j++) {
00552 
00553         double som1 = 0 ;
00554         double som2 = 0 ;
00555         
00556         xco[nr-1] = 0 ;
00557         for (int i = nr-2 ; i >= 0 ; i-- ) {
00558             som1 += cx1[i] * xci[i+1] ;
00559             som2 += cx2[i] * xci[i+1] ;
00560             xco[i] = som1 - cx3[i] * som2 ;
00561         }   // Fin de la boucle sur r
00562         xco[0] *= .5 ;  // normalisation        
00563         xci += nr ; // au
00564         xco += nr ; //  suivant
00565         }   // Fin de la boucle sur theta
00566 
00567     // k = 1
00568     xci += nr*nt ;
00569     xco += nr*nt ;
00570     
00571     // k >= 2
00572     for (int k=2 ; k<np+1 ; k++) {
00573     int m = (k/2) % 2 ;
00574     cx1 = cx1t[m] ;
00575     cx2 = cx2t[m] ;
00576     cx3 = cx3t[m] ;
00577         for (int j=0 ; j<nt ; j++) {
00578 
00579         double som1 = 0 ;
00580         double som2 = 0 ;
00581         
00582         xco[nr-1] = 0 ;
00583         for (int i = nr-2 ; i >= 0 ; i-- ) {
00584             som1 += cx1[i] * xci[i+1] ;
00585             som2 += cx2[i] * xci[i+1] ;
00586             xco[i] = som1 - cx3[i] * som2 ;
00587         }   // Fin de la boucle sur r
00588         if (m == 0) xco[0] *= .5 ;  // normalisation eventuelle
00589         xci += nr ; // au
00590         xco += nr ; //  suivant
00591         }   // Fin de la boucle sur theta
00592     }   // Fin de la boucle sur phi
00593 
00594     // On remet les choses la ou il faut
00595     delete [] tb->t ;
00596     tb->t = xo ;
00597     
00598     // base de developpement
00599     // inchangee
00600 }
00601 
00602 // cas R_CHEBPIM_I
00603 //----------------
00604 void _d2sdx2_r_chebpim_i(Tbl *tb, int & )
00605 {
00606 
00607     // Peut-etre rien a faire ?
00608     if (tb->get_etat() == ETATZERO) {
00609     return ;
00610     }
00611     
00612     // Protection
00613     assert(tb->get_etat() == ETATQCQ) ;
00614     
00615     // Pour le confort
00616     int nr = (tb->dim).dim[0] ;     // Nombre
00617     int nt = (tb->dim).dim[1] ;     //   de points
00618     int np = (tb->dim).dim[2] ;     //      physiques REELS
00619     np = np - 2 ;           // Nombre de points physiques
00620     
00621      // Variables statiques
00622     static double* cx1p = 0x0 ;
00623     static double* cx2p = 0x0 ;
00624     static double* cx3p = 0x0 ;
00625     static double* cx1i = 0x0 ;
00626     static double* cx2i = 0x0 ;
00627     static double* cx3i = 0x0 ;
00628     static int nr_pre = 0 ;
00629 
00630     // Test sur np pour initialisation eventuelle
00631     if (nr > nr_pre) {
00632     nr_pre = nr ;
00633     if (cx1p != 0x0) delete [] cx1p ;
00634     if (cx2p != 0x0) delete [] cx2p ;
00635     if (cx3p != 0x0) delete [] cx3p ;
00636     if (cx1i != 0x0) delete [] cx1i ;
00637     if (cx2i != 0x0) delete [] cx2i ;
00638     if (cx3i != 0x0) delete [] cx3i ;
00639     cx1p = new double[nr] ;
00640     cx2p = new double[nr] ;
00641     cx3p = new double[nr] ;
00642     cx1i = new double[nr] ;
00643     cx2i = new double[nr] ;
00644     cx3i = new double[nr] ;
00645     for (int i=0 ; i<nr ; i++) {
00646         cx1p[i] =  8*(i+1)*(i+1)*(i+1) ;
00647         cx2p[i] =  2*(i+1) ;
00648         cx3p[i] =  4*i*i ;
00649 
00650         cx1i[i] =  (2*i+3)*(2*i+3)*(2*i+3) ;
00651         cx2i[i] =  (2*i+3) ;
00652         cx3i[i] =  (2*i+1)*(2*i+1) ;
00653     }
00654     }
00655 
00656     double* cx1t[2] ;
00657     double* cx2t[2] ;
00658     double* cx3t[2] ;
00659     cx1t[1] = cx1p ; cx1t[0] = cx1i ;
00660     cx2t[1] = cx2p ; cx2t[0] = cx2i ;
00661     cx3t[1] = cx3p ; cx3t[0] = cx3i ;
00662 
00663     // pt. sur le tableau de double resultat
00664     double* xo = new double[(tb->dim).taille] ;
00665     
00666     // Initialisation a zero :
00667     for (int i=0; i<(tb->dim).taille; i++) {
00668     xo[i] = 0 ; 
00669     }
00670     
00671     // On y va...
00672     double* xi = tb->t ;
00673     double* xci ;   // Pointeurs
00674     double* xco ;   //  courants
00675     
00676     // Position depart
00677     xci = xi ;
00678     xco = xo ;
00679 
00680     double *cx1, *cx2, *cx3 ;
00681 
00682     // k = 0
00683     cx1 = cx1t[0] ;
00684     cx2 = cx2t[0] ;
00685     cx3 = cx3t[0] ;
00686         for (int j=0 ; j<nt ; j++) {
00687 
00688         double som1 = 0 ;
00689         double som2 = 0 ;
00690         
00691         xco[nr-1] = 0 ;
00692         for (int i = nr-2 ; i >= 0 ; i-- ) {
00693             som1 += cx1[i] * xci[i+1] ;
00694             som2 += cx2[i] * xci[i+1] ;
00695             xco[i] = som1 - cx3[i] * som2 ;
00696         }   // Fin de la boucle sur r       
00697         xci += nr ;
00698         xco += nr ;
00699         }   // Fin de la boucle sur theta
00700 
00701     // k = 1
00702     xci += nr*nt ;
00703     xco += nr*nt ;
00704     
00705     // k >= 2
00706     for (int k=2 ; k<np+1 ; k++) {
00707     int m = (k/2) % 2 ;
00708     cx1 = cx1t[m] ;
00709     cx2 = cx2t[m] ;
00710     cx3 = cx3t[m] ;
00711         for (int j=0 ; j<nt ; j++) {
00712 
00713         double som1 = 0 ;
00714         double som2 = 0 ;
00715         
00716         xco[nr-1] = 0 ;
00717         for (int i = nr-2 ; i >= 0 ; i-- ) {
00718             som1 += cx1[i] * xci[i+1] ;
00719             som2 += cx2[i] * xci[i+1] ;
00720             xco[i] = som1 - cx3[i] * som2 ;
00721         }   // Fin de la boucle sur r
00722         if (m == 1) xco[0] *= .5 ;
00723         xci += nr ;
00724         xco += nr ;
00725         }   // Fin de la boucle sur theta
00726     }   // Fin de la boucle sur phi
00727 
00728     // On remet les choses la ou il faut
00729     delete [] tb->t ;
00730     tb->t = xo ;
00731     
00732     // base de developpement
00733     // inchangee
00734 }
00735 
00736 // cas R_CHEBPI_P
00737 //----------------
00738 void _d2sdx2_r_chebpi_p(Tbl *tb, int & )
00739 {
00740 
00741     // Peut-etre rien a faire ?
00742     if (tb->get_etat() == ETATZERO) {
00743     return ;
00744     }
00745     
00746     // Protection
00747     assert(tb->get_etat() == ETATQCQ) ;
00748     
00749     // Pour le confort
00750     int nr = (tb->dim).dim[0] ;     // Nombre
00751     int nt = (tb->dim).dim[1] ;     //   de points
00752     int np = (tb->dim).dim[2] ;     //      physiques REELS
00753     np = np - 2 ;           // Nombre de points physiques
00754     
00755     // Variables statiques
00756     static double* cx1p = 0x0 ;
00757     static double* cx2p = 0x0 ;
00758     static double* cx3p = 0x0 ;
00759     static double* cx1i = 0x0 ;
00760     static double* cx2i = 0x0 ;
00761     static double* cx3i = 0x0 ;
00762     static int nr_pre = 0 ;
00763 
00764     // Test sur np pour initialisation eventuelle
00765     if (nr > nr_pre) {
00766     nr_pre = nr ;
00767     if (cx1p != 0x0) delete [] cx1p ;
00768     if (cx2p != 0x0) delete [] cx2p ;
00769     if (cx3p != 0x0) delete [] cx3p ;
00770     if (cx1i != 0x0) delete [] cx1i ;
00771     if (cx2i != 0x0) delete [] cx2i ;
00772     if (cx3i != 0x0) delete [] cx3i ;
00773     cx1p = new double[nr] ;
00774     cx2p = new double[nr] ;
00775     cx3p = new double[nr] ;
00776     cx1i = new double[nr] ;
00777     cx2i = new double[nr] ;
00778     cx3i = new double[nr] ;
00779     for (int i=0 ; i<nr ; i++) {
00780         cx1p[i] =  8*(i+1)*(i+1)*(i+1) ;
00781         cx2p[i] =  2*(i+1) ;
00782         cx3p[i] =  4*i*i ;
00783 
00784         cx1i[i] =  (2*i+3)*(2*i+3)*(2*i+3) ;
00785         cx2i[i] =  (2*i+3) ;
00786         cx3i[i] =  (2*i+1)*(2*i+1) ;
00787     }
00788     }
00789     
00790     double* cx1t[2] ;
00791     double* cx2t[2] ;
00792     double* cx3t[2] ;
00793     cx1t[0] = cx1p ; cx1t[1] = cx1i ;
00794     cx2t[0] = cx2p ; cx2t[1] = cx2i ;
00795     cx3t[0] = cx3p ; cx3t[1] = cx3i ;
00796 
00797     // pt. sur le tableau de double resultat
00798     double* xo = new double[(tb->dim).taille] ;
00799     
00800     // Initialisation a zero :
00801     for (int i=0; i<(tb->dim).taille; i++) {
00802     xo[i] = 0 ; 
00803     }
00804     
00805     // On y va...
00806     double* xi = tb->t ;
00807     double* xci ;   // Pointeurs
00808     double* xco ;   //  courants
00809     
00810     // Position depart
00811     xci = xi ;
00812     xco = xo ;
00813     
00814     double *cx1, *cx2, *cx3 ;
00815 
00816     // k = 0
00817     for (int j=0 ; j<nt ; j++) {
00818       int l = j % 2 ;
00819       cx1 = cx1t[l] ;
00820       cx2 = cx2t[l] ;
00821       cx3 = cx3t[l] ;
00822       double som1 = 0 ;
00823       double som2 = 0 ;
00824       
00825       xco[nr-1] = 0 ;
00826       for (int i = nr-2 ; i >= 0 ; i-- ) {
00827     som1 += cx1[i] * xci[i+1] ;
00828     som2 += cx2[i] * xci[i+1] ;
00829     xco[i] = som1 - cx3[i] * som2 ;
00830       } // Fin de la boucle sur r
00831       if (l == 0) xco[0] *= .5 ;    // normalisation        
00832       xci += nr ;   // au
00833       xco += nr ;   //  suivant
00834     }   // Fin de la boucle sur theta
00835     
00836     // k = 1
00837     xci += nr*nt ;
00838     xco += nr*nt ;
00839     
00840     // k >= 2
00841     for (int k=2 ; k<np+1 ; k++) {
00842       for (int j=0 ; j<nt ; j++) {
00843     int l = j % 2 ;
00844     cx1 = cx1t[l] ;
00845     cx2 = cx2t[l] ;
00846     cx3 = cx3t[l] ;
00847     double som1 = 0 ;
00848     double som2 = 0 ;
00849     
00850     xco[nr-1] = 0 ;
00851     for (int i = nr-2 ; i >= 0 ; i-- ) {
00852       som1 += cx1[i] * xci[i+1] ;
00853       som2 += cx2[i] * xci[i+1] ;
00854       xco[i] = som1 - cx3[i] * som2 ;
00855     }   // Fin de la boucle sur r
00856     if (l == 0) xco[0] *= .5 ;  // normalisation eventuelle
00857     xci += nr ; // au
00858     xco += nr ; //  suivant
00859       }   // Fin de la boucle sur theta
00860     }   // Fin de la boucle sur phi
00861     
00862     // On remet les choses la ou il faut
00863     delete [] tb->t ;
00864     tb->t = xo ;
00865     
00866     // base de developpement
00867     // inchangee
00868 }
00869 
00870 // cas R_CHEBPI_I
00871 //----------------
00872 void _d2sdx2_r_chebpi_i(Tbl *tb, int & )
00873 {
00874 
00875     // Peut-etre rien a faire ?
00876     if (tb->get_etat() == ETATZERO) {
00877     return ;
00878     }
00879     
00880     // Protection
00881     assert(tb->get_etat() == ETATQCQ) ;
00882     
00883     // Pour le confort
00884     int nr = (tb->dim).dim[0] ;     // Nombre
00885     int nt = (tb->dim).dim[1] ;     //   de points
00886     int np = (tb->dim).dim[2] ;     //      physiques REELS
00887     np = np - 2 ;           // Nombre de points physiques
00888     
00889      // Variables statiques
00890     static double* cx1p = 0x0 ;
00891     static double* cx2p = 0x0 ;
00892     static double* cx3p = 0x0 ;
00893     static double* cx1i = 0x0 ;
00894     static double* cx2i = 0x0 ;
00895     static double* cx3i = 0x0 ;
00896     static int nr_pre = 0 ;
00897 
00898     // Test sur np pour initialisation eventuelle
00899     if (nr > nr_pre) {
00900     nr_pre = nr ;
00901     if (cx1p != 0x0) delete [] cx1p ;
00902     if (cx2p != 0x0) delete [] cx2p ;
00903     if (cx3p != 0x0) delete [] cx3p ;
00904     if (cx1i != 0x0) delete [] cx1i ;
00905     if (cx2i != 0x0) delete [] cx2i ;
00906     if (cx3i != 0x0) delete [] cx3i ;
00907     cx1p = new double[nr] ;
00908     cx2p = new double[nr] ;
00909     cx3p = new double[nr] ;
00910     cx1i = new double[nr] ;
00911     cx2i = new double[nr] ;
00912     cx3i = new double[nr] ;
00913     for (int i=0 ; i<nr ; i++) {
00914         cx1p[i] =  8*(i+1)*(i+1)*(i+1) ;
00915         cx2p[i] =  2*(i+1) ;
00916         cx3p[i] =  4*i*i ;
00917 
00918         cx1i[i] =  (2*i+3)*(2*i+3)*(2*i+3) ;
00919         cx2i[i] =  (2*i+3) ;
00920         cx3i[i] =  (2*i+1)*(2*i+1) ;
00921     }
00922     }
00923 
00924     double* cx1t[2] ;
00925     double* cx2t[2] ;
00926     double* cx3t[2] ;
00927     cx1t[1] = cx1p ; cx1t[0] = cx1i ;
00928     cx2t[1] = cx2p ; cx2t[0] = cx2i ;
00929     cx3t[1] = cx3p ; cx3t[0] = cx3i ;
00930 
00931     // pt. sur le tableau de double resultat
00932     double* xo = new double[(tb->dim).taille] ;
00933     
00934     // Initialisation a zero :
00935     for (int i=0; i<(tb->dim).taille; i++) {
00936     xo[i] = 0 ; 
00937     }
00938     
00939     // On y va...
00940     double* xi = tb->t ;
00941     double* xci ;   // Pointeurs
00942     double* xco ;   //  courants
00943     
00944     // Position depart
00945     xci = xi ;
00946     xco = xo ;
00947 
00948     double *cx1, *cx2, *cx3 ;
00949 
00950     // k = 0
00951     for (int j=0 ; j<nt ; j++) {
00952       int l = j % 2 ;
00953       cx1 = cx1t[l] ;
00954       cx2 = cx2t[l] ;
00955       cx3 = cx3t[l] ;
00956       double som1 = 0 ;
00957       double som2 = 0 ;
00958       
00959       xco[nr-1] = 0 ;
00960       for (int i = nr-2 ; i >= 0 ; i-- ) {
00961     som1 += cx1[i] * xci[i+1] ;
00962     som2 += cx2[i] * xci[i+1] ;
00963     xco[i] = som1 - cx3[i] * som2 ;
00964       } // Fin de la boucle sur r
00965       if (l == 1) xco[0] *= .5 ;
00966       xci += nr ;
00967       xco += nr ;
00968     }   // Fin de la boucle sur theta
00969     
00970     // k = 1
00971     xci += nr*nt ;
00972     xco += nr*nt ;
00973     
00974     // k >= 2
00975     for (int k=2 ; k<np+1 ; k++) {
00976       for (int j=0 ; j<nt ; j++) {
00977     int l = j % 2 ;
00978     cx1 = cx1t[l] ;
00979     cx2 = cx2t[l] ;
00980     cx3 = cx3t[l] ;
00981     double som1 = 0 ;
00982     double som2 = 0 ;
00983     
00984     xco[nr-1] = 0 ;
00985     for (int i = nr-2 ; i >= 0 ; i-- ) {
00986       som1 += cx1[i] * xci[i+1] ;
00987       som2 += cx2[i] * xci[i+1] ;
00988       xco[i] = som1 - cx3[i] * som2 ;
00989     }   // Fin de la boucle sur r
00990     if (l == 1) xco[0] *= .5 ;
00991     xci += nr ;
00992     xco += nr ;
00993       }   // Fin de la boucle sur theta
00994     }   // Fin de la boucle sur phi
00995     
00996     // On remet les choses la ou il faut
00997     delete [] tb->t ;
00998     tb->t = xo ;
00999     
01000     // base de developpement
01001     // inchangee
01002 }
01003 
01004 
01005 
01006 
01007 // cas R_JACO02
01008 //-----------
01009 void _d2sdx2_r_jaco02(Tbl *tb, int & )
01010 {
01011 
01012     // Peut-etre rien a faire ?
01013     if (tb->get_etat() == ETATZERO) {
01014     return ;
01015     }
01016     
01017     // Protection
01018     assert(tb->get_etat() == ETATQCQ) ;
01019     
01020     // Pour le confort
01021     int nr = (tb->dim).dim[0] ;     // Nombre
01022     int nt = (tb->dim).dim[1] ;     //   de points
01023     int np = (tb->dim).dim[2] ;     //      physiques REELS
01024     np = np - 2 ;           // Nombre de points physiques
01025     
01026     // Variables statiques
01027     static double* cx1 = 0x0 ;
01028     static double* cx2 = 0x0 ;
01029     static double* cx3 = 0x0 ;
01030     static int nr_pre = 0 ;
01031 
01032     // Test sur np pour initialisation eventuelle
01033     if (nr > nr_pre) {
01034     nr_pre = nr ;
01035     if (cx1 != 0x0) delete [] cx1 ;
01036     if (cx2 != 0x0) delete [] cx2 ;
01037     if (cx3 != 0x0) delete [] cx3 ;
01038     cx1 = new double [nr] ;
01039     cx2 = new double [nr] ;
01040     cx3 = new double [nr] ;  
01041     
01042     for (int i=0 ; i<nr ; i++) {
01043         cx1[i] =  (i+2)*(i+2)*(i+2) ;
01044         cx2[i] =  (i+2) ;
01045         cx3[i] =  i*i ;
01046     }
01047     }
01048     
01049     // pt. sur le tableau de double resultat
01050     double* xo = new double[(tb->dim).taille] ;
01051     
01052     // Initialisation a zero :
01053     for (int i=0; i<(tb->dim).taille; i++) {
01054     xo[i] = 0 ; 
01055     }
01056     
01057     // On y va...
01058     double* xi = tb->t ;
01059     double* xci = xi ;  // Pointeurs
01060     double* xco = xo ;  //  courants
01061     
01062     for (int k=0 ; k<np+1 ; k++) 
01063     if (k == 1)  {
01064         xci += nr*nt ;
01065         xco += nr*nt ;
01066         }
01067     else {
01068     for (int j=0 ; j<nt ; j++) {
01069 
01070         double* tb1 = new double[nr] ;
01071         for (int m =0 ; m<nr ; m++) { 
01072             tb1[m]=xci[m]; 
01073         }
01074         double* res = new double[nr] ;
01075         _d2sdx2_1d_r_jaco02(nr,tb1,res) ;
01076         for (int i = 0 ; i<nr ; i++ ) {
01077         xco[i] = res[i] ;
01078         }
01079         delete [] res ;
01080         delete [] tb1 ;     
01081         
01082         xci += nr ;
01083         xco += nr ;
01084     }   // Fin de la boucle sur theta
01085     }   // Fin de la boucle sur phi
01086     
01087     // On remet les choses la ou il faut
01088     delete [] tb->t ;
01089     tb->t = xo ;
01090     
01091     // base de developpement
01092     // inchangee
01093 }
01094 
01095
Generated on Tue Feb 7 01:35:18 2012 for LORENE by  doxygen 1.4.6