LORENE
op_mult_cp.C
1 /*
2  * Sets of routines for multiplication by cos(phi)
3  * (internal use)
4  *
5  * void _mult_cp_XXXX(Tbl * t, int & b)
6  * t pointer on the Tbl containing the spectral coefficients
7  * b spectral base
8  *
9  */
10 
11 /*
12  * Copyright (c) 2000-2001 Eric Gourgoulhon
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation; either version 2 of the License, or
19  * (at your option) any later version.
20  *
21  * LORENE is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with LORENE; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  *
30  */
31 
32 
33 char op_mult_cp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_cp.C,v 1.5 2014/10/13 08:53:25 j_novak Exp $" ;
34 
35 /*
36  * $Id: op_mult_cp.C,v 1.5 2014/10/13 08:53:25 j_novak Exp $
37  * $Log: op_mult_cp.C,v $
38  * Revision 1.5 2014/10/13 08:53:25 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.4 2007/12/14 10:19:33 jl_cornou
42  * *** empty log message ***
43  *
44  * Revision 1.3 2007/10/05 12:37:20 j_novak
45  * Corrected a few errors in the theta-nonsymmetric case (bases T_COSSIN_C and
46  * T_COSSIN_S).
47  *
48  * Revision 1.2 2004/11/23 15:16:01 m_forot
49  *
50  * Added the bases for the cases without any equatorial symmetry
51  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
52  *
53  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
54  * LORENE
55  *
56  * Revision 2.4 2000/11/14 15:11:45 eric
57  * Traitement du cas np=1
58  *
59  * Revision 2.3 2000/09/18 10:14:42 eric
60  * Ajout des bases P_COSSIN_P et P_COSSIN_I
61  *
62  * Revision 2.2 2000/09/12 13:36:11 eric
63  * Met les bonnes bases meme dans le cas ETATZERO
64  *
65  * Revision 2.1 2000/09/12 08:28:47 eric
66  * Traitement des bases qui dependent de la parite de m
67  *
68  * Revision 2.0 2000/09/11 13:53:32 eric
69  * *** empty log message ***
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_cp.C,v 1.5 2014/10/13 08:53:25 j_novak Exp $
73  *
74  */
75 
76 // Headers Lorene
77 #include "tbl.h"
78 
79  //-----------------------------------
80  // Routine for unknown cases
81  //-----------------------------------
82 
83 namespace Lorene {
84 void _mult_cp_pas_prevu(Tbl* , int& base) {
85  cout << "mult_cp() is not not implemented for the basis " << base << " !"
86  << endl ;
87  abort() ;
88 }
89 
90  //----------------
91  // basis P_COSSIN
92  //----------------
93 
94 void _mult_cp_p_cossin(Tbl* tb, int& base) {
95 
96  assert(tb->get_etat() != ETATNONDEF) ; // Protection
97 
98  // The spectral bases in r and theta which depend on the parity of m
99  // are changed
100  // -----------------------------------------------------------------
101 
102  int base_r = base & MSQ_R ;
103  int base_t = base & MSQ_T ;
104  const int base_p = base & MSQ_P ;
105 
106  switch (base_r) {
107 
108  case R_CHEBPIM_P : {
109  base_r = R_CHEBPIM_I ;
110  break ;
111  }
112 
113  case R_CHEBPIM_I : {
114  base_r = R_CHEBPIM_P ;
115  break ;
116  }
117  case R_CHEBPI_P : {
118  break ;
119  }
120 
121  case R_CHEBPI_I : {
122  break ;
123  }
124 
125 
126  case R_CHEB : {
127  break ;
128  }
129 
130  case R_JACO02 : {
131  break ;
132  }
133 
134  case R_CHEBU : {
135  break ;
136  }
137 
138  default : {
139  cout << "_mult_cp_p_cossin : unknown basis in r !" << endl ;
140  cout << " base_r = " << hex << base_r << endl ;
141  abort() ;
142  }
143  }
144 
145  switch (base_t) {
146 
147  case T_COSSIN_CP : {
148  base_t = T_COSSIN_SI ;
149  break ;
150  }
151 
152  case T_COSSIN_SI : {
153  base_t = T_COSSIN_CP ;
154  break ;
155  }
156 
157  case T_COSSIN_CI : {
158  base_t = T_COSSIN_SP ;
159  break ;
160  }
161 
162  case T_COSSIN_SP : {
163  base_t = T_COSSIN_CI ;
164  break ;
165  }
166 
167  case T_COSSIN_S : {
168  base_t = T_COSSIN_C ;
169  break ;
170  }
171  case T_COSSIN_C : {
172  base_t = T_COSSIN_S ;
173  break ;
174  }
175 
176  default : {
177  cout << "_mult_cp_p_cossin : unknown basis in theta !" << endl ;
178  cout << " base_t = " << hex << base_t << endl ;
179  abort() ;
180  }
181  }
182 
183  base = base_r | base_t | base_p ;
184 
185  //----------------------------------
186  // Start of computation
187  //----------------------------------
188 
189  // Nothing to do ?
190  if (tb->get_etat() == ETATZERO) {
191  return ;
192  }
193 
194  assert(tb->get_etat() == ETATQCQ) ; // Protection
195 
196  // Number of degrees of freedom
197  int nr = tb->get_dim(0) ;
198  int nt = tb->get_dim(1) ;
199  int np = tb->get_dim(2) - 2 ;
200 
201  // Special case np = 1 (axisymmetry) --> identity
202  // ---------------------------------
203 
204  if (np==1) {
205  return ;
206  }
207 
208  assert(np >= 4) ;
209 
210  int ntnr = nt*nr ;
211 
212  double* const cf = tb->t ; // input coefficients
213  double* const resu = new double[ tb->get_taille() ] ; // final result
214  double* co = resu ; // initial position
215 
216  // Case k=0 (m=0)
217  // --------------
218 
219  int q = 2 * ntnr ;
220  for (int i=0; i<ntnr; i++) {
221  co[i] = 0.5 * cf[q + i] ;
222  }
223  co += ntnr ;
224 
225  // Case k = 1
226  // ----------
227 
228  for (int i=0; i<ntnr; i++) {
229  co[i] = 0 ;
230  }
231  co += ntnr ;
232 
233  // Case k = 2
234  // ----------
235 
236  q = 4*ntnr ;
237  for (int i=0; i<ntnr; i++) {
238  co[i] = cf[i] + 0.5 * cf[q+i] ;
239  }
240  co += ntnr ;
241 
242  if (np==4) {
243 
244  // Case k = 3 for np=4
245  // ----------
246 
247  for (int i=0; i<ntnr; i++) {
248  co[i] = 0 ;
249  }
250  co += ntnr ;
251  }
252  else {
253 
254  // Case k = 3 for np>4
255  // ----------
256 
257  q = 5*ntnr ;
258  for (int i=0; i<ntnr; i++) {
259  co[i] = 0.5 * cf[q+i] ;
260  }
261  co += ntnr ;
262 
263  // Cases 4 <= k <= np-2
264  // --------------------
265 
266  for (int k=4; k<np-1; k++) {
267  int q1 = (k-2)*ntnr ;
268  int q2 = (k+2)*ntnr ;
269  for (int i=0; i<ntnr; i++) {
270  co[i] = 0.5 * ( cf[q1+i] + cf[q2+i] ) ;
271  }
272  co += ntnr ;
273  }
274 
275  // Case k = np - 1 for np>4
276  // ---------------
277 
278  q = (np-3)*ntnr ;
279  for (int i=0; i<ntnr; i++) {
280  co[i] = 0.5 * cf[q+i] ;
281  }
282  co += ntnr ;
283 
284  } // End of case np > 4
285 
286 
287  // Case k = np
288  // -----------
289 
290  q = (np-2)*ntnr ;
291  for (int i=0; i<ntnr; i++) {
292  co[i] = 0.5 * cf[q+i] ;
293  }
294  co += ntnr ;
295 
296  // Case k = np+1
297  // -------------
298 
299  for (int i=0; i<ntnr; i++) {
300  co[i] = 0 ;
301  }
302 
303  //## verif :
304  co += ntnr ;
305  assert( co == resu + (np+2)*ntnr ) ;
306 
307 
308  // The result is set to tb :
309  // -----------------------
310  delete [] tb->t ;
311  tb->t = resu ;
312 
313  return ;
314 }
315 
316 
317  //-----------------
318  // basis P_COSSIN_P
319  //-----------------
320 
321 void _mult_cp_p_cossin_p(Tbl* tb, int& base) {
322 
323  assert(tb->get_etat() != ETATNONDEF) ; // Protection
324 
325  // New spectral bases
326  // ------------------
327  int base_r = base & MSQ_R ;
328  int base_t = base & MSQ_T ;
329  base = base_r | base_t | P_COSSIN_I ;
330 
331  //----------------------------------
332  // Start of computation
333  //----------------------------------
334 
335  // Nothing to do ?
336  if (tb->get_etat() == ETATZERO) {
337  return ;
338  }
339 
340  assert(tb->get_etat() == ETATQCQ) ; // Protection
341 
342  // Number of degrees of freedom
343  int nr = tb->get_dim(0) ;
344  int nt = tb->get_dim(1) ;
345  int np = tb->get_dim(2) - 2 ;
346 
347  // Special case np = 1 (axisymmetry) --> identity
348  // ---------------------------------
349 
350  if (np==1) {
351  return ;
352  }
353 
354  assert(np >= 4) ;
355 
356  int ntnr = nt*nr ;
357 
358  double* const cf = tb->t ; // input coefficients
359  double* const resu = new double[ tb->get_taille() ] ; // final result
360  double* co = resu ; // initial position
361 
362  // Case k=0
363  // --------
364 
365  int q = 2 * ntnr ;
366  for (int i=0; i<ntnr; i++) {
367  co[i] = cf[i] + 0.5 * cf[q + i] ;
368  }
369  co += ntnr ;
370 
371  // Case k = 1
372  // ----------
373 
374  for (int i=0; i<ntnr; i++) {
375  co[i] = 0 ;
376  }
377  co += ntnr ;
378 
379  // Case k = 2
380  // ----------
381 
382  q = 3*ntnr ;
383  for (int i=0; i<ntnr; i++) {
384  co[i] = 0.5 * cf[q+i] ;
385  }
386  co += ntnr ;
387 
388 
389  // Cases 3 <= k <= np-1
390  // --------------------
391 
392  for (int k=3; k<np; k++) {
393  int q1 = (k-1)*ntnr ;
394  int q2 = (k+1)*ntnr ;
395  for (int i=0; i<ntnr; i++) {
396  co[i] = 0.5 * ( cf[q1+i] + cf[q2+i] ) ;
397  }
398  co += ntnr ;
399  }
400 
401  // Case k = np
402  // -----------
403 
404  q = (np-1)*ntnr ;
405  for (int i=0; i<ntnr; i++) {
406  co[i] = 0.5 * cf[q+i] ;
407  }
408  co += ntnr ;
409 
410  // Case k = np+1
411  // -------------
412 
413  for (int i=0; i<ntnr; i++) {
414  co[i] = 0 ;
415  }
416 
417  //## verif :
418  co += ntnr ;
419  assert( co == resu + (np+2)*ntnr ) ;
420 
421  // The result is set to tb :
422  // -----------------------
423  delete [] tb->t ;
424  tb->t = resu ;
425 
426  return ;
427 }
428 
429 
430 
431  //-----------------
432  // basis P_COSSIN_I
433  //-----------------
434 
435 void _mult_cp_p_cossin_i(Tbl* tb, int& base) {
436 
437  assert(tb->get_etat() != ETATNONDEF) ; // Protection
438 
439  // New spectral bases
440  // ------------------
441  int base_r = base & MSQ_R ;
442  int base_t = base & MSQ_T ;
443  base = base_r | base_t | P_COSSIN_P ;
444 
445  //----------------------------------
446  // Start of computation
447  //----------------------------------
448 
449  // Nothing to do ?
450  if (tb->get_etat() == ETATZERO) {
451  return ;
452  }
453 
454  assert(tb->get_etat() == ETATQCQ) ; // Protection
455 
456  // Number of degrees of freedom
457  int nr = tb->get_dim(0) ;
458  int nt = tb->get_dim(1) ;
459  int np = tb->get_dim(2) - 2 ;
460 
461  // Special case np = 1 (axisymmetry) --> identity
462  // ---------------------------------
463 
464  if (np==1) {
465  return ;
466  }
467 
468  assert(np >= 4) ;
469 
470  int ntnr = nt*nr ;
471 
472  double* const cf = tb->t ; // input coefficients
473  double* const resu = new double[ tb->get_taille() ] ; // final result
474  double* co = resu ; // initial position
475 
476  // Case k=0
477  // --------
478 
479  for (int i=0; i<ntnr; i++) {
480  co[i] = 0.5 * cf[i] ;
481  }
482  co += ntnr ;
483 
484  // Case k = 1
485  // ----------
486 
487  for (int i=0; i<ntnr; i++) {
488  co[i] = 0 ;
489  }
490  co += ntnr ;
491 
492  // Case k = 2
493  // ----------
494 
495  int q = 3*ntnr ;
496  for (int i=0; i<ntnr; i++) {
497  co[i] = 0.5 * ( cf[i] + cf[q+i] ) ;
498  }
499  co += ntnr ;
500 
501  // Cases 3 <= k <= np-1
502  // --------------------
503 
504  for (int k=3; k<np; k++) {
505  int q1 = (k-1)*ntnr ;
506  int q2 = (k+1)*ntnr ;
507  for (int i=0; i<ntnr; i++) {
508  co[i] = 0.5 * ( cf[q1+i] + cf[q2+i] ) ;
509  }
510  co += ntnr ;
511  }
512 
513  // Case k = np
514  // -----------
515 
516  q = (np-1)*ntnr ;
517  for (int i=0; i<ntnr; i++) {
518  co[i] = 0.5 * cf[q+i] ;
519  }
520  co += ntnr ;
521 
522  // Case k = np+1
523  // -------------
524 
525  for (int i=0; i<ntnr; i++) {
526  co[i] = 0 ;
527  }
528 
529  //## verif :
530  co += ntnr ;
531  assert( co == resu + (np+2)*ntnr ) ;
532 
533  // The result is set to tb :
534  // -----------------------
535  delete [] tb->t ;
536  tb->t = resu ;
537 
538  return ;
539 }
540 
541 
542 
543 
544 
545 
546 
547 
548 
549 
550 
551 
552 
553 
554 
555 
556 
557 
558 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition: type_parite.h:154
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
P_COSSIN_I
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180