LORENE
map_et.C
1/*
2 * Methods of class Map_et
3 */
4
5/*
6 * Copyright (c) 1999-2001 Eric Gourgoulhon
7 *
8 * This file is part of LORENE.
9 *
10 * LORENE is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
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
27char map_et_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et.C,v 1.16 2014/10/13 08:53:03 j_novak Exp $" ;
28
29/*
30 * $Id: map_et.C,v 1.16 2014/10/13 08:53:03 j_novak Exp $
31 * $Log: map_et.C,v $
32 * Revision 1.16 2014/10/13 08:53:03 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.15 2014/10/06 15:13:13 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.14 2013/06/05 15:10:42 j_novak
39 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
40 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
41 *
42 * Revision 1.13 2008/09/29 13:23:51 j_novak
43 * Implementation of the angular mapping associated with an affine
44 * mapping. Things must be improved to take into account the domain index.
45 *
46 * Revision 1.12 2008/08/27 08:48:26 jl_cornou
47 * Added_R_JACO02 case
48 *
49 * Revision 1.11 2005/11/30 11:09:07 p_grandclement
50 * Changes for the Bin_ns_bh project
51 *
52 * Revision 1.10 2004/03/25 10:29:23 j_novak
53 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
54 *
55 * Revision 1.9 2004/01/29 08:50:03 p_grandclement
56 * Modification of Map::operator==(const Map&) and addition of the surface
57 * integrales using Scalar.
58 *
59 * Revision 1.8 2003/10/15 10:36:52 e_gourgoulhon
60 * In method fait_poly(): changed local variable name x to x1, not to shadow
61 * Coord's x.
62 *
63 * Revision 1.7 2003/07/07 20:01:43 p_grandclement
64 * change assert in constructor for map_et from a surface
65 *
66 * Revision 1.6 2003/06/04 21:11:55 p_grandclement
67 * Correction of separation in odd-even harmonics
68 *
69 * Revision 1.5 2002/10/16 14:36:41 j_novak
70 * Reorganization of #include instructions of standard C++, in order to
71 * use experimental version 3 of gcc.
72 *
73 * Revision 1.4 2002/05/07 07:10:44 e_gourgoulhon
74 * Compatibilty with xlC compiler on IBM SP2:
75 * suppressed the parenthesis around argument of instruction new:
76 * e.g. aa = new (Tbl*[nzone]) ---> aa = new Tbl*[nzone]
77 * result = new (Param) ---> result = new Param
78 *
79 * Revision 1.3 2002/01/15 15:53:06 p_grandclement
80 * I have had a constructor fot map_et using the equation of the surface
81 * of the star.
82 *
83 * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
84 *
85 * All writing/reading to a binary file are now performed according to
86 * the big endian convention, whatever the system is big endian or
87 * small endian, thanks to the functions fwrite_be and fread_be
88 *
89 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
90 * LORENE
91 *
92 * Revision 1.11 2001/02/28 11:04:20 eric
93 * 1ere version testee de resize.
94 *
95 * Revision 1.10 2001/02/26 17:29:42 eric
96 * Ajout de la fonction resize.
97 *
98 * Revision 1.9 2000/08/18 11:10:48 eric
99 * Ajout de l'operateur d'affectation a un autre Map_et.
100 *
101 * Revision 1.8 2000/01/24 16:42:36 eric
102 * Ajout de la fonction virtuelle operator=(const Map_af& ).
103 *
104 * Revision 1.7 2000/01/24 11:03:28 eric
105 * Correction d'une erreur dans le constructeur par lecture de fichier:
106 * ff et gg doivent etre construits sur mgi.get_angu() et non sur mgi.
107 *
108 * Revision 1.6 1999/12/20 10:24:49 eric
109 * Ajout des fonctions de lecture des parametres de Map_et:
110 * get_alpha(), get_beta(), get_ff(), get_gg().
111 *
112 * Revision 1.5 1999/12/17 11:20:08 eric
113 * Ajout de la fonction homothetie.
114 *
115 * Revision 1.4 1999/12/17 09:14:30 eric
116 * Amelioration de l'affichage.
117 *
118 * Revision 1.3 1999/11/24 16:31:41 eric
119 * Ajout des fonctions set_ff et set_gg.
120 *
121 * Revision 1.2 1999/11/24 11:22:44 eric
122 * Map_et : fonctions de constructions amies.
123 *
124 * Revision 1.1 1999/11/22 10:37:36 eric
125 * Initial revision
126 *
127 *
128 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et.C,v 1.16 2014/10/13 08:53:03 j_novak Exp $
129 *
130 */
131
132// headers C
133#include <cmath>
134
135// headers Lorene
136#include "proto.h"
137#include "map.h"
138#include "utilitaires.h"
139#include "unites.h"
140
141 //--------------//
142 // Constructors //
143 //--------------//
144
145// -----------------------
146// Constructor from a grid
147// -----------------------
148namespace Lorene {
149Map_et::Map_et(const Mg3d& mgrille, const double* bornes)
150 : Map_radial(mgrille),
151 aasx( mgrille.get_nr(0) ),
152 aasx2( mgrille.get_nr(0) ),
153 zaasx( mgrille.get_nr(mgrille.get_nzone()-1) ),
154 zaasx2( mgrille.get_nr(mgrille.get_nzone()-1) ),
155 bbsx( mgrille.get_nr(0) ),
156 bbsx2( mgrille.get_nr(0) ),
157 ff(mgrille.get_angu()) ,
158 gg(mgrille.get_angu())
159{
160 // The Coord rsxdxdr and rsx2drdx are constructed by the default Coord
161 // constructor
162
163 // Assignement of the building functions of the Coord's
164 // ----------------------------------------------------
165 set_coord() ;
166
167
168 // alpha and beta
169 // --------------
170 int nzone = mg->get_nzone() ;
171
172 alpha = new double[nzone] ;
173 beta = new double[nzone] ;
174
175 for (int l=0 ; l<nzone ; l++) {
176 switch (mg->get_type_r(l)) {
177 case RARE: {
178 alpha[l] = bornes[l+1] - bornes[l] ;
179 beta[l] = bornes[l] ;
180 break ;
181 }
182
183 case FIN: {
184 alpha[l] = (bornes[l+1] - bornes[l]) * .5 ;
185 beta[l] = (bornes[l+1] + bornes[l]) * .5 ;
186 break ;
187 }
188
189 case UNSURR: {
190 double umax = double(1) / bornes[l] ;
191 double umin = double(1) /bornes[l+1] ;
192 alpha[l] = (umin - umax) * double(0.5) ; // u est une fonction decroissante
193 beta[l] = (umin + umax) * double(0.5) ; // de l'indice i en r
194 break ;
195 }
196
197 default: {
198 cout << "Map_et::Map_et: unkown type_r ! " << endl ;
199 abort () ;
200 break ;
201 }
202
203 }
204 } // End of the loop onto the domains
205
206
207 // Radial polynomials A(x) and B(x)
208 // --------------------------------
209
210 fait_poly() ;
211
212 // Initialisation at zero of the functions F(theta',phi') and G(theta',phi')
213 // -------------------------------------------------------------------------
214
215 ff.set_etat_zero() ;
216 gg.set_etat_zero() ;
217
218 ff.std_base_scal() ; // Standard spectral bases for F
219 gg.std_base_scal() ; // Standard spectral bases for G
220
221}
222
223Map_et::Map_et(const Mg3d& grille, const double* r_lim, const Tbl& S_0) :
224 Map_radial(grille),
225 aasx(grille.get_nr(0) ),
226 aasx2(grille.get_nr(0) ),
227 zaasx(grille.get_nr(grille.get_nzone()-1) ),
228 zaasx2(grille.get_nr(grille.get_nzone()-1) ),
229 bbsx(grille.get_nr(0) ),
230 bbsx2(grille.get_nr(0) ),
231 ff(grille.get_angu()) ,
232 gg(grille.get_angu()) {
233
234 assert (S_0.get_ndim() == 2) ;
235 assert (S_0.get_dim(0) == grille.get_nt(0)) ;
236 assert (S_0.get_dim(1) == grille.get_np(0)) ;
237
238 Map_et mapping (grille, r_lim) ;
239
240 int nz = grille.get_nzone() ;
241 assert (nz >2) ;
242
243 // Le noyau :
244 int np = grille.get_np(0) ;
245 int nt = grille.get_nt(0) ;
246
247 double * cf = new double [nt*(np+2)] ;
248 for (int k=0 ; k<np ; k++)
249 for (int j=0 ; j<nt ; j++)
250 cf[k*nt+j] = S_0(k,j) - S_0(0,0) ;
251
252 int* deg = new int [3] ;
253 deg[0] = np ;
254 deg[1] = nt ;
255 deg[2] = 1 ;
256
257 int* dim = new int [3] ;
258 dim[0] = np+2 ;
259 dim[1] = nt ;
260 dim[2] = 1 ;
261
262 Tbl ff_nucleus (np,nt) ;
263 ff_nucleus.set_etat_qcq() ;
264
265 Tbl gg_nucleus (np,nt) ;
266 gg_nucleus.set_etat_qcq() ;
267
268 // On recupere la base en phi :
269 int base_p = grille.std_base_scal().get_base_p(0) ;
270 // Selon les cas (pas tres propre mais bon ...)
271 double * odd ;
272 double * even ;
273 double * coloc_odd ;
274 double * coloc_even ;
275
276 switch (base_p) {
277 case P_COSSIN:
278 cfpcossin (deg,dim,cf) ;
279
280 // Separation des harmoniques paires et impaires :
281 odd = new double [nt*(np+2)] ;
282 even = new double [nt*(np+2)] ;
283
284
285 for (int k=0 ; k<np+2 ; k++)
286 if ((k%4 == 0) || (k%4==1))
287 for (int j=0 ; j<nt ; j++) {
288 odd[k*nt+j] = 0 ;
289 even[k*nt+j] = cf[k*nt+j] ;
290 }
291 else
292 if ((k%4 == 2) || (k%4 == 3))
293 for (int j=0 ; j<nt ; j++) {
294 even[k*nt+j] = 0 ;
295 odd[k*nt+j] = cf[k*nt+j] ;
296 }
297
298 else {
299 cout << "Erreur bizzare..." << endl ;
300 abort() ;
301 }
302
303 coloc_odd = new double [nt*np] ;
304 coloc_even = new double [nt*np] ;
305
306 cipcossin (deg,dim,deg,odd,coloc_odd) ;
307 cipcossin (deg,dim,deg,even,coloc_even) ;
308 for (int k=0 ; k<np ; k++)
309 for (int j=0 ; j<nt ; j++) {
310 gg_nucleus.set(k,j) = coloc_even[k*nt+j] ;
311 ff_nucleus.set(k,j) = coloc_odd[k*nt+j] ;
312 }
313
314 delete [] even ;
315 delete [] odd ;
316 delete [] coloc_even ;
317 delete [] coloc_odd ;
318 delete[] dim ;
319 delete [] deg ;
320 delete [] cf ;
321
322 break;
323 default:
324 cout << "Base_p != P_COSSIN not implemented in Map_et constructor" <<
325 endl ;
326 abort() ;
327 }
328
329 double mu_nucleus = - min(gg_nucleus) ;
330 double alpha_nucleus = S_0(0,0)-mu_nucleus ;
331
332 ff_nucleus /= alpha_nucleus ;
333 gg_nucleus += mu_nucleus ;
334 gg_nucleus /= alpha_nucleus ;
335
336 // First shell : much simpler no ?
337 Tbl ff_shell (np,nt) ;
338 ff_shell.set_etat_qcq() ;
339 ff_shell = S_0 - S_0(0,0) ;
340
341 double lambda_shell = -max(ff_shell) ;
342
343 double R_ext = r_lim[2] ;
344
345 double beta_shell = (R_ext+S_0(0,0)-lambda_shell)/2. ;
346 double alpha_shell = (R_ext-S_0(0,0)+lambda_shell)/2. ;
347
348 ff_shell += lambda_shell ;
349 ff_shell /= alpha_shell ;
350
351 ff.annule_hard() ;
352 gg.annule_hard() ;
353
356
357 for (int k=0 ; k<np ; k++)
358 for (int j=0 ; j<nt ; j++) {
359 ff.set(0,k,j,0) = ff_nucleus(k,j) ;
360 gg.set(0,k,j,0) = gg_nucleus(k,j) ;
361 ff.set(1,k,j,0) = ff_shell(k,j) ;
362 }
363
364 gg.annule(1,nz-1) ;
365 ff.annule(2,nz-1) ;
366
367 ff.std_base_scal() ;
368 gg.std_base_scal() ;
369
370 alpha = new double[nz] ;
371 alpha[0] = alpha_nucleus ;
372 alpha[1] = alpha_shell ;
373
374 beta = new double[nz] ;
375 beta[0] = 0 ;
376 beta[1] = beta_shell ;
377 for (int i=2 ; i<nz ; i++) {
378 alpha[i] = mapping.get_alpha()[i] ;
379 beta[i] = mapping.get_beta()[i] ;
380 }
381
382 fait_poly() ;
383 set_coord() ;
384}
385// ------------------
386// Copy constructor
387// ------------------
388Map_et::Map_et(const Map_et& mpi) : Map_radial(mpi) ,
389 aasx( mpi.aasx ),
390 aasx2( mpi.aasx2 ),
391 zaasx( mpi.zaasx ),
392 zaasx2( mpi.zaasx2 ),
393 bbsx( mpi.bbsx ),
394 bbsx2( mpi.bbsx2 ),
395 ff(mpi.ff) ,
396 gg(mpi.gg)
397{
398 // Assignement of the building functions of the Coord's
399 // ----------------------------------------------------
400 set_coord() ;
401
402 // alpha and beta
403 // --------------
404 int nzone = mg->get_nzone() ;
405
406 alpha = new double[nzone] ;
407 beta = new double[nzone] ;
408
409 for (int l=0 ; l<nzone ; l++) {
410 alpha[l] = mpi.alpha[l] ;
411 beta[l] = mpi.beta[l] ;
412 }
413
414 // Radial polynomials A(x) and B(x)
415 // --------------------------------
416
417 fait_poly() ;
418
419}
420 //------------------------------------------//
421 // Modification of the mapping parameters //
422 //------------------------------------------//
423
424void Map_et::set_alpha(double alpha0, int l) {
425
426 assert(l>=0) ;
427 assert(l<mg->get_nzone()) ;
428
429 alpha[l] = alpha0 ;
430
431 reset_coord() ;
432
433}
434
435void Map_et::set_beta(double beta0, int l) {
436
437 assert(l>=0) ;
438 assert(l<mg->get_nzone()) ;
439
440 beta[l] = beta0 ;
441
442 reset_coord() ;
443
444}
445
446// ---------------------
447// Constructor from file
448// ---------------------
449Map_et::Map_et(const Mg3d& mgi, FILE* fich)
450 : Map_radial(mgi, fich),
451 aasx( mgi.get_nr(0) ),
452 aasx2( mgi.get_nr(0) ),
453 zaasx( mgi.get_nr(mgi.get_nzone()-1) ),
454 zaasx2( mgi.get_nr(mgi.get_nzone()-1) ),
455 bbsx( mgi.get_nr(0) ),
456 bbsx2( mgi.get_nr(0) ),
457 ff(*(mgi.get_angu()), fich) ,
458 gg(*(mgi.get_angu()), fich)
459{
460 // The Coord rsxdxdr and rsx2drdx are constructed by the default Coord
461 // constructor
462
463 // alpha and beta
464 // --------------
465 int nz = mg->get_nzone() ;
466 alpha = new double[nz] ;
467 beta = new double[nz] ;
468 fread_be(alpha, sizeof(double), nz, fich) ;
469 fread_be(beta, sizeof(double), nz, fich) ;
470
471 // Assignement of the building functions of the Coord's
472 // ----------------------------------------------------
473 set_coord() ;
474
475 // Radial polynomials A(x) and B(x)
476 // --------------------------------
477
478 fait_poly() ;
479
480}
481
482 //------------//
483 // Destructor //
484 //------------//
485
487
488 delete [] alpha ;
489 delete [] beta ;
490
491 for (int l=0 ; l<mg->get_nzone(); l++) {
492 delete aa[l] ;
493 delete daa[l] ;
494 delete ddaa[l] ;
495 delete bb[l] ;
496 delete dbb[l] ;
497 delete ddbb[l] ;
498 }
499 delete [] aa ;
500 delete [] daa ;
501 delete [] ddaa ;
502 delete [] bb ;
503 delete [] dbb ;
504 delete [] ddbb ;
505
506}
507
508
509 //------------//
510 // Assignment //
511 //------------//
512
513void Map_et::operator=(const Map_et& mpi) {
514
515 assert(mpi.get_mg() == mg) ;
516
517 set_ori( mpi.get_ori_x(), mpi.get_ori_y(), mpi.get_ori_z() ) ;
518
519 set_rot_phi( mpi.get_rot_phi() ) ;
520
521 // The members bvect_spher and bvect_cart are treated by the functions
522 // set_ori and set_rot_phi.
523
524 for (int l=0; l<mg->get_nzone(); l++){
525 alpha[l] = mpi.get_alpha()[l] ;
526 beta[l] = mpi.get_beta()[l] ;
527 }
528
529 ff = mpi.ff ;
530 gg = mpi.gg ;
531
532 reset_coord() ; // update of all the Coords
533
534}
535
536
537
538
539void Map_et::operator=(const Map_af& mpi) {
540
541 assert(mpi.get_mg() == mg) ;
542
543 set_ori( mpi.get_ori_x(), mpi.get_ori_y(), mpi.get_ori_z() ) ;
544
545 set_rot_phi( mpi.get_rot_phi() ) ;
546
547 // The members bvect_spher and bvect_cart are treated by the functions
548 // set_ori and set_rot_phi.
549
550 for (int l=0; l<mg->get_nzone(); l++){
551 alpha[l] = mpi.get_alpha()[l] ;
552 beta[l] = mpi.get_beta()[l] ;
553 }
554
555 ff = 0 ;
556 gg = 0 ;
557
558 reset_coord() ; // update of all the Coords
559
560}
561
562void Map_et::set_ff(const Valeur& ffi) {
563
564 ff = ffi ;
565
566 reset_coord() ; // update of all the Coords
567
568}
569
570void Map_et::set_gg(const Valeur& ggi) {
571
572 gg = ggi ;
573
574 reset_coord() ; // update of all the Coords
575
576}
577
578
579
580 //-------------------------------------------------//
581 // Assignment of the Coord building functions //
582 //-------------------------------------------------//
583
585
586 // ... Coord's introduced by the base class Map :
587 r.set(this, map_et_fait_r) ;
588 tet.set(this, map_et_fait_tet) ;
589 phi.set(this, map_et_fait_phi) ;
590 sint.set(this, map_et_fait_sint) ;
591 cost.set(this, map_et_fait_cost) ;
592 sinp.set(this, map_et_fait_sinp) ;
593 cosp.set(this, map_et_fait_cosp) ;
594
595 x.set(this, map_et_fait_x) ;
596 y.set(this, map_et_fait_y) ;
597 z.set(this, map_et_fait_z) ;
598
599 xa.set(this, map_et_fait_xa) ;
600 ya.set(this, map_et_fait_ya) ;
601 za.set(this, map_et_fait_za) ;
602
603 // ... Coord's introduced by the base class Map_radial :
604 xsr.set(this, map_et_fait_xsr) ;
605 dxdr.set(this, map_et_fait_dxdr) ;
606 drdt.set(this, map_et_fait_drdt) ;
607 stdrdp.set(this, map_et_fait_stdrdp) ;
608 srdrdt.set(this, map_et_fait_srdrdt) ;
609 srstdrdp.set(this, map_et_fait_srstdrdp) ;
610 sr2drdt.set(this, map_et_fait_sr2drdt) ;
611 sr2stdrdp.set(this, map_et_fait_sr2stdrdp) ;
612 d2rdx2.set(this, map_et_fait_d2rdx2) ;
613 lapr_tp.set(this, map_et_fait_lapr_tp) ;
614 d2rdtdx.set(this, map_et_fait_d2rdtdx) ;
615 sstd2rdpdx.set(this, map_et_fait_sstd2rdpdx) ;
616 sr2d2rdt2.set(this, map_et_fait_sr2d2rdt2) ;
617
618 //... Coord's which belong to the class Map_et only :
619 rsxdxdr.set(this, map_et_fait_rsxdxdr) ;
620 rsx2drdx.set(this, map_et_fait_rsx2drdx) ;
621
622}
623
624 //--------------------------//
625 // Reset of the Coord's //
626 //--------------------------//
627
629
630 // Coord's of all the class derived from Map_radial:
631
633
634 // Coord's specific to Map_et
635
636 rsxdxdr.del_t() ;
637 rsx2drdx.del_t() ;
638
639}
640
641 //------------------------------------------------------//
642 // Construction of the radial polynomials A(x) and B(x) //
643 //------------------------------------------------------//
644
646
647 int nzone = mg->get_nzone() ;
648
649 aa = new Tbl*[nzone] ;
650 daa = new Tbl*[nzone] ;
651 ddaa = new Tbl*[nzone] ;
652 bb = new Tbl*[nzone] ;
653 dbb = new Tbl*[nzone] ;
654 ddbb = new Tbl*[nzone] ;
655
656 for (int l=0; l<nzone; l++) {
657 int nr = mg->get_nr(l) ;
658 aa[l] = new Tbl(nr) ;
659 daa[l] = new Tbl(nr) ;
660 ddaa[l] = new Tbl(nr) ;
661 bb[l] = new Tbl(nr) ;
662 dbb[l] = new Tbl(nr) ;
663 ddbb[l] = new Tbl(nr) ;
664 }
665
666 // Values in the nucleus
667 // ---------------------
668 assert( mg->get_type_r(0) == RARE || mg->get_type_r(0) == FIN ) ;
669
670 aa[0]->set_etat_qcq() ; // Memory allocation for the Tbl
671 daa[0]->set_etat_qcq() ;
672 ddaa[0]->set_etat_qcq() ;
673 aasx.set_etat_qcq() ;
675
676 bb[0]->set_etat_qcq() ;
677 dbb[0]->set_etat_qcq() ;
678 ddbb[0]->set_etat_qcq() ;
679 bbsx.set_etat_qcq() ;
681
682 for (int i=0; i<mg->get_nr(0); i++) {
683
684 double x1 = (mg->get_grille3d(0))->x[i] ;
685 double x2 = x1 * x1 ;
686 double x3 = x1 * x2 ;
687
688 //##...... A(x) = 2 x^2 - x^4 :
689 // (aa[0])->t[i] = x2 * (2. - x2) ;
690 // (daa[0])->t[i] = 4. * x * (1. + x) * (1. - x) ;
691 // (ddaa[0])->t[i] = 4. *(1. - 3.* x2) ;
692 // aasx->t[i] = x * (2. - x2) ;
693 // aasx2->t[i] = 2. - x2 ;
694
695 //...... A(x) = 3 x^4 - 2 x^6 :
696
697 aa[0]->set(i) = x2 * x2 * (3. - 2.*x2) ;
698 daa[0]->set(i) = 12. * x3 * (1. + x1) * (1. - x1) ;
699 ddaa[0]->set(i) = 12. *x2 *(3. - 5.* x2) ;
700 aasx.set(i) = x3 * (3. - 2.*x2) ;
701 aasx2.set(i) = x2 * (3. - 2.*x2) ;
702
703 //... B(x) = 5/2 x^3 - 3/2 x^5 :
704
705 bb[0]->set(i) = 0.5 * x3 * (5. - 3.* x2) ;
706 dbb[0]->set(i) = 7.5 * x2 * (1. + x1) * (1. - x1) ;
707 ddbb[0]->set(i) = 15. * x1 * ( 1. - 2.*x2 ) ;
708 bbsx.set(i) = 0.5 * x2 * (5. - 3.* x2) ;
709 bbsx2.set(i) = 0.5 * x1 * (5. - 3.* x2) ;
710 }
711
712 // Values in the shells and the outermost domain
713 // ---------------------------------------------
714
715 for (int l=1; l<nzone; l++) {
716
717 assert( (mg->get_type_r(l) == FIN)|| (mg->get_type_r(l) == UNSURR) ) ;
718
719 aa[l]->set_etat_qcq() ; // Memory allocation for the Tbl
720 daa[l]->set_etat_qcq() ;
721 ddaa[l]->set_etat_qcq() ;
722
723 bb[l]->set_etat_qcq() ;
724 dbb[l]->set_etat_qcq() ;
725 ddbb[l]->set_etat_qcq() ;
726
727 for (int i=0; i<mg->get_nr(l); i++) {
728
729 double x1 = (mg->get_grille3d(l))->x[i] ;
730 double xm1 = x1 - 1. ;
731 double xp1 = x1 + 1. ;
732
733 //... A(x) = 1/4 (x-1)^2 (x+2) = 1/4(x^3 -3x +2) :
734
735 aa[l]->set(i) = 0.25* xm1 * xm1 * (x1 + 2.) ;
736 daa[l]->set(i) = 0.75* xm1 * xp1 ;
737 ddaa[l]->set(i) = 1.5* x1 ;
738
739 //... B(x) = 1/4 (x+1)^2 (-x+2) = 1/4(-x^3 +3x +2) :
740
741 bb[l]->set(i) = 0.25* xp1 * xp1 * (2. - x1) ;
742 dbb[l]->set(i) = - 0.75* xm1 * xp1 ;
743 ddbb[l]->set(i) = - 1.5* x1 ;
744 }
745
746 } // End of the loop onto the domains
747
748 // Special case of a compactified outermost domain
749 // -----------------------------------------------
750
751 int nzm1 = nzone - 1 ;
752 if ( mg->get_type_r(nzm1) == UNSURR ) {
753
754 zaasx.set_etat_qcq() ; // Memory allocation for the Tbl
756
757 for (int i=0; i<mg->get_nr(nzm1); i++) {
758
759 double x1 = (mg->get_grille3d(nzm1))->x[i] ;
760 zaasx.set(i) = 0.25 * (x1 - 1.) * (2. + x1) ; // A(x)/(x-1)
761 zaasx2.set(i) = 0.25 * (2. + x1) ; // A(x)/(x-1)^2
762
763 }
764
765 bb[nzm1]->set_etat_zero() ;
766 dbb[nzm1]->set_etat_zero() ;
767 ddbb[nzm1]->set_etat_zero() ;
768
769 }
770
771}
772
773
774
775 //----------------//
776 // Save in a file //
777 //----------------//
778
779void Map_et::sauve(FILE* fich) const {
780
781 Map_radial::sauve(fich) ; // Write of the elements common to all the
782 // classes derived from Map_radial
783
784 ff.sauve(fich) ; // Write of F(theta',phi')
785 gg.sauve(fich) ; // Write of G(theta',phi')
786
787 // Write of alpha and beta :
788 int nz = mg->get_nzone() ;
789 fwrite_be(alpha, sizeof(double), nz, fich) ;
790 fwrite_be(beta, sizeof(double), nz, fich) ;
791
792}
793
794 //---------------------------//
795 // Printing //
796 //---------------------------//
797
798ostream & Map_et::operator>>(ostream & ost) const {
799
800 using namespace Unites ;
801
802 ost <<
803 "Radial mapping of form r = xi + A(xi)F(t,p) + B(xi)G(t,p) (class Map_et)"
804 << endl ;
805 int nz = mg->get_nzone() ;
806 for (int l=0; l<nz; l++) {
807 ost << " Domain #" << l << " : alpha_l = " << alpha[l]
808 << " , beta_l = " << beta[l] << endl ;
809 }
810 ost << endl << "Function F(theta', phi') : " << endl ;
811 ost << "------------------------- " << endl ;
812 ff.affiche_seuil(ost) ;
813 ost << endl <<"Function G(theta', phi') : " << endl ;
814 ost << "------------------------- " << endl ;
815 gg.affiche_seuil(ost) ;
816
817 int type_t = mg->get_type_t() ;
818 int type_p = mg->get_type_p() ;
819
820 ost << endl
821 << "Values of r at the outer boundary of each domain [km] :" << endl ;
822 ost << "------------------------------------------------------" << endl ;
823 ost << " 1/ for theta = Pi/2 and phi = 0 : " << endl ;
824 ost << " val_r : " ;
825 for (int l=0; l<nz; l++) {
826 ost << " " << val_r(l, 1., M_PI/2, 0) / km ;
827 }
828 ost << endl ;
829
830 if ( type_t == SYM ) {
831 assert( (type_p == SYM) || (type_p == NONSYM) ) ;
832 ost << " Coord r : " ;
833 for (int l=0; l<nz; l++) {
834 int nrm1 = mg->get_nr(l) - 1 ;
835 int ntm1 = mg->get_nt(l) - 1 ;
836 ost << " " << (+r)(l, 0, ntm1, nrm1) / km ;
837 }
838 ost << endl ;
839 }
840
841 ost << " 2/ for theta = Pi/2 and phi = Pi/2 : " << endl ;
842 ost << " val_r : " ;
843 for (int l=0; l<nz; l++) {
844 ost << " " << val_r(l, 1., M_PI/2, M_PI/2) / km ;
845 }
846 ost << endl ;
847
848 if ( type_t == SYM ) {
849 ost << " Coord r : " ;
850 for (int l=0; l<nz; l++) {
851 int nrm1 = mg->get_nr(l) - 1 ;
852 int ntm1 = mg->get_nt(l) - 1 ;
853 int np = mg->get_np(l) ;
854 if ( (type_p == NONSYM) && (np % 4 == 0) ) {
855 ost << " " << (+r)(l, np/4, ntm1, nrm1) / km ;
856 }
857 if ( type_p == SYM ) {
858 ost << " " << (+r)(l, np/2, ntm1, nrm1) / km ;
859 }
860 }
861 ost << endl ;
862 }
863
864 ost << " 3/ for theta = Pi/2 and phi = Pi : " << endl ;
865 ost << " val_r : " ;
866 for (int l=0; l<nz; l++) {
867 ost << " " << val_r(l, 1., M_PI/2, M_PI) / km ;
868 }
869 ost << endl ;
870
871 if ( (type_t == SYM) && (type_p == NONSYM) ) {
872 ost << " Coord r : " ;
873 for (int l=0; l<nz; l++) {
874 int nrm1 = mg->get_nr(l) - 1 ;
875 int ntm1 = mg->get_nt(l) - 1 ;
876 int np = mg->get_np(l) ;
877 ost << " " << (+r)(l, np/2, ntm1, nrm1) / km ;
878 }
879 ost << endl ;
880 }
881
882 ost << " 4/ for theta = 0 : " << endl ;
883 ost << " val_r : " ;
884 for (int l=0; l<nz; l++) {
885 ost << " " << val_r(l, 1., 0., 0.) / km ;
886 }
887 ost << endl ;
888
889 ost << " Coord r : " ;
890 for (int l=0; l<nz; l++) {
891 int nrm1 = mg->get_nr(l) - 1 ;
892 ost << " " << (+r)(l, 0, 0, nrm1) / km ;
893 }
894 ost << endl ;
895
896 return ost ;
897
898}
899
900 //------------------//
901 // Homothetie //
902 //------------------//
903
904
905void Map_et::homothetie(double fact) {
906
907 int nz = mg->get_nzone() ;
908
909 for (int l=0; l<nz; l++) {
910 if (mg->get_type_r(l) == UNSURR) {
911 alpha[l] /= fact ;
912 beta[l] /= fact ;
913 }
914 else {
915 alpha[l] *= fact ;
916 beta[l] *= fact ;
917 }
918 }
919
920 reset_coord() ;
921
922}
923
924 //----------------------------//
925 // Rescaling of one domain //
926 //----------------------------//
927
928void Map_et::resize(int l, double lambda) {
929
930 // Protections
931 // -----------
932 if (mg->get_type_r(l) != FIN) {
933 cout << "Map_et::resize can be applied only to a shell !" << endl ;
934 abort() ;
935 }
936
937 // New values of alpha, beta, F and G in domain l :
938 // ----------------------------------------------
939 double n_alpha = 0.5 * ( (lambda + 1.) * alpha[l] +
940 (lambda - 1.) * beta[l] ) ;
941
942 double n_beta = 0.5 * ( (lambda - 1.) * alpha[l] +
943 (lambda + 1.) * beta[l] ) ;
944
945 ff.set(l) = alpha[l] / n_alpha * ff(l) ;
946 gg.set(l) = lambda * alpha[l] / n_alpha * gg(l) ;
947
948 alpha[l] = n_alpha ;
949 beta[l] = n_beta ;
950
951 // New values of alpha, beta, F and G in in domain l+1 :
952 // ----------------------------------------------------
953 assert(l<mg->get_nzone()-1) ;
954 int lp1 = l + 1 ;
955
956 if (mg->get_type_r(lp1) == UNSURR) { // compactified case
957
958 assert(ff(lp1).get_etat() == ETATZERO ) ;
959 assert(gg(lp1).get_etat() == ETATZERO ) ;
960
961 alpha[lp1] = - 0.5 / ( alpha[l] + beta[l] ) ;
962 beta[lp1] = - alpha[lp1] ;
963
964 }
965 else{ // non-compactified case
966
967 assert( mg->get_type_r(lp1) == FIN ) ;
968 n_alpha = 0.5 * ( alpha[lp1] - alpha[l] + beta[lp1] - beta[l] ) ;
969 n_beta = 0.5 * ( alpha[lp1] + alpha[l] + beta[lp1] + beta[l] ) ;
970
971 ff.set(lp1) = alpha[l] / n_alpha * gg(l) ;
972 gg.set(lp1) = alpha[lp1] / n_alpha * gg(lp1) ;
973
974 alpha[lp1] = n_alpha ;
975 beta[lp1] = n_beta ;
976 }
977
978 // The coords are no longer up to date :
979 reset_coord() ;
980}
981
982
983// Comparison operator :
984bool Map_et::operator==(const Map& mpi) const {
985
986 // Precision of the comparison
987 double precis = 1e-10 ;
988 bool resu = true ;
989
990 // Dynamic cast pour etre sur meme Map...
991 const Map_et* mp0 = dynamic_cast<const Map_et*>(&mpi) ;
992 if (mp0 == 0x0)
993 resu = false ;
994 else {
995 if (*mg != *(mpi.get_mg()))
996 resu = false ;
997
998 if (fabs(ori_x-mpi.get_ori_x()) > precis) resu = false ;
999 if (fabs(ori_y-mpi.get_ori_y()) > precis) resu = false ;
1000 if (fabs(ori_z-mpi.get_ori_z()) > precis) resu = false ;
1001
1002 if (bvect_spher != mpi.get_bvect_spher()) resu = false ;
1003 if (bvect_cart != mpi.get_bvect_cart()) resu = false ;
1004
1005 int nz = mg->get_nzone() ;
1006 for (int i=0 ; i<nz ; i++) {
1007 if (fabs(alpha[i]-mp0->alpha[i])/fabs(alpha[i]) > precis)
1008 resu = false ;
1009 if ((i!=0) && (i!=nz-1))
1010 if (fabs(beta[i]-mp0->beta[i])/fabs(beta[i]) > precis)
1011 resu = false ;
1012 }
1013
1014 if (max(diffrelmax(ff, mp0->ff)) > precis)
1015 resu = false ;
1016 if (max(diffrelmax(gg, mp0->gg)) > precis)
1017 resu = false ;
1018 }
1019
1020 return resu ;
1021}
1022 //--------------------------------------//
1023 // Extraction of the mapping parameters //
1024 //--------------------------------------//
1025
1026const double* Map_et::get_alpha() const {
1027 return alpha ;
1028}
1029
1030const double* Map_et::get_beta() const {
1031 return beta ;
1032}
1033
1034const Valeur& Map_et::get_ff() const {
1035 return ff ;
1036}
1037
1038const Valeur& Map_et::get_gg() const {
1039 return gg ;
1040}
1041
1042
1043// To be done
1044//-----------
1045const Map_af& Map_et::mp_angu(int) const {
1046 const char* f = __FILE__ ;
1047 c_est_pas_fait(f) ;
1048 p_mp_angu = new Map_af(*this) ;
1049 return *p_mp_angu ;
1050}
1051
1052}
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
void del_t() const
Logical destructor (deletes the Mtbl member *c ).
Definition coord.C:125
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition coord.C:134
Affine radial mapping.
Definition map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Radial mapping of rather general form.
Definition map.h:2752
const Valeur & get_ff() const
Returns a (constant) reference to the function .
Definition map_et.C:1034
const double * get_alpha() const
Returns a pointer on the array alpha (values of in each domain)
Definition map_et.C:1026
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2760
Tbl ** ddaa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2775
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
virtual void operator=(const Map_et &mp)
Assignment to another Map_et
Definition map_et.C:513
virtual const Map_af & mp_angu(int) const
Returns the "angular" mapping for the outside of domain l_zone.
Definition map_et.C:1045
virtual ostream & operator>>(ostream &) const
Operator >>
Definition map_et.C:798
Tbl ** ddbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2806
Tbl zaasx2
Values at the nr collocation points of in the outermost compactified domain.
Definition map.h:2791
Coord rsx2drdx
in the nucleus and the shells; \ in the outermost compactified domain.
Definition map.h:2841
virtual bool operator==(const Map &) const
Comparison operator (egality)
Definition map_et.C:984
Tbl bbsx
Values at the nr collocation points of in the nucleus.
Definition map.h:2809
void set_gg(const Valeur &)
Assigns a given value to the function .
Definition map_et.C:570
const Valeur & get_gg() const
Returns a (constant) reference to the function .
Definition map_et.C:1038
Map_et(const Mg3d &mgrille, const double *r_limits)
Standard Constructor.
Definition map_et.C:149
Tbl ** dbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2801
void fait_poly()
Construction of the polynomials and .
Definition map_et.C:645
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
Tbl aasx
Values at the nr collocation points of in the nucleus.
Definition map.h:2778
const double * get_beta() const
Returns a pointer on the array beta (values of in each domain)
Definition map_et.C:1030
Tbl aasx2
Values at the nr collocation points of in the nucleus.
Definition map.h:2781
Tbl zaasx
Values at the nr collocation points of in the outermost compactified domain.
Definition map.h:2786
Tbl ** daa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2770
void set_coord()
Assignement of the building functions to the member Coords.
Definition map_et.C:584
Coord rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
Definition map.h:2834
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition map_et.C:435
Tbl ** aa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2765
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2819
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition map_et.C:424
Tbl ** bb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Definition map.h:2796
virtual ~Map_et()
Destructor.
Definition map_et.C:486
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2758
Tbl bbsx2
Values at the nr collocation points of in the nucleus.
Definition map.h:2812
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2826
virtual void sauve(FILE *) const
Save in a file.
Definition map_et.C:779
void set_ff(const Valeur &)
Assigns a given value to the function .
Definition map_et.C:562
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Definition map_et.C:928
virtual void reset_coord()
Resets all the member Coords.
Definition map_et.C:628
Base class for pure radial mappings.
Definition map.h:1536
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1619
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1600
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1592
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1640
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1648
virtual void reset_coord()
Resets all the member Coords.
Definition map_radial.C:126
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1631
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1608
Coord drdt
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
Definition map.h:1568
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1584
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
virtual void sauve(FILE *) const
Save in a file.
Definition map_radial.C:116
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1560
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1657
Coord stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
Definition map.h:1576
Base class for coordinate mappings.
Definition map.h:670
double get_ori_z() const
Returns the z coordinate of the origin.
Definition map.h:772
Coord cosp
Definition map.h:724
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord sint
Definition map.h:721
double ori_x
Absolute coordinate x of the origin.
Definition map.h:678
Coord ya
Absolute y coordinate.
Definition map.h:731
void set_rot_phi(double phi0)
Sets a new rotation angle.
Definition map.C:263
double ori_y
Absolute coordinate y of the origin.
Definition map.h:679
Coord r
r coordinate centered on the grid
Definition map.h:718
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:689
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Definition map.C:253
double get_ori_y() const
Returns the y coordinate of the origin.
Definition map.h:770
Coord za
Absolute z coordinate.
Definition map.h:732
Coord tet
coordinate centered on the grid
Definition map.h:719
Coord sinp
Definition map.h:723
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Base_vect_cart bvect_cart
Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition map.h:697
Coord z
z coordinate centered on the grid
Definition map.h:728
double get_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
Map_af * p_mp_angu
Pointer on the "angular" mapping.
Definition map.h:715
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord phi
coordinate centered on the grid
Definition map.h:720
double ori_z
Absolute coordinate z of the origin.
Definition map.h:680
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
Coord cost
Definition map.h:722
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
Multi-domain grid.
Definition grilles.h:273
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition grilles.h:485
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
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_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
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
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
Basic array class.
Definition tbl.h:161
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition tbl.h:400
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition valeur.C:689
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void affiche_seuil(ostream &ostr, int type=0, int precision=4, double threshold=1.e-7) const
Prints only the values greater than a given threshold.
Definition valeur.C:558
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
void sauve(FILE *) const
Save in a file.
Definition valeur.C:475
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
#define P_COSSIN
dev. standart
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.