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