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
23char 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
57namespace Lorene {
58void 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
173void 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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
void fait() const
Computes, at each point of the grid, the value of the coordinate or mapping derivative represented by...
Definition coord.C:116
Mtbl * c
The coordinate values at each grid point.
Definition coord.h:97
Base class for pure radial mappings.
Definition map.h:1536
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.
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...
Base class for coordinate mappings.
Definition map.h:670
Coord r
r coordinate centered on the grid
Definition map.h:718
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
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...
Multi-domain array.
Definition mtbl.h:118
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Parameter storage.
Definition param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
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
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64