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 //--------------------------------
101 namespace Lorene {
102 void _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 
112 void _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 
215 void _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 
322 void _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
Lorene prototypes.
Definition: app_hor.h:67