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