LORENE
chb_legmp_cos.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_legmp_cos_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_legmp_cos.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $" ;
25 
26 /*
27  * Calcule les coefficients du developpement (suivant theta)
28  * en cos(j*theta)
29  * a partir des coefficients du developpement en fonctions
30  * associees de Legendre P_l^m(cos(theta)) (m pair)
31  * pour une une fonction 3-D symetrique par le retournement (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 a_l du develop. en fonctions de
42  * Legendre associees P_n^m:
43  *
44  * f(theta) = som_{l=m}^{nt-1} a_l P_l^m( cos(theta) )
45  *
46  * (m pair)
47  *
48  * ou P_n^m(x) represente la fonction de Legendre associee
49  * de degre n et d'ordre m normalisee de facon a ce que
50  *
51  * int_0^pi [ P_n^m(cos(theta)) ]^2 sin(theta) dtheta = 1
52  *
53  * L'espace memoire correspondant au pointeur cfi doit etre
54  * nr*nt*(np+2) et doit avoir ete alloue avant
55  * l'appel a la routine.
56  * Le coefficient a_l (0 <= l <= nt-1) doit etre stoke dans le
57  * tableau cfi comme suit
58  * a_l = cfi[ nr*nt* k + i + nr* l ]
59  * ou k et i sont les indices correspondant a phi et r
60  * respectivement: m = 2 (k/2).
61  * NB: pour l < m, a_l = 0
62  *
63  * Sortie:
64  * -------
65  * double* cfo : tableau des coefficients c_j du develop. en cos/sin definis
66  * comme suit (a r et phi fixes) :
67  *
68  * f(theta) = som_{j=0}^{nt-1} c_j cos( j theta )
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 c_j (0 <= j <= nt-1) est stoke dans le
74  * tableau cfo comme suit
75  * c_j = cfo[ nr*nt* k + i + nr* j ]
76  * ou k et i sont les indices correspondant a
77  * phi et r respectivement: m = 2 (k/2).
78  * Pour m impair, c_0 = c_{nt-1} = 0.
79  *
80  *
81  * NB:
82  * ---
83  * Il n'est pas possible d'avoir le pointeur cfo egal a cfi.
84  */
85 
86 /*
87  * $Id: chb_legmp_cos.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
88  * $Log: chb_legmp_cos.C,v $
89  * Revision 1.3 2014/10/13 08:53:11 j_novak
90  * Lorene classes and functions now belong to the namespace Lorene.
91  *
92  * Revision 1.2 2014/10/06 15:16:00 j_novak
93  * Modified #include directives to use c++ syntax.
94  *
95  * Revision 1.1 2009/10/13 13:49:36 j_novak
96  * New base T_LEG_MP.
97  *
98  *
99  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_legmp_cos.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
100  *
101  */
102 
103 
104 
105 // headers du C
106 #include <cstdlib>
107 #include <cassert>
108 
109 // Prototypage
110 #include "headcpp.h"
111 #include "proto.h"
112 
113 namespace Lorene {
114 //******************************************************************************
115 
116 void chb_legmp_cos(const int* deg , const double* cfi, double* cfo) {
117 
118 int k2, l, j, i, m ;
119 
120 // Nombres de degres de liberte en phi et theta :
121  int np = deg[0] ;
122  int nt = deg[1] ;
123  int nr = deg[2] ;
124 
125  assert(np < 4*nt) ;
126 
127  // Tableau de travail
128  double* som = new double[nr] ;
129 
130 // Recherche de la matrice de passage Legendre --> cos/sin
131  double* bb = mat_legmp_cos(np, nt) ;
132 
133 // Increment en m pour la matrice bb :
134  int mbb = nt * nt ;
135 
136 // Pointeurs de travail :
137  double* resu = cfo ;
138  const double* cc = cfi ;
139 
140 // Increment en phi :
141  int ntnr = nt * nr ;
142 
143 // Indice courant en phi :
144  int k = 0 ;
145 
146 //----------------------------------------------------------------
147 // Cas axisymetrique
148 //----------------------------------------------------------------
149 
150  if (np == 1) {
151 
152  m = 0 ;
153 
154 // Boucle sur l'indice j du developpement en cos(j theta)
155 
156  for (j=0; j<nt; j++) {
157 
158 // ... produit matriciel (parallelise sur r)
159  for (i=0; i<nr; i++) {
160  som[i] = 0 ;
161  }
162 
163  for (l=m; l<nt; l++) {
164 
165  double bmjl = bb[nt*j + l] ;
166  for (i=0; i<nr; i++) {
167  som[i] += bmjl * cc[nr*l + 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 j
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  // On sort
186  delete [] som ;
187  return ;
188 
189  } // fin du cas np=1
190 
191 
192 //----------------------------------------------------------------
193 // Cas 3-D
194 //----------------------------------------------------------------
195 
196 
197 // Boucle sur phi :
198 
199  for (m=0; m < np + 1 ; m+=2) {
200 
201  for (k2=0; k2 < 2; k2++) { // k2=0 : cos(m phi) ; k2=1 : sin(m phi)
202 
203  if ( (k == 1) || (k == np+1) ) { // On met les coef de sin(0 phi)
204  // et sin( np phi) a zero
205  for (j=0; j<nt; j++) {
206  for (i=0; i<nr; i++) {
207  *resu = 0 ;
208  resu++ ;
209  }
210  }
211  }
212  else {
213 
214 // Boucle sur l'indice j du developpement en cos(2 j theta)
215 
216  for (j=0; j<nt; j++) {
217 
218 // ... produit matriciel (parallelise sur r)
219  for (i=0; i<nr; i++) {
220  som[i] = 0 ;
221  }
222 
223  for (l=m; l<nt; l++) {
224 
225  double bmjl = bb[nt*j + l] ;
226  for (i=0; i<nr; i++) {
227  som[i] += bmjl * cc[nr*l + i] ;
228  }
229  }
230 
231  for (i=0; i<nr; i++) {
232  *resu = som[i] ;
233  resu++ ;
234  }
235 
236  } // fin de la boucle sur j
237 
238  } // fin du cas k != 1
239 
240 // On passe au phi suivant :
241  cc = cc + ntnr ;
242  k++ ;
243 
244  } // fin de la boucle sur k2
245 
246 // On passe a l'harmonique en phi suivante :
247 
248  bb += mbb ; // pointeur sur la nouvelle matrice de passage
249 
250  } // fin de la boucle (m) sur phi
251 
252 //## verif :
253  assert(resu == cfo + (np+2)*ntnr) ;
254 
255  // Menage
256  delete [] som ;
257 
258 }
259 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64