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 
24 
25 /*
26  * $Id: map_radial_reevaluate.C,v 1.6 2021/10/25 14:30:27 j_novak Exp $
27  * $Log: map_radial_reevaluate.C,v $
28  * Revision 1.6 2021/10/25 14:30:27 j_novak
29  * The cancelling of the non reevaluated domains is not done if nezt = nz
30  *
31  * Revision 1.5 2016/12/05 16:17:58 j_novak
32  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
33  *
34  * Revision 1.4 2014/10/13 08:53:06 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.3 2014/10/06 15:13:14 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.2 2007/05/15 12:43:57 p_grandclement
41  * Scalar version of reevaluate
42  *
43  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
44  * LORENE
45  *
46  * Revision 2.0 2000/01/04 15:24:00 eric
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reevaluate.C,v 1.6 2021/10/25 14:30:27 j_novak Exp $
51  *
52  */
53 
54 // Headers C
55 #include <cstdlib>
56 
57 // Headers Lorene
58 #include "map.h"
59 #include "cmp.h"
60 #include "param.h"
61 #include "scalar.h"
62 
63 namespace Lorene {
64 void Map_radial::reevaluate(const Map* mp_prev0, int nzet, Cmp& uu) const {
65 
66  const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
67 
68  if (mp_prev == 0x0) {
69  cout <<
70  "Map_radial::reevaluate : the mapping mp_prev does not belong"
71  << endl ;
72  cout << " to the class Map_radial !" << endl ;
73  abort() ;
74  }
75 
76  int nz = mg->get_nzone() ;
77 
78  // Protections
79  // -----------
80  assert(uu.get_mp() == this) ;
81  assert(uu.get_dzpuis() == 0) ;
82  assert(uu.get_etat() != ETATNONDEF) ;
83  assert(mp_prev->mg == mg) ;
84  assert(nzet <= nz) ;
85 
86 
87  // Maybe nothing to do ?
88  if ( uu.get_etat() == ETATZERO ) {
89  return ;
90  }
91 
92  assert(uu.get_etat() == ETATQCQ) ;
93  (uu.va).coef() ; // the coefficients of uu are required
94 
95  // Copy of the coefficients
96  Mtbl_cf va_cf = *((uu.va).c_cf) ;
97 
98  // Preparation of uu.va for the storage of the new values.
99  // ------------------------------------------------------
100 
101  (uu.va).set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
102  // if it does not exist already
103  // Destroys the coefficients
104 
105  Mtbl& va_c = *((uu.va).c) ; // Reference to the Mtbl uu.va.c
106 
107  va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
108  // domain if they do not exist already
109 
110 
111  // Values of the Coord r
112  // ---------------------
113 
114  if ( r.c == 0x0 ) r.fait() ;
115  const Mtbl& rc = *(r.c) ;
116 
117  // Precision for val_lx_jk :
118  // ------------------------
119  int nitermax = 100 ; // Maximum number of iterations in the secant method
120  int niter ; // Number of iterations effectively used
121  double precis = 1.e-15 ; // Absolute precision in the secant method
122  Param par ;
123  par.add_int(nitermax) ;
124  par.add_int_mod(niter) ;
125  par.add_double(precis) ;
126 
127  // Loop on the domains where the evaluation is to be performed
128  // -----------------------------------------------------------
129 
130  for (int l=0; l<nzet; l++) {
131  int nr = mg->get_nr(l) ;
132  int nt = mg->get_nt(l) ;
133  int np = mg->get_np(l) ;
134 
135  va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
136  // store the result
137 
138  double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
139 
140  double* pr = (rc.t[l])->t ; // Pointer on the values of r
141 
142  // Loop on all the grid points in the considered domain
143 
144  for (int k=0; k<np; k++) {
145  for (int j=0; j<nt; j++) {
146  for (int i=0; i<nr; i++) {
147 
148  // domain l0 and value of the coordinate xi0 of the
149  // considered point in the previous mapping
150 
151  int l0 ;
152  double xi0 ;
153  mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
154 
155  // Value of uu at this point
156 
157  *ptx = va_cf.val_point_jk(l0, xi0, j, k) ;
158 
159  // next point
160  pr++ ;
161  ptx++ ;
162  }
163  }
164  }
165 
166 
167  } // End of the loop on the domains where the evaluation had to be performed
168 
169  // In the remaining domains, uu is set to zero:
170  // -------------------------------------------
171 
172  uu.annule(nzet, nz - 1) ;
173 
174 
175 
176 
177 }
178 
179 void Map_radial::reevaluate(const Map* mp_prev0, int nzet, Scalar& uu) const {
180 
181  const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
182 
183  if (mp_prev == 0x0) {
184  cout <<
185  "Map_radial::reevaluate : the mapping mp_prev does not belong"
186  << endl ;
187  cout << " to the class Map_radial !" << endl ;
188  abort() ;
189  }
190 
191  int nz = mg->get_nzone() ;
192 
193  // Protections
194  // -----------
195  assert(uu.get_mp() == *this) ;
196  assert(uu.get_dzpuis() == 0) ;
197  assert(uu.get_etat() != ETATNONDEF) ;
198  assert(mp_prev->mg == mg) ;
199  assert(nzet <= nz) ;
200 
201 
202  // Maybe nothing to do ?
203  if ( uu.get_etat() == ETATZERO ) {
204  return ;
205  }
206 
207  assert(uu.get_etat() == ETATQCQ) ;
208  uu.set_spectral_va().coef() ; // the coefficients of uu are required
209 
210  // Copy of the coefficients
211  Mtbl_cf va_cf = *(uu.set_spectral_va().c_cf) ;
212 
213  // Preparation of uu.va for the storage of the new values.
214  // ------------------------------------------------------
215 
216  uu.set_spectral_va().set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
217  // if it does not exist already
218  // Destroys the coefficients
219 
220  Mtbl& va_c = *(uu.set_spectral_va().c) ; // Reference to the Mtbl uu.va.c
221 
222  va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
223  // domain if they do not exist already
224 
225 
226  // Values of the Coord r
227  // ---------------------
228 
229  if ( r.c == 0x0 ) r.fait() ;
230  const Mtbl& rc = *(r.c) ;
231 
232  // Precision for val_lx_jk :
233  // ------------------------
234  int nitermax = 100 ; // Maximum number of iterations in the secant method
235  int niter ; // Number of iterations effectively used
236  double precis = 1.e-15 ; // Absolute precision in the secant method
237  Param par ;
238  par.add_int(nitermax) ;
239  par.add_int_mod(niter) ;
240  par.add_double(precis) ;
241 
242  // Loop on the domains where the evaluation is to be performed
243  // -----------------------------------------------------------
244 
245  for (int l=0; l<nzet; l++) {
246  int nr = mg->get_nr(l) ;
247  int nt = mg->get_nt(l) ;
248  int np = mg->get_np(l) ;
249 
250  va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
251  // store the result
252 
253  double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
254 
255  double* pr = (rc.t[l])->t ; // Pointer on the values of r
256 
257  // Loop on all the grid points in the considered domain
258 
259  for (int k=0; k<np; k++) {
260  for (int j=0; j<nt; j++) {
261  for (int i=0; i<nr; i++) {
262 
263  // domain l0 and value of the coordinate xi0 of the
264  // considered point in the previous mapping
265 
266  int l0 ;
267  double xi0 ;
268  mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
269 
270  // Value of uu at this point
271 
272  *ptx = va_cf.val_point_jk(l0, xi0, j, k) ;
273 
274  // next point
275  pr++ ;
276  ptx++ ;
277  }
278  }
279  }
280 
281 
282  } // End of the loop on the domains where the evaluation had to be performed
283 
284  // In the remaining domains, uu is set to zero:
285  // -------------------------------------------
286  if (nzet < nz)
287  uu.annule(nzet, nz - 1) ;
288 
289 
290 
291 
292 }
293 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:67
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
Mtbl * c
The coordinate values at each grid point.
Definition: coord.h:97
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
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...
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.
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Parameter storage.
Definition: param.h:125
Base class for pure radial mappings.
Definition: map.h:1557
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
void fait() const
Computes, at each point of the grid, the value of the coordinate or mapping derivative represented by...
Definition: coord.C:119
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...
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:694
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:704
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
Coord r
r coordinate centered on the grid
Definition: map.h:736