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
30char 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
78namespace Lorene {
79double valeur_equipot_fonc(double, const Param&) ;
80
81//****************************************************************************
82
83void 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Parameter storage.
Definition param.h:125
Basic array class.
Definition tbl.h:161
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).
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
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 prototypes.
Definition app_hor.h:64