LORENE
chb_cos_legmp.C
1 /*
2  * Copyright (c) 1999-2001 Eric Gourgoulhon
3  * 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_cos_legmp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_cos_legmp.C,v 1.4 2014/10/13 08:53:10 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 cos(j*theta)
30  * representant une fonction 3-D symetrique 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 cos( 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-1) 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_l du develop. en fonctions de
58  * Legendre associees P_n^m:
59  *
60  * f(theta) =
61  * som_{l=m}^{nt-1} a_l P_l^m( cos(theta) )
62  *
63  * avec m pair : m = 0, 2, ..., np.
64  *
65  * P_n^m(x) represente la fonction de Legendre associee
66  * de degre n et d'ordre m normalisee de facon a ce que
67  *
68  * int_0^pi [ P_n^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_l (0 <= l <= nt-1) est stoke dans le
74  * tableau cfo comme suit
75  * a_l = cfo[ nr*nt* k + i + nr* l ]
76  * ou k et i sont les indices correspondant a phi et r
77  * respectivement: m = 2( k/2 ).
78  * NB: pour l < m, a_l = 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_cos_legmp.C,v 1.4 2014/10/13 08:53:10 j_novak Exp $
87  * $Log: chb_cos_legmp.C,v $
88  * Revision 1.4 2014/10/13 08:53:10 j_novak
89  * Lorene classes and functions now belong to the namespace Lorene.
90  *
91  * Revision 1.3 2014/10/06 15:15:59 j_novak
92  * Modified #include directives to use c++ syntax.
93  *
94  * Revision 1.2 2009/10/23 12:54:47 j_novak
95  * New base T_LEG_MI
96  *
97  * Revision 1.1 2009/10/13 13:49:36 j_novak
98  * New base T_LEG_MP.
99  *
100  *
101  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_cos_legmp.C,v 1.4 2014/10/13 08:53:10 j_novak Exp $
102  *
103  */
104 
105 
106 // headers du C
107 #include <cassert>
108 #include <cstdlib>
109 
110 // Prototypage
111 #include "headcpp.h"
112 #include "proto.h"
113 
114 namespace Lorene {
115 //******************************************************************************
116 
117 void chb_cos_legmp(const int* deg , const double* cfi, double* cfo) {
118 
119 int k2, l, jmin, j, i, m ;
120 
121 // Nombres de degres de liberte en phi et theta :
122  int np = deg[0] ;
123  int nt = deg[1] ;
124  int nr = deg[2] ;
125 
126  assert(np < 4*nt) ;
127 
128  // Tableau de travail
129  double* som = new double[nr] ;
130 
131 // Recherche de la matrice de passage cos --> Legendre
132  double* aa = mat_cos_legmp(np, nt) ;
133 
134 // Increment en m pour la matrice aa :
135  int maa = nt * nt ;
136 
137 // Pointeurs de travail :
138  double* resu = cfo ;
139  const double* cc = cfi ;
140 
141 // Increment en phi :
142  int ntnr = nt * nr ;
143 
144 // Indice courant en phi :
145  int k = 0 ;
146 
147 //----------------------------------------------------------------
148 // Cas axisymetrique
149 //----------------------------------------------------------------
150 
151  if (np == 1) {
152 
153  m = 0 ;
154 
155 // Boucle sur l'indice l du developpement en Legendre
156 
157 // ... produit matriciel (parallelise sur r)
158  for (l=m; l<nt; l++) {
159  for (i=0; i<nr; i++) {
160  som[i] = 0 ;
161  }
162 
163  jmin = l ; // pour m=0, aa_lj = 0 pour j<l
164  for (j=jmin; j<nt; j++) {
165  double amlj = aa[nt*l + j] ;
166  for (i=0; i<nr; i++) {
167  som[i] += amlj * cc[nr*j + i] ;
168  }
169  }
170 
171  for (i=0; i<nr; i++) {
172  *resu = som[i] ;
173  resu++ ;
174  }
175 
176  } // fin de la boucle sur l
177 
178  // Mise a zero des coefficients k=1 et k=2 :
179  // ---------------------------------------
180 
181  for (i=ntnr; i<3*ntnr; i++) {
182  cfo[i] = 0 ;
183  }
184 
185 
186  // on sort
187  delete [] som ;
188  return ;
189 
190  } // fin du cas np=1
191 
192 
193 //----------------------------------------------------------------
194 // Cas 3-D
195 //----------------------------------------------------------------
196 
197 
198 // Boucle sur phi :
199 
200 
201  for (m=0; m < np + 1 ; m+=2) {
202 
203  for (k2=0; k2 < 2; k2++) { // k2=0 : cos(m phi) ; k2=1 : sin(m phi)
204 
205  if ( (k == 1) || (k == np+1) ) { // On met les coef de sin(0 phi)
206  // et sin( np phi) a zero
207  for (l=0; l<nt; l++) {
208  for (i=0; i<nr; i++) {
209  *resu = 0 ;
210  resu++ ;
211  }
212  }
213  }
214  else {
215 
216 // Boucle sur l'indice l du developpement en Legendre
217 
218  int lmax = (m<nt-1 ? m : nt-1) ;
219  for (l=0; l<lmax; l++) {
220  for (i=0; i<nr; i++) {
221  *resu = 0 ;
222  resu++ ;
223  }
224  }
225 // ... produit matriciel (parallelise sur r)
226  for (l=m; l<nt; l++) {
227  for (i=0; i<nr; i++) {
228  som[i] = 0 ;
229  }
230 
231  jmin = ( m == 0 ) ? l : 0 ; // pour m=0, aa_lj = 0 pour j<l
232  for (j=jmin; j<nt; j++) {
233  double amlj = aa[nt*l + j] ;
234  for (i=0; i<nr; i++) {
235  som[i] += amlj * cc[nr*j + i] ;
236  }
237  }
238 
239  for (i=0; i<nr; i++) {
240  *resu = som[i] ;
241  resu++ ;
242  }
243 
244  } // fin de la boucle sur l
245 
246  } // fin du cas k != 1 et k!=np+1
247 
248 // On passe au phi suivant :
249  cc = cc + ntnr ;
250  k++ ;
251 
252  } // fin de la boucle sur k2
253 
254 // On passe a l'harmonique en phi suivante :
255 
256  aa += maa ; // pointeur sur la nouvelle matrice de passage
257 
258  } // fin de la boucle (m) sur phi
259 
260 //## verif :
261  assert(resu == cfo + (np+2)*ntnr) ;
262 
263  // Menage
264  delete [] som ;
265 
266 }
267 }
Lorene prototypes.
Definition: app_hor.h:64