LORENE
map_radial_reevaluate.C
1 /*
2  * Copyright (c) 2000-2001 Eric Gourgoulhon
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char map_radial_reevaluate_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reevaluate.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $" ;
24 
25 /*
26  * $Id: map_radial_reevaluate.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
27  * $Log: map_radial_reevaluate.C,v $
28  * Revision 1.4 2014/10/13 08:53:06 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.3 2014/10/06 15:13:14 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.2 2007/05/15 12:43:57 p_grandclement
35  * Scalar version of reevaluate
36  *
37  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
38  * LORENE
39  *
40  * Revision 2.0 2000/01/04 15:24:00 eric
41  * *** empty log message ***
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reevaluate.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
45  *
46  */
47 
48 // Headers C
49 #include <cstdlib>
50 
51 // Headers Lorene
52 #include "map.h"
53 #include "cmp.h"
54 #include "param.h"
55 #include "scalar.h"
56 
57 namespace Lorene {
58 void Map_radial::reevaluate(const Map* mp_prev0, int nzet, Cmp& uu) const {
59 
60  const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
61 
62  if (mp_prev == 0x0) {
63  cout <<
64  "Map_radial::reevaluate : the mapping mp_prev does not belong"
65  << endl ;
66  cout << " to the class Map_radial !" << endl ;
67  abort() ;
68  }
69 
70  int nz = mg->get_nzone() ;
71 
72  // Protections
73  // -----------
74  assert(uu.get_mp() == this) ;
75  assert(uu.get_dzpuis() == 0) ;
76  assert(uu.get_etat() != ETATNONDEF) ;
77  assert(mp_prev->mg == mg) ;
78  assert(nzet <= nz) ;
79 
80 
81  // Maybe nothing to do ?
82  if ( uu.get_etat() == ETATZERO ) {
83  return ;
84  }
85 
86  assert(uu.get_etat() == ETATQCQ) ;
87  (uu.va).coef() ; // the coefficients of uu are required
88 
89  // Copy of the coefficients
90  Mtbl_cf va_cf = *((uu.va).c_cf) ;
91 
92  // Preparation of uu.va for the storage of the new values.
93  // ------------------------------------------------------
94 
95  (uu.va).set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
96  // if it does not exist already
97  // Destroys the coefficients
98 
99  Mtbl& va_c = *((uu.va).c) ; // Reference to the Mtbl uu.va.c
100 
101  va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
102  // domain if they do not exist already
103 
104 
105  // Values of the Coord r
106  // ---------------------
107 
108  if ( r.c == 0x0 ) r.fait() ;
109  const Mtbl& rc = *(r.c) ;
110 
111  // Precision for val_lx_jk :
112  // ------------------------
113  int nitermax = 100 ; // Maximum number of iterations in the secant method
114  int niter ; // Number of iterations effectively used
115  double precis = 1.e-15 ; // Absolute precision in the secant method
116  Param par ;
117  par.add_int(nitermax) ;
118  par.add_int_mod(niter) ;
119  par.add_double(precis) ;
120 
121  // Loop on the domains where the evaluation is to be performed
122  // -----------------------------------------------------------
123 
124  for (int l=0; l<nzet; l++) {
125  int nr = mg->get_nr(l) ;
126  int nt = mg->get_nt(l) ;
127  int np = mg->get_np(l) ;
128 
129  va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
130  // store the result
131 
132  double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
133 
134  double* pr = (rc.t[l])->t ; // Pointer on the values of r
135 
136  // Loop on all the grid points in the considered domain
137 
138  for (int k=0; k<np; k++) {
139  for (int j=0; j<nt; j++) {
140  for (int i=0; i<nr; i++) {
141 
142  // domain l0 and value of the coordinate xi0 of the
143  // considered point in the previous mapping
144 
145  int l0 ;
146  double xi0 ;
147  mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
148 
149  // Value of uu at this point
150 
151  *ptx = va_cf.val_point_jk(l0, xi0, j, k) ;
152 
153  // next point
154  pr++ ;
155  ptx++ ;
156  }
157  }
158  }
159 
160 
161  } // End of the loop on the domains where the evaluation had to be performed
162 
163  // In the remaining domains, uu is set to zero:
164  // -------------------------------------------
165 
166  uu.annule(nzet, nz - 1) ;
167 
168 
169 
170 
171 }
172 
173 void Map_radial::reevaluate(const Map* mp_prev0, int nzet, Scalar& uu) const {
174 
175  const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
176 
177  if (mp_prev == 0x0) {
178  cout <<
179  "Map_radial::reevaluate : the mapping mp_prev does not belong"
180  << endl ;
181  cout << " to the class Map_radial !" << endl ;
182  abort() ;
183  }
184 
185  int nz = mg->get_nzone() ;
186 
187  // Protections
188  // -----------
189  assert(uu.get_mp() == *this) ;
190  assert(uu.get_dzpuis() == 0) ;
191  assert(uu.get_etat() != ETATNONDEF) ;
192  assert(mp_prev->mg == mg) ;
193  assert(nzet <= nz) ;
194 
195 
196  // Maybe nothing to do ?
197  if ( uu.get_etat() == ETATZERO ) {
198  return ;
199  }
200 
201  assert(uu.get_etat() == ETATQCQ) ;
202  uu.set_spectral_va().coef() ; // the coefficients of uu are required
203 
204  // Copy of the coefficients
205  Mtbl_cf va_cf = *(uu.set_spectral_va().c_cf) ;
206 
207  // Preparation of uu.va for the storage of the new values.
208  // ------------------------------------------------------
209 
210  uu.set_spectral_va().set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
211  // if it does not exist already
212  // Destroys the coefficients
213 
214  Mtbl& va_c = *(uu.set_spectral_va().c) ; // Reference to the Mtbl uu.va.c
215 
216  va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
217  // domain if they do not exist already
218 
219 
220  // Values of the Coord r
221  // ---------------------
222 
223  if ( r.c == 0x0 ) r.fait() ;
224  const Mtbl& rc = *(r.c) ;
225 
226  // Precision for val_lx_jk :
227  // ------------------------
228  int nitermax = 100 ; // Maximum number of iterations in the secant method
229  int niter ; // Number of iterations effectively used
230  double precis = 1.e-15 ; // Absolute precision in the secant method
231  Param par ;
232  par.add_int(nitermax) ;
233  par.add_int_mod(niter) ;
234  par.add_double(precis) ;
235 
236  // Loop on the domains where the evaluation is to be performed
237  // -----------------------------------------------------------
238 
239  for (int l=0; l<nzet; l++) {
240  int nr = mg->get_nr(l) ;
241  int nt = mg->get_nt(l) ;
242  int np = mg->get_np(l) ;
243 
244  va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
245  // store the result
246 
247  double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
248 
249  double* pr = (rc.t[l])->t ; // Pointer on the values of r
250 
251  // Loop on all the grid points in the considered domain
252 
253  for (int k=0; k<np; k++) {
254  for (int j=0; j<nt; j++) {
255  for (int i=0; i<nr; i++) {
256 
257  // domain l0 and value of the coordinate xi0 of the
258  // considered point in the previous mapping
259 
260  int l0 ;
261  double xi0 ;
262  mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
263 
264  // Value of uu at this point
265 
266  *ptx = va_cf.val_point_jk(l0, xi0, j, k) ;
267 
268  // next point
269  pr++ ;
270  ptx++ ;
271  }
272  }
273  }
274 
275 
276  } // End of the loop on the domains where the evaluation had to be performed
277 
278  // In the remaining domains, uu is set to zero:
279  // -------------------------------------------
280 
281  uu.annule(nzet, nz - 1) ;
282 
283 
284 
285 
286 }
287 }
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::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::Coord::c
Mtbl * c
The coordinate values at each grid point.
Definition: coord.h:97
Lorene::Map_radial::reevaluate
virtual void reevaluate(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_reevaluate.C:58
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::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::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::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Map_radial
Base class for pure radial mappings.
Definition: map.h:1536
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid
Definition: valeur.h:299
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
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::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Cmp::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
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::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::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition: mtbl.h:132
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
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::Coord::fait
void fait() const
Computes, at each point of the grid, the value of the coordinate or mapping derivative represented by...
Definition: coord.C:116
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Lorene::Map::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676