LORENE
gval_from_spectral.C
1 /*
2  * Functions for spectral summation to a Valencia-type grid (see grille_val.h)
3  *
4  */
5 
6 /*
7  * Copyright (c) 2001 and 2004 Jerome Novak
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 char gval_from_spectral_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valencia/gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $" ;
27 
28 
29 /*
30  * $Id: gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $
31  * $Log: gval_from_spectral.C,v $
32  * Revision 1.14 2014/10/13 08:53:48 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.13 2014/10/06 15:13:22 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.12 2009/10/28 13:40:23 j_novak
39  * General case for the theta symmetry (now should work).
40  *
41  * Revision 1.11 2009/10/21 13:19:04 j_novak
42  * Going back (temporary) to previous version.
43  *
44  * Revision 1.9 2007/12/21 10:46:29 j_novak
45  * In "from_spectral..." functions: better treatment of ETATZERO case.
46  *
47  * Revision 1.8 2007/11/02 16:49:12 j_novak
48  * Suppression of intermediate array for spectral summation.
49  *
50  * Revision 1.7 2006/10/02 07:41:03 j_novak
51  * Corrected an error in the case r=0, when exporting to a cartesian grid.
52  *
53  * Revision 1.6 2005/06/23 13:44:18 j_novak
54  * Removed some old comments.
55  *
56  * Revision 1.5 2005/06/23 13:40:08 j_novak
57  * The tests on the number of dimensions have been changed to handle better the
58  * axisymmetric case.
59  *
60  * Revision 1.4 2005/06/22 09:11:17 lm_lin
61  *
62  * Grid wedding: convert from the old C++ object "Cmp" to "Scalar".
63  *
64  * Revision 1.3 2004/12/17 13:35:04 m_forot
65  * Add the case T_LEG
66  *
67  * Revision 1.2 2004/05/07 13:19:24 j_novak
68  * Prevention of warnings
69  *
70  * Revision 1.1 2004/05/07 12:32:13 j_novak
71  * New summation from spectral to FD grid. Much faster!
72  *
73  *
74  * $Header: /cvsroot/Lorene/C++/Source/Valencia/gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $
75  *
76  */
77 
78 #include <cmath>
79 
80 // Lorene headers
81 #include "grille_val.h"
82 #include "proto_f77.h"
83 
84 
85  //--------------------------------------
86  // Sommation depuis une grille spectrale
87  //--------------------------------------
88 
89 namespace Lorene {
90 void Grille_val::somme_spectrale1(const Scalar& meudon, double* resu, int taille_in) const {
91 
92  int taille = dim.dim[0]+2*nfantome ;
93  if (taille != taille_in) {
94  cout << "Gval_spher::somme_spectral2():\n" ;
95  cout << "grid size incompatible with array size... exiting!" << endl ;
96  abort() ;
97  }
98  int nrv = dim.dim[0]+nfantome ;
99  const Map& mp = meudon.get_mp() ;
100  int l ;
101  double xi ;
102  for (int i=0; i<nfantome; i++) resu[i] = 0 ;
103  for (int i=nfantome; i<nrv; i++) {
104  mp.val_lx(zr->t[i],0.,0.,l,xi) ;
105  resu[i] = meudon.get_spectral_va().val_point_jk(l, xi, 0, 0) ;
106  }
107  for (int i=nrv; i<taille; i++) resu[i] = 0 ;
108 }
109 
110 void Gval_cart::somme_spectrale2(const Scalar& meudon, double* resu, int taille_in) const {
111  int nzv = dim.dim[0] + nfantome ;
112  int nxv = dim.dim[1] + nfantome ;
113  int nzv2 = dim.dim[0] + 2*nfantome ;
114  int nxv2 = dim.dim[1] + 2*nfantome ;
115  int taille = nxv2*nzv2 ;
116  if (taille != taille_in) {
117  cout << "Gval_spher::somme_spectral2():\n" ;
118  cout << "grid size incompatible with array size... exiting!" << endl ;
119  abort() ;
120  }
121  const Map& mp = meudon.get_mp() ;
122  int l ;
123  double xi0, rr, theta ;
124  double phi = 0 ;
125  int inum = 0 ;
126  for (int ix=0; ix<nfantome; ix++) {
127  for (int iz=0; iz<nzv2; iz++) {
128  resu[inum] = 0. ;
129  inum++ ;
130  }
131  }
132  for (int ix=nfantome; ix<nxv; ix++) {
133  for (int iz=0; iz<nfantome; iz++) {
134  resu[inum] = 0. ;
135  inum++ ;
136  }
137  double xx2 = (x->t[ix])*(x->t[ix]) ;
138  for (int iz=nfantome; iz<nzv; iz++) {
139  rr = sqrt((zr->t[iz])*(zr->t[iz]) + xx2) ;
140  theta = (rr != 0. ? acos((zr->t[iz])/rr) : 0) ;
141  mp.val_lx(rr, theta, phi, l, xi0) ;
142  resu[inum] = meudon.get_spectral_va().val_point(l, xi0, theta, phi) ;
143  inum++ ;
144  }
145  for (int iz=nzv; iz<nzv2; iz++) {
146  resu[inum] = 0. ;
147  inum++ ;
148  }
149  }
150  for (int ix=nxv; ix<nxv2; ix++) {
151  for (int iz=0; iz<nzv2; iz++) {
152  resu[inum] = 0. ;
153  inum++ ;
154  }
155  }
156 }
157 
158 void Gval_cart::somme_spectrale3(const Scalar& meudon, double* resu, int taille_in) const{
159  int nzv = dim.dim[0] + nfantome ;
160  int nxv = dim.dim[1] + nfantome ;
161  int nyv = dim.dim[2] + nfantome ;
162  int nzv2 = dim.dim[0] + 2*nfantome ;
163  int nxv2 = dim.dim[1] + 2*nfantome ;
164  int nyv2 = dim.dim[2] + 2*nfantome ;
165  int taille = nyv2*nxv2*nzv2 ;
166  if (taille != taille_in) {
167  cout << "Gval_spher::somme_spectral2():\n" ;
168  cout << "grid size incompatible with array size... exiting!" << endl ;
169  abort() ;
170  }
171  const Map& mp = meudon.get_mp() ;
172  int l ;
173  double xi0, rr, theta, phi ;
174  int inum = 0 ;
175  for (int iy=0; iy<nfantome; iy++) {
176  for (int ix=0; ix<nxv2; ix++) {
177  for (int iz=0; iz<nzv2; iz++){
178  resu[inum] = 0. ;
179  inum++ ;
180  }
181  }
182  }
183  for (int iy=nfantome; iy<nyv; iy++) {
184  double yy = x->t[iy] ;
185  double yy2 = yy*yy ;
186  for (int ix=0; ix<nfantome; ix++) {
187  for (int iz=0; iz<nzv2; iz++) {
188  resu[inum] = 0. ;
189  inum++ ;
190  }
191  }
192  for (int ix=nfantome; ix<nxv; ix++) {
193  for (int iz=0; iz<nfantome; iz++) {
194  resu[inum] = 0. ;
195  inum++ ;
196  }
197  double xx = x->t[ix] ;
198  double xx2 = xx*xx ;
199  for (int iz=nfantome; iz<nzv; iz++) {
200  rr = sqrt((zr->t[iz])*(zr->t[iz]) + xx2 + yy2) ;
201  theta = (rr != 0. ? acos((zr->t[iz])/rr) : 0. );
202  phi = (rr != 0. ? atan2(yy, xx) : 0. ) ; // return value in [-M_PI,M_PI], should work
203  mp.val_lx(rr, theta, phi, l, xi0) ;
204  resu[inum] = meudon.get_spectral_va().val_point(l, xi0, theta, phi) ;
205  inum++ ;
206  }
207  for (int iz=nzv; iz<nzv2; iz++) {
208  resu[inum] = 0. ;
209  inum++ ;
210  }
211  }
212  for (int ix=nxv; ix<nxv2; ix++) {
213  for (int iz=0; iz<nzv2; iz++) {
214  resu[inum] = 0. ;
215  inum++ ;
216  }
217  }
218  }
219  for (int iy=nyv; iy<nyv2; iy++) {
220  for (int ix=0; ix<nxv2; ix++) {
221  for (int iz=0; iz<nzv2; iz++){
222  resu[inum] = 0. ;
223  inum++ ;
224  }
225  }
226  }
227 }
228 
229 void Gval_spher::somme_spectrale2(const Scalar& meudon, double* resu, int taille_in) const {
230 
231  assert (dim.ndim >=2) ;
232  int nrv = dim.dim[0] + nfantome ;
233  int ntv = dim.dim[1] + nfantome ;
234  int nrv2 = dim.dim[0] + 2*nfantome ;
235  int ntv2 = dim.dim[1] + 2*nfantome ;
236  int taille = ntv2*nrv2 ;
237  if (taille != taille_in) {
238  cout << "Gval_spher::somme_spectral2():\n" ;
239  cout << "grid size incompatible with array size... exiting!" << endl ;
240  abort() ;
241  }
242  const Map& mp = meudon.get_mp() ;
243  int l ;
244  double xi, rr, theta ;
245  double phi0 = 0 ;
246  int inum = 0 ;
247  for (int it=0; it<nfantome; it++) {
248  for (int ir=0; ir<nrv2; ir++) {
249  resu[inum] = 0. ;
250  inum++ ;
251  }
252  }
253  for (int it=nfantome; it<ntv; it++) {
254  for (int ir=0; ir<nfantome; ir++) {
255  resu[inum] = 0. ;
256  inum++ ;
257  }
258  theta = tet->t[it] ;
259  for (int ir=nfantome; ir<nrv; ir++) {
260  rr = zr->t[ir] ;
261  mp.val_lx(rr, theta, phi0, l, xi) ;
262  resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
263  inum++ ;
264  }
265  for (int ir=nrv; ir<nrv2; ir++) {
266  resu[inum] = 0. ;
267  inum++ ;
268  }
269  }
270  for (int it=ntv; it<ntv2; it++) {
271  for (int ir=0; ir<nrv2; ir++) {
272  resu[inum] = 0. ;
273  inum++ ;
274  }
275  }
276 }
277 
278 double* Gval_spher::somme_spectrale2ri(const Scalar& meudon) const {
279  int nrv = dim.dim[0] + 1 + nfantome ;
280  int ntv = dim.dim[1] + nfantome ;
281  int nrv2 = dim.dim[0] + 1 + 2*nfantome ;
282  int ntv2 = dim.dim[1] + 2*nfantome ;
283  int taille = ntv2*nrv2 ;
284  const Map& mp = meudon.get_mp() ;
285  double* resu = new double[taille] ;
286  int l ;
287  double xi, rr, theta ;
288  double phi0 = 0 ;
289  int inum = 0 ;
290  for (int it=0; it<nfantome; it++) {
291  for (int ir=0; ir<nrv2; ir++) {
292  resu[inum] = 0. ;
293  inum++ ;
294  }
295  }
296  for (int it=nfantome; it<ntv; it++) {
297  for (int ir=0; ir<nfantome; ir++) {
298  resu[inum] = 0. ;
299  inum++ ;
300  }
301  theta = tet->t[it] ;
302  for (int ir=nfantome; ir<nrv; ir++) {
303  rr = zri->t[ir] ;
304  mp.val_lx(rr, theta, phi0, l, xi) ;
305  resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
306  inum++ ;
307  }
308  for (int ir=nrv; ir<nrv2; ir++) {
309  resu[inum] = 0. ;
310  inum++ ;
311  }
312  }
313  for (int it=ntv; it<ntv2; it++) {
314  for (int ir=0; ir<nrv2; ir++) {
315  resu[inum] = 0. ;
316  inum++ ;
317  }
318  }
319  return resu ;
320 }
321 
322 double* Gval_spher::somme_spectrale2ti(const Scalar& meudon) const {
323  int nrv = dim.dim[0] + nfantome ;
324  int ntv = dim.dim[1] + 1 + nfantome ;
325  int nrv2 = dim.dim[0] + 2*nfantome ;
326  int ntv2 = dim.dim[1] + 1 + 2*nfantome ;
327  int taille = ntv2*nrv2 ;
328  const Map& mp = meudon.get_mp() ;
329  double* resu = new double[taille] ;
330  int l ;
331  double xi, rr, theta ;
332  double phi0 = 0 ;
333  int inum = 0 ;
334  for (int it=0; it<nfantome; it++) {
335  for (int ir=0; ir<nrv2; ir++) {
336  resu[inum] = 0. ;
337  inum++ ;
338  }
339  }
340  for (int it=nfantome; it<ntv; it++) {
341  for (int ir=0; ir<nfantome; ir++) {
342  resu[inum] = 0. ;
343  inum++ ;
344  }
345  theta = teti->t[it] ;
346  for (int ir=nfantome; ir<nrv; ir++) {
347  rr = zr->t[ir] ;
348  mp.val_lx(rr, theta, phi0, l, xi) ;
349  resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
350  inum++ ;
351  }
352  for (int ir=nrv; ir<nrv2; ir++) {
353  resu[inum] = 0. ;
354  inum++ ;
355  }
356  }
357  for (int it=ntv; it<ntv2; it++) {
358  for (int ir=0; ir<nrv2; ir++) {
359  resu[inum] = 0. ;
360  inum++ ;
361  }
362  }
363  return resu ;
364 }
365 
366 void Gval_spher::somme_spectrale3(const Scalar& meudon, double* resu, int taille_in) const{
367 
368  assert(meudon.get_etat() == ETATQCQ) ;
369  meudon.get_spectral_va().coef() ;
370 
371  //Sizes of both grids
372  //-------------------
373  int nrv0 = dim.dim[0] ;
374  int ntv0 = dim.dim[1] ;
375  int nrv = dim.dim[0] + nfantome ;
376  int ntv = dim.dim[1] + nfantome ;
377  int npv = dim.dim[2] + nfantome ;
378  int nrv2 = dim.dim[0] + 2*nfantome ;
379  int ntv2 = dim.dim[1] + 2*nfantome ;
380  int npv2 = dim.dim[2] + 2*nfantome ;
381  int taille = npv2*ntv2*nrv2 ;
382  if (taille != taille_in) {
383  cout << "Gval_spher::somme_spectral3():\n" ;
384  cout << "grid size incompatible with array size... exiting!" << endl ;
385  abort() ;
386  }
387  const Map& mp = meudon.get_mp() ;
388 #ifndef NDEBUG
389  const Map_af* mpaff = dynamic_cast<const Map_af*>(&mp) ;
390  assert(mpaff != 0x0) ;
391 #endif
392  const Mg3d* mg = mp.get_mg() ;
393  int ntm = mg->get_nt(0) ;
394  int npm = mg->get_np(0) ;
395  int nz = mg->get_nzone() ;
396 #ifndef NDEBUG
397  for (int lz=1; lz<nz; lz++) {
398  assert (ntm == mg->get_nt(lz)) ; //Same angular grids in all domains...
399  assert (npm == mg->get_np(lz)) ;
400  }
401 #endif
402 
403  //Intermediate quantities
404  //-----------------------
405  double* alpha = new double[nrv0*(npm+2)*ntm] ;
406  double* p_coef = alpha ;
407  double* chebnri = 0x0 ; //size ~ nrv0 * (npm+2) * nr ...
408  int* idom = 0x0 ;
409  initialize_spectral_r(mp, meudon.get_spectral_va().get_base(), idom, chebnri) ;
410  double* p_func = chebnri ;
411  Mtbl_cf& mtbcf = *meudon.get_spectral_va().c_cf ;
412  double** coefm = new double*[nz] ;
413  for (int lz=0; lz<nz; lz++) {
414  assert((mtbcf.t[lz])->get_etat() != ETATNONDEF) ;
415  coefm[lz] = (mtbcf.t[lz])->t ;
416  if (coefm[lz] == 0x0) {
417  int sizem = mg->get_nr(lz)*ntm*(npm+2) ;
418  coefm[lz] = new double[sizem] ;
419  double* pcf = coefm[lz] ;
420  for (int i=0; i<sizem; i++)
421  pcf[i] = 0. ;
422  }
423  }
424 
425  //First partial summation
426  //-----------------------
427  for (int irv=0; irv<nrv0; irv++) {
428  int lz = idom[irv] ;
429  double* tbcf = coefm[lz] ;
430  int nrm = mg->get_nr(lz) ;
431  for (int mpm=0; mpm<npm+2; mpm++) {
432  for (int ltm=0; ltm<ntm; ltm++) {
433  *p_coef = 0 ;
434  for (int irm=0; irm<nrm; irm++) {
435  *p_coef += (*tbcf)*(*p_func) ;
436  tbcf++ ;
437  p_func++ ;
438 // cout << *p_func << ", " << *tbcf << ", " << *p_coef << endl ;
439  }
440  p_coef++ ;
441  }
442  }
443  }
444 
445  for (int lz=0; lz<nz; lz++) {
446  if ((mtbcf.t[lz])->t == 0x0) delete [] coefm[lz] ;
447  }
448  delete [] coefm ;
449  delete [] chebnri ;
450  delete [] idom ;
451 
452  double* beta = new double[ntv0*nrv0*(npm+2)] ;
453  p_coef = beta ;
454  double* tetlj = 0x0 ;
455  initialize_spectral_theta(mp, meudon.get_spectral_va().get_base(), tetlj) ;
456  p_func = tetlj ;
457  double* p_interm = alpha ;
458 
459  //Second partial summation
460  //------------------------
461  for (int jtv=0; jtv<ntv0; jtv++) {
462  for (int irv=0; irv<nrv0; irv++) {
463  for (int mpm=0; mpm<npm+2; mpm++) {
464  *p_coef = 0 ;
465  for (int ltm=0; ltm<ntm; ltm++) {
466  *p_coef += (*p_interm) * (*p_func) ;
467  p_interm++ ;
468  p_func++ ;
469  }
470  p_coef++ ;
471  } // Loop on m
472  p_func -= (npm+2)*ntm ;
473  } //Loop on irv
474  p_interm = alpha ;
475  p_func += (npm+2)*ntm ;
476  } //Loop on jtv
477 
478  delete [] alpha ;
479  delete [] tetlj ;
480 
481 
482 
483  // Final summation
484  //----------------
485  p_interm = beta ;
486  double* expmk = 0x0 ;
487  initialize_spectral_phi(mp, meudon.get_spectral_va().get_base(), expmk) ;
488  p_func = expmk ;
489  p_coef = resu ;
490  for (int ip=0; ip<nfantome; ip++) {
491  for (int it=0; it<ntv2; it++) {
492  for (int ir=0; ir<nrv2; ir++){
493  *p_coef = 0. ;
494  p_coef++ ;
495  }
496  }
497  }
498  for (int ip=nfantome; ip<npv; ip++) {
499  for (int it=0; it<nfantome; it++) {
500  for (int ir=0; ir<nrv2; ir++) {
501  *p_coef = 0. ;
502  p_coef++ ;
503  }
504  }
505  for (int it=nfantome; it<ntv; it++) {
506  for (int ir=0; ir<nfantome; ir++) {
507  *p_coef = 0. ;
508  p_coef++ ;
509  }
510  for (int ir=nfantome; ir<nrv; ir++) {
511  *p_coef = 0. ;
512  for (int mpm=0; mpm<npm+2; mpm++) {
513  *p_coef += (*p_interm) * (*p_func) ;
514  p_interm++ ;
515  p_func++ ;
516  }
517  p_coef++ ;
518  p_func -= (npm+2) ;
519  }
520  for (int ir=nrv; ir<nrv2; ir++) {
521  *p_coef = 0. ;
522  p_coef++ ;
523  }
524  }
525  for (int it=ntv; it<ntv2; it++) {
526  for (int ir=0; ir<nrv2; ir++) {
527  *p_coef = 0. ;
528  p_coef++ ;
529  }
530  }
531  p_func += npm+2 ; //Next point in phi
532  p_interm = beta ;
533  }
534  for (int ip=npv; ip<npv2; ip++) {
535  for (int it=0; it<ntv2; it++) {
536  for (int ir=0; ir<nrv2; ir++){
537  *p_coef = 0. ;
538  p_coef++ ;
539  }
540  }
541  }
542  delete [] expmk ;
543  delete [] beta ;
544 }
545 
546 
547 void Gval_spher::initialize_spectral_r(const Map& mp, const Base_val& base,
548  int*& idom, double*& chebnri) const {
549 
550  int nrv0 = dim.dim[0] ;
551  const Mg3d* mg = mp.get_mg() ;
552  int npm = mg->get_np(0) ;
553  int ntm = mg->get_nt(0) ;
554 
555  assert (idom == 0x0) ;
556  idom = new int[nrv0] ;
557  double* xi = new double[nrv0] ;
558  int nrmax = 0 ;
559 
560  for (int i=0; i<nrv0; i++) {
561  mp.val_lx(zr->t[i+nfantome], 0., 0., idom[i], xi[i]) ;
562  nrmax += mg->get_nr(idom[i]) ;
563  }
564 
565  assert (chebnri == 0x0) ;
566  chebnri = new double[(npm+2)*ntm*nrmax] ;
567  double* p_out = chebnri ;
568  for (int irv=0; irv<nrv0; irv++) {
569  bool nucleus = (mg->get_type_r(idom[irv]) == RARE) ;
570  int nmax = (nucleus ? 2*mg->get_nr(idom[irv]) + 1
571  : mg->get_nr(idom[irv])) ;
572  double* cheb = new double[nmax] ;
573  cheb[0] = 1. ;
574  cheb[1] = xi[irv] ;
575  for (int ir=2; ir<nmax; ir++) {
576  cheb[ir] = 2*xi[irv]*cheb[ir-1] - cheb[ir-2] ;
577  }
578 
579  int base_r = base.get_base_r(idom[irv]) ;
580 
581  for (int ip=0; ip<npm+2; ip++) {
582  for (int it=0; it<ntm; it++) {
583  int fact = 1 ;
584  int par = 0 ;
585  if (nucleus) {
586  fact = 2 ;
587  switch (base_r) {
588 
589  case R_CHEBP : {
590  break ;
591  }
592 
593  case R_CHEBI : {
594  par = 1 ;
595  break ;
596  }
597 
598  case R_CHEBPI_P : {
599  par = it % 2 ;
600  break ;
601  }
602 
603  case R_CHEBPI_I : {
604  par = 1 - (it % 2) ;
605  break ;
606  }
607  case R_CHEBPIM_P : {
608  par = (ip/2) % 2 ;
609  break ;
610  }
611 
612  case R_CHEBPIM_I : {
613  par = 1 - ((ip/2) % 2) ;
614  break ;
615  }
616 
617  default : {
618  cout << "Gval_spher::initialize_spectral_r : " << '\n'
619  << "Unexpected radial base !" << '\n'
620  << "Base : " << base_r << endl ;
621  abort() ;
622  break ;
623  }
624  }
625  }
626  for (int ir=0; ir<mg->get_nr(idom[irv]); ir++) {
627  *p_out = cheb[fact*ir+par] ;
628  p_out++ ;
629  }
630 
631  } // Loop on it
632  } // Loop on ip
633  delete [] cheb ;
634 
635  }// Loop on irv
636 
637  delete [] xi ;
638 
639 }
640 
641 void Gval_spher::initialize_spectral_theta(const Map& mp, const Base_val& base,
642  double*& tetlj) const {
643 
644  int ntv0 = dim.dim[1] ;
645  const Mg3d* mg = mp.get_mg() ;
646  int npm = mg->get_np(0) ;
647  int ntm = mg->get_nt(0) ;
648  int base_t = base.get_base_t(0) ;
649 
650  assert (tetlj == 0x0) ;
651  tetlj = new double[(npm+2)*ntv0*ntm] ;
652  double* p_out = tetlj ;
653 
654  for (int jtv=0; jtv<ntv0; jtv++) {
655  double teta = tet->t[jtv+nfantome] ;
656  for (int mpm=0; mpm < npm+2; mpm++) {
657  for (int ltm=0; ltm<ntm; ltm++) {
658  switch (base_t) { //## One should use array of functions...
659  case T_COS : {
660  *p_out = cos(ltm*teta) ;
661  break ;
662  }
663  case T_SIN : {
664  *p_out = sin(ltm*teta) ;
665  break ;
666  }
667  case T_COS_P : {
668  *p_out = cos(2*ltm*teta) ;
669  break ;
670  }
671  case T_COS_I : {
672  *p_out = cos((2*ltm+1)*teta) ;
673  break ;
674  }
675  case T_SIN_P : {
676  *p_out = sin(2*ltm*teta) ;
677  break ;
678  }
679  case T_SIN_I : {
680  *p_out = sin((2*ltm+1)*teta) ;
681  break ;
682  }
683  case T_COSSIN_CP : {
684  *p_out = ( ((mpm/2) % 2 == 0) ? cos(2*ltm*teta)
685  : sin((2*ltm+1)*teta)) ;
686  break ;
687  }
688  case T_COSSIN_CI : {
689  *p_out = ( ((mpm/2) % 2 == 0) ? cos((2*ltm+1)*teta)
690  : sin(2*ltm*teta)) ;
691  break ;
692  }
693  case T_COSSIN_SP : {
694  *p_out = ( ((mpm/2) % 2 == 0) ? sin(2*ltm*teta)
695  : cos((2*ltm+1)*teta)) ;
696  break ;
697  }
698  case T_COSSIN_SI : {
699  *p_out = ( ((mpm/2) % 2 == 0) ? sin((2*ltm+1)*teta)
700  : cos(2*ltm*teta)) ;
701  break ;
702  }
703  case T_COSSIN_C : {
704  *p_out = ( ((mpm/2) % 2 == 0) ? cos(ltm*teta)
705  : sin(ltm*teta)) ;
706  break ;
707  }
708  case T_COSSIN_S : {
709  *p_out = ( ((mpm/2) % 2 == 0) ? sin(ltm*teta)
710  : cos(ltm*teta)) ;
711  break ;
712  }
713  default : {
714  cout << "Gval_spher::initialize_spectral_theta : " << '\n'
715  << "Unexpected theta base !" << '\n'
716  << "Base : " << base_t << endl ;
717  abort() ;
718  break ;
719  }
720  }
721  p_out++ ;
722  }
723  if ( (base_t == T_COS_I) || (base_t == T_SIN_P) || (base_t == T_SIN_I) )
724  {
725  p_out-- ;
726  *p_out = 0. ;
727  p_out++ ;
728  }
729  } //Loop on mpm
730  } //Loop on jtv
731 
732 }
733 
734 
735 void Gval_spher::initialize_spectral_phi(const Map& mp, const Base_val& base,
736  double*& expmk) const {
737 
738  int npv0 = dim.dim[2] ;
739  const Mg3d* mg = mp.get_mg() ;
740  int npm = mg->get_np(0) ;
741  int base_p = base.get_base_p(0) ;
742 
743  assert (expmk == 0x0) ;
744  expmk = new double[(npm+2)*npv0] ;
745  double* p_out = expmk ;
746 
747  for (int kpv=0; kpv<npv0; kpv++) {
748  double fi = phi->t[kpv+nfantome] ;
749  for (int mpm=0; mpm < npm+2; mpm++) {
750  switch (base_p) { //## One should use array of functions...
751  case P_COSSIN : {
752  int m = mpm / 2 ;
753  *p_out = ( (mpm%2 == 0) ? cos(m*fi) : sin(m*fi) ) ;
754  break ;
755  }
756  case P_COSSIN_P : {
757  int m = mpm / 2 ;
758  *p_out = ( (mpm%2 == 0) ? cos(2*m*fi) : sin(2*m*fi) ) ;
759  break ;
760  }
761  case P_COSSIN_I : {
762  int m = mpm / 2 ;
763  *p_out = ( (mpm%2 == 0) ? cos((2*m+1)*fi) : sin((2*m+1)*fi) ) ;
764  break ;
765  }
766  default : {
767  cout << "Gval_spher::initialize_spectral_phi : " << '\n'
768  << "Unexpected phi base !" << '\n'
769  << "Base : " << base_p << endl ;
770  abort() ;
771  break ;
772  }
773  }
774  p_out++ ;
775  } //Loop on mpm
776  } //Loop on kpv
777 
778 }
779 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
virtual void somme_spectrale2(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
Tbl * zri
Arrays containing the values of coordinate z (or r) on the interfaces.
Definition: grille_val.h:126
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
#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
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:411
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Base class for coordinate mappings.
Definition: map.h:670
double * somme_spectrale2ri(const Scalar &meudon) const
Same as before but at radial grid interfaces.
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition: base_val.h:400
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
Definition: valeur.C:882
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
double * t
The array of double.
Definition: tbl.h:173
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes.
Definition: grille_val.h:124
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:422
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Definition: valeur.C:900
double * somme_spectrale2ti(const Scalar &meudon) const
Same as before but at angular grid interfaces.
virtual void somme_spectrale2(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
#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
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
Definition: dim_tbl.h:101
int nfantome
The number of hidden cells (same on each side)
Definition: grille_val.h:104
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
Bases of the spectral expansions.
Definition: base_val.h:322
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
Affine radial mapping.
Definition: map.h:2027
Dim_tbl dim
The dimensions of the grid.
Definition: grille_val.h:102
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:169
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
int * dim
Array of dimensions (size: ndim).
Definition: dim_tbl.h:102
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point given by its physical coordina...
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
void somme_spectrale1(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...