LORENE
valeur_equipot.C
1 /*
2  * Method of the class Valeur to compute an equipotential surface.
3  *
4  * (see file valeur.h for the documentation).
5  */
6 
7 /*
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: valeur_equipot.C,v 1.6 2016/12/05 16:18:20 j_novak Exp $
33  * $Log: valeur_equipot.C,v $
34  * Revision 1.6 2016/12/05 16:18:20 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.5 2014/10/13 08:53:50 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.4 2014/10/06 15:13:23 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.3 2003/12/19 15:05:56 j_novak
44  * Added some initializations
45  *
46  * Revision 1.2 2002/10/16 14:37:15 j_novak
47  * Reorganization of #include instructions of standard C++, in order to
48  * use experimental version 3 of gcc.
49  *
50  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
51  * LORENE
52  *
53  * Revision 1.3 2000/01/28 17:06:37 eric
54  * Changement du cas l2<nz_search-1.
55  *
56  * Revision 1.2 2000/01/03 15:07:50 eric
57  * *** empty log message ***
58  *
59  * Revision 1.1 1999/12/29 13:12:08 eric
60  * Initial revision
61  *
62  *
63  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_equipot.C,v 1.6 2016/12/05 16:18:20 j_novak Exp $
64  *
65  */
66 
67 // Headers C
68 #include <cstdlib>
69 #include <cmath>
70 
71 // Headers Lorene
72 #include "valeur.h"
73 #include "itbl.h"
74 #include "param.h"
75 #include "nbr_spx.h"
76 #include "utilitaires.h"
77 
78 // Local prototypes
79 namespace Lorene {
80 double valeur_equipot_fonc(double, const Param&) ;
81 
82 //****************************************************************************
83 
84 void Valeur::equipot(double uu0, int nz_search, double precis, int nitermax,
85  int& niter, Itbl& l_iso, Tbl& xi_iso) const {
86 
87 
88  int nz = mg->get_nzone() ;
89  int nt = mg->get_nt(0) ;
90  int np = mg->get_np(0) ;
91 
92  // Protections
93  // -----------
94  assert(etat == ETATQCQ) ;
95  assert(nz_search > 0) ;
96  assert(nz_search <= nz) ;
97  for (int l=1; l<nz_search; l++) {
98  assert( mg->get_nt(l) == nt ) ;
99  assert( mg->get_np(l) == np ) ;
100  }
101 
102  coef() ; // the spectral coefficients are required
103 
104  l_iso.set_etat_qcq() ;
105  xi_iso.set_etat_qcq() ;
106 
107  Param parf ;
108  int j, k, l2 ;
109  parf.add_int_mod(j, 0) ;
110  parf.add_int_mod(k, 1) ;
111  parf.add_int_mod(l2, 2) ;
112  parf.add_double(uu0) ;
113  parf.add_mtbl_cf(*c_cf) ;
114 
115 
116  // Loop of phi and theta
117  // ---------------------
118 
119  niter = 0 ;
120 
121  for (k=0; k<np; k++) {
122 
123  for (j=0; j<nt; j++) {
124 
125 
126 // 1. Recherche du point de collocation en r (l2,i2) egal a ou situe
127 // immediatement avant r_iso(phi,theta)
128 
129  int i2 = 0;
130  l2 = nz ;
131  for (int l=nz_search-1; l>= 0; l--) { // On part de la zone la plus extreme
132  int nr = mg->get_nr(l) ;
133 
134  for (int i=nr-1; i >= 0; i--) {
135  double uux = (*this)(l, k, j, i) ;
136  if ( ( (uux > uu0) || ( fabs(uux-uu0) < 1.e-14 ) ) &&
137  (uux != __infinity) ) {
138  l2 = l ; // on a trouve le point
139  i2 = i ; //
140  break ;
141  }
142  }
143 
144  if (l2 == l) break ;
145  } // fin de la boucle sur les zones
146 
147  if (l2==nz) {
148  cout << "Valeur::equipot: the point uu >= uu0 has not been found" << endl ;
149  cout << " for the phi index " << k << endl ;
150  cout << " and the theta index " << j << endl ;
151  cout << " uu0 = " << uu0 << endl ;
152  abort () ;
153  }
154 
155 // 2. Point suivant (l2,i3)
156  int i3 = 0 ;
157  if (i2 < mg->get_nr(l2) -1) { // on reste dans la meme zone
158  i3 = i2 + 1 ;
159  }
160  else {
161  if (l2<nz_search-1) { // on change de zone
162 
163  double uux = (*this)(l2, k, j, i2) ;
164  if ( ( fabs(uux-uu0) < 1.e-14 ) ) {
165  // it is OK
166  cout <<
167  "Valeur::equipot: WARNING : fabs(uux-uu0) < 1.e-14" << endl ;
168  l_iso.set(k, j) = l2 ;
169  xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
170  }
171  else{
172  cout << "Valeur::equipot: PROBLEM !!!" << endl ;
173  cout << " k, j : " << k << " " << j << endl ;
174  cout << " uu0 : " << uu0 << endl ;
175  abort () ;
176  }
177  }
178  else { // on est a l'extremite de la grille
179  l_iso.set(k, j) = nz_search-1 ;
180  xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
181  continue ; // on passe au theta suivant
182  }
183  }
184 
185  l_iso.set(k, j) = l2 ;
186 
187  double x2 = (mg->get_grille3d(l2))->x[i2] ;
188  double x3 = (mg->get_grille3d(l2))->x[i3] ;
189 
190  double y2 = (*this)(l2, k, j, i2) - uu0 ;
191 
192 // 3. Recherche de x0
193 
194  int niter0 = 0 ;
195  if (fabs(y2) < 1.e-14) {
196  xi_iso.set(k, j) = x2 ;
197  }
198  else {
199  xi_iso.set(k, j) = zerosec(valeur_equipot_fonc, parf, x2, x3, precis,
200  nitermax, niter0 ) ;
201  }
202 
203  niter = ( niter0 > niter ) ? niter0 : niter ;
204 
205  } // fin de la boucle sur theta
206  } // fin de la boucle sur phi
207 
208 }
209 
210 
211 
212 //****************************************************************************
213 
214 double valeur_equipot_fonc(double xi, const Param& par) {
215 
216  int j = par.get_int_mod(0) ;
217  int k = par.get_int_mod(1) ;
218  int l = par.get_int_mod(2) ;
219  double uu0 = par.get_double() ;
220  const Mtbl_cf& cuu = par.get_mtbl_cf() ;
221 
222  return cuu.val_point_jk(l, xi, j, k) - uu0 ;
223 }
224 
225 }
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:517
void equipot(double uu0, int nz_search, double precis, int nitermax, int &niter, Itbl &l_iso, Tbl &xi_iso) const
Determines an equipotential surface of the field represented by *this (inward search).
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Lorene prototypes.
Definition: app_hor.h:67
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Basic integer array class.
Definition: itbl.h:122
void add_mtbl_cf(const Mtbl_cf &mi, int position=0)
Adds the address of a new Mtbl_cf to the list.
Definition: param.C:1282
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition: zerosec.C:92
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition: valeur.h:302
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...
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:433
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: itbl.C:264
const Mtbl_cf & get_mtbl_cf(int position=0) const
Returns the reference of a Mtbl_cf stored in the list.
Definition: param.C:1328
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: valeur.h:305
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388