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