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
92namespace Lorene {
93void _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
105void _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
158void _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
211void _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
225void _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
268void _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
312void _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
367void _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
409void _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
455void 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}
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition type_parite.h:184
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_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_LEG
#define R_LEG
base de Legendre ordinaire (fin)
Definition type_parite.h:182
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67