LORENE
d2sdx2_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 
24 
25 /*
26  * $Id: d2sdx2_1d.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
27  * $Log: d2sdx2_1d.C,v $
28  * Revision 1.8 2016/12/05 16:18:07 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.7 2015/03/05 08:49:31 j_novak
32  * Implemented operators with Legendre bases.
33  *
34  * Revision 1.6 2014/10/13 08:53:23 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.5 2014/10/06 15:16:06 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.4 2007/12/11 15:28:18 jl_cornou
41  * Jacobi(0,2) polynomials partially implemented
42  *
43  * Revision 1.3 2002/10/16 15:05:54 j_novak
44  * *** empty log message ***
45  *
46  * Revision 1.2 2002/10/16 14:36:58 j_novak
47  * Reorganization of #include instructions of standard C++, in order to
48  * use experimental version 3 of gcc.
49  *
50  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
51  * LORENE
52  *
53  * Revision 2.4 1999/10/11 15:18:14 phil
54  * *** empty log message ***
55  *
56  * Revision 2.3 1999/10/11 14:27:06 phil
57  * initialisation des variables statiques
58  *
59  * Revision 2.2 1999/09/03 14:15:56 phil
60  * Correction termes 0 (/2)
61  *
62  * Revision 2.1 1999/07/08 09:53:13 phil
63  * correction gestion memoire
64  *
65  * Revision 2.0 1999/07/07 10:15:26 phil
66  * *** empty log message ***
67  *
68  *
69  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/d2sdx2_1d.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
70  *
71  */
72 
73 #include <cstdlib>
74 #include "type_parite.h"
75 #include "headcpp.h"
76 #include "proto.h"
77 
78 /*
79  * Routine appliquant l'operateur d2sdx2.
80  *
81  * Entree : tb contient les coefficients du developpement
82  * int nr : nombre de points en r.
83  *
84  * Sortie : tb contient d2sdx2
85  *
86  */
87 
88  //----------------------------------
89  // Routine pour les cas non prevus --
90  //----------------------------------
91 
92 namespace Lorene {
93 void _d2sdx2_1d_pas_prevu(int nr, double* tb, double *xo) {
94  cout << "d2sdx2 pas prevu..." << endl ;
95  cout << "Nombre de points : " << nr << endl ;
96  cout << "Valeurs : " << tb << " " << xo <<endl ;
97  abort() ;
98  exit(-1) ;
99 }
100 
101  //----------------
102  // cas R_CHEBU ---
103  //----------------
104 
105 void _d2sdx2_1d_r_chebu(int nr, double* tb, double *xo)
106 {
107  // Variables statiques
108  static double* cx1 = 0x0 ;
109  static double* cx2 = 0x0 ;
110  static double* cx3 = 0x0 ;
111  static int nr_pre = 0 ;
112 
113  // Test sur np pour initialisation eventuelle
114  if (nr > nr_pre) {
115  nr_pre = nr ;
116  if (cx1 != 0x0) delete [] cx1 ;
117  if (cx2 != 0x0) delete [] cx2 ;
118  if (cx3 != 0x0) delete [] cx3 ;
119  cx1 = new double [nr] ;
120  cx2 = new double [nr] ;
121  cx3 = new double [nr] ;
122  for (int i=0 ; i<nr ; i++) {
123  cx1[i] = (i+2)*(i+2)*(i+2) ;
124  cx2[i] = (i+2) ;
125  cx3[i] = i*i ;
126  }
127  }
128 
129  double som1, som2 ;
130 
131  xo[nr-1] = 0 ;
132  som1 = (nr-1)*(nr-1)*(nr-1) * tb[nr-1] ;
133  som2 = (nr-1) * tb[nr-1] ;
134  xo[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
135  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
136  som1 += cx1[i] * tb[i+2] ;
137  som2 += cx2[i] * tb[i+2] ;
138  xo[i] = som1 - cx3[i] * som2 ;
139  } // Fin de la premiere boucle sur r
140 
141  xo[nr-2] = 0 ;
142  som1 = (nr-2)*(nr-2)*(nr-2) * tb[nr-2] ;
143  som2 = (nr-2) * tb[nr-2] ;
144  xo[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
145  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
146  som1 += cx1[i] * tb[i+2] ;
147  som2 += cx2[i] * tb[i+2] ;
148  xo[i] = som1 - cx3[i] * som2 ;
149  } // Fin de la deuxieme boucle sur r
150  xo[0] *= 0.5 ;
151 
152 }
153 
154  //---------------
155  // cas R_CHEB ---
156  //---------------
157 
158 void _d2sdx2_1d_r_cheb(int nr, double* tb, double *xo)
159 {
160 
161  // Variables statiques
162  static double* cx1 = 0x0 ;
163  static double* cx2 = 0x0 ;
164  static double* cx3 = 0x0 ;
165  static int nr_pre = 0 ;
166 
167  // Test sur np pour initialisation eventuelle
168  if (nr > nr_pre) {
169  nr_pre = nr ;
170  if (cx1 != 0x0) delete [] cx1 ;
171  if (cx2 != 0x0) delete [] cx2 ;
172  if (cx3 != 0x0) delete [] cx3 ;
173  cx1 = new double [nr] ;
174  cx2 = new double [nr] ;
175  cx3 = new double [nr] ;
176  for (int i=0 ; i<nr ; i++) {
177  cx1[i] = (i+2)*(i+2)*(i+2) ;
178  cx2[i] = (i+2) ;
179  cx3[i] = i*i ;
180  }
181  }
182 
183  double som1, som2 ;
184 
185  xo[nr-1] = 0 ;
186  som1 = (nr-1)*(nr-1)*(nr-1) * tb[nr-1] ;
187  som2 = (nr-1) * tb[nr-1] ;
188  xo[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
189  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
190  som1 += cx1[i] * tb[i+2] ;
191  som2 += cx2[i] * tb[i+2] ;
192  xo[i] = som1 - cx3[i] * som2 ;
193  } // Fin de la premiere boucle sur r
194  xo[nr-2] = 0 ;
195  som1 = (nr-2)*(nr-2)*(nr-2) * tb[nr-2] ;
196  som2 = (nr-2) * tb[nr-2] ;
197  xo[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
198  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
199  som1 += cx1[i] * tb[i+2] ;
200  som2 += cx2[i] * tb[i+2] ;
201  xo[i] = som1 - cx3[i] * som2 ;
202  } // Fin de la deuxieme boucle sur r
203  xo[0] *= .5 ;
204 
205 }
206 
207  //----------------
208  // cas R_JACO02 --
209  //----------------
210 
211 void _d2sdx2_1d_r_jaco02(int nr, double* tb, double *xo)
212 {
213  dsdx_1d(nr, &tb, R_JACO02 >> TRA_R) ;
214  dsdx_1d(nr, &tb, R_JACO02 >> TRA_R) ;
215  for (int i = 0 ; i<nr ; i++) {
216  xo[i] = tb[i] ;
217  }
218 }
219 
220 
221  //-----------------
222  // cas R_CHEBP ---
223  //----------------
224 
225 void _d2sdx2_1d_r_chebp(int nr, double* tb, double *xo)
226 {
227  // Variables statiques
228  static double* cx1 = 0x0 ;
229  static double* cx2 = 0x0 ;
230  static double* cx3 = 0x0 ;
231  static int nr_pre = 0 ;
232 
233  // Test sur np pour initialisation eventuelle
234  if (nr > nr_pre) {
235  nr_pre = nr ;
236  if (cx1 != 0x0) delete [] cx1 ;
237  if (cx2 != 0x0) delete [] cx2 ;
238  if (cx3 != 0x0) delete [] cx3 ;
239  cx1 = new double [nr] ;
240  cx2 = new double [nr] ;
241  cx3 = new double [nr] ;
242  for (int i=0 ; i<nr ; i++) {
243  cx1[i] = 8*(i+1)*(i+1)*(i+1) ;
244  cx2[i] = 2*(i+1) ;
245  cx3[i] = 4*i*i ;
246  }
247  }
248 
249  double som1, som2 ;
250 
251  xo[nr-1] = 0 ;
252  som1 = 8*(nr-1)*(nr-1)*(nr-1) * tb[nr-1] ;
253  som2 = 2*(nr-1) * tb[nr-1] ;
254  xo[nr-2] = som1 - 4*(nr-2)*(nr-2)*som2 ;
255 
256  for (int i = nr-3 ; i >= 0 ; i-- ) {
257  som1 += cx1[i] * tb[i+1] ;
258  som2 += cx2[i] * tb[i+1] ;
259  xo[i] = som1 - cx3[i] * som2 ;
260  } // Fin de la boucle sur r
261  xo[0] *= .5 ;
262 }
263 
264  //---------------
265  // cas R_CHEBI --
266  //---------------
267 
268 void _d2sdx2_1d_r_chebi(int nr, double* tb, double *xo)
269 {
270  // Variables statiques
271  static double* cx1 = 0x0 ;
272  static double* cx2 = 0x0 ;
273  static double* cx3 = 0x0 ;
274  static int nr_pre = 0 ;
275 
276  // Test sur np pour initialisation eventuelle
277 
278  if (nr > nr_pre) {
279  nr_pre = nr ;
280  if (cx1 != 0x0) delete [] cx1 ;
281  if (cx2 != 0x0) delete [] cx2 ;
282  if (cx3 != 0x0) delete [] cx3 ;
283  cx1 = new double [nr] ;
284  cx2 = new double [nr] ;
285  cx3 = new double [nr] ;
286  for (int i=0 ; i<nr ; i++) {
287  cx1[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
288  cx2[i] = (2*i+3) ;
289  cx3[i] = (2*i+1)*(2*i+1) ;
290  }
291  }
292 
293  // pt. sur le tableau de double resultat
294  double som1, som2 ;
295 
296  xo[nr-1] = 0 ;
297  som1 = (2*nr-1)*(2*nr-1)*(2*nr-1) * tb[nr-1] ;
298  som2 = (2*nr-1) * tb[nr-1] ;
299  xo[nr-2] = som1 - (2*nr-3)*(2*nr-3)*som2 ;
300  for (int i = nr-3 ; i >= 0 ; i-- ) {
301  som1 += cx1[i] * tb[i+1] ;
302  som2 += cx2[i] * tb[i+1] ;
303  xo[i] = som1 - cx3[i] * som2 ;
304  } // Fin de la boucle su r
305 
306 }
307 
308  //--------------
309  // cas R_LEG ---
310  //--------------
311 
312 void _d2sdx2_1d_r_leg(int nr, double* tb, double *xo)
313 {
314 
315  // Variables statiques
316  static double* cx1 = 0x0 ;
317  static double* cx2 = 0x0 ;
318  static double* cx3 = 0x0 ;
319  static int nr_pre = 0 ;
320 
321  // Test sur np pour initialisation eventuelle
322  if (nr > nr_pre) {
323  nr_pre = nr ;
324  if (cx1 != 0x0) delete [] cx1 ;
325  if (cx2 != 0x0) delete [] cx2 ;
326  if (cx3 != 0x0) delete [] cx3 ;
327  cx1 = new double [nr] ;
328  cx2 = new double [nr] ;
329  cx3 = new double [nr] ;
330  for (int i=0 ; i<nr ; i++) {
331  cx1[i] = (i+2)*(i+3) ;
332  cx2[i] = i*(i+1) ;
333  cx3[i] = double(i) + 0.5 ;
334  }
335  }
336 
337  double som1, som2 ;
338 
339  xo[nr-1] = 0 ;
340  som1 = (nr-1)* nr * tb[nr-1] ;
341  som2 = tb[nr-1] ;
342  if (nr > 2)
343  xo[nr-3] = (double(nr) - 2.5) * (som1 - (nr-3)*(nr-2)*som2) ;
344  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
345  som1 += cx1[i] * tb[i+2] ;
346  som2 += tb[i+2] ;
347  xo[i] = cx3[i]*(som1 - cx2[i] * som2) ;
348  } // Fin de la premiere boucle sur r
349  if (nr > 1) xo[nr-2] = 0 ;
350  if (nr > 3) {
351  som1 = (nr-2)*(nr-1) * tb[nr-2] ;
352  som2 = tb[nr-2] ;
353  xo[nr-4] = (double(nr) - 3.5) * (som1 - (nr-4)*(nr-3)*som2) ;
354  }
355  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
356  som1 += cx1[i] * tb[i+2] ;
357  som2 += tb[i+2] ;
358  xo[i] = cx3[i]*(som1 - cx2[i] * som2) ;
359  } // Fin de la deuxieme boucle sur r
360 
361 }
362 
363  //---------------
364  // cas R_LEGP ---
365  //---------------
366 
367 void _d2sdx2_1d_r_legp(int nr, double* tb, double *xo)
368 {
369  // Variables statiques
370  static double* cx1 = 0x0 ;
371  static double* cx2 = 0x0 ;
372  static double* cx3 = 0x0 ;
373  static int nr_pre = 0 ;
374 
375  // Test sur np pour initialisation eventuelle
376  if (nr > nr_pre) {
377  nr_pre = nr ;
378  if (cx1 != 0x0) delete [] cx1 ;
379  if (cx2 != 0x0) delete [] cx2 ;
380  if (cx3 != 0x0) delete [] cx3 ;
381  cx1 = new double [nr] ;
382  cx2 = new double [nr] ;
383  cx3 = new double [nr] ;
384  for (int i=0 ; i<nr ; i++) {
385  cx1[i] = (2*i+2)*(2*i+3) ;
386  cx2[i] = 2*i*2*(i+1) ;
387  cx3[i] = double(2*i)+ 0.5 ;
388  }
389  }
390 
391  double som1, som2 ;
392 
393  xo[nr-1] = 0 ;
394  som1 = (2*nr-2)*(2*nr-1)* tb[nr-1] ;
395  som2 = tb[nr-1] ;
396  if (nr > 1)
397  xo[nr-2] = (double(2*nr) - 1.5)*(som1 - 2*(nr-2)*(2*nr-1)*som2) ;
398  for (int i = nr-3 ; i >= 0 ; i-- ) {
399  som1 += cx1[i] * tb[i+1] ;
400  som2 += tb[i+1] ;
401  xo[i] = cx3[i]*(som1 - cx2[i]*som2) ;
402  } // Fin de la boucle sur r
403 }
404 
405  //---------------
406  // cas R_LEGI --
407  //---------------
408 
409 void _d2sdx2_1d_r_legi(int nr, double* tb, double *xo)
410 {
411  // Variables statiques
412  static double* cx1 = 0x0 ;
413  static double* cx2 = 0x0 ;
414  static double* cx3 = 0x0 ;
415  static int nr_pre = 0 ;
416 
417  // Test sur np pour initialisation eventuelle
418 
419  if (nr > nr_pre) {
420  nr_pre = nr ;
421  if (cx1 != 0x0) delete [] cx1 ;
422  if (cx2 != 0x0) delete [] cx2 ;
423  if (cx3 != 0x0) delete [] cx3 ;
424  cx1 = new double [nr] ;
425  cx2 = new double [nr] ;
426  cx3 = new double [nr] ;
427  for (int i=0 ; i<nr ; i++) {
428  cx1[i] = (2*i+3)*(2*i+4) ;
429  cx2[i] = (2*i+1)*(2*i+2) ;
430  cx3[i] = double(2*i) + 1.5 ;
431  }
432  }
433 
434  // pt. sur le tableau de double resultat
435  double som1, som2 ;
436 
437  xo[nr-1] = 0 ;
438  som1 = (2*nr-1)*(2*nr) * tb[nr-1] ;
439  som2 = tb[nr-1] ;
440  if (nr > 1)
441  xo[nr-2] = (double(nr) - 1.5)*(som1 - (2*nr-3)*(2*nr-2)*som2) ;
442  for (int i = nr-3 ; i >= 0 ; i-- ) {
443  som1 += cx1[i] * tb[i+1] ;
444  som2 += tb[i+1] ;
445  xo[i] = cx3[i]*(som1 - cx2[i] * som2) ;
446  } // Fin de la boucle sur r
447 }
448 
449 
450  // ---------------------
451  // La routine a appeler
452  //----------------------
453 
454 
455 void d2sdx2_1d(int nr, double** tb, int base_r)
456 {
457 
458  // Routines de derivation
459  static void (*d2sdx2_1d[MAX_BASE])(int, double*, double *) ;
460  static int nap = 0 ;
461 
462  // Premier appel
463  if (nap==0) {
464  nap = 1 ;
465  for (int i=0 ; i<MAX_BASE ; i++) {
466  d2sdx2_1d[i] = _d2sdx2_1d_pas_prevu ;
467  }
468  // Les routines existantes
469  d2sdx2_1d[R_CHEB >> TRA_R] = _d2sdx2_1d_r_cheb ;
470  d2sdx2_1d[R_CHEBU >> TRA_R] = _d2sdx2_1d_r_chebu ;
471  d2sdx2_1d[R_CHEBP >> TRA_R] = _d2sdx2_1d_r_chebp ;
472  d2sdx2_1d[R_CHEBI >> TRA_R] = _d2sdx2_1d_r_chebi ;
473  d2sdx2_1d[R_LEG >> TRA_R] = _d2sdx2_1d_r_leg ;
474  d2sdx2_1d[R_LEGP >> TRA_R] = _d2sdx2_1d_r_legp ;
475  d2sdx2_1d[R_LEGI >> TRA_R] = _d2sdx2_1d_r_legi ;
476  d2sdx2_1d[R_JACO02 >> TRA_R] = _d2sdx2_1d_r_jaco02 ;
477  }
478 
479  double *result = new double[nr] ;
480 
481  d2sdx2_1d[base_r](nr, *tb, result) ;
482 
483  delete [] (*tb) ;
484  (*tb) = result ;
485 }
486 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67
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
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_LEG
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166