map_et_adapt.C
/*
* Method of the class Map_et to compute the mapping parameters to let the
* domain boundaries coincide with isosurfaces of a given scalar field.
*
* (see file map.h for documentation).
*/
/*
* Copyright (c) 2000-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 map_et_adapt_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $" ;
/*
* $Id: map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $
* $Log: map_et_adapt.C,v $
* Revision 1.11 2014/10/13 08:53:03 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.10 2014/10/06 15:13:12 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.9 2010/01/31 16:58:39 e_gourgoulhon
* Back to rev. 1.7 (1.8 has been committed by mistake).
*
* Revision 1.8 2010/01/31 16:51:47 e_gourgoulhon
* File local_settings_linux is now called local_settings_linux_old (for
* it is quite old and not compatible with the gcc 4.x series).
*
* Revision 1.7 2006/09/05 13:39:45 p_grandclement
* update of the bin_ns_bh project
*
* Revision 1.6 2006/05/17 13:27:12 j_novak
* Removed strange character on line 534.
*
* Revision 1.5 2006/05/03 11:15:25 p_grandclement
* modif filtering
*
* Revision 1.4 2006/05/03 07:07:52 p_grandclement
* Petite correction
*
* Revision 1.3 2006/04/25 07:21:59 p_grandclement
* Various changes for the NS_BH project
*
* Revision 1.2 2003/10/03 15:58:48 j_novak
* Cleaning of some headers
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 2.6 2000/11/10 15:19:300 eric
* Appel de Valeur::equipot_outward plutot que Valeur::equipot
* dans la determination des isopotentielles.
*
* Revision 2.5 2000/10/23 14:01:17 eric
* Changement des arguments (en autre, ajout de nz_search, qui est
* desormais a priori independant de nzadapt).
*
* Revision 2.4 2000/10/20 13:13:01 eric
* nzet --> nzadapt.
* nz_search = nzadapt --> nz_search = nzadapt + 1
*
* Revision 2.3 2000/02/17 19:00:53 eric
* nz_search = nzet et non plus nz_search = nz-1.
*
* Revision 2.2 2000/02/15 15:09:12 eric
* Changement du Param : fact_echelle est desormais passe en double_mod.
*
* Revision 2.1 2000/01/03 15:46:59 eric
* *** empty log message ***
*
* Revision 2.0 2000/01/03 13:30:59 eric
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Map/map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $
*
*/
// Headers C
#include <cassert>
// Headers Lorene
#include "cmp.h"
#include "itbl.h"
#include "param.h"
#include "proto.h"
namespace Lorene {
void Map_et::adapt(const Cmp& ent, const Param& par, int nbr_filtre) {
// Parameters of the computation
// -----------------------------
int nitermax = par.get_int(0) ;
int& niter = par.get_int_mod(0) ;
int nzadapt = par.get_int(1) ;
int nz_search = par.get_int(2) ;
int adapt_flag = par.get_int(3) ;
int j_bord = par.get_int(4) ; // index of theta_*
int k_bord = par.get_int(5) ; // index of phi_*
double precis = par.get_double(0) ;
double fact_lamu = par.get_double(1) ;
double fact_echelle = par.get_double(2) ;
const Tbl& ent_limit = par.get_tbl(0) ;
int nz = mg->get_nzone() ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
// Protections
// -----------
assert(ent.get_mp() == this) ;
assert(nzadapt < nz) ;
assert(nzadapt > 0) ;
assert(nz_search >= nzadapt) ;
for (int l=1; l<nz; l++) {
assert( mg->get_nt(l) == nt ) ;
assert( mg->get_np(l) == np ) ;
}
assert(ent_limit.get_ndim() == 1) ;
assert(ent_limit.get_taille() >= nzadapt) ;
assert(ent_limit.get_etat() == ETATQCQ) ;
// Initializations
// ---------------
niter = 0 ;
const double xi_max = 1 ;
//=======================================================================
// Special case of a fixed mapping : only a rescaling is performed
//=======================================================================
if ( (adapt_flag == 0) || (nzadapt == 0) ) {
homothetie(fact_echelle) ;
return ;
}
//=======================================================================
// General case
//=======================================================================
// New mapping parameters
// ----------------------
double* nalpha = new double[nz] ;
double* nbeta = new double[nz] ;
Itbl l_ext(np, nt) ;
Tbl x_ext(np, nt) ;
Tbl xtgg(np, nt, 1) ; // working array constructed from the isosurface
// equation
Tbl xtff(np, nt, 1) ; // idem
//----------------------------------------------------------------
// Adaptation of the nucleus
//----------------------------------------------------------------
double ent0 = ent_limit(0) ;
// Search for the equipotential surface ent = ent0
// ---> l = l_ext(theta, phi)
// xi = x_ext(theta, phi)
// -----------------------------------------------
int niter0 ;
(ent.va).equipot_outward(ent0, nz_search, precis, nitermax, niter0,
l_ext, x_ext) ;
niter = ( niter0 > niter ) ? niter0 : niter ;
// The new mapping must be such that the found isosurface corresponds
// to xi = 1
// Computation of r_ext(theta, phi) - r_eq ---> xtgg
// -------------------------------------------------
xtgg.set_etat_qcq() ;
assert(l_ext.get_etat() == ETATQCQ) ;
double r_eq = val_r_jk(0, xi_max, j_bord, k_bord) ;
double* pxtgg = xtgg.t ;
int* pl_ext = l_ext.t ;
double* px_ext = x_ext.t ;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
*pxtgg = val_r_jk(*pl_ext, *px_ext, j, k) - r_eq ;
// ... next one
pl_ext++ ;
px_ext++ ;
pxtgg++ ;
}
}
// Decomposition into even and odd harmonics in phi of xtgg = r_ext - r_eq :
// even phi harmonics --> xtgg
// odd phi harmonics --> xtff
// ------------------------------------------------------------------------
int base_p = ( ((ent.va).base).b[0] ) & MSQ_P ;
switch( base_p ) {
case P_COSSIN_P : { // case of only even harmonics in phi
xtff.set_etat_zero() ;
break ;
}
case P_COSSIN : { // general case: a Fourier transform must be
// done
xtff.set_etat_qcq() ;
double* pxtff = xtff.t ;
pxtgg = xtgg.t ;
assert(np>=4) ;
int deg[3] ;
int dimc[3] ;
deg[0] = np ;
deg[1] = nt ;
deg[2] = 1;
dimc[0] = np + 2 ;
dimc[1] = nt ;
dimc[2] = 1 ;
double* cf = new double[(np+2)*nt] ;
double* cf0 = new double[(np+2)*nt] ;
double* ff0 = new double[np*nt] ;
for (int i=0; i < np*nt; i++) {
cf[i] = *pxtgg ;
pxtgg++ ;
}
cfpcossin(deg, dimc, cf) ; // Fourier transform
// Even harmonics
// --------------
double* pcf0 = cf0 ;
double* pcf = cf ;
for (int k=0; k<np-1; k += 4) {
for(int j=0; j<2*nt; j++) {
*pcf0 = *pcf ;
pcf0++ ;
pcf++ ;
}
for(int j=0; j<2*nt; j++) {
*pcf0 = 0 ;
pcf0++ ;
pcf++ ;
}
}
if (np%4 == 0) { // particular case of np multiple of 4
for(int j=0; j<2*nt; j++) {
*pcf0 = *pcf ;
pcf0++ ;
pcf++ ;
}
}
cipcossin(deg, dimc, deg, cf0, ff0) ; // Inverse Fourier transform
pxtgg = xtgg.t ;
for (int i=0; i < np*nt; i++) {
*pxtgg = ff0[i] ;
pxtgg++ ;
}
// Odd harmonics
// -------------
pcf0 = cf0 ;
pcf = cf ;
for (int k=0; k<np-1; k += 4) {
for(int j=0; j<2*nt; j++) {
*pcf0 = 0 ;
pcf0++ ;
pcf++ ;
}
for(int j=0; j<2*nt; j++) {
*pcf0 = *pcf ;
pcf0++ ;
pcf++ ;
}
}
if (np%4 == 0) { // particular case of np multiple of 4
for(int j=0; j<2*nt; j++) {
*pcf0 = 0 ;
pcf0++ ;
}
}
cipcossin(deg, dimc, deg, cf0, ff0) ; // TF inverse
pxtff = xtff.t ;
for (int i=0; i < np*nt; i++) {
*pxtff = ff0[i] ;
pxtff++ ;
}
delete [] cf ;
delete [] cf0 ;
delete [] ff0 ;
break ;
}
default : {
cout << "Map_et::adapt: unknown phi basis !" << endl ;
cout << " base_p = " << base_p << endl ;
abort() ;
}
}
// Computation of lambda and mu in the nucleus :
// --------------------------------------------
double lambda = 0 ; // lambda is set to zero because F(theta, phi)
// must not have constant terms in the nucleus
double mu = - fact_lamu * min(xtgg) ;
// Computation of alpha and beta in the nucleus :
// --------------------------------------------
nalpha[0] = r_eq - lambda - mu ;
nbeta[0] = 0 ;
// Computation of F_0, G_0 and {\tilde F_1} :
// ------------------------------------------
Mtbl nff(ff.get_mg()) ;
Mtbl ngg(gg.get_mg()) ;
nff.set_etat_qcq() ;
ngg.set_etat_qcq() ;
*(nff.t[0]) = ( xtff + lambda ) / nalpha[0] ;
*(ngg.t[0]) = ( xtgg + mu ) / nalpha[0] ;
xtff += xtgg ;
//----------------------------------------------------------------
// Adaptation of the shells
//----------------------------------------------------------------
double r_eqlm1 = r_eq ;
// Loop on the shells where the adaptation must be performed
for (int l=1; l<nzadapt; l++) {
ent0 = ent_limit(l) ;
// Search for the equipotential surface ent = ent0
// ---> l = l_ext(theta, phi)
// xi = x_ext(theta, phi)
// -----------------------------------------------
(ent.va).equipot_outward(ent0, nz_search, precis, nitermax, niter0,
l_ext, x_ext) ;
niter = ( niter0 > niter ) ? niter0 : niter ;
// The new mapping must be such that the found isosurface corresponds
// to xi = 1
// Computation of r_ext(theta, phi) - r_eq ---> xtgg
// -------------------------------------------------
xtgg.set_etat_qcq() ;
assert(l_ext.get_etat() == ETATQCQ) ;
r_eq = val_r_jk(l, xi_max, j_bord, k_bord) ;
pxtgg = xtgg.t ;
pl_ext = l_ext.t ;
px_ext = x_ext.t ;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
*pxtgg = val_r_jk(*pl_ext, *px_ext, j, k) - r_eq ;
// ... next one
pl_ext++ ;
px_ext++ ;
pxtgg++ ;
}
}
// Computation of lambda and mu in domain no. l :
// --------------------------------------------
lambda = - fact_lamu * max(xtff) ;
mu = - fact_lamu * min(xtgg) ;
// Computation of alpha and beta in domain no. l :
// --------------------------------------------
nalpha[l] = .5 * ( r_eq - r_eqlm1 + lambda - mu ) ;
nbeta[l] = .5 * ( r_eq + r_eqlm1 - lambda - mu ) ;
// Computation of F_l, G_l and {\tilde F_(l+1)} :
// ------------------------------------------
*(nff.t[l]) = ( xtff + lambda ) / nalpha[l] ;
*(ngg.t[l]) = ( xtgg + mu ) / nalpha[l] ;
xtff = xtgg ;
r_eqlm1 = r_eq ;
} // end of the loop on the shells where the adaptation must be performed
//----------------------------------------------------------------
// Adaptation of the domain of index nzadapt
//----------------------------------------------------------------
if (mg->get_type_r(nzadapt) == UNSURR) {
// Case of a compactified domain
// -----------------------------
xtff = 1 / (xtff + r_eqlm1) - double(1) / r_eqlm1 ;
lambda = - fact_lamu * min(xtff) ;
nalpha[nzadapt] = .5 * ( lambda - double(1) / r_eqlm1 ) ;
nbeta[nzadapt] = - nalpha[nzadapt] ;
// Computation of F_nzadapt :
*(nff.t[nzadapt]) = ( xtff + lambda ) / nalpha[nzadapt] ;
}
else {
// Case of a non-compactified domain
// ----------------------------------
r_eq = val_r_jk(nzadapt, xi_max, j_bord, k_bord) ;
lambda = - fact_lamu * max(xtff) ;
nalpha[nzadapt] = .5 * ( r_eq - r_eqlm1 + lambda ) ;
nbeta[nzadapt] = .5 * ( r_eq + r_eqlm1 - lambda ) ;
// Computation of F_l :
*(nff.t[nzadapt]) = ( xtff + lambda ) / nalpha[nzadapt] ;
}
// In all cases, G_nzadapt is set to zero :
ngg.t[nzadapt]->set_etat_zero() ;
//----------------------------------------------------------------
// Values of alpha, beta, F and G in the most external domains
// where the mapping is unchanged
//----------------------------------------------------------------
for (int l=nzadapt+1; l<nz; l++) {
nalpha[l] = alpha[l] ;
nbeta[l] = beta[l] ;
}
if (ff.get_etat() == ETATZERO) {
for (int l=nzadapt+1; l<nz; l++) {
nff.t[l]->set_etat_zero() ;
}
}
else {
assert(ff.get_etat() == ETATQCQ) ;
assert(ff.c != 0x0) ;
assert(ff.c->get_etat() == ETATQCQ) ;
for (int l=nzadapt+1; l<nz; l++) {
*(nff.t[l]) = *(ff.c->t[l]) ;
}
}
if (gg.get_etat() == ETATZERO) {
for (int l=nzadapt+1; l<nz; l++) {
ngg.t[l]->set_etat_zero() ;
}
}
else {
assert(gg.get_etat() == ETATQCQ) ;
assert(gg.c != 0x0) ;
assert(gg.c->get_etat() == ETATQCQ) ;
for (int l=nzadapt+1; l<nz; l++) {
*(ngg.t[l]) = *(gg.c->t[l]) ;
}
}
//=============================================================================
// Assignment of the mapping parameters alpha, beta, ff and gg to
// the newly computed values
//=============================================================================
for (int l=0; l<nz; l++) {
if (mg->get_type_r(l) == UNSURR) {
alpha[l] = nalpha[l] / fact_echelle ;
beta[l] = nbeta[l] / fact_echelle ;
}
else {
alpha[l] = fact_echelle * nalpha[l] ;
beta[l] = fact_echelle * nbeta[l] ;
}
}
ff = nff ;
gg = ngg ;
ff.std_base_scal() ; // Standard spectral bases for F
gg.std_base_scal() ; // Standard spectral bases for G
if (nbr_filtre !=0) {
ff.coef() ;
gg.coef() ;
ff.set_etat_cf_qcq() ;
gg.set_etat_cf_qcq() ;
for (int l=0 ; l<nzadapt+1 ; l++)
for (int k=np-nbr_filtre ; k<np ; k++)
for (int j=0 ; j<nt ; j++) {
if (ff.c_cf->t[l]->get_etat() != ETATZERO)
ff.c_cf->set(l, k,j,0) = 0 ;
if (gg.c_cf->t[l]->get_etat() != ETATZERO)
gg.c_cf->set(l,k,j,0) = 0 ;
}
}
// The derived quantities must be reset
// ------------------------------------
reset_coord() ;
// Clean exit
// ----------
delete [] nalpha ;
delete [] nbeta ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::get_etat
int get_etat() const
Gives the logical state.
Definition itbl.h:317
Lorene::Itbl::t
int * t
The array of int 's.
Definition itbl.h:132
Lorene::Map_et::beta
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2760
Lorene::Map_et::homothetie
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
Lorene::Map_et::adapt
virtual void adapt(const Cmp &ent, const Param &par, int nbr_filtre=0)
Adaptation of the mapping to a given scalar field.
Definition map_et_adapt.C:108
Lorene::Map_et::ff
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2819
Lorene::Map_et::alpha
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2758
Lorene::Map_et::gg
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2826
Lorene::Map_et::val_r_jk
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...
Definition map_et_radius.C:367
Lorene::Map_et::reset_coord
virtual void reset_coord()
Resets all the member Coords.
Definition map_et.C:628
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
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::get_tbl
const Tbl & get_tbl(int position=0) const
Returns the reference of a Tbl stored in the list.
Definition param.C:567
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::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim).
Definition tbl.h:400
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition type_parite.h:247
P_COSSIN
#define P_COSSIN
dev. standart
Definition type_parite.h:245
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
Lorene
Lorene prototypes.
Definition app_hor.h:64