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
33char 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
83namespace Lorene {
84void _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
94void _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
357void _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
492void _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 P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
#define MSQ_R
Extraction de l'info sur R.
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define MSQ_T
Extraction de l'info sur Theta.
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
#define R_CHEB
base de Chebychev ordinaire (fin)
#define MSQ_P
Extraction de l'info sur Phi.
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64