LORENE
map_af.C
1 /*
2  * Methods of class Map_af
3  *
4  * (see file map.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 1999-2003 Eric Gourgoulhon
10  * Copyright (c) 1999-2001 Philippe Grandclement
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: map_af.C,v 1.26 2024/06/24 14:11:43 j_novak Exp $
35  * $Log: map_af.C,v $
36  * Revision 1.26 2024/06/24 14:11:43 j_novak
37  * Removed parameter names to avoid compilation warnings.
38  *
39  * Revision 1.25 2023/05/26 15:39:31 g_servignat
40  * Added c_est_pas_fait() to poisson_angu(Cmp)
41  *
42  * Revision 1.24 2023/01/04 10:23:37 j_novak
43  * Removed spurious output.
44  *
45  * Revision 1.23 2022/09/30 14:44:03 j_novak
46  * Improved treatment of outer boundary in constructor from a file
47  *
48  * Revision 1.22 2020/01/27 11:00:19 j_novak
49  * New include <stdexcept> to be compatible with older versions of g++
50  *
51  * Revision 1.21 2018/12/05 15:43:45 j_novak
52  * New Map_af constructor from a formatted file.
53  *
54  * Revision 1.20 2016/12/05 16:17:56 j_novak
55  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
56  *
57  * Revision 1.19 2014/10/13 08:53:02 j_novak
58  * Lorene classes and functions now belong to the namespace Lorene.
59  *
60  * Revision 1.18 2014/10/06 15:13:11 j_novak
61  * Modified #include directives to use c++ syntax.
62  *
63  * Revision 1.17 2013/06/05 15:10:42 j_novak
64  * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
65  * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
66  *
67  * Revision 1.16 2012/01/17 15:34:35 j_penner
68  * *** empty log message ***
69  *
70  * Revision 1.15 2012/01/17 10:31:32 j_penner
71  * added access to computational coordinate xi
72  *
73  * Revision 1.14 2008/09/29 13:23:51 j_novak
74  * Implementation of the angular mapping associated with an affine
75  * mapping. Things must be improved to take into account the domain index.
76  *
77  * Revision 1.13 2008/08/19 06:42:00 j_novak
78  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
79  * cast-type operations, and constant strings that must be defined as const char*
80  *
81  * Revision 1.12 2007/12/21 11:05:33 j_novak
82  * Put back the constructor from a Mg3d and a Tbl (it had disappeared?).
83  *
84  * Revision 1.11 2007/12/20 09:11:04 jl_cornou
85  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
86  *
87  * Revision 1.10 2007/12/11 15:28:14 jl_cornou
88  * Jacobi(0,2) polynomials partially implemented
89  *
90  * Revision 1.9 2006/04/25 07:21:59 p_grandclement
91  * Various changes for the NS_BH project
92  *
93  * Revision 1.8 2005/08/29 15:10:18 p_grandclement
94  * Addition of things needed :
95  * 1) For BBH with different masses
96  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
97  * WORKING YET !!!
98  *
99  * Revision 1.7 2004/12/02 09:33:06 p_grandclement
100  * *** empty log message ***
101  *
102  * Revision 1.6 2004/03/25 10:29:23 j_novak
103  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
104  *
105  * Revision 1.5 2004/01/29 08:50:03 p_grandclement
106  * Modification of Map::operator==(const Map&) and addition of the surface
107  * integrales using Scalar.
108  *
109  * Revision 1.4 2003/10/15 10:33:11 e_gourgoulhon
110  * Added new Coord's : drdt, srdrdp.
111  *
112  * Revision 1.3 2002/10/16 14:36:41 j_novak
113  * Reorganization of #include instructions of standard C++, in order to
114  * use experimental version 3 of gcc.
115  *
116  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
117  *
118  * All writing/reading to a binary file are now performed according to
119  * the big endian convention, whatever the system is big endian or
120  * small endian, thanks to the functions fwrite_be and fread_be
121  *
122  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
123  * LORENE
124  *
125  * Revision 2.23 2001/02/28 11:04:05 eric
126  * 1ere version testee de resize.
127  *
128  * Revision 2.22 2001/02/26 17:29:22 eric
129  * Ajout de la fonction resize.
130  *
131  * Revision 2.21 2001/01/10 11:03:13 phil
132  * ajout de homothetie interne
133  *
134  * Revision 2.20 2000/01/24 17:09:13 eric
135  * suppression de la fonction convert.
136  * suppression du constructeur par convertion d'un Map_et.
137  * ajout du constructeur par conversion d'un Map.
138  *
139  * Revision 2.19 2000/01/24 16:42:48 eric
140  * Dans operator=(const Map_af& ), appel de set_rot_phi.
141  * Ajout de la fonction convert(const Map& ).
142  *
143  * Revision 2.18 1999/12/21 16:27:26 eric
144  * Ajout du constructeur par conversion Map_af::Map_af(const Map_et&).
145  * Ajout des fonctions Map_af::set_alpha et Map_af::set_beta.
146  *
147  * Revision 2.17 1999/12/17 09:14:43 eric
148  * Amelioration de l'affichage.
149  *
150  * Revision 2.16 1999/12/06 13:12:23 eric
151  * Introduction de la fonction homothetie.
152  *
153  * Revision 2.15 1999/11/25 16:28:53 eric
154  * Le calcul des derivees partielles est transfere dans le fichier
155  * map_af_deriv.C.
156  *
157  * Revision 2.14 1999/11/24 14:32:54 eric
158  * Les prototypes des fonctions de constructions des coords sont desormais
159  * dans map.h.
160  * Introduction des fonctions get_alpha() et get_beta().
161  * /
162  *
163  * Revision 2.13 1999/11/22 10:36:14 eric
164  * Introduction de la fonction set_coord().
165  * Fonction del_coord() rebaptisee reset_coord().
166  *
167  * Revision 2.12 1999/10/27 08:46:29 eric
168  * Introduction de Cmp::va a la place de *(Cmp::c).
169  *
170  * Revision 2.11 1999/10/15 09:23:10 eric
171  * *** empty log message ***
172  *
173  * Revision 2.10 1999/10/15 09:16:23 eric
174  * Changement prototypes: const.
175  *
176  * Revision 2.9 1999/10/14 14:27:05 eric
177  * Depoussierage.
178  *
179  * Revision 2.8 1999/04/12 12:54:05 phil
180  * *** empty log message ***
181  *
182  * Revision 2.7 1999/04/12 12:09:03 phil
183  * Mise a jour des bases dans dsdr()
184  *
185  * Revision 2.6 1999/03/04 13:11:48 eric
186  * Ajout des Coord representant les derivees du changement de variable.
187  *
188  * Revision 2.5 1999/03/03 11:19:08 hyc
189  * *** empty log message ***
190  *
191  *
192  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af.C,v 1.26 2024/06/24 14:11:43 j_novak Exp $
193  *
194  */
195 
196 // headers C++
197 #include <stdexcept>
198 
199 // headers C
200 #include <cmath>
201 
202 // headers Lorene
203 #include "cmp.h"
204 #include "nbr_spx.h"
205 #include "utilitaires.h"
206 #include "proto.h"
207 #include "unites.h"
208 
209  //---------------//
210  // Constructeurs //
211  //---------------//
212 
213 // Constructor from a grid
214 // -----------------------
215 namespace Lorene {
216 Map_af::Map_af(const Mg3d& mgrille, const double* bornes) : Map_radial(mgrille)
217 {
218  // Les coordonnees et les derivees du changement de variable
219  set_coord() ;
220 
221  // Les bornes
222  int nzone = mg->get_nzone() ;
223 
224  alpha = new double[nzone] ;
225  beta = new double[nzone] ;
226 
227  for (int l=0 ; l<nzone ; l++) {
228  switch (mg->get_type_r(l)) {
229  case RARE: {
230  alpha[l] = bornes[l+1] - bornes[l] ;
231  beta[l] = bornes[l] ;
232  break ;
233  }
234 
235  case FIN: {
236  alpha[l] = (bornes[l+1] - bornes[l]) * .5 ;
237  beta[l] = (bornes[l+1] + bornes[l]) * .5 ;
238  break ;
239  }
240 
241  case UNSURR: {
242  double umax = 1./bornes[l] ;
243  double umin = 1./bornes[l+1] ;
244  alpha[l] = (umin - umax) * .5 ; // u est une fonction decroissante
245  beta[l] = (umin + umax) * .5 ; // de l'indice i en r
246  break ;
247  }
248 
249  default: {
250  cout << "Map_af::Map_af: unkown type_r ! " << endl ;
251  abort () ;
252  break ;
253  }
254 
255  }
256  } // Fin de la boucle sur zone
257 }
258 
259 // Constructor from a grid
260 // -----------------------
261 Map_af::Map_af(const Mg3d& mgrille, const Tbl& bornes) : Map_radial(mgrille)
262 {
263  // Les coordonnees et les derivees du changement de variable
264  set_coord() ;
265 
266  // Les bornes
267  int nzone = mg->get_nzone() ;
268 
269  alpha = new double[nzone] ;
270  beta = new double[nzone] ;
271 
272  for (int l=0 ; l<nzone ; l++) {
273  switch (mg->get_type_r(l)) {
274  case RARE: {
275  alpha[l] = bornes(l+1) - bornes(l) ;
276  beta[l] = bornes(l) ;
277  break ;
278  }
279 
280  case FIN: {
281  alpha[l] = (bornes(l+1) - bornes(l)) * .5 ;
282  beta[l] = (bornes(l+1) + bornes(l)) * .5 ;
283  break ;
284  }
285 
286  case UNSURR: {
287  assert (l==nzone-1) ;
288  double umax = 1./bornes(l) ;
289  double umin = 0 ;
290  alpha[l] = (umin - umax) * .5 ; // u est une fonction decroissante
291  beta[l] = (umin + umax) * .5 ; // de l'indice i en r
292  break ;
293  }
294 
295  default: {
296  cout << "Map_af::Map_af: unkown type_r ! " << endl ;
297  abort () ;
298  break ;
299  }
300 
301  }
302  } // Fin de la boucle sur zone
303 }
304 
305 // Copy constructor
306 // ----------------
308 {
309  // Les coordonnees et les derivees du changement de variable
310  set_coord() ;
311 
312  // Les bornes
313  int nzone = mg->get_nzone() ;
314 
315  alpha = new double[nzone] ;
316  beta = new double[nzone] ;
317 
318  for (int l=0; l<nzone; l++){
319  alpha[l] = mp.alpha[l] ;
320  beta[l] = mp.beta[l] ;
321  }
322 }
323 
324  // Constructor from a formatted file
325  //-----------------------------------
326  Map_af::Map_af(const Mg3d& mgrille, const string& filename) : Map_radial(mgrille)
327 {
328  // Opening of the file
329  //---------------------
330  ifstream parfile(filename.c_str()) ;
331  if (!parfile) {
332  string message = "Unable to open file " ;
333  message += filename ;
334  throw ios_base::failure(message);
335  };
336 
337  string tmp_str = "Definition of the Map_af" ;
338  if (!search_file(parfile, tmp_str)) {
339  string mesg = "No data found to contruct an object Map_af in " ;
340  mesg += filename ;
341  throw invalid_argument(mesg) ;
342  }
343  parfile.ignore(1000, '\n') ;
344 
345  // Number of domains
346  //---------------------
347  int nz ;
348  parfile >> nz ; parfile.ignore(1000, '\n') ; // Number of domains
349  if (mg->get_nzone() != nz) {
350  string mesg = "Wrong number of domains for Map_af in " ;
351  mesg += filename ;
352  throw(invalid_argument(mesg)) ;
353  }
354  Tbl bornes(nz+1) ;
355  bornes.set_etat_qcq() ;
356  for (int i=0; i<nz; i++) {
357  parfile >> bornes.set(i) ;
358  parfile.ignore(1000, '\n') ;
359  }
360 
361  string last_boundary ;
362  parfile >> last_boundary ;
363 
364  if (last_boundary.compare("infinity") == 0) {
365  bornes.set(nz) = __infinity ;
366  if (mgrille.get_type_r(nz-1) != UNSURR) {
367  string mesg = filename
368  + ": infinite outer boundary can be set only in the ced case." ;
369  throw(invalid_argument(mesg)) ;
370  }
371  }
372  else {
373  double rout = 0. ;
374  try {
375  rout = stod(last_boundary) ;
376  }
377  catch(invalid_argument& ia) {
378  cerr << "Map_af constructor from a file, invalid argument in file\n"
379  << last_boundary << endl ;
380  abort() ;
381  }
382  bornes.set(nz) = rout ;
383  if (mgrille.get_type_r(nz-1) == UNSURR) {
384  string mesg = filename
385  + ": the ced grid type cannot accept a finite outer radius." ;
386  throw(invalid_argument(mesg)) ;
387  }
388  }
389 
390  set_coord() ;
391 
392  alpha = new double[nz] ;
393  beta = new double[nz] ;
394 
395  for (int l=0 ; l<nz ; l++) {
396  switch (mg->get_type_r(l)) {
397  case RARE: {
398  alpha[l] = bornes(l+1) - bornes(l) ;
399  beta[l] = bornes(l) ;
400  break ;
401  }
402 
403  case FIN: {
404  alpha[l] = (bornes(l+1) - bornes(l)) * .5 ;
405  beta[l] = (bornes(l+1) + bornes(l)) * .5 ;
406  break ;
407  }
408 
409  case UNSURR: {
410  assert (l==nz-1) ;
411  double umax = 1./bornes(l) ;
412  double umin = 0 ;
413  alpha[l] = (umin - umax) * .5 ; // u est une fonction decroissante
414  beta[l] = (umin + umax) * .5 ; // de l'indice i en r
415  break ;
416  }
417 
418  default: {
419  cout << "Map_af::Map_af: unkown type_r ! " << endl ;
420  abort () ;
421  break ;
422  }
423 
424  }
425  } // Fin de la boucle sur zone
426 }
427 
428 
429 // Constructor from file
430 // ---------------------
431 Map_af::Map_af(const Mg3d& mgi, FILE* fd) : Map_radial(mgi, fd)
432 {
433  int nz = mg->get_nzone() ;
434  alpha = new double[nz] ;
435  beta = new double[nz] ;
436  fread_be(alpha, sizeof(double), nz, fd) ;
437  fread_be(beta, sizeof(double), nz, fd) ;
438 
439  // Les coordonnees et les derivees du changement de variable
440  set_coord() ;
441 }
442 
443 
444 
445 // Constructor from a Map
446 // -----------------------
447 Map_af::Map_af(const Map& mpi) : Map_radial(*(mpi.get_mg()))
448 {
449  // Les coordonnees et les derivees du changement de variable
450  set_coord() ;
451 
452  // Les bornes
453  int nz = mg->get_nzone() ;
454 
455  alpha = new double[nz] ;
456  beta = new double[nz] ;
457 
458  const Map_af* mp0 = dynamic_cast<const Map_af*>(&mpi) ;
459  const Map_et* mp1 = dynamic_cast<const Map_et*>(&mpi) ;
460 
461  if( (mp0 == 0x0) && (mp1 == 0x0) ) {
462  cout << "Map_af::Map_af(const Map& ) : unkown mapping type !"
463  << endl ;
464  abort() ;
465  }
466 
467  if (mp0 != 0x0) {
468  assert( mp1 == 0x0 ) ;
469  for (int l=0; l<nz; l++){
470  alpha[l] = mp0->get_alpha()[l] ;
471  beta[l] = mp0->get_beta()[l] ;
472  }
473  }
474 
475 
476  if (mp1 != 0x0) {
477  assert( mp0 == 0x0 ) ;
478  for (int l=0; l<nz; l++){
479  alpha[l] = mp1->get_alpha()[l] ;
480  beta[l] = mp1->get_beta()[l] ;
481  }
482  }
483 
484 
485  set_ori( mpi.get_ori_x(), mpi.get_ori_y(), mpi.get_ori_z() ) ;
486 
487  set_rot_phi( mpi.get_rot_phi() ) ;
488 
489 }
490 
491 
492 
493 
494  //--------------//
495  // Destructeurs //
496  //--------------//
497 
499  delete [] alpha ;
500  delete [] beta ;
501 }
502 
503  //-------------//
504  // Assignement //
505  //-------------//
506 
507 // From another Map_af
508 // -------------------
509 
510 void Map_af::operator=(const Map_af & mpi) {
511 
512  assert(mpi.mg == mg) ;
513 
514  set_ori( mpi.ori_x, mpi.ori_y, mpi.ori_z ) ;
515 
516  set_rot_phi( mpi.rot_phi ) ;
517 
518  for (int l = 0; l<mg->get_nzone(); l++) {
519  alpha[l] = mpi.alpha[l] ;
520  beta[l] = mpi.beta[l] ;
521  }
522 
523  reset_coord() ;
524 }
525 
526 
527 
528 
529  //-------------------------------------------------//
530  // Assignement of the Coord building functions //
531  //-------------------------------------------------//
532 
534 
535  // ... Coord's introduced by the base class Map :
536  r.set(this, map_af_fait_r) ;
537  tet.set(this, map_af_fait_tet) ;
538  phi.set(this, map_af_fait_phi) ;
539  sint.set(this, map_af_fait_sint) ;
540  cost.set(this, map_af_fait_cost) ;
541  sinp.set(this, map_af_fait_sinp) ;
542  cosp.set(this, map_af_fait_cosp) ;
543 
544  x.set(this, map_af_fait_x) ;
545  y.set(this, map_af_fait_y) ;
546  z.set(this, map_af_fait_z) ;
547 
548  xa.set(this, map_af_fait_xa) ;
549  ya.set(this, map_af_fait_ya) ;
550  za.set(this, map_af_fait_za) ;
551 
552  // ... Coord's introduced by the base class Map_radial :
553  xsr.set(this, map_af_fait_xsr) ;
554  dxdr.set(this, map_af_fait_dxdr) ;
555  drdt.set(this, map_af_fait_drdt) ;
556  stdrdp.set(this, map_af_fait_stdrdp) ;
557  srdrdt.set(this, map_af_fait_srdrdt) ;
558  srstdrdp.set(this, map_af_fait_srstdrdp) ;
559  sr2drdt.set(this, map_af_fait_sr2drdt) ;
560  sr2stdrdp.set(this, map_af_fait_sr2stdrdp) ;
561  d2rdx2.set(this, map_af_fait_d2rdx2) ;
562  lapr_tp.set(this, map_af_fait_lapr_tp) ;
563  d2rdtdx.set(this, map_af_fait_d2rdtdx) ;
564  sstd2rdpdx.set(this, map_af_fait_sstd2rdpdx) ;
565  sr2d2rdt2.set(this, map_af_fait_sr2d2rdt2) ;
566 
567 }
568 // Comparison operator :
569 bool Map_af::operator==(const Map& mpi) const {
570 
571  // Precision of the comparison
572  double precis = 1e-10 ;
573  bool resu = true ;
574 
575  // Dynamic cast pour etre sur meme Map...
576  const Map_af* mp0 = dynamic_cast<const Map_af*>(&mpi) ;
577  if (mp0 == 0x0)
578  resu = false ;
579  else {
580  if (*mg != *(mpi.get_mg()))
581  resu = false ;
582 
583  if (fabs(ori_x-mpi.get_ori_x()) > precis) resu = false ;
584  if (fabs(ori_y-mpi.get_ori_y()) > precis) resu = false ;
585  if (fabs(ori_z-mpi.get_ori_z()) > precis) resu = false ;
586 
587  if (bvect_spher != mpi.get_bvect_spher()) resu = false ;
588  if (bvect_cart != mpi.get_bvect_cart()) resu = false ;
589 
590  int nz = mg->get_nzone() ;
591  for (int i=0 ; i<nz ; i++) {
592  if (fabs(alpha[i]-mp0->alpha[i])/fabs(alpha[i]) > precis)
593  resu = false ;
594  if ((i!=0) && (i!=nz-1))
595  if (fabs(beta[i]-mp0->beta[i])/fabs(beta[i]) > precis)
596  resu = false ;
597  }
598  }
599 
600  return resu ;
601 }
602 
603  //--------------------------------------//
604  // Extraction of the mapping parameters //
605  //--------------------------------------//
606 
607 const double* Map_af::get_alpha() const {
608  return alpha ;
609 }
610 
611 const double* Map_af::get_beta() const {
612  return beta ;
613 }
614 
615  //------------//
616  // Sauvegarde //
617  //------------//
618 
619 void Map_af::sauve(FILE* fd) const {
620 
621  Map_radial::sauve(fd) ;
622 
623  int nz = mg->get_nzone() ;
624  fwrite_be(alpha, sizeof(double), nz, fd) ;
625  fwrite_be(beta, sizeof(double), nz, fd) ;
626 
627 }
628 
629  //------------//
630  // Impression //
631  //------------//
632 
633 ostream & Map_af::operator>>(ostream & ost) const {
634 
635  using namespace Unites ;
636 
637  ost << "Affine mapping (class Map_af)" << endl ;
638  int nz = mg->get_nzone() ;
639  for (int l=0; l<nz; l++) {
640  ost << " Domain #" << l << " : alpha_l = " << alpha[l]
641  << " , beta_l = " << beta[l] << endl ;
642  }
643 
644  ost << endl << " Values of r at the outer boundary of each domain [km] :"
645  << endl ;
646  ost << " val_r : " ;
647  for (int l=0; l<nz; l++) {
648  ost << " " << val_r(l, 1., 0., 0.) / km ;
649  }
650  ost << endl ;
651 
652  ost << " Coord r : " ;
653  for (int l=0; l<nz; l++) {
654  int nrm1 = mg->get_nr(l) - 1 ;
655  ost << " " << (+r)(l, 0, 0, nrm1) / km ;
656  }
657  ost << endl ;
658 
659  return ost ;
660 }
661 
662  //------------------//
663  // Homothetie //
664  //------------------//
665 
666 
667 void Map_af::homothetie(double fact) {
668 
669  int nz = mg->get_nzone() ;
670 
671  for (int l=0; l<nz; l++) {
672  if (mg->get_type_r(l) == UNSURR) {
673  alpha[l] /= fact ;
674  beta[l] /= fact ;
675  }
676  else {
677  alpha[l] *= fact ;
678  beta[l] *= fact ;
679  }
680  }
681 
682  reset_coord() ;
683 
684 }
685 
686  //----------------------------//
687  // Rescaling of one domain //
688  //----------------------------//
689 
690 void Map_af::resize(int l, double lambda) {
691 
692  // Protections
693  // -----------
694  if (mg->get_type_r(l) != FIN) {
695  cout << "Map_af::resize can be applied only to a shell !" << endl ;
696  abort() ;
697  }
698 
699  // New values of alpha and beta in domain l :
700  // ----------------------------------------
701  double n_alpha = 0.5 * ( (lambda + 1.) * alpha[l] +
702  (lambda - 1.) * beta[l] ) ;
703 
704  double n_beta = 0.5 * ( (lambda - 1.) * alpha[l] +
705  (lambda + 1.) * beta[l] ) ;
706 
707  alpha[l] = n_alpha ;
708  beta[l] = n_beta ;
709 
710  // New values of alpha and beta in domain l+1 :
711  // ------------------------------------------
712  assert(l<mg->get_nzone()-1) ;
713  int lp1 = l + 1 ;
714 
715  if (mg->get_type_r(lp1) == UNSURR) { // compactified case
716 
717  alpha[lp1] = - 0.5 / ( alpha[l] + beta[l] ) ;
718  beta[lp1] = - alpha[lp1] ;
719 
720  }
721  else{ // non-compactified case
722 
723  assert( mg->get_type_r(lp1) == FIN ) ;
724  n_alpha = 0.5 * ( alpha[lp1] - alpha[l] + beta[lp1] - beta[l] ) ;
725  n_beta = 0.5 * ( alpha[lp1] + alpha[l] + beta[lp1] + beta[l] ) ;
726  alpha[lp1] = n_alpha ;
727  beta[lp1] = n_beta ;
728  }
729 
730  // The coords are no longer up to date :
731  reset_coord() ;
732 
733 }
734 
735 
736 
737 
738 
739  //---------------------------//
740  // Homothetie partielle //
741  //-------------------------//
742 
743 
744 void Map_af::homothetie_interne(double fact) {
745 
746  // Dans le noyau
747  alpha[0] *= fact ;
748 
749  // Dans la premiere coquille :
750  double asauve = alpha[1] ;
751  alpha[1] = (1-fact)/2.*beta[1] + (1+fact)/2. * alpha[1] ;
752  beta[1] = (1+fact)/2.*beta[1]+ (1-fact)/2. * asauve ;
753 
754  reset_coord() ;
755 }
756  //------------------------------------------//
757  // Modification of the mapping parameters //
758  //------------------------------------------//
759 
760 void Map_af::set_alpha(double alpha0, int l) {
761 
762  assert(l>=0) ;
763  assert(l<mg->get_nzone()) ;
764 
765  alpha[l] = alpha0 ;
766 
767  reset_coord() ;
768 
769 }
770 
771 void Map_af::set_beta(double beta0, int l) {
772 
773  assert(l>=0) ;
774  assert(l<mg->get_nzone()) ;
775 
776  beta[l] = beta0 ;
777 
778  reset_coord() ;
779 
780 }
781 
782  //------------------------------------//
783  // Angular part of the mapping //
784  //------------------------------------//
785 
786 const Map_af& Map_af::mp_angu(int l_zone) const {
787 //## the handling of l_zone must be improved
788  if (p_mp_angu == 0x0) {
789  const Mg3d& g_angu = (*get_mg()->get_angu_1dom()) ;
790  double Rb = val_r_jk(l_zone, 1., 0, 0) ;
791  Tbl rlim(2) ;
792  rlim.set_etat_qcq() ;
793  rlim.set(0) = Rb ;
794  rlim.set(1) = Rb ;
795  p_mp_angu = new Map_af(g_angu, rlim) ;
796  }
797  return *p_mp_angu ;
798 }
799 
800 // To be done
801 //-----------
802 
803 void Map_af::adapt(const Cmp&, const Param&, int) {
804  const char* f = __FILE__ ;
805  c_est_pas_fait(f) ;
806 }
807 
808 void Map_af::poisson_angu(const Cmp&, Param&, Cmp&, double) const
809 {
810  c_est_pas_fait("Map_af::poisson_angu(Cmp)") ;
811 }
812 
813 
814 
815 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
void set_coord()
Assignment of the building functions to the member Coords.
Definition: map_af.C:533
virtual void sauve(FILE *) const
Save in a file.
Definition: map_af.C:619
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1640
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1678
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1629
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
Radial mapping of rather general form.
Definition: map.h:2783
const double * get_alpha() const
Returns a pointer on the array alpha (values of in each domain)
Definition: map_et.C:1049
const double * get_beta() const
Returns a pointer on the array beta (values of in each domain)
Definition: map_et.C:1053
void set_rot_phi(double phi0)
Sets a new rotation angle.
Definition: map.C:266
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:667
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2054
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:788
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1621
Base class for coordinate mappings.
Definition: map.h:688
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
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.
Definition: map_af_radius.C:99
virtual const Map_af & mp_angu(int) const
Returns the "angular" mapping for the outside of domain l_zone.
Definition: map_af.C:786
virtual void operator=(const Map_af &)
Assignment to another affine mapping.
Definition: map_af.C:510
Map_af * p_mp_angu
Pointer on the "angular" mapping.
Definition: map.h:733
virtual bool operator==(const Map &) const
Comparison operator (egality)
Definition: map_af.C:569
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1613
virtual double val_r_jk(int l, double xi, int j, int k) const
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
Coord tet
coordinate centered on the grid
Definition: map.h:737
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: map_af.C:633
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Definition: map.C:256
Coord phi
coordinate centered on the grid
Definition: map.h:738
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:771
Coord sint
Definition: map.h:739
virtual ~Map_af()
Destructor.
Definition: map_af.C:498
Map_af(const Mg3d &mgrille, const double *r_limits)
Standard Constructor.
Definition: map_af.C:216
Coord stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED)...
Definition: map.h:1597
bool search_file(ifstream &infile, const string &pattern)
A function that searches for a pattern in a file and places the file stream after the found pattern...
Definition: misc.C:53
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1581
double ori_y
Absolute coordinate y of the origin.
Definition: map.h:697
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
double ori_z
Absolute coordinate z of the origin.
Definition: map.h:698
Parameter storage.
Definition: param.h:125
Base class for pure radial mappings.
Definition: map.h:1557
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2056
Coord sinp
Definition: map.h:741
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1669
virtual void sauve(FILE *) const
Save in a file.
Definition: map_radial.C:119
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:760
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1570
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:73
Coord drdt
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED)...
Definition: map.h:1589
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:72
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Coord ya
Absolute y coordinate.
Definition: map.h:749
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:694
Multi-domain grid.
Definition: grilles.h:279
double ori_x
Absolute coordinate x of the origin.
Definition: map.h:696
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
virtual void poisson_angu(const Scalar &source, Param &par, Scalar &uu, double lambda=0) const
Computes the solution of the generalized angular Poisson equation.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord za
Absolute z coordinate.
Definition: map.h:750
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1652
Coord cosp
Definition: map.h:742
Coord x
x coordinate centered on the grid
Definition: map.h:744
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:707
double get_ori_z() const
Returns the z coordinate of the origin.
Definition: map.h:790
virtual void reset_coord()
Resets all the member Coords.
Definition: map_radial.C:129
Basic array class.
Definition: tbl.h:164
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
Definition: map_af.C:744
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1605
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:625
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)
Adaptation of the mapping to a given scalar field.
Definition: map_af.C:803
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Definition: map_af.C:690
Coord z
z coordinate centered on the grid
Definition: map.h:746
double rot_phi
Angle between the x –axis and X –axis.
Definition: map.h:699
Base_vect_cart bvect_cart
Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:715
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1661
Coord cost
Definition: map.h:740