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
26char 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
89namespace Lorene {
90void 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
110void 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
158void 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
229void 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
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
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
366void 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
547void 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
641void 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
735void 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}
Bases of the spectral expansions.
Definition base_val.h:322
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
int * dim
Array of dimensions (size: ndim).
Definition dim_tbl.h:102
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
Definition dim_tbl.h:101
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
int nfantome
The number of hidden cells (same on each side)
Definition grille_val.h:104
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...
Dim_tbl dim
The dimensions of the grid.
Definition grille_val.h:102
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
Definition grille_val.h:124
Tbl * zri
Arrays containing the values of coordinate z (or r) on the interfaces.
Definition grille_val.h:126
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...
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
Tbl * x
Arrays containing the values of coordinate x on the nodes.
Definition grille_val.h:343
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
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 * phi
Arrays containing the values of coordinate on the nodes.
Definition grille_val.h:591
Tbl * teti
Arrays containing the values of coordinate on the interfaces.
Definition grille_val.h:589
double * somme_spectrale2ri(const Scalar &meudon) const
Same as before but at radial grid interfaces.
double * somme_spectrale2ti(const Scalar &meudon) const
Same as before but at angular grid interfaces.
Tbl * tet
Arrays containing the values of coordinate on the nodes.
Definition grille_val.h:587
Affine radial mapping.
Definition map.h:2027
Base class for coordinate mappings.
Definition map.h:670
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 Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
double * t
The array of double.
Definition tbl.h:173
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:169
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define T_SIN_P
dev. sin seulement, harmoniques paires
#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_COS_P
dev. cos seulement, harmoniques paires
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
#define P_COSSIN
dev. standart
#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 T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
#define T_SIN_I
dev. sin seulement, harmoniques impaires
#define T_COS
dev. cos seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define T_SIN
dev. sin seulement
#define T_COS_I
dev. cos seulement, harmoniques impaires
#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