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
63namespace Lorene {
64//*****************************************************************************
65
66void 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
157void 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
255void 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