LORENE
chb_sin_legmi.C
1 /*
2  * Copyright (c) 1999-2001 Eric Gourgoulhon
3  * Copyright (c) 2009 Jerome Novak
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 char chb_sin_legmi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $" ;
25 
26 /*
27  * Calcule les coefficients du developpement (suivant theta) en fonctions
28  * associees de Legendre P_l^m(cos(theta)) a partir des coefficients du
29  * developpement en sin(j*theta)
30  * representant une fonction 3-D antisymetrique par le retournement
31  * (x, y, z) --> (-x, -y, z).
32  *
33  * Entree:
34  * -------
35  * const int* deg : tableau du nombre effectif de degres de liberte dans chacune
36  * des 3 dimensions:
37  * deg[0] = np : nombre de points de collocation en phi
38  * deg[1] = nt : nombre de points de collocation en theta
39  * deg[2] = nr : nombre de points de collocation en r
40  *
41  * const double* cfi : tableau des coefficients c_j du develop. en cos defini
42  * comme suit (a r et phi fixes)
43  *
44  * f(theta) = som_{j=0}^{nt-1} c_j sin( j theta )
45  *
46  * L'espace memoire correspondant au pointeur cfi doit etre
47  * nr*nt*(np+2) et doit avoir ete alloue avant
48  * l'appel a la routine.
49  * Le coefficient c_j (0 <= j <= nt-2) doit etre stoke dans le
50  * tableau cfi comme suit
51  * c_j = cfi[ nr*nt* k + i + nr* j ]
52  * ou k et i sont les indices correspondant a
53  * phi et r respectivement.
54  *
55  * Sortie:
56  * -------
57  * double* cfo : tableau des coefficients a_j du develop. en fonctions de
58  * Legendre associees P_l^m (m impair)
59  *
60  * f(theta) =
61  * som_{l=m}^{nt-1} a_j P_j^m( cos(theta) )
62  *
63  * avec m impair.
64  *
65  * P_l^m(x) represente la fonction de Legendre associee
66  * de degre l et d'ordre m normalisee de facon a ce que
67  *
68  * int_0^pi [ P_l^m(cos(theta)) ]^2 sin(theta) dtheta = 1
69  *
70  * L'espace memoire correspondant au pointeur cfo doit etre
71  * nr*nt*(np+2) et doit avoir ete alloue avant
72  * l'appel a la routine.
73  * Le coefficient a_j (0 <= j <= nt-1) est stoke dans le
74  * tableau cfo comme suit
75  * a_j = cfo[ nr*nt* k + i + nr* j ]
76  * ou k et i sont les indices correspondant a phi et r
77  * respectivement: m = 2( k/2 ).
78  * NB: pour j< m, a_j = 0
79  *
80  * NB:
81  * ---
82  * Il n'est pas possible d'avoir le pointeur cfo egal a cfi.
83  */
84 
85 /*
86  * $Id: chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
87  * $Log: chb_sin_legmi.C,v $
88  * Revision 1.3 2014/10/13 08:53:11 j_novak
89  * Lorene classes and functions now belong to the namespace Lorene.
90  *
91  * Revision 1.2 2014/10/06 15:16:01 j_novak
92  * Modified #include directives to use c++ syntax.
93  *
94  * Revision 1.1 2009/10/23 12:54:47 j_novak
95  * New base T_LEG_MI
96  *
97  *
98  *
99  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
100  *
101  */
102 
103 // headers du C
104 #include <cassert>
105 #include <cstdlib>
106 
107 // Prototypage
108 #include "headcpp.h"
109 #include "proto.h"
110 
111 namespace Lorene {
112 //******************************************************************************
113 
114 void chb_sin_legmi(const int* deg , const double* cfi, double* cfo) {
115 
116 int k2, l, jmin, j, i, m ;
117 
118  // Nombres de degres de liberte en phi et theta :
119  int np = deg[0] ;
120  int nt = deg[1] ;
121  int nr = deg[2] ;
122 
123  assert(np < 4*nt) ;
124  assert( cfi != cfo ) ;
125 
126  // Tableau de travail
127  double* som = new double[nr] ;
128 
129  // Recherche de la matrice de passage sin --> Legendre
130  double* aa = mat_sin_legmi(np, nt) ;
131 
132  // Increment en m pour la matrice aa :
133  int maa = nt * nt ;
134 
135  // Pointeurs de travail :
136  double* resu = cfo ;
137  const double* cc = cfi ;
138 
139  // Increment en phi :
140  int ntnr = nt * nr ;
141 
142  // Indice courant en phi :
143  int k = 0 ;
144 
145  // Cas k=0 (m=1 : cos(phi))
146  // ------------------------
147 
148  // Cas l=0 : a_l = 0
149  for (i=0; i<nr; i++) {
150  *resu = 0. ;
151  resu++ ;
152  }
153 
154  // ... produit matriciel
155  for (l=1; l<nt-1; l++) {
156  for (i=0; i<nr; i++) {
157  som[i] = 0 ;
158  }
159 
160  jmin = l ; // pour m=1, aa_lj = 0 pour j<l
161 
162  for (j=jmin; j<nt-1; j++) {
163  double amlj = aa[nt*l + j] ;
164  for (i=0; i<nr; i++) {
165  som[i] += amlj * cc[nr*j + i] ;
166  }
167  }
168 
169  for (i=0; i<nr; i++) {
170  *resu = som[i] ;
171  resu++ ;
172  }
173 
174  } // fin de la boucle sur l
175 
176  // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
177  for (i=0; i<nr; i++) {
178  *resu = 0 ;
179  resu++ ;
180  }
181 
182  // Special case np=1 (axisymmetry)
183  // -------------------------------
184  if (np==1) {
185  for (i=0; i<2*ntnr; i++) {
186  *resu = 0 ;
187  resu++ ;
188  }
189  delete [] som ;
190  return ;
191  }
192 
193 
194  // On passe au phi suivant :
195  cc = cc + ntnr ;
196  k++ ;
197 
198  // Cas k=1 : tout est mis a zero
199  // -----------------------------
200 
201  for (l=0; l<nt; l++) {
202  for (i=0; i<nr; i++) {
203  *resu = 0 ;
204  resu++ ;
205  }
206  }
207 
208  // On passe au phi suivant :
209  cc = cc + ntnr ;
210  k++ ;
211 
212  // Cas k=2 (m=1 : sin(phi))
213  // ------------------------
214 
215  // Cas l=0 : a_l = 0
216  for (i=0; i<nr; i++) {
217  *resu = 0. ;
218  resu++ ;
219  }
220 
221  // ... produit matriciel
222  for (l=1; l<nt-1; l++) {
223  for (i=0; i<nr; i++) {
224  som[i] = 0 ;
225  }
226 
227  jmin = l ; // pour m=1, aa_lj = 0 pour j<l
228 
229  for (j=jmin; j<nt-1; j++) {
230  double amlj = aa[nt*l + j] ;
231  for (i=0; i<nr; i++) {
232  som[i] += amlj * cc[nr*j + i] ;
233  }
234  }
235 
236  for (i=0; i<nr; i++) {
237  *resu = som[i] ;
238  resu++ ;
239  }
240 
241  } // fin de la boucle sur l
242 
243  // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
244  for (i=0; i<nr; i++) {
245  *resu = 0 ;
246  resu++ ;
247  }
248 
249  // On passe au phi suivant :
250  cc = cc + ntnr ;
251  k++ ;
252 
253  // On passe au m suivant
254  aa += maa ; // pointeur sur la nouvelle matrice de passage
255 
256  // Cas k >= 3
257  // ----------
258 
259  for (m=3; m < np ; m+=2) {
260 
261  for (k2=0; k2 < 2; k2++) { // k2=0 : cos(m phi) ; k2=1 : sin(m phi)
262  int lmax = (m<nt-1 ? m : nt -1) ;
263  for (l=0; l<lmax; l++) {
264  for (i=0; i<nr; i++) {
265  *resu = 0 ;
266  resu++ ;
267  }
268  }
269 
270  // ... produit matriciel
271  for (l=m; l<nt-1; l++) {
272  for (i=0; i<nr; i++) {
273  som[i] = 0 ;
274  }
275 
276  jmin = 1 ;
277 
278  for (j=jmin; j<nt-1; j++) {
279  double amlj = aa[nt*l + j] ;
280  for (i=0; i<nr; i++) {
281  som[i] += amlj * cc[nr*j + i] ;
282  }
283  }
284 
285  for (i=0; i<nr; i++) {
286  *resu = som[i] ;
287  resu++ ;
288  }
289 
290  } // fin de la boucle sur l
291 
292  // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
293  for (i=0; i<nr; i++) {
294  *resu = 0 ;
295  resu++ ;
296  }
297 
298  // On passe au phi suivant :
299  cc = cc + ntnr ;
300  k++ ;
301 
302  } // fin de la boucle sur k2
303 
304  // On passe a l'harmonique en phi suivante :
305 
306  aa += maa ; // pointeur sur la nouvelle matrice de passage
307 
308  } // fin de la boucle (m) sur phi
309 
310  // Cas k=np+1 : tout est mis a zero
311  // --------------------------------
312 
313  for (l=0; l<nt; l++) {
314  for (i=0; i<nr; i++) {
315  *resu = 0 ;
316  resu++ ;
317  }
318  }
319 
320 
321 //## verif :
322  assert(resu == cfo + (np+2)*ntnr) ;
323 
324  // Menage
325  delete [] som ;
326 
327 }
328 }
Lorene prototypes.
Definition: app_hor.h:64