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