LORENE
sx2_1d.C
1 /*
2  * Copyright (c) 1999-2001 Philippe Grandclement
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char sx2_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx2_1d.C,v 1.6 2015/03/05 08:49:32 j_novak Exp $" ;
24 
25 /*
26  * $Id: sx2_1d.C,v 1.6 2015/03/05 08:49:32 j_novak Exp $
27  * $Log: sx2_1d.C,v $
28  * Revision 1.6 2015/03/05 08:49:32 j_novak
29  * Implemented operators with Legendre bases.
30  *
31  * Revision 1.5 2014/10/13 08:53:27 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.4 2014/10/06 15:16:07 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.3 2007/12/11 15:28:18 jl_cornou
38  * Jacobi(0,2) polynomials partially implemented
39  *
40  * Revision 1.2 2002/10/16 14:36:59 j_novak
41  * Reorganization of #include instructions of standard C++, in order to
42  * use experimental version 3 of gcc.
43  *
44  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
45  * LORENE
46  *
47  * Revision 2.1 1999/07/08 09:55:21 phil
48  * correction gestion memoire
49  *
50  * Revision 2.0 1999/07/07 10:15:03 phil
51  * *** empty log message ***
52  *
53  *
54  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx2_1d.C,v 1.6 2015/03/05 08:49:32 j_novak Exp $
55  *
56  */
57 
58 // Includes
59 #include <cstdlib>
60 
61 
62 #include "headcpp.h"
63 #include "type_parite.h"
64 #include "proto.h"
65 
66 /*
67  * Operateurs non singuliers :
68  * R_CHEB : Identite
69  * R_CHEBP, R_CHEBI : (f-f(0)-f'(0)x)/x^2
70  * R_CHEBU : (f-f(1)-(x-1)*f'(1))/(x-1)^2
71  *
72  *Entree : tb contient coefficients de developpement en r
73  * nr : nombre de points en r.
74  *Sortie : tb contient les coefficients du resultat
75  *
76  */
77 namespace Lorene {
78 
79  void _sx_1d_r_legp(int, double* , double *) ;
80  void _sx_1d_r_legi(int, double* , double *) ;
81 
82 
83  //-----------------------------------
84  // Routine pour les cas non prevus --
85  //-----------------------------------
86 
87 void _sx2_1d_pas_prevu(int nr, double* tb, double *res) {
88  cout << "sx2 pas prevu..." << tb << " " << res << endl ;
89  cout << "nr : " << nr << endl ;
90  abort() ;
91  exit(-1) ;
92 }
93 
94 
95  //-------------
96  // Identite --
97  //------------
98 
99 void _sx2_1d_identite(int nr, double*tb ,double *res) {
100  for (int i=0 ; i<nr ; i++)
101  res[i] = tb[i] ;
102 }
103 
104 
105  //---------------
106  // cas R_CHEBP --
107  //---------------
108 
109 void _sx2_1d_r_chebp(int nr, double* tb, double *xo)
110 {
111 
112 
113 
114  double somp, somn ;
115  int sgn = 1 ;
116 
117  xo[nr-1] = 0 ;
118  somp = 4 * sgn * (nr-1) * tb[nr-1] ;
119  somn = 2 * sgn * tb[nr-1] ;
120  xo[nr-2] = somp - 2*(nr-2)*somn ;
121  for (int i = nr-3 ; i >= 0 ; i-- ) {
122  sgn = - sgn ;
123  somp += 4 * (i+1) * sgn * tb[i+1] ;
124  somn += 2 * sgn * tb[i+1] ;
125  xo[i] = somp - 2*i * somn ;
126  } // Fin de la premiere boucle sur r
127  for (int i=0 ; i<nr ; i+=2) {
128  xo[i] = -xo[i] ;
129  } // Fin de la deuxieme boucle sur r
130  xo[0] *= .5 ;
131 
132 }
133 
134 
135  //---------------
136  // cas R_CHEBI --
137  //---------------
138 
139 void _sx2_1d_r_chebi(int nr, double* tb, double *xo)
140 {
141  double somp, somn ;
142  int sgn = 1 ;
143 
144  xo[nr-1] = 0 ;
145  somp = 2 * sgn * (2*(nr-1)+1) * tb[nr-1] ;
146  somn = 2 * sgn * tb[nr-1] ;
147  xo[nr-2] = somp - (2*(nr-2)+1)*somn ;
148  for (int i = nr-3 ; i >= 0 ; i-- ) {
149  sgn = - sgn ;
150  somp += 2 * (2*(i+1)+1) * sgn * tb[i+1] ;
151  somn += 2 * sgn * tb[i+1] ;
152  xo[i] = somp - (2*i+1) * somn ;
153  } // Fin de la premiere boucle sur r
154  for (int i=0 ; i<nr ; i+=2) {
155  xo[i] = -xo[i] ;
156  } // Fin de la deuxieme boucle sur r
157 
158 }
159  //----------------
160  // cas R_CHEBU --
161  //---------------
162 
163 void _sxm12_1d_r_chebu(int nr, double* tb, double *xo) {
164 
165  // On appelle juste deux fois la routine existante :
166  sxm1_1d_cheb(nr, tb) ;
167  sxm1_1d_cheb(nr, tb) ;
168  for (int i=0 ; i<nr ; i++)
169  xo[i] = tb[i] ;
170 }
171 
172  //--------------
173  // cas R_LEGP --
174  //--------------
175 
176 void _sx2_1d_r_legp(int nr, double* tb, double *xo)
177 {
178  // On appelle juste deux fois la routine existante :
179  double* interm = new double[nr] ;
180  _sx_1d_r_legp(nr, tb, interm) ;
181  _sx_1d_r_legi(nr, interm, xo) ;
182  delete [] interm ;
183 
184 }
185 
186 
187  //--------------
188  // cas R_LEGI --
189  //--------------
190 
191 void _sx2_1d_r_legi(int nr, double* tb, double *xo)
192 {
193  // On appelle juste deux fois la routine existante :
194  double* interm = new double[nr] ;
195  _sx_1d_r_legi(nr, tb, interm) ;
196  _sx_1d_r_legp(nr, interm, xo) ;
197  delete [] interm ;
198 }
199 
200  //----------------------
201  // La routine a appeler
202  //----------------------
203 
204 void sx2_1d(int nr, double **tb, int base_r) // Version appliquee a this
205 {
206 
207 // Routines de derivation
208 static void (*sx2_1d[MAX_BASE])(int, double *, double*) ;
209 static int nap = 0 ;
210 
211  // Premier appel
212  if (nap==0) {
213  nap = 1 ;
214  for (int i=0 ; i<MAX_BASE ; i++) {
215  sx2_1d[i] = _sx2_1d_pas_prevu ;
216  }
217  // Les routines existantes
218  sx2_1d[R_CHEB >> TRA_R] = _sx2_1d_identite ;
219  sx2_1d[R_CHEBP >> TRA_R] = _sx2_1d_r_chebp ;
220  sx2_1d[R_CHEBI >> TRA_R] = _sx2_1d_r_chebi ;
221  sx2_1d[R_LEG >> TRA_R] = _sx2_1d_identite ;
222  sx2_1d[R_LEGP >> TRA_R] = _sx2_1d_r_legp ;
223  sx2_1d[R_LEGI >> TRA_R] = _sx2_1d_r_legi ;
224  sx2_1d[R_CHEBU >> TRA_R] = _sxm12_1d_r_chebu ;
225  }
226 
227  double *result = new double[nr] ;
228  sx2_1d[base_r](nr, *tb, result) ;
229 
230  delete [] (*tb) ;
231  (*tb) = result ;
232 }
233 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166