LORENE
cirleg.C
1 /*
2  * Copyright (c) 2013 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 
24 
25 
26 /*
27  * Transformation de Legendre inverse sur le troisieme indice
28  * (indice correspondant a r) d'un tableau 3-D.
29  *
30  *
31  */
32 
33 /*
34  * $Id: cirleg.C,v 1.4 2016/12/05 16:18:01 j_novak Exp $
35  * $Log: cirleg.C,v $
36  * Revision 1.4 2016/12/05 16:18:01 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.3 2014/10/13 08:53:11 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.2 2013/06/13 14:17:47 j_novak
43  * Implementation of Legendre inverse coefficient transform.
44  *
45  * Revision 1.1 2013/06/06 15:31:33 j_novak
46  * Functions to compute Legendre coefficients (not fully tested yet).
47  *
48  *
49  * $Header $
50  *
51  */
52 
53 // headers du C
54 #include <cassert>
55 #include <cstdlib>
56 
57 //Lorene prototypes
58 #include "tbl.h"
59 #include "proto.h"
60 #include "utilitaires.h"
61 
62 
63 namespace Lorene {
64 //*****************************************************************************
65 
66 void cirleg(const int* deg, const int* dimc, double* cf, const int* dimf,
67  double* ff)
68 
69 {
70 
71 // Dimensions des tableaux ff et cf :
72  int n1f = dimf[0] ;
73  int n2f = dimf[1] ;
74  int n3f = dimf[2] ;
75  int n1c = dimc[0] ;
76  int n2c = dimc[1] ;
77  int n3c = dimc[2] ;
78 
79 // Nombres de degres de liberte en r :
80  int nr = deg[2] ;
81 
82 // Tests de dimension:
83  if (nr > n3c) {
84  cout << "cirleg: nr > n3c : nr = " << nr << " , n3c = "
85  << n3c << endl ;
86  abort () ;
87  exit(-1) ;
88  }
89  if (nr > n3f) {
90  cout << "cirleg: nr > n3f : nr = " << nr << " , n3f = "
91  << n3f << endl ;
92  abort () ;
93  exit(-1) ;
94  }
95  if (n1c > n1f) {
96  cout << "cirleg: n1c > n1f : n1c = " << n1c << " , n1f = "
97  << n1f << endl ;
98  abort () ;
99  exit(-1) ;
100  }
101  if (n2c > n2f) {
102  cout << "cirleg: n2c > n2f : n2c = " << n2c << " , n2f = "
103  << n2f << endl ;
104  abort () ;
105  exit(-1) ;
106  }
107 
108 // boucle sur phi et theta
109 
110  int n2n3f = n2f * n3f ;
111  int n2n3c = n2c * n3c ;
112 
113 /*
114  * Borne de la boucle sur phi:
115  * si n1c = 1, on effectue la boucle une fois seulement.
116  * si n1c > 1, on va jusqu'a j = n1c-2 en sautant j = 1 (les coefficients
117  * j=n1c-1 et j=0 ne sont pas consideres car nuls).
118  */
119  int borne_phi = ( n1c > 1 ) ? n1c-1 : 1 ;
120 
121  double* colloc = new double[nr] ;
122  legendre_collocation_points(nr, colloc) ;
123 
124  for (int j=0; j< borne_phi; j++) {
125 
126  if (j==1) continue ; // on ne traite pas le terme en sin(0 phi)
127 
128  for (int k=0; k<n2c; k++) {
129 
130  int i0 = n2n3c * j + n3c * k ; // indice de depart
131  double* cf0 = cf + i0 ; // tableau des donnees a transformer
132 
133  i0 = n2n3f * j + n3f * k ; // indice de depart
134  double* ff0 = ff + i0 ; // tableau resultat
135 
136  for (int i = 0; i<nr; i++) {
137  double x0 = colloc[i] ;
138  double Pi = 1. ;
139  double Pip1 = x0 ;
140  double som = cf0[0] + cf0[1]*x0 ;
141  for (int h=2; h<nr; h++) {
142  double Pip2 = (2. - 1./double(h))*x0*Pip1
143  - (1. - 1./double(h))*Pi ;
144  som += cf0[h]*Pip2 ;
145  Pi = Pip1 ;
146  Pip1 = Pip2 ;
147  }
148  ff0[i] = som ;
149  }
150  } // fin de la boucle sur theta
151  } // fin de la boucle sur phi
152  delete [] colloc ;
153 }
154 
155 //*****************************************************************************
156 
157 void cirlegp(const int* deg, const int* dimc, double* cf,
158  const int* dimf, double* ff)
159 
160 {
161 
162 // Dimensions des tableaux ff et cf :
163  int n1f = dimf[0] ;
164  int n2f = dimf[1] ;
165  int n3f = dimf[2] ;
166  int n1c = dimc[0] ;
167  int n2c = dimc[1] ;
168  int n3c = dimc[2] ;
169 
170 // Nombres de degres de liberte en r :
171  int nr = deg[2] ;
172 
173 // Tests de dimension:
174  if (nr > n3c) {
175  cout << "cirlegp: nr > n3c : nr = " << nr << " , n3c = "
176  << n3c << endl ;
177  abort () ;
178  exit(-1) ;
179  }
180  if (nr > n3f) {
181  cout << "cirlegp: nr > n3f : nr = " << nr << " , n3f = "
182  << n3f << endl ;
183  abort () ;
184  exit(-1) ;
185  }
186  if (n1c > n1f) {
187  cout << "cirlegp: n1c > n1f : n1c = " << n1c << " , n1f = "
188  << n1f << endl ;
189  abort () ;
190  exit(-1) ;
191  }
192  if (n2c > n2f) {
193  cout << "cirlegp: n2c > n2f : n2c = " << n2c << " , n2f = "
194  << n2f << endl ;
195  abort () ;
196  exit(-1) ;
197  }
198 
199 // Nombre de points
200  int nm1 = nr - 1;
201  int dnm1 = 2*nr - 1 ;
202 
203 // boucle sur phi et theta
204 
205  int n2n3f = n2f * n3f ;
206  int n2n3c = n2c * n3c ;
207 
208 /*
209  * Borne de la boucle sur phi:
210  * si n1c = 1, on effectue la boucle une fois seulement.
211  * si n1c > 1, on va jusqu'a j = n1c-2 en sautant j = 1 (les coefficients
212  * j=n1c-1 et j=0 ne sont pas consideres car nuls).
213  */
214  int borne_phi = ( n1c > 1 ) ? n1c-1 : 1 ;
215 
216  double* colloc = new double[dnm1] ;
217  legendre_collocation_points(dnm1, colloc) ;
218 
219  for (int j=0; j< borne_phi; j++) {
220 
221  if (j==1) continue ; // on ne traite pas le terme en sin(0 phi)
222 
223  for (int k=0; k<n2c; k++) {
224 
225  int i0 = n2n3c * j + n3c * k ; // indice de depart
226  double* cf0 = cf + i0 ; // tableau des donnees a transformer
227 
228  i0 = n2n3f * j + n3f * k ; // indice de depart
229  double* ff0 = ff + i0 ; // tableau resultat
230 
231  for (int i = 0; i<nr; i++) {
232  double x0 = colloc[nm1+i] ;
233  double Pi = 1. ;
234  double Pip1 = x0 ;
235  double som = cf0[0] ;
236  for (int h=2; h<dnm1; h++) {
237  double Pip2 = (2. - 1./double(h))*x0*Pip1
238  - (1. - 1./double(h))*Pi ;
239  if (h%2 == 0) som += cf0[h/2]*Pip2 ;
240  Pi = Pip1 ;
241  Pip1 = Pip2 ;
242  }
243  ff0[i] = som ;
244  }
245 
246  } // fin de la boucle sur theta
247  } // fin de la boucle sur phi
248 
249  delete [] colloc ;
250 
251 }
252 
253 //*****************************************************************************
254 
255 void cirlegi(const int* deg, const int* dimc, double* cf,
256  const int* dimf, double* ff)
257 
258 {
259 
260 // Dimensions des tableaux ff et cf :
261  int n1f = dimf[0] ;
262  int n2f = dimf[1] ;
263  int n3f = dimf[2] ;
264  int n1c = dimc[0] ;
265  int n2c = dimc[1] ;
266  int n3c = dimc[2] ;
267 
268 // Nombres de degres de liberte en r :
269  int nr = deg[2] ;
270 
271 // Tests de dimension:
272  if (nr > n3c) {
273  cout << "cirlegi: nr > n3c : nr = " << nr << " , n3c = "
274  << n3c << endl ;
275  abort () ;
276  exit(-1) ;
277  }
278  if (nr > n3f) {
279  cout << "cirlegi: nr > n3f : nr = " << nr << " , n3f = "
280  << n3f << endl ;
281  abort () ;
282  exit(-1) ;
283  }
284  if (n1c > n1f) {
285  cout << "cirlegi: n1c > n1f : n1c = " << n1c << " , n1f = "
286  << n1f << endl ;
287  abort () ;
288  exit(-1) ;
289  }
290  if (n2c > n2f) {
291  cout << "cirlegi: n2c > n2f : n2c = " << n2c << " , n2f = "
292  << n2f << endl ;
293  abort () ;
294  exit(-1) ;
295  }
296 
297 // Nombre de points
298  int nm1 = nr - 1;
299  int dnm1 = 2*nr - 1 ;
300 
301 // boucle sur phi et theta
302 
303  int n2n3f = n2f * n3f ;
304  int n2n3c = n2c * n3c ;
305 
306 /*
307  * Borne de la boucle sur phi:
308  * si n1c = 1, on effectue la boucle une fois seulement.
309  * si n1c > 1, on va jusqu'a j = n1c-2 en sautant j = 1 (les coefficients
310  * j=n1c-1 et j=0 ne sont pas consideres car nuls).
311  */
312  int borne_phi = ( n1c > 1 ) ? n1c-1 : 1 ;
313 
314  double* colloc = new double[dnm1] ;
315  legendre_collocation_points(dnm1, colloc) ;
316 
317  for (int j=0; j< borne_phi; j++) {
318 
319  if (j==1) continue ; // on ne traite pas le terme en sin(0 phi)
320 
321  for (int k=0; k<n2c; k++) {
322 
323  int i0 = n2n3c * j + n3c * k ; // indice de depart
324  double* cf0 = cf + i0 ; // tableau des donnees a transformer
325 
326  i0 = n2n3f * j + n3f * k ; // indice de depart
327  double* ff0 = ff + i0 ; // tableau resultat
328 
329  for (int i = 0; i<nr; i++) {
330  double x0 = colloc[nm1+i] ;
331  double Pi = 1. ;
332  double Pip1 = x0 ;
333  double som = cf0[0]*x0 ;
334  for (int h=2; h<dnm1; h++) {
335  double Pip2 = (2. - 1./double(h))*x0*Pip1
336  - (1. - 1./double(h))*Pi ;
337  if (h%2 == 1) som += cf0[h/2]*Pip2 ;
338  Pi = Pip1 ;
339  Pip1 = Pip2 ;
340  }
341  ff0[i] = som ;
342  }
343 
344  } // fin de la boucle sur theta
345  } // fin de la boucle sur phi
346 
347  delete [] colloc ;
348 
349 }
350 
351 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67