LORENE
dsdx_1d.C
1  /*
2  * Copyright (c) 2004 Jerome Novak
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 dsdx_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/dsdx_1d.C,v 1.9 2015/03/25 15:03:00 j_novak Exp $" ;
24 
25 /*
26  * $Id: dsdx_1d.C,v 1.9 2015/03/25 15:03:00 j_novak Exp $
27  * $Log: dsdx_1d.C,v $
28  * Revision 1.9 2015/03/25 15:03:00 j_novak
29  * Correction of a bug...
30  *
31  * Revision 1.8 2015/03/05 08:49:31 j_novak
32  * Implemented operators with Legendre bases.
33  *
34  * Revision 1.7 2014/10/13 08:53:23 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.6 2014/10/06 15:16:06 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.5 2007/12/12 12:30:48 jl_cornou
41  * *** empty log message ***
42  *
43  * Revision 1.4 2007/12/11 15:28:18 jl_cornou
44  * Jacobi(0,2) polynomials partially implemented
45  *
46  * Revision 1.3 2006/04/10 15:19:20 j_novak
47  * New definition of 1D operators dsdx and sx in the nucleus (bases R_CHEBP and
48  * R_CHEBI).
49  *
50  * Revision 1.2 2005/01/10 16:34:53 j_novak
51  * New class for 1D mono-domain differential operators.
52  *
53  * Revision 1.1 2004/02/06 10:53:53 j_novak
54  * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
55  *
56  *
57  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/dsdx_1d.C,v 1.9 2015/03/25 15:03:00 j_novak Exp $
58  *
59  */
60 
61 #include <cmath>
62 #include <cstdlib>
63 #include "type_parite.h"
64 #include "headcpp.h"
65 #include "proto.h"
66 
67 /*
68  * Routine appliquant l'operateur dsdx.
69  *
70  * Entree : tb contient les coefficients du developpement
71  * int nr : nombre de points en r.
72  *
73  * Sortie : tb contient dsdx
74  *
75  */
76 
77  //----------------------------------
78  // Routine pour les cas non prevus --
79  //----------------------------------
80 
81 namespace Lorene {
82 void _dsdx_1d_pas_prevu(int nr, double* tb, double *xo) {
83  cout << "dsdx pas prevu..." << endl ;
84  cout << "Nombre de points : " << nr << endl ;
85  cout << "Valeurs : " << tb << " " << xo <<endl ;
86  abort() ;
87  exit(-1) ;
88 }
89 
90  //----------------
91  // cas R_CHEBU ---
92  //----------------
93 
94 void _dsdx_1d_r_chebu(int nr, double* tb, double *xo)
95 {
96 
97  double som ;
98 
99  xo[nr-1] = 0 ;
100  som = 2*(nr-1) * tb[nr-1] ;
101  xo[nr-2] = som ;
102  for (int i = nr-4 ; i >= 0 ; i -= 2 ) {
103  som += 2*(i+1) * tb[i+1] ;
104  xo[i] = som ;
105  } // Fin de la premiere boucle sur r
106  som = 2*(nr-2) * tb[nr-2] ;
107  xo[nr-3] = som ;
108  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
109  som += 2*(i+1) * tb[i+1] ;
110  xo[i] = som ;
111  } // Fin de la deuxieme boucle sur r
112  xo[0] *= .5 ;
113 
114 }
115 
116 
117  //----------------
118  // cas R_JACO02 --
119  //----------------
120 
121 void _dsdx_1d_r_jaco02(int nr, double* tb, double *xo)
122 {
123 
124  double som ;
125 
126  xo[nr-1] = 0 ;
127  for (int i = 0 ; i < nr-1 ; i++ ) {
128  som = 0 ;
129  for ( int j = i+1 ; j < nr ; j++ ) {
130  som += (1 - pow(double(-1),(j-i))*(i+1)*(i+2)/double((j+1)*(j+2)))*tb[j] ;
131  } // Fin de la boucle auxiliaire
132  xo[i] = (i+1.5)*som ;
133  } // Fin de la boucle sur R
134 
135 }
136 
137  //----------------
138  // cas R_CHEBI ---
139  //----------------
140 
141 void _dsdx_1d_r_chebi(int nr, double* tb, double *xo)
142 {
143 
144  double som ;
145 
146  xo[nr-1] = 0 ;
147  som = 2*(2*nr-3) * tb[nr-2] ;
148  xo[nr-2] = som ;
149  for (int i = nr-3 ; i >= 0 ; i -- ) {
150  som += 2*(2*i+1) * tb[i] ;
151  xo[i] = som ;
152  }
153  xo[0] *= .5 ;
154 }
155 
156  //----------------
157  // cas R_CHEBP ---
158  //----------------
159 
160 void _dsdx_1d_r_chebp(int nr, double* tb, double *xo)
161 {
162 
163  double som ;
164 
165  xo[nr-1] = 0 ;
166  som = 4*(nr-1) * tb[nr-1] ;
167  xo[nr-2] = som ;
168  for (int i = nr-3 ; i >= 0 ; i --) {
169  som += 4*(i+1) * tb[i+1] ;
170  xo[i] = som ;
171  }
172 }
173 
174  //--------------
175  // cas R_LEG ---
176  //--------------
177 
178 void _dsdx_1d_r_leg(int nr, double* tb, double *xo)
179 {
180  double som ;
181 
182  xo[nr-1] = 0 ;
183  som = tb[nr-1] ;
184  xo[nr-2] = double(2*nr-3)*som ;
185  for (int i = nr-4 ; i >= 0 ; i -= 2 ) {
186  som += tb[i+1] ;
187  xo[i] = double(2*i+1)*som ;
188  } // Fin de la premiere boucle sur r
189  som = tb[nr-2] ;
190  if (nr > 2) xo[nr-3] = double(2*nr-5)*som ;
191  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
192  som += tb[i+1] ;
193  xo[i] = double(2*i+1)*som ;
194  } // Fin de la deuxieme boucle sur r
195 }
196 
197  //---------------
198  // cas R_LEGI ---
199  //---------------
200 
201 void _dsdx_1d_r_legi(int nr, double* tb, double *xo)
202 {
203 
204  double som ;
205 
206  xo[nr-1] = 0 ;
207  som = tb[nr-2] ;
208  if (nr > 1) xo[nr-2] = double(4*nr-7)*som ;
209  for (int i = nr-3 ; i >= 0 ; i -- ) {
210  som += tb[i] ;
211  xo[i] = double(4*i+1)*som ;
212  }
213 }
214 
215  //---------------
216  // cas R_LEGP ---
217  //---------------
218 
219 void _dsdx_1d_r_legp(int nr, double* tb, double *xo)
220 {
221 
222  double som ;
223 
224  xo[nr-1] = 0 ;
225  som = tb[nr-1] ;
226  if (nr > 1) xo[nr-2] = double(4*nr-5)*som ;
227  for (int i = nr-3 ; i >= 0 ; i --) {
228  som += tb[i+1] ;
229  xo[i] = double(4*i+3)*som ;
230  }
231 }
232 
233  // ---------------------
234  // La routine a appeler
235  //----------------------
236 
237 
238 void dsdx_1d(int nr, double** tb, int base_r)
239 {
240 
241  // Routines de derivation
242  static void (*dsdx_1d[MAX_BASE])(int, double*, double *) ;
243  static int nap = 0 ;
244 
245  // Premier appel
246  if (nap==0) {
247  nap = 1 ;
248  for (int i=0 ; i<MAX_BASE ; i++) {
249  dsdx_1d[i] = _dsdx_1d_pas_prevu ;
250  }
251  // Les routines existantes
252  dsdx_1d[R_CHEBU >> TRA_R] = _dsdx_1d_r_chebu ;
253  dsdx_1d[R_CHEBP >> TRA_R] = _dsdx_1d_r_chebp ;
254  dsdx_1d[R_CHEBI >> TRA_R] = _dsdx_1d_r_chebi ;
255  dsdx_1d[R_CHEB >> TRA_R] = _dsdx_1d_r_cheb ;
256  dsdx_1d[R_LEGP >> TRA_R] = _dsdx_1d_r_legp ;
257  dsdx_1d[R_LEGI >> TRA_R] = _dsdx_1d_r_legi ;
258  dsdx_1d[R_LEG >> TRA_R] = _dsdx_1d_r_leg ;
259  dsdx_1d[R_JACO02 >> TRA_R] = _dsdx_1d_r_jaco02 ;
260 
261  }
262 
263  double *result = new double[nr] ;
264 
265  dsdx_1d[base_r](nr, *tb, result) ;
266 
267  delete [] (*tb) ;
268  (*tb) = result ;
269 }
270 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
R_LEG
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182