LORENE
op_mult_sp.C
1 /*
2  * Sets of routines for multiplication by sin(phi)
3  * (internal use)
4  *
5  * void _mult_sp_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_sp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_sp.C,v 1.5 2014/10/13 08:53:25 j_novak Exp $" ;
34 
35 /*
36  * $Id: op_mult_sp.C,v 1.5 2014/10/13 08:53:25 j_novak Exp $
37  * $Log: op_mult_sp.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:12:01 eric
57  * Traitement du cas np=1
58  *
59  * Revision 2.3 2000/09/18 10:14:59 eric
60  * Ajout des bases P_COSSIN_P et P_COSSIN_I
61  *
62  * Revision 2.2 2000/09/12 13:36:38 eric
63  * Met les bonnes bases meme dans le cas ETATZERO
64  *
65  * Revision 2.1 2000/09/12 08:29:11 eric
66  * Traitement des bases qui dependent de la parite de m
67  *
68  * Revision 2.0 2000/09/11 13:53:42 eric
69  * *** empty log message ***
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_sp.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_sp_pas_prevu(Tbl* , int& base) {
85  cout << "mult_sp() is not not implemented for the basis " << base << " !"
86  << endl ;
87  abort() ;
88 }
89 
90  //----------------
91  // basis P_COSSIN
92  //----------------
93 
94 void _mult_sp_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 
118  case R_CHEBPI_P : {
119  break ;
120  }
121 
122  case R_CHEBPI_I : {
123  break ;
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 
172  case T_COSSIN_C : {
173  base_t = T_COSSIN_S ;
174  break ;
175  }
176 
177  default : {
178  cout << "_mult_cp_p_cossin : unknown basis in theta !" << endl ;
179  cout << " base_t = " << hex << base_t << endl ;
180  abort() ;
181  }
182  }
183 
184  base = base_r | base_t | base_p ;
185 
186 
187  //----------------------------------
188  // Start of computation
189  //----------------------------------
190 
191  // Nothing to do ?
192  if (tb->get_etat() == ETATZERO) {
193  return ;
194  }
195 
196  assert(tb->get_etat() == ETATQCQ) ; // Protection
197 
198  // Number of degrees of freedom
199  int nr = tb->get_dim(0) ;
200  int nt = tb->get_dim(1) ;
201  int np = tb->get_dim(2) - 2 ;
202 
203 
204  // Special case np = 1 (axisymmetry) --> zero is returned
205  // ---------------------------------
206 
207  if (np==1) {
208  tb->set_etat_zero() ;
209  return ;
210  }
211 
212 
213  assert(np >= 4) ;
214 
215  int ntnr = nt*nr ;
216 
217  double* const cf = tb->t ; // input coefficients
218  double* const resu = new double[ tb->get_taille() ] ; // final result
219  double* co = resu ; // initial position
220 
221  // Case k=0 (m=0)
222  // --------------
223 
224  int q = 3 * ntnr ;
225  for (int i=0; i<ntnr; i++) {
226  co[i] = 0.5 * cf[q + i] ;
227  }
228  co += ntnr ;
229 
230  // Case k = 1
231  // ----------
232 
233  for (int i=0; i<ntnr; i++) {
234  co[i] = 0 ;
235  }
236  co += ntnr ;
237 
238  if (np==4) {
239 
240  // Case k = 2 for np=4
241  // ----------
242 
243  for (int i=0; i<ntnr; i++) {
244  co[i] = 0 ;
245  }
246  co += ntnr ;
247 
248  // Case k = 3 for np=4
249  // ----------
250 
251  q = 4*ntnr ;
252  for (int i=0; i<ntnr; i++) {
253  co[i] = cf[i] - 0.5 * cf[q+i] ;
254  }
255  co += ntnr ;
256 
257  }
258  else{
259 
260  // Case k = 2 for np>4
261  // ----------
262 
263  q = 5*ntnr ;
264  for (int i=0; i<ntnr; i++) {
265  co[i] = 0.5 * cf[q+i] ;
266  }
267  co += ntnr ;
268 
269  // Case k = 3 for np>4
270  // ----------
271 
272  q = 4*ntnr ;
273  for (int i=0; i<ntnr; i++) {
274  co[i] = cf[i] - 0.5 * cf[q+i] ;
275  }
276  co += ntnr ;
277 
278 
279  // Cases 4 <= k <= np-3
280  // --------------------
281 
282  for (int k=4; k<np-2; k+=2) {
283 
284  int k1 = k ; // k1 even
285 
286  int q1 = (k1-1)*ntnr ;
287  int q2 = (k1+3)*ntnr ;
288  for (int i=0; i<ntnr; i++) {
289  co[i] = 0.5 * ( cf[q2+i] - cf[q1+i] ) ;
290  }
291  co += ntnr ;
292  k1++ ; // k1 odd
293 
294  q1 = (k1-3)*ntnr ;
295  q2 = (k1+1)*ntnr ;
296  for (int i=0; i<ntnr; i++) {
297  co[i] = 0.5 * ( cf[q1+i] - cf[q2+i] ) ;
298  }
299  co += ntnr ;
300  }
301 
302  // Case k = np - 2 for np>4
303  // ---------------
304 
305  q = (np-3)*ntnr ;
306  for (int i=0; i<ntnr; i++) {
307  co[i] = - 0.5 * cf[q + i] ;
308  }
309  co += ntnr ;
310 
311  // Case k = np - 1 for np>4
312  // ---------------
313 
314  int q1 = (np-4)*ntnr ;
315  int q2 = np*ntnr ;
316  for (int i=0; i<ntnr; i++) {
317  co[i] = 0.5 * ( cf[q1+i] - cf[q2+i] ) ;
318  }
319  co += ntnr ;
320 
321  } // End of case np > 4
322 
323 
324  // Case k = np
325  // -----------
326 
327  q = (np-1)*ntnr ;
328  for (int i=0; i<ntnr; i++) {
329  co[i] = - 0.5 * cf[q+i] ;
330  }
331  co += ntnr ;
332 
333  // Case k = np+1
334  // -------------
335 
336  for (int i=0; i<ntnr; i++) {
337  co[i] = 0 ;
338  }
339 
340  //## verif :
341  co += ntnr ;
342  assert( co == resu + (np+2)*ntnr ) ;
343 
344 
345  // The result is set to tb :
346  // -----------------------
347  delete [] tb->t ;
348  tb->t = resu ;
349 
350  return ;
351 }
352 
353  //------------------
354  // basis P_COSSIN_P
355  //------------------
356 
357 void _mult_sp_p_cossin_p(Tbl* tb, int& base) {
358 
359  assert(tb->get_etat() != ETATNONDEF) ; // Protection
360 
361  // New spectral bases
362  // ------------------
363  int base_r = base & MSQ_R ;
364  int base_t = base & MSQ_T ;
365  base = base_r | base_t | P_COSSIN_I ;
366 
367  //----------------------------------
368  // Start of computation
369  //----------------------------------
370 
371  // Nothing to do ?
372  if (tb->get_etat() == ETATZERO) {
373  return ;
374  }
375 
376  assert(tb->get_etat() == ETATQCQ) ; // Protection
377 
378  // Number of degrees of freedom
379  int nr = tb->get_dim(0) ;
380  int nt = tb->get_dim(1) ;
381  int np = tb->get_dim(2) - 2 ;
382 
383  // Special case np = 1 (axisymmetry) --> zero is returned
384  // ---------------------------------
385 
386  if (np==1) {
387  tb->set_etat_zero() ;
388  return ;
389  }
390 
391  assert(np >= 4) ;
392 
393  int ntnr = nt*nr ;
394 
395  double* const cf = tb->t ; // input coefficients
396  double* const resu = new double[ tb->get_taille() ] ; // final result
397  double* co = resu ; // initial position
398 
399  // Case k=0
400  // --------
401 
402  int q = 3 * ntnr ;
403  for (int i=0; i<ntnr; i++) {
404  co[i] = 0.5 * cf[q + i] ;
405  }
406  co += ntnr ;
407 
408  // Case k = 1
409  // ----------
410 
411  for (int i=0; i<ntnr; i++) {
412  co[i] = 0 ;
413  }
414  co += ntnr ;
415 
416  // Case k = 2
417  // ----------
418 
419  q = 2*ntnr ;
420  for (int i=0; i<ntnr; i++) {
421  co[i] = cf[i] - 0.5 * cf[q+i] ;
422  }
423  co += ntnr ;
424 
425  // Cases 3 <= k <= np-2
426  // --------------------
427 
428  for (int k=3; k<np-1; k+=2) {
429 
430  int k1 = k ; // k1 odd
431 
432  int q1 = k1*ntnr ;
433  int q2 = (k1+2)*ntnr ;
434  for (int i=0; i<ntnr; i++) {
435  co[i] = 0.5 * ( cf[q2+i] - cf[q1+i] ) ;
436  }
437  co += ntnr ;
438  k1++ ; // k1 even
439 
440  q1 = (k1-2)*ntnr ;
441  q2 = k1*ntnr ;
442  for (int i=0; i<ntnr; i++) {
443  co[i] = 0.5 * ( cf[q1+i] - cf[q2+i] ) ;
444  }
445  co += ntnr ;
446  }
447 
448  // Case k = np - 1
449  // ---------------
450 
451  q = (np-1)*ntnr ;
452  for (int i=0; i<ntnr; i++) {
453  co[i] = - 0.5 * cf[q+i] ;
454  }
455  co += ntnr ;
456 
457  // Case k = np
458  // -----------
459 
460  int q1 = (np-2)*ntnr ;
461  int q2 = np*ntnr ;
462  for (int i=0; i<ntnr; i++) {
463  co[i] = 0.5 * ( cf[q1+i] - cf[q2+i] ) ;
464  }
465  co += ntnr ;
466 
467  // Case k = np+1
468  // -------------
469 
470  for (int i=0; i<ntnr; i++) {
471  co[i] = 0 ;
472  }
473 
474  //## verif :
475  co += ntnr ;
476  assert( co == resu + (np+2)*ntnr ) ;
477 
478 
479  // The result is set to tb :
480  // -----------------------
481  delete [] tb->t ;
482  tb->t = resu ;
483 
484  return ;
485 }
486 
487 
488  //------------------
489  // basis P_COSSIN_I
490  //------------------
491 
492 void _mult_sp_p_cossin_i(Tbl* tb, int& base) {
493 
494  assert(tb->get_etat() != ETATNONDEF) ; // Protection
495 
496  // New spectral bases
497  // ------------------
498  int base_r = base & MSQ_R ;
499  int base_t = base & MSQ_T ;
500  base = base_r | base_t | P_COSSIN_P ;
501 
502  //----------------------------------
503  // Start of computation
504  //----------------------------------
505 
506  // Nothing to do ?
507  if (tb->get_etat() == ETATZERO) {
508  return ;
509  }
510 
511  assert(tb->get_etat() == ETATQCQ) ; // Protection
512 
513  // Number of degrees of freedom
514  int nr = tb->get_dim(0) ;
515  int nt = tb->get_dim(1) ;
516  int np = tb->get_dim(2) - 2 ;
517 
518  // Special case np = 1 (axisymmetry) --> zero is returned
519  // ---------------------------------
520 
521  if (np==1) {
522  tb->set_etat_zero() ;
523  return ;
524  }
525 
526  assert(np >= 4) ;
527 
528  int ntnr = nt*nr ;
529 
530  double* const cf = tb->t ; // input coefficients
531  double* const resu = new double[ tb->get_taille() ] ; // final result
532  double* co = resu ; // initial position
533 
534  // Case k=0
535  // --------
536 
537  int q = 2 * ntnr ;
538  for (int i=0; i<ntnr; i++) {
539  co[i] = 0.5 * cf[q + i] ;
540  }
541  co += ntnr ;
542 
543  // Case k = 1
544  // ----------
545 
546  for (int i=0; i<ntnr; i++) {
547  co[i] = 0 ;
548  }
549  co += ntnr ;
550 
551  // Case k = 2
552  // ----------
553 
554  int q1 = 2*ntnr ;
555  int q2 = 4*ntnr ;
556  for (int i=0; i<ntnr; i++) {
557  co[i] = 0.5 * ( cf[q2+i] - cf[q1+i] ) ;
558  }
559  co += ntnr ;
560 
561  // Case k = 3
562  // ----------
563 
564  q = 3*ntnr ;
565  for (int i=0; i<ntnr; i++) {
566  co[i] = 0.5 * ( cf[i] - cf[q+i] ) ;
567  }
568  co += ntnr ;
569 
570  // Cases 4 <= k <= np-1
571  // --------------------
572 
573  for (int k=4; k<np; k+=2) {
574 
575  int k1 = k ; // k1 even
576 
577  q1 = k1*ntnr ;
578  q2 = (k1+2)*ntnr ;
579  for (int i=0; i<ntnr; i++) {
580  co[i] = 0.5 * ( cf[q2+i] - cf[q1+i] ) ;
581  }
582  co += ntnr ;
583  k1++ ; // k1 odd
584 
585  q1 = (k1-2)*ntnr ;
586  q2 = k1*ntnr ;
587  for (int i=0; i<ntnr; i++) {
588  co[i] = 0.5 * ( cf[q1+i] - cf[q2+i] ) ;
589  }
590  co += ntnr ;
591  }
592 
593  // Case k = np
594  // -----------
595 
596  q = np*ntnr ;
597  for (int i=0; i<ntnr; i++) {
598  co[i] = - 0.5 * cf[q+i] ;
599  }
600  co += ntnr ;
601 
602  // Case k = np+1
603  // -------------
604 
605  for (int i=0; i<ntnr; i++) {
606  co[i] = 0 ;
607  }
608 
609  //## verif :
610  co += ntnr ;
611  assert( co == resu + (np+2)*ntnr ) ;
612 
613 
614  // The result is set to tb :
615  // -----------------------
616  delete [] tb->t ;
617  tb->t = resu ;
618 
619  return ;
620 }
621 
622 }
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene prototypes.
Definition: app_hor.h:64
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166