LORENE
op_primr.C
1/*
2 * Computation of primitive in a single domain
3 *
4 *
5 */
6
7/*
8 * Copyright (c) 2004 Eric Gourgoulhon.
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27char op_primr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_primr.C,v 1.10 2014/10/13 08:53:26 j_novak Exp $" ;
28
29/*
30 * $Id: op_primr.C,v 1.10 2014/10/13 08:53:26 j_novak Exp $
31 * $Log: op_primr.C,v $
32 * Revision 1.10 2014/10/13 08:53:26 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.9 2014/10/06 15:16:06 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.8 2013/04/25 15:46:06 j_novak
39 * Added special treatment in the case np = 1, for type_p = NONSYM.
40 *
41 * Revision 1.7 2007/12/21 13:59:02 j_novak
42 * Suppression of call to pow(-1, something).
43 *
44 * Revision 1.6 2007/12/11 15:28:18 jl_cornou
45 * Jacobi(0,2) polynomials partially implemented
46 *
47 * Revision 1.5 2006/05/17 15:01:16 j_novak
48 * Treatment of the case nr = 1 and R_CHEB
49 *
50 * Revision 1.4 2004/11/23 15:16:01 m_forot
51 *
52 * Added the bases for the cases without any equatorial symmetry
53 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
54 *
55 * Revision 1.3 2004/10/12 09:58:24 j_novak
56 * Better memory management.
57 *
58 * Revision 1.2 2004/06/14 15:24:57 e_gourgoulhon
59 * First operationnal version (tested).
60 *
61 * Revision 1.1 2004/06/13 21:33:13 e_gourgoulhon
62 * First version.
63 *
64 *
65 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_primr.C,v 1.10 2014/10/13 08:53:26 j_novak Exp $
66 *
67 */
68
69
70// C headers
71#include <cstdlib>
72#include <cmath>
73
74// Lorene headers
75#include "tbl.h"
76
77// Unexpected case
78//----------------
79namespace Lorene {
80void _primr_pas_prevu(const Tbl&, int bin, const Tbl&, Tbl&, int&, Tbl& ) {
81
82 cout << "Unexpected basis in primr : basis = " << hex << bin << endl ;
83 abort() ;
84
85}
86
87// case R_CHEB
88//------------
89void _primr_r_cheb(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
90 int& bout, Tbl& valp1) {
91
92 assert(tin.dim == tout.dim) ;
93
94 // Output spectral basis
95 bout = bin ;
96
97 // Number of coefficients
98 int nr = tin.get_dim(0) ;
99 int nt = tin.get_dim(1) ;
100 int np = tin.get_dim(2) - 2 ;
101 int borne_phi = np + 1 ;
102 if (np == 1) borne_phi = 1 ;
103
104 // Case of a zero input or pure angular grid
105 // -----------------------------------------
106 if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
107 if (valm1.get_etat() == ETATZERO) {
108 tout.set_etat_zero() ;
109 valp1.set_etat_zero() ;
110 return ;
111 }
112 else {
113 assert(valm1.get_etat() == ETATQCQ) ;
114 tout.set_etat_qcq() ;
115 valp1.set_etat_qcq() ;
116 double* xco = tout.t ;
117 for (int k=0 ; k< borne_phi ; k++) {
118 if (k==1) { // jump over the coefficient of sin(0*phi)
119 xco += nr*nt ;
120 }
121 else {
122 for (int j=0 ; j<nt ; j++) {
123 xco[0] = valm1(k,j) ; // constant value = boundary value
124 for (int i=1; i<nr; i++) xco[i] = 0 ;
125 valp1.set(k,j) = xco[0] ;
126 xco += nr ;
127 }
128 }
129 }
130 return ;
131 }
132 }
133
134 // Case of a non-zero input
135 // ------------------------
136
137 assert(tin.get_etat() == ETATQCQ ) ;
138 tout.annule_hard() ;
139 valp1.annule_hard() ;
140
141 const double* xci = tin.t ;
142 double* xco = tout.t ;
143
144 for (int k=0 ; k< borne_phi ; k++) {
145 if (k==1) { // jump over the coefficient of sin(0*phi)
146 xci += nr*nt ;
147 xco += nr*nt ;
148 }
149 else {
150 for (int j=0 ; j<nt ; j++) {
151
152 xco[1] = xci[0] - 0.5 * xci[2] ; // special case i = 1
153
154 for (int i=2; i<nr-2; i++) {
155 xco[i] = (xci[i-1] - xci[i+1]) / double(2*i) ;
156 }
157
158 xco[nr-2] = xci[nr-3] / double(2*nr - 4) ;
159 xco[nr-1] = xci[nr-2] / double(2*nr - 2) ;
160
161 // Determination of the T_0 coefficient by matching with
162 // provided value at xi = - 1 :
163 double som = - xco[1] ;
164 for (int i=2; i<nr; i+=2) som += xco[i] ;
165 for (int i=3; i<nr; i+=2) som -= xco[i] ;
166 xco[0] = valm1(k,j) - som ;
167
168 // Value of primitive at xi = + 1 :
169 som = xco[0] ;
170 for (int i=1; i<nr; i++) som += xco[i] ;
171 valp1.set(k,j) = som ;
172
173 xci += nr ;
174 xco += nr ;
175 } // end of theta loop
176 }
177 } // end of phi loop
178
179}
180
181
182
183// case R_CHEBP
184//-------------
185void _primr_r_chebp(const Tbl& tin, int bin, const Tbl&, Tbl& tout,
186 int& bout, Tbl& valp1) {
187
188 assert(tin.dim == tout.dim) ;
189
190 // Output spectral basis
191 int base_t = bin & MSQ_T ;
192 int base_p = bin & MSQ_P ;
193 bout = base_p | base_t | R_CHEBI ;
194
195 // Number of coefficients
196 int nr = tin.get_dim(0) ;
197 int nt = tin.get_dim(1) ;
198 int np = tin.get_dim(2) - 2 ;
199 int borne_phi = np + 1 ;
200 if (np == 1) borne_phi = 1 ;
201
202 // Case of a zero input
203 // --------------------
204 if (tin.get_etat() == ETATZERO) {
205 tout.set_etat_zero() ;
206 valp1.set_etat_zero() ;
207 return ;
208 }
209
210 // Case of a non-zero input
211 // ------------------------
212
213 assert(tin.get_etat() == ETATQCQ ) ;
214 tout.annule_hard() ;
215 valp1.annule_hard() ;
216
217 const double* xci = tin.t ;
218 double* xco = tout.t ;
219
220 for (int k=0 ; k< borne_phi ; k++) {
221 if (k==1) { // jump over the coefficient of sin(0*phi)
222 xci += nr*nt ;
223 xco += nr*nt ;
224 }
225 else {
226 for (int j=0 ; j<nt ; j++) {
227
228 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
229
230 for (int i=1; i<nr-2; i++) {
231 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
232 }
233
234 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
235 xco[nr-1] = 0 ;
236
237 // Value of primitive at xi = + 1 :
238 double som = xco[0] ;
239 for (int i=1; i<nr; i++) som += xco[i] ;
240 valp1.set(k,j) = som ;
241
242 xci += nr ;
243 xco += nr ;
244 } // end of theta loop
245 }
246 } // end of phi loop
247
248}
249
250
251// case R_CHEBI
252//-------------
253void _primr_r_chebi(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
254 int& bout, Tbl& valp1) {
255
256 assert(tin.dim == tout.dim) ;
257
258 // Output spectral basis
259 int base_t = bin & MSQ_T ;
260 int base_p = bin & MSQ_P ;
261 bout = base_p | base_t | R_CHEBP ;
262
263 // Number of coefficients
264 int nr = tin.get_dim(0) ;
265 int nt = tin.get_dim(1) ;
266 int np = tin.get_dim(2) - 2 ;
267 int borne_phi = np + 1 ;
268 if (np == 1) borne_phi = 1 ;
269
270
271 // Case of a zero input
272 // --------------------
273 if (tin.get_etat() == ETATZERO) {
274 if (val0.get_etat() == ETATZERO) {
275 tout.set_etat_zero() ;
276 valp1.set_etat_zero() ;
277 return ;
278 }
279 else {
280 assert(val0.get_etat() == ETATQCQ) ;
281 tout.annule_hard() ;
282 valp1.annule_hard() ;
283 double* xco = tout.t ;
284 for (int k=0 ; k< borne_phi ; k++) {
285 if (k==1) { // jump over the coefficient of sin(0*phi)
286 xco += nr*nt ;
287 }
288 else {
289 for (int j=0 ; j<nt ; j++) {
290 xco[0] = val0(k,j) ; // constant value = boundary value
291 for (int i=1; i<nr; i++) xco[i] = 0 ;
292 valp1.set(k,j) = xco[0] ;
293 xco += nr ;
294 }
295 }
296 }
297 return ;
298 }
299 }
300
301 // Case of a non-zero input
302 // ------------------------
303
304 assert(tin.get_etat() == ETATQCQ ) ;
305 tout.annule_hard() ;
306 valp1.annule_hard() ;
307
308 const double* xci = tin.t ;
309 double* xco = tout.t ;
310
311 for (int k=0 ; k< borne_phi ; k++) {
312 if (k==1) { // jump over the coefficient of sin(0*phi)
313 xci += nr*nt ;
314 xco += nr*nt ;
315 }
316 else {
317 for (int j=0 ; j<nt ; j++) {
318
319 for (int i=1; i<nr-1; i++) {
320 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
321 }
322
323 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
324
325 // Determination of the T_0 coefficient by maching with
326 // provided value at xi = 0 :
327 double som = - xco[1] ;
328 for (int i=2; i<nr; i+=2) som += xco[i] ;
329 for (int i=3; i<nr; i+=2) som -= xco[i] ;
330 xco[0] = val0(k,j) - som ;
331
332 // Value of primitive at xi = + 1 :
333 som = xco[0] ;
334 for (int i=1; i<nr; i++) som += xco[i] ;
335 valp1.set(k,j) = som ;
336
337 xci += nr ;
338 xco += nr ;
339 } // end of theta loop
340 }
341 } // end of phi loop
342
343}
344
345
346
347// case R_CHEBPIM_P
348//-----------------
349void _primr_r_chebpim_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
350 int& bout, Tbl& valp1) {
351
352 assert(tin.dim == tout.dim) ;
353
354 // Output spectral basis
355 int base_t = bin & MSQ_T ;
356 int base_p = bin & MSQ_P ;
357 bout = base_p | base_t | R_CHEBPIM_I ;
358
359 // Number of coefficients
360 int nr = tin.get_dim(0) ;
361 int nt = tin.get_dim(1) ;
362 int np = tin.get_dim(2) - 2 ;
363 int borne_phi = np + 1 ;
364 if (np == 1) borne_phi = 1 ;
365
366
367 // Case of a zero input
368 // --------------------
369 if (tin.get_etat() == ETATZERO) {
370 if (val0.get_etat() == ETATZERO) {
371 tout.set_etat_zero() ;
372 valp1.set_etat_zero() ;
373 return ;
374 }
375 else {
376 assert(val0.get_etat() == ETATQCQ) ;
377 tout.annule_hard() ;
378 valp1.annule_hard() ;
379 double* xco = tout.t ;
380
381 // m even part
382 for (int k=0 ; k<borne_phi ; k += 4) {
383 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
384 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
385 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
386 xco += nr*nt ;
387 }
388 else {
389 for (int j=0 ; j<nt ; j++) {
390 assert( val0(k+kmod,j) == double(0) ) ;
391 for (int i=0; i<nr; i++) xco[i] = 0 ;
392 valp1.set(k+kmod,j) = 0. ;
393 xco += nr ;
394 }
395 }
396 }
397 xco += 2*nr*nt ; // next even m
398 }
399
400 // m odd part
401 xco = tout.t + 2*nr*nt ;
402 for (int k=2 ; k<borne_phi ; k += 4) {
403 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
404 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
405 for (int j=0 ; j<nt ; j++) {
406 xco[0] = val0(k+kmod,j) ; // constant value = boundary value
407 for (int i=1; i<nr; i++) xco[i] = 0 ;
408 valp1.set(k+kmod,j) = xco[0] ;
409 xco += nr ;
410 }
411 }
412 xco += 2*nr*nt ; // next odd m
413 }
414 return ;
415 }
416 }
417
418 // Case of a non-zero input
419 // ------------------------
420
421 assert(tin.get_etat() == ETATQCQ ) ;
422 tout.annule_hard() ;
423 valp1.annule_hard() ;
424
425 const double* xci = tin.t ;
426 double* xco = tout.t ;
427
428 // m even part
429 // -----------
430 for (int k=0 ; k<borne_phi ; k += 4) {
431 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
432 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
433 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
434 xci += nr*nt ;
435 xco += nr*nt ;
436 }
437 else {
438 for (int j=0 ; j<nt ; j++) {
439 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
440
441 for (int i=1; i<nr-2; i++) {
442 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
443 }
444
445 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
446 xco[nr-1] = 0 ;
447
448 // Value of primitive at xi = + 1 :
449 double som = xco[0] ;
450 for (int i=1; i<nr; i++) som += xco[i] ;
451 valp1.set(k+kmod,j) = som ;
452
453 xci += nr ;
454 xco += nr ;
455 } // end of theta loop
456
457 }
458 }
459 xci += 2*nr*nt ; // next even m
460 xco += 2*nr*nt ; //
461 }
462
463 // m odd part
464 // ----------
465 xci = tin.t + 2*nr*nt ;
466 xco = tout.t + 2*nr*nt ;
467 for (int k=2 ; k<borne_phi ; k += 4) {
468 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
469 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
470 for (int j=0 ; j<nt ; j++) {
471
472 for (int i=1; i<nr-1; i++) {
473 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
474 }
475
476 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
477
478 // Determination of the T_0 coefficient by maching with
479 // provided value at xi = 0 :
480 double som = - xco[1] ;
481 for (int i=2; i<nr; i+=2) som += xco[i] ;
482 for (int i=3; i<nr; i+=2) som -= xco[i] ;
483 xco[0] = val0(k+kmod,j) - som ;
484
485 // Value of primitive at xi = + 1 :
486 som = xco[0] ;
487 for (int i=1; i<nr; i++) som += xco[i] ;
488 valp1.set(k+kmod,j) = som ;
489
490 xci += nr ;
491 xco += nr ;
492 } // end of theta loop
493 }
494 xci += 2*nr*nt ; // next odd m
495 xco += 2*nr*nt ; //
496 }
497
498
499}
500
501
502
503// case R_CHEBPIM_I
504//-----------------
505void _primr_r_chebpim_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
506 int& bout, Tbl& valp1) {
507
508 assert(tin.dim == tout.dim) ;
509
510 // Output spectral basis
511 int base_t = bin & MSQ_T ;
512 int base_p = bin & MSQ_P ;
513 bout = base_p | base_t | R_CHEBPIM_P ;
514
515 // Number of coefficients
516 int nr = tin.get_dim(0) ;
517 int nt = tin.get_dim(1) ;
518 int np = tin.get_dim(2) - 2 ;
519 int borne_phi = np + 1 ;
520 if (np == 1) borne_phi = 1 ;
521
522 // Case of a zero input
523 // --------------------
524 if (tin.get_etat() == ETATZERO) {
525 if (val0.get_etat() == ETATZERO) {
526 tout.set_etat_zero() ;
527 valp1.set_etat_zero() ;
528 return ;
529 }
530 else {
531 assert(val0.get_etat() == ETATQCQ) ;
532 tout.annule_hard() ;
533 valp1.annule_hard() ;
534 double* xco = tout.t ;
535
536 // m odd part
537 for (int k=0 ; k<borne_phi ; k += 4) {
538 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
539 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
540 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
541 xco += nr*nt ;
542 }
543 else {
544 for (int j=0 ; j<nt ; j++) {
545 xco[0] = val0(k+kmod,j) ; // constant value = boundary value
546 for (int i=1; i<nr; i++) xco[i] = 0 ;
547 valp1.set(k+kmod,j) = xco[0] ;
548 xco += nr ;
549 }
550 }
551 }
552 xco += 2*nr*nt ; // next odd m
553 }
554
555 // m even part
556 xco = tout.t + 2*nr*nt ;
557 for (int k=2 ; k<borne_phi ; k += 4) {
558 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
559 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
560 for (int j=0 ; j<nt ; j++) {
561 assert( val0(k+kmod,j) == double(0) ) ;
562 for (int i=0; i<nr; i++) xco[i] = 0 ;
563 valp1.set(k+kmod,j) = 0. ;
564 xco += nr ;
565 }
566 }
567 xco += 2*nr*nt ; // next even m
568 }
569 return ;
570 }
571 }
572
573 // Case of a non-zero input
574 // ------------------------
575
576 assert(tin.get_etat() == ETATQCQ ) ;
577 tout.annule_hard() ;
578 valp1.annule_hard() ;
579
580 const double* xci = tin.t ;
581 double* xco = tout.t ;
582
583 // m odd part
584 // ----------
585 for (int k=0 ; k<borne_phi ; k += 4) {
586 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
587 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
588 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
589 xci += nr*nt ;
590 xco += nr*nt ;
591 }
592 else {
593 for (int j=0 ; j<nt ; j++) {
594
595 for (int i=1; i<nr-1; i++) {
596 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
597 }
598
599 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
600
601 // Determination of the T_0 coefficient by maching with
602 // provided value at xi = 0 :
603 double som = - xco[1] ;
604 for (int i=2; i<nr; i+=2) som += xco[i] ;
605 for (int i=3; i<nr; i+=2) som -= xco[i] ;
606 xco[0] = val0(k+kmod,j) - som ;
607
608 // Value of primitive at xi = + 1 :
609 som = xco[0] ;
610 for (int i=1; i<nr; i++) som += xco[i] ;
611 valp1.set(k+kmod,j) = som ;
612
613 xci += nr ;
614 xco += nr ;
615 } // end of theta loop
616 }
617 }
618 xci += 2*nr*nt ; // next odd m
619 xco += 2*nr*nt ; //
620 }
621
622 // m even part
623 // -----------
624 xci = tin.t + 2*nr*nt ;
625 xco = tout.t + 2*nr*nt ;
626 for (int k=2 ; k<borne_phi ; k += 4) {
627 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
628 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
629
630 for (int j=0 ; j<nt ; j++) {
631 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
632
633 for (int i=1; i<nr-2; i++) {
634 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
635 }
636
637 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
638 xco[nr-1] = 0 ;
639
640 // Value of primitive at xi = + 1 :
641 double som = xco[0] ;
642 for (int i=1; i<nr; i++) som += xco[i] ;
643 valp1.set(k+kmod,j) = som ;
644
645 xci += nr ;
646 xco += nr ;
647 } // end of theta loop
648
649 }
650 xci += 2*nr*nt ; // next even m
651 xco += 2*nr*nt ; //
652 }
653
654
655}
656
657// case R_CHEBPI_P
658//-------------
659void _primr_r_chebpi_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
660 int& bout, Tbl& valp1) {
661
662 assert(tin.dim == tout.dim) ;
663
664 // Output spectral basis
665 int base_t = bin & MSQ_T ;
666 int base_p = bin & MSQ_P ;
667 bout = base_p | base_t | R_CHEBPI_I ;
668
669 // Number of coefficients
670 int nr = tin.get_dim(0) ;
671 int nt = tin.get_dim(1) ;
672 int np = tin.get_dim(2) - 2 ;
673 int borne_phi = np + 1 ;
674 if (np == 1) borne_phi = 1 ;
675
676
677 // Case of a zero input
678 // --------------------
679 if (tin.get_etat() == ETATZERO) {
680 if (val0.get_etat() == ETATZERO) {
681 tout.set_etat_zero() ;
682 valp1.set_etat_zero() ;
683 return ;
684 }
685 else {
686 assert(val0.get_etat() == ETATQCQ) ;
687 tout.annule_hard() ;
688 valp1.annule_hard() ;
689 double* xco = tout.t ;
690 for (int k=0 ; k< borne_phi ; k++) {
691 if (k==1) { // jump over the coefficient of sin(0*phi)
692 xco += nr*nt ;
693 }
694 else {
695 for (int j=0 ; j<nt ; j++) {
696 int l = j%2;
697 if(l==0){
698 for (int i=0; i<nr; i++) xco[i] = 0 ;
699 valp1.set(k,j) = 0. ;
700 } else {
701 xco[0] = val0(k,j) ; // constant value = boundary value
702 for (int i=1; i<nr; i++) xco[i] = 0 ;
703 valp1.set(k,j) = xco[0] ;
704 }
705 xco += nr ;
706 }
707 }
708 }
709 return ;
710 }
711 }
712
713 // Case of a non-zero input
714 // ------------------------
715
716 assert(tin.get_etat() == ETATQCQ ) ;
717 tout.annule_hard() ;
718 valp1.annule_hard() ;
719
720 const double* xci = tin.t ;
721 double* xco = tout.t ;
722
723 for (int k=0 ; k< borne_phi ; k++) {
724 if (k==1) { // jump over the coefficient of sin(0*phi)
725 xci += nr*nt ;
726 xco += nr*nt ;
727 }
728 else {
729 for (int j=0 ; j<nt ; j++) {
730
731 int l = j%2;
732 if(l==0){
733 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
734
735 for (int i=1; i<nr-2; i++) {
736 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
737 }
738
739 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
740 xco[nr-1] = 0 ;
741
742 // Value of primitive at xi = + 1 :
743 double som = xco[0] ;
744 for (int i=1; i<nr; i++) som += xco[i] ;
745 valp1.set(k,j) = som ;
746 } else {
747 for (int i=1; i<nr-1; i++) {
748 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
749 }
750
751 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
752
753 // Determination of the T_0 coefficient by maching with
754 // provided value at xi = 0 :
755 double som = - xco[1] ;
756 for (int i=2; i<nr; i+=2) som += xco[i] ;
757 for (int i=3; i<nr; i+=2) som -= xco[i] ;
758 xco[0] = val0(k,j) - som ;
759
760 // Value of primitive at xi = + 1 :
761 som = xco[0] ;
762 for (int i=1; i<nr; i++) som += xco[i] ;
763 valp1.set(k,j) = som ;
764 }
765 xci += nr ;
766 xco += nr ;
767 } // end of theta loop
768 }
769 } // end of phi loop
770
771}
772// case R_CHEBPI_I
773//-------------
774void _primr_r_chebpi_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
775 int& bout, Tbl& valp1) {
776
777 assert(tin.dim == tout.dim) ;
778
779 // Output spectral basis
780 int base_t = bin & MSQ_T ;
781 int base_p = bin & MSQ_P ;
782 bout = base_p | base_t | R_CHEBPI_P ;
783
784 // Number of coefficients
785 int nr = tin.get_dim(0) ;
786 int nt = tin.get_dim(1) ;
787 int np = tin.get_dim(2) - 2 ;
788 int borne_phi = np + 1 ;
789 if (np == 1) borne_phi = 1 ;
790
791
792 // Case of a zero input
793 // --------------------
794 if (tin.get_etat() == ETATZERO) {
795 if (val0.get_etat() == ETATZERO) {
796 tout.set_etat_zero() ;
797 valp1.set_etat_zero() ;
798 return ;
799 }
800 else {
801 assert(val0.get_etat() == ETATQCQ) ;
802 tout.annule_hard() ;
803 valp1.annule_hard() ;
804 double* xco = tout.t ;
805 for (int k=0 ; k< borne_phi ; k++) {
806 if (k==1) { // jump over the coefficient of sin(0*phi)
807 xco += nr*nt ;
808 }
809 else {
810 for (int j=0 ; j<nt ; j++) {
811 int l = j%2;
812 if(l==1){
813 for (int i=0; i<nr; i++) xco[i] = 0 ;
814 valp1.set(k,j) = 0. ;
815 } else {
816 xco[0] = val0(k,j) ; // constant value = boundary value
817 for (int i=1; i<nr; i++) xco[i] = 0 ;
818 valp1.set(k,j) = xco[0] ;
819 }
820 xco += nr ;
821 }
822 }
823 }
824 return ;
825 }
826 }
827
828 // Case of a non-zero input
829 // ------------------------
830
831 assert(tin.get_etat() == ETATQCQ ) ;
832 tout.annule_hard() ;
833 valp1.annule_hard() ;
834
835 const double* xci = tin.t ;
836 double* xco = tout.t ;
837
838 for (int k=0 ; k< borne_phi ; k++) {
839 if (k==1) { // jump over the coefficient of sin(0*phi)
840 xci += nr*nt ;
841 xco += nr*nt ;
842 }
843 else {
844 for (int j=0 ; j<nt ; j++) {
845
846 int l = j%2;
847 if(l==1){
848 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
849
850 for (int i=1; i<nr-2; i++) {
851 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
852 }
853
854 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
855 xco[nr-1] = 0 ;
856
857 // Value of primitive at xi = + 1 :
858 double som = xco[0] ;
859 for (int i=1; i<nr; i++) som += xco[i] ;
860 valp1.set(k,j) = som ;
861 } else {
862 for (int i=1; i<nr-1; i++) {
863 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
864 }
865
866 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
867
868 // Determination of the T_0 coefficient by maching with
869 // provided value at xi = 0 :
870 double som = - xco[1] ;
871 for (int i=2; i<nr; i+=2) som += xco[i] ;
872 for (int i=3; i<nr; i+=2) som -= xco[i] ;
873 xco[0] = val0(k,j) - som ;
874
875 // Value of primitive at xi = + 1 :
876 som = xco[0] ;
877 for (int i=1; i<nr; i++) som += xco[i] ;
878 valp1.set(k,j) = som ;
879 }
880 xci += nr ;
881 xco += nr ;
882 } // end of theta loop
883 }
884 } // end of phi loop
885
886}
887
888
889// case R_JACO02
890//------------
891void _primr_r_jaco02(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
892 int& bout, Tbl& valp1) {
893
894 assert(tin.dim == tout.dim) ;
895
896 // Output spectral basis
897 bout = bin ;
898
899 // Number of coefficients
900 int nr = tin.get_dim(0) ;
901 int nt = tin.get_dim(1) ;
902 int np = tin.get_dim(2) - 2 ;
903 int borne_phi = np + 1 ;
904 if (np == 1) borne_phi = 1 ;
905
906 // Case of a zero input or pure angular grid
907 // -----------------------------------------
908 if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
909 if (valm1.get_etat() == ETATZERO) {
910 tout.set_etat_zero() ;
911 valp1.set_etat_zero() ;
912 return ;
913 }
914 else {
915 assert(valm1.get_etat() == ETATQCQ) ;
916 tout.set_etat_qcq() ;
917 valp1.set_etat_qcq() ;
918 double* xco = tout.t ;
919 for (int k=0 ; k< borne_phi ; k++) {
920 if (k==1) { // jump over the coefficient of sin(0*phi)
921 xco += nr*nt ;
922 }
923 else {
924 for (int j=0 ; j<nt ; j++) {
925 xco[0] = valm1(k,j) ; // constant value = boundary value
926 for (int i=1; i<nr; i++) xco[i] = 0 ;
927 valp1.set(k,j) = xco[0] ;
928 xco += nr ;
929 }
930 }
931 }
932 return ;
933 }
934 }
935
936 // Case of a non-zero input
937 // ------------------------
938
939 assert(tin.get_etat() == ETATQCQ ) ;
940 tout.annule_hard() ;
941 valp1.annule_hard() ;
942
943 const double* xci = tin.t ;
944 double* xco = tout.t ;
945
946 for (int k=0 ; k< borne_phi ; k++) {
947 if (k==1) { // jump over the coefficient of sin(0*phi)
948 xci += nr*nt ;
949 xco += nr*nt ;
950 }
951 else {
952 for (int j=0 ; j<nt ; j++) {
953
954 for (int i=1; i<nr-1; i++) {
955 xco[i] = (i+2)/double((i+1)*(2*i+1))*xci[i-1] - xci[i]/double((i+1)*(i+2)) - (i+1)/double((i+2)*(2*i+5))*xci[i+1] ;
956 }
957 xco[nr-1] = (nr+1)/double((nr)*(2*nr-1))*xci[nr-2] - xci[nr-1]/double((nr)*(nr+1));
958
959 // Determination of the J_0 coefficient by matching with
960 // provided value at xi = - 1 :
961
962 double som = -3*xco[1] ;
963 for (int i=2; i<nr; i++) {
964 int signe = (i%2 == 0 ? 1 : -1) ;
965 som += xco[i]*signe*(i+1)*(i+2)/double(2) ;
966 }
967 xco[0] = valm1(k,j) - som ;
968
969 // Value of primitive at xi = + 1 :
970 som = xco[0] ;
971 for (int i=1; i<nr; i++) som += xco[i] ;
972 valp1.set(k,j) = som ;
973
974 xci += nr ;
975 xco += nr ;
976 } // end of theta loop
977 }
978 } // end of phi loop
979
980}
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001}
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define MSQ_T
Extraction de l'info sur Theta.
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define MSQ_P
Extraction de l'info sur Phi.
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64