map_radial_reeval_symy.C
/*
* 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_radial_reeval_symy_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reeval_symy.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $" ;
/*
* $Id: map_radial_reeval_symy.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
* $Log: map_radial_reeval_symy.C,v $
* Revision 1.4 2014/10/13 08:53:06 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2014/10/06 15:13:14 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.2 2007/05/15 12:43:57 p_grandclement
* Scalar version of reevaluate
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 2.0 2000/03/06 11:30:19 eric
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reeval_symy.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
// Headers Lorene
#include "map.h"
#include "cmp.h"
#include "param.h"
#include "scalar.h"
namespace Lorene {
void Map_radial::reevaluate_symy(const Map* mp_prev0, int nzet, Cmp& uu) const {
const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
if (mp_prev == 0x0) {
cout <<
"Map_radial::reevaluate_symy : the mapping mp_prev does not belong"
<< endl ;
cout << " to the class Map_radial !" << endl ;
abort() ;
}
int nz = mg->get_nzone() ;
// Protections
// -----------
assert(uu.get_mp() == this) ;
assert(uu.get_dzpuis() == 0) ;
assert(uu.get_etat() != ETATNONDEF) ;
assert(mp_prev->mg == mg) ;
assert(nzet <= nz) ;
assert(mg->get_type_p() == NONSYM) ;
// Maybe nothing to do ?
if ( uu.get_etat() == ETATZERO ) {
return ;
}
assert(uu.get_etat() == ETATQCQ) ;
(uu.va).coef() ; // the coefficients of uu are required
// Copy of the coefficients
Mtbl_cf va_cf = *((uu.va).c_cf) ;
// Preparation of uu.va for the storage of the new values.
// ------------------------------------------------------
(uu.va).set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
// if it does not exist already
// Destroys the coefficients
Mtbl& va_c = *((uu.va).c) ; // Reference to the Mtbl uu.va.c
va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
// domain if they do not exist already
// Values of the Coord r
// ---------------------
if ( r.c == 0x0 ) r.fait() ;
const Mtbl& rc = *(r.c) ;
// Precision for val_lx_jk :
// ------------------------
int nitermax = 100 ; // Maximum number of iterations in the secant method
int niter ; // Number of iterations effectively used
double precis = 1.e-15 ; // Absolute precision in the secant method
Param par ;
par.add_int(nitermax) ;
par.add_int_mod(niter) ;
par.add_double(precis) ;
// Loop on the domains where the evaluation is to be performed
// -----------------------------------------------------------
for (int l=0; l<nzet; l++) {
int nr = mg->get_nr(l) ;
int nt = mg->get_nt(l) ;
int np = mg->get_np(l) ;
va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
// store the result
double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
double* pr = (rc.t[l])->t ; // Pointer on the values of r
// Loop on half the grid points in the considered arrival domain
// (the other half will be obtained by symmetry with respect to
// the y=0 plane).
for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
// domain l0 and value of the coordinate xi0 of the
// considered point in the previous mapping
int l0 ;
double xi0 ;
mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
// Value of uu at this point
*ptx = va_cf.val_point_jk_symy(l0, xi0, j, k) ;
// next point
pr++ ;
ptx++ ;
}
}
}
// The remaining points are obtained by symmetry with rspect to the
// y=0 plane
for (int k=np/2+1 ; k<np ; k++) {
// pointer on the value (already computed) at the point symmetric
// with respect to the plane y=0
double* ptx_symy = (va_c.t[l])->t + (np-k)*nt*nr ;
// copy :
for (int j=0 ; j<nt ; j++) {
for (int i=0 ; i<nr ; i++) {
*ptx = *ptx_symy ;
ptx++ ;
ptx_symy++ ;
}
}
}
} // End of the loop on the domains where the evaluation had to be performed
// In the remaining domains, uu is set to zero:
// -------------------------------------------
uu.annule(nzet, nz - 1) ;
}
void Map_radial::reevaluate_symy(const Map* mp_prev0, int nzet, Scalar& uu) const {
const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
if (mp_prev == 0x0) {
cout <<
"Map_radial::reevaluate_symy : the mapping mp_prev does not belong"
<< endl ;
cout << " to the class Map_radial !" << endl ;
abort() ;
}
int nz = mg->get_nzone() ;
// Protections
// -----------
assert(uu.get_mp() == *this) ;
assert(uu.get_dzpuis() == 0) ;
assert(uu.get_etat() != ETATNONDEF) ;
assert(mp_prev->mg == mg) ;
assert(nzet <= nz) ;
assert(mg->get_type_p() == NONSYM) ;
// Maybe nothing to do ?
if ( uu.get_etat() == ETATZERO ) {
return ;
}
assert(uu.get_etat() == ETATQCQ) ;
uu.set_spectral_va().coef() ; // the coefficients of uu are required
// Copy of the coefficients
Mtbl_cf va_cf = *(uu.set_spectral_va().c_cf) ;
// Preparation of uu.va for the storage of the new values.
// ------------------------------------------------------
uu.set_spectral_va().set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
// if it does not exist already
// Destroys the coefficients
Mtbl& va_c = *(uu.set_spectral_va().c) ; // Reference to the Mtbl uu.va.c
va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
// domain if they do not exist already
// Values of the Coord r
// ---------------------
if ( r.c == 0x0 ) r.fait() ;
const Mtbl& rc = *(r.c) ;
// Precision for val_lx_jk :
// ------------------------
int nitermax = 100 ; // Maximum number of iterations in the secant method
int niter ; // Number of iterations effectively used
double precis = 1.e-15 ; // Absolute precision in the secant method
Param par ;
par.add_int(nitermax) ;
par.add_int_mod(niter) ;
par.add_double(precis) ;
// Loop on the domains where the evaluation is to be performed
// -----------------------------------------------------------
for (int l=0; l<nzet; l++) {
int nr = mg->get_nr(l) ;
int nt = mg->get_nt(l) ;
int np = mg->get_np(l) ;
va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
// store the result
double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
double* pr = (rc.t[l])->t ; // Pointer on the values of r
// Loop on half the grid points in the considered arrival domain
// (the other half will be obtained by symmetry with respect to
// the y=0 plane).
for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
// domain l0 and value of the coordinate xi0 of the
// considered point in the previous mapping
int l0 ;
double xi0 ;
mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
// Value of uu at this point
*ptx = va_cf.val_point_jk_symy(l0, xi0, j, k) ;
// next point
pr++ ;
ptx++ ;
}
}
}
// The remaining points are obtained by symmetry with rspect to the
// y=0 plane
for (int k=np/2+1 ; k<np ; k++) {
// pointer on the value (already computed) at the point symmetric
// with respect to the plane y=0
double* ptx_symy = (va_c.t[l])->t + (np-k)*nt*nr ;
// copy :
for (int j=0 ; j<nt ; j++) {
for (int i=0 ; i<nr ; i++) {
*ptx = *ptx_symy ;
ptx++ ;
ptx_symy++ ;
}
}
}
} // End of the loop on the domains where the evaluation had to be performed
// In the remaining domains, uu is set to zero:
// -------------------------------------------
uu.annule(nzet, nz - 1) ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Lorene::Cmp::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Map_radial::reevaluate_symy
virtual void reevaluate_symy(const Map *mp_prev, int nzet, Cmp &uu) const
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
Definition map_radial_reeval_symy.C:59
Lorene::Map_radial::val_lx_jk
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
Lorene::Map_radial::Map_radial
Map_radial(const Mg3d &mgrid)
Constructor from a grid (protected to make Map_radial an abstract class).
Definition map_radial.C:89
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::val_point_jk_symy
double val_point_jk_symy(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_vp_symy.C:201
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::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::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::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Valeur::set_etat_c_qcq
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:718