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 
27 char 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 //----------------
79 namespace Lorene {
80 void _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 //------------
89 void _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 //-------------
185 void _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 //-------------
253 void _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 //-----------------
349 void _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 //-----------------
505 void _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 //-------------
659 void _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 //-------------
774 void _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 //------------
891 void _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_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 R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
#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