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
33char 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
83namespace Lorene {
84void _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
94void _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
321void _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
435void _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}
#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