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