valeur_equipot.C
/*
* Method of the class Valeur to compute an equipotential surface.
*
* (see file valeur.h for the documentation).
*/
/*
* Copyright (c) 1999-2001 Eric Gourgoulhon
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char valeur_equipot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $" ;
/*
* $Id: valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
* $Log: valeur_equipot.C,v $
* Revision 1.5 2014/10/13 08:53:50 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2014/10/06 15:13:23 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.3 2003/12/19 15:05:56 j_novak
* Added some initializations
*
* Revision 1.2 2002/10/16 14:37:15 j_novak
* Reorganization of #include instructions of standard C++, in order to
* use experimental version 3 of gcc.
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 1.3 2000/01/28 17:06:37 eric
* Changement du cas l2<nz_search-1.
*
* Revision 1.2 2000/01/03 15:07:50 eric
* *** empty log message ***
*
* Revision 1.1 1999/12/29 13:12:08 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cmath>
// Headers Lorene
#include "valeur.h"
#include "itbl.h"
#include "param.h"
#include "nbr_spx.h"
#include "utilitaires.h"
// Local prototypes
namespace Lorene {
double valeur_equipot_fonc(double, const Param&) ;
//****************************************************************************
void Valeur::equipot(double uu0, int nz_search, double precis, int nitermax,
int& niter, Itbl& l_iso, Tbl& xi_iso) const {
int nz = mg->get_nzone() ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
// Protections
// -----------
assert(etat == ETATQCQ) ;
assert(nz_search > 0) ;
assert(nz_search <= nz) ;
for (int l=1; l<nz_search; l++) {
assert( mg->get_nt(l) == nt ) ;
assert( mg->get_np(l) == np ) ;
}
coef() ; // the spectral coefficients are required
l_iso.set_etat_qcq() ;
xi_iso.set_etat_qcq() ;
Param parf ;
int j, k, l2 ;
parf.add_int_mod(j, 0) ;
parf.add_int_mod(k, 1) ;
parf.add_int_mod(l2, 2) ;
parf.add_double(uu0) ;
parf.add_mtbl_cf(*c_cf) ;
// Loop of phi and theta
// ---------------------
niter = 0 ;
for (k=0; k<np; k++) {
for (j=0; j<nt; j++) {
// 1. Recherche du point de collocation en r (l2,i2) egal a ou situe
// immediatement avant r_iso(phi,theta)
int i2 = 0;
l2 = nz ;
for (int l=nz_search-1; l>= 0; l--) { // On part de la zone la plus extreme
int nr = mg->get_nr(l) ;
for (int i=nr-1; i >= 0; i--) {
double uux = (*this)(l, k, j, i) ;
if ( ( (uux > uu0) || ( fabs(uux-uu0) < 1.e-14 ) ) &&
(uux != __infinity) ) {
l2 = l ; // on a trouve le point
i2 = i ; //
break ;
}
}
if (l2 == l) break ;
} // fin de la boucle sur les zones
if (l2==nz) {
cout << "Valeur::equipot: the point uu >= uu0 has not been found" << endl ;
cout << " for the phi index " << k << endl ;
cout << " and the theta index " << j << endl ;
cout << " uu0 = " << uu0 << endl ;
abort () ;
}
// 2. Point suivant (l2,i3)
int i3 = 0 ;
if (i2 < mg->get_nr(l2) -1) { // on reste dans la meme zone
i3 = i2 + 1 ;
}
else {
if (l2<nz_search-1) { // on change de zone
double uux = (*this)(l2, k, j, i2) ;
if ( ( fabs(uux-uu0) < 1.e-14 ) ) {
// it is OK
cout <<
"Valeur::equipot: WARNING : fabs(uux-uu0) < 1.e-14" << endl ;
l_iso.set(k, j) = l2 ;
xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
}
else{
cout << "Valeur::equipot: PROBLEM !!!" << endl ;
cout << " k, j : " << k << " " << j << endl ;
cout << " uu0 : " << uu0 << endl ;
abort () ;
}
}
else { // on est a l'extremite de la grille
l_iso.set(k, j) = nz_search-1 ;
xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
continue ; // on passe au theta suivant
}
}
l_iso.set(k, j) = l2 ;
double x2 = (mg->get_grille3d(l2))->x[i2] ;
double x3 = (mg->get_grille3d(l2))->x[i3] ;
double y2 = (*this)(l2, k, j, i2) - uu0 ;
// 3. Recherche de x0
int niter0 = 0 ;
if (fabs(y2) < 1.e-14) {
xi_iso.set(k, j) = x2 ;
}
else {
xi_iso.set(k, j) = zerosec(valeur_equipot_fonc, parf, x2, x3, precis,
nitermax, niter0 ) ;
}
niter = ( niter0 > niter ) ? niter0 : niter ;
} // fin de la boucle sur theta
} // fin de la boucle sur phi
}
//****************************************************************************
double valeur_equipot_fonc(double xi, const Param& par) {
int j = par.get_int_mod(0) ;
int k = par.get_int_mod(1) ;
int l = par.get_int_mod(2) ;
double uu0 = par.get_double() ;
const Mtbl_cf& cuu = par.get_mtbl_cf() ;
return cuu.val_point_jk(l, xi, j, k) - uu0 ;
}
}
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition itbl.C:261
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::val_point_jk
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...
Definition mtbl_cf_val_point.C:258
Lorene::Param
Parameter storage.
Definition param.h:125
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::Param::get_double
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
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::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::Param::get_int_mod
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition param.C:430
Lorene::Param::get_mtbl_cf
const Mtbl_cf & get_mtbl_cf(int position=0) const
Returns the reference of a Mtbl_cf stored in the list.
Definition param.C:1325
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Valeur::equipot
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).
Definition valeur_equipot.C:81
Lorene::Valeur::mg
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Valeur::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
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
Lorene prototypes.
Definition app_hor.h:64