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