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
23char 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
60namespace Lorene {
61//*****************************************************************************
62
63void 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
154void 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
252void 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 prototypes.
Definition app_hor.h:64