LORENE
op_dsdphi.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24
25
26/*
27 * Ensemble des routines de base pour la derivation par rapport a phi
28 * (Utilisation interne)
29 *
30 * void _dsdphi_XXXX(Tbl * t, int & b)
31 * t pointeur sur le Tbl d'entree-sortie
32 * b base spectrale
33 *
34 */
35
36/*
37 * $Id: op_dsdphi.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
38 * $Log: op_dsdphi.C,v $
39 * Revision 1.7 2016/12/05 16:18:07 j_novak
40 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
41 *
42 * Revision 1.6 2014/10/13 08:53:25 j_novak
43 * Lorene classes and functions now belong to the namespace Lorene.
44 *
45 * Revision 1.5 2013/04/25 15:46:06 j_novak
46 * Added special treatment in the case np = 1, for type_p = NONSYM.
47 *
48 * Revision 1.4 2006/03/10 12:45:38 j_novak
49 * Use of C++-style cast.
50 *
51 * Revision 1.3 2003/12/19 16:21:48 j_novak
52 * Shadow hunt
53 *
54 * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
55 * Modification of declaration of Fortran 77 prototypes for
56 * a better portability (in particular on IBM AIX systems):
57 * All Fortran subroutine names are now written F77_* and are
58 * defined in the new file C++/Include/proto_f77.h.
59 *
60 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
61 * LORENE
62 *
63 * Revision 2.8 2000/11/14 15:08:16 eric
64 * Traitement du cas np=1 dans P_COSSIN_I.
65 *
66 * Revision 2.7 2000/10/04 15:54:09 eric
67 * *** empty log message ***
68 *
69 * Revision 2.6 2000/10/04 12:50:38 eric
70 * Ajout de la base P_COSSIN_I.
71 *
72 * Revision 2.5 2000/10/04 11:50:32 eric
73 * Ajout d' abort() dans le cas non prevu.
74 *
75 * Revision 2.4 2000/08/18 13:20:01 eric
76 * Traitement du cas np = 1.
77 *
78 * Revision 2.3 2000/02/24 16:41:03 eric
79 * Initialisation a zero du tableau xo avant le calcul.
80 *
81 * Revision 2.2 1999/11/15 16:38:03 eric
82 * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
83 *
84 * Revision 2.1 1999/02/23 11:36:34 hyc
85 * *** empty log message ***
86 *
87 * Revision 2.0 1999/02/22 17:24:59 hyc
88 * *** empty log message ***
89 *
90 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdphi.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
91 *
92 */
93
94// Headers Lorene
95#include "tbl.h"
96#include "proto_f77.h"
97
98
99// Routine pour les cas non prevus
100//--------------------------------
101namespace Lorene {
102void _dsdphi_pas_prevu(Tbl* , int & b) {
103 cout << "Unknown phi basis in Mtbl_cf::dsdp() !" << endl ;
104 cout << " basis: " << hex << b << endl ;
105 abort() ;
106}
107
108 //------------------//
109 // cas P_COSSIN //
110 //------------------//
111
112void _dsdphi_p_cossin(Tbl* tb, int & )
113{
114
115 // Peut-etre rien a faire ?
116 if (tb->get_etat() == ETATZERO) {
117 return ;
118 }
119
120 // Protection
121 assert(tb->get_etat() == ETATQCQ) ;
122
123 // Pour le confort
124 int nr = (tb->dim).dim[0] ; // Nombre
125 int nt = (tb->dim).dim[1] ; // de points
126 int np = (tb->dim).dim[2] ; // physiques REELS
127 np = np - 2 ; // Nombre de points physiques
128
129 // Cas particulier de la symetrie axiale :
130 if (np == 1) {
131 tb->set_etat_zero() ;
132 return ;
133 }
134
135 // Variables statiques
136 static double* cx = 0 ;
137 static int np_pre =0 ;
138
139 // Test sur np pour initialisation eventuelle
140 {
141 if (np > np_pre) {
142 np_pre = np ;
143 cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
144 for (int i=0 ; i<np+2 ; i += 2) {
145 cx[i] = - (i/2) ;
146 cx[i+1] = (i/2) ;
147 }
148 }
149 } // Fin de region critique
150
151 // pt. sur le tableau de double resultat
152 double* xo = new double[(tb->dim).taille] ;
153
154 // Initialisation a zero :
155 for (int i=0; i<(tb->dim).taille; i++) {
156 xo[i] = 0 ;
157 }
158
159 // On y va...
160 double* xi = tb->t ;
161 double* xci = xi ; // Pointeurs
162 double* xco = xo ; // courants
163 int k, j ;
164
165 // k = 0: inutile, deviendra k=1
166 xci += nr*nt ; // Positionnement
167 xco += nr*nt ;
168
169 // k = 1: 0
170 for (j=0 ; j<nr*nt ; j++) {
171 xco[j] = 0 ;
172 } // Fin de la boucle sur r * theta
173 xci += nr*nt ; // Positionnement
174 xco += nr*nt ;
175
176 // 2 =< k <= np-1
177 for (k=2 ; k<np ; k++) {
178 for (j=0 ; j<nr*nt ; j++) {
179 xco[j] = cx[k] * xci[j] ;
180 } // Fin de la boucle sur r * theta
181 xci += nr*nt ; // Positionnement
182 xco += nr*nt ;
183 } // Fin de la boucle sur phi
184
185 // k = np: inutile, deviendra k=np+1
186 xci += nr*nt ; // Positionnement
187 xco += nr*nt ;
188
189 // k = np+1
190 for (j=0 ; j<nr*nt ; j++) {
191 xco[j] = 0 ;
192 } // Fin de la boucle sur r * theta
193
194 // Inversion des sinus et cosinus (appel a dswap de BLAS)
195 int nbr = nr*nt ;
196 int inc = 1 ;
197 double* p1 = xo ;
198 double* p2 = p1 + nr*nt ;
199 for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
200 F77_dswap(&nbr, p1, &inc, p2, &inc) ;
201 }
202
203 // On remet les choses la ou il faut
204 delete [] tb->t ;
205 tb->t = xo ;
206
207 // base de developpement
208 // inchangee
209}
210
211 //------------------//
212 // cas P_COSSIN_P //
213 //------------------//
214
215void _dsdphi_p_cossin_p(Tbl* tb, int & )
216{
217
218 // Peut-etre rien a faire ?
219 if (tb->get_etat() == ETATZERO) {
220 return ;
221 }
222
223 // Protection
224 assert(tb->get_etat() == ETATQCQ) ;
225
226 // Pour le confort
227 int nr = (tb->dim).dim[0] ; // Nombre
228 int nt = (tb->dim).dim[1] ; // de points
229 int np = (tb->dim).dim[2] ; // physiques REELS
230 np = np - 2 ; // Nombre de points physiques
231
232 // Cas particulier de la symetrie axiale :
233 if (np == 1) {
234 tb->set_etat_zero() ;
235 return ;
236 }
237
238 // Variables statiques
239 static double* cx = 0 ;
240 static int np_pre =0 ;
241
242 // Test sur np pour initialisation eventuelle
243 {
244 if (np > np_pre) {
245 np_pre = np ;
246 cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
247 int i ;
248 for (i=0 ; i<np+2 ; i += 2) {
249 cx[i] = - (i/2) ;
250 cx[i+1] = (i/2) ;
251 }
252 for (i=0 ; i<np+2 ; i++) {
253 cx[i] = 2 * cx[i] ;
254 }
255 }
256 } // Fin de region critique
257
258 // pt. sur le tableau de double resultat
259 double* xo = new double[(tb->dim).taille] ;
260
261 // Initialisation a zero :
262 for (int i=0; i<(tb->dim).taille; i++) {
263 xo[i] = 0 ;
264 }
265
266 // On y va...
267 double* xi = tb->t ;
268 double* xci = xi ; // Pointeurs
269 double* xco = xo ; // courants
270 int k, j ;
271
272 // k = 0: inutile, deviendra k=1
273 xci += nr*nt ; // Positionnement
274 xco += nr*nt ;
275
276 // k = 1: 0
277 for (j=0 ; j<nr*nt ; j++) {
278 xco[j] = 0 ;
279 } // Fin de la boucle sur r * theta
280 xci += nr*nt ; // Positionnement
281 xco += nr*nt ;
282
283 // 2 =< k <= np-1
284 for (k=2 ; k<np ; k++) {
285 for (j=0 ; j<nr*nt ; j++) {
286 xco[j] = cx[k] * xci[j] ;
287 } // Fin de la boucle sur r * theta
288 xci += nr*nt ; // Positionnement
289 xco += nr*nt ;
290 } // Fin de la boucle sur phi
291
292 // k = np: inutile, deviendra k=np+1
293 xci += nr*nt ; // Positionnement
294 xco += nr*nt ;
295
296 // k = np+1
297 for (j=0 ; j<nr*nt ; j++) {
298 xco[j] = 0 ;
299 } // Fin de la boucle sur r * theta
300
301 // Inversion des sinus et cosinus (appel a dswap de BLAS)
302 int nbr = nr*nt ;
303 int inc = 1 ;
304 double* p1 = xo ;
305 double* p2 = p1 + nr*nt ;
306 for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
307 F77_dswap(&nbr, p1, &inc, p2, &inc) ;
308 }
309
310 // On remet les choses la ou il faut
311 delete [] tb->t ;
312 tb->t = xo ;
313
314 // base de developpement
315 // inchangee
316}
317
318 //------------------//
319 // cas P_COSSIN_I //
320 //------------------//
321
322void _dsdphi_p_cossin_i(Tbl* tb, int & )
323{
324
325 // Peut-etre rien a faire ?
326 if (tb->get_etat() == ETATZERO) {
327 return ;
328 }
329
330 // Protection
331 assert(tb->get_etat() == ETATQCQ) ;
332
333 // Pour le confort
334 int nr = (tb->dim).dim[0] ; // Nombre
335 int nt = (tb->dim).dim[1] ; // de points
336 int np = (tb->dim).dim[2] ; // physiques REELS
337 np = np - 2 ; // Nombre de points physiques
338
339 int ntnr = nt*nr ; // increment en phi
340
341 // Cas particulier de la symetrie axiale :
342 if (np == 1) {
343 tb->set_etat_zero() ;
344 return ;
345 }
346
347 // pt. sur le tableau de double resultat
348 double* xo = new double[(tb->dim).taille] ;
349
350 // pt. sur le tableau de double entree
351 const double* xi = tb->t ;
352
353 // k = 0 : resultat sur cos(phi)
354 // ------------------------------
355
356 double* xco = xo ; // coef de cos(phi)
357 const double* xci = xi + 2*ntnr ; // coef de sin(phi)
358 for (int i=0; i<ntnr; i++) {
359 xco[i] = xci[i] ;
360 }
361
362 // k = 1 : mise a zero
363 // --------------------
364
365 xco = xo + ntnr ;
366 for (int i=0; i<ntnr; i++) {
367 xco[i] = 0 ;
368 }
369
370 // k = 2 : resultat sur sin(phi)
371 // ------------------------------
372
373 xco = xo + 2*ntnr ; // coef de sin(phi)
374 xci = xi ; // coef de cos(phi)
375 for (int i=0; i<ntnr; i++) {
376 xco[i] = - xci[i] ;
377 }
378
379 // k >= 3
380 // ------
381
382 for (int k=3; k<np; k+=2) {
383
384 // resultat sur cos(k phi)
385 // -----------------------
386
387 xco = xo + k*ntnr ; // coef de cos(k phi)
388 xci = xi + (k+1)*ntnr ; // coef de sin(k phi)
389
390 for (int i=0; i<ntnr; i++) {
391 xco[i] = k * xci[i] ;
392 }
393
394 // resultat sur sin(k phi)
395 // -----------------------
396
397 xco = xo + (k+1)*ntnr ; // coef de sin(k phi)
398 xci = xi + k*ntnr ; // coef de cos(k phi)
399
400 for (int i=0; i<ntnr; i++) {
401 xco[i] = - k * xci[i] ;
402 }
403
404 }
405
406 // k = np+1 : mise a zero
407 // -----------------------
408
409 xco = xo + (np+1)*ntnr ;
410 for (int i=0; i<ntnr; i++) {
411 xco[i] = 0 ;
412 }
413
414 // On remet les choses la ou il faut
415 // ---------------------------------
416 delete [] tb->t ;
417 tb->t = xo ;
418
419 // base de developpement
420 // inchangee
421}
422}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene
Lorene prototypes.
Definition app_hor.h:67