LORENE
map_radial_reeval_symy.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_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 $" ;
24
25/*
26 * $Id: map_radial_reeval_symy.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
27 * $Log: map_radial_reeval_symy.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/03/06 11:30:19 eric
41 * *** empty log message ***
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Map/map_radial_reeval_symy.C,v 1.4 2014/10/13 08:53:06 j_novak Exp $
45 *
46 */
47
48
49// Headers C
50#include <cstdlib>
51
52// Headers Lorene
53#include "map.h"
54#include "cmp.h"
55#include "param.h"
56#include "scalar.h"
57
58namespace Lorene {
59void Map_radial::reevaluate_symy(const Map* mp_prev0, int nzet, Cmp& uu) const {
60
61 const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
62
63 if (mp_prev == 0x0) {
64 cout <<
65 "Map_radial::reevaluate_symy : the mapping mp_prev does not belong"
66 << endl ;
67 cout << " to the class Map_radial !" << endl ;
68 abort() ;
69 }
70
71 int nz = mg->get_nzone() ;
72
73 // Protections
74 // -----------
75 assert(uu.get_mp() == this) ;
76 assert(uu.get_dzpuis() == 0) ;
77 assert(uu.get_etat() != ETATNONDEF) ;
78 assert(mp_prev->mg == mg) ;
79 assert(nzet <= nz) ;
80 assert(mg->get_type_p() == NONSYM) ;
81
82
83 // Maybe nothing to do ?
84 if ( uu.get_etat() == ETATZERO ) {
85 return ;
86 }
87
88 assert(uu.get_etat() == ETATQCQ) ;
89 (uu.va).coef() ; // the coefficients of uu are required
90
91 // Copy of the coefficients
92 Mtbl_cf va_cf = *((uu.va).c_cf) ;
93
94 // Preparation of uu.va for the storage of the new values.
95 // ------------------------------------------------------
96
97 (uu.va).set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
98 // if it does not exist already
99 // Destroys the coefficients
100
101 Mtbl& va_c = *((uu.va).c) ; // Reference to the Mtbl uu.va.c
102
103 va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
104 // domain if they do not exist already
105
106
107 // Values of the Coord r
108 // ---------------------
109
110 if ( r.c == 0x0 ) r.fait() ;
111 const Mtbl& rc = *(r.c) ;
112
113 // Precision for val_lx_jk :
114 // ------------------------
115 int nitermax = 100 ; // Maximum number of iterations in the secant method
116 int niter ; // Number of iterations effectively used
117 double precis = 1.e-15 ; // Absolute precision in the secant method
118 Param par ;
119 par.add_int(nitermax) ;
120 par.add_int_mod(niter) ;
121 par.add_double(precis) ;
122
123 // Loop on the domains where the evaluation is to be performed
124 // -----------------------------------------------------------
125
126 for (int l=0; l<nzet; l++) {
127 int nr = mg->get_nr(l) ;
128 int nt = mg->get_nt(l) ;
129 int np = mg->get_np(l) ;
130
131 va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
132 // store the result
133
134 double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
135
136 double* pr = (rc.t[l])->t ; // Pointer on the values of r
137
138 // Loop on half the grid points in the considered arrival domain
139 // (the other half will be obtained by symmetry with respect to
140 // the y=0 plane).
141
142 for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
143 for (int j=0; j<nt; j++) {
144 for (int i=0; i<nr; i++) {
145
146 // domain l0 and value of the coordinate xi0 of the
147 // considered point in the previous mapping
148
149 int l0 ;
150 double xi0 ;
151 mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
152
153 // Value of uu at this point
154
155 *ptx = va_cf.val_point_jk_symy(l0, xi0, j, k) ;
156
157 // next point
158 pr++ ;
159 ptx++ ;
160 }
161 }
162 }
163
164 // The remaining points are obtained by symmetry with rspect to the
165 // y=0 plane
166
167 for (int k=np/2+1 ; k<np ; k++) {
168
169 // pointer on the value (already computed) at the point symmetric
170 // with respect to the plane y=0
171 double* ptx_symy = (va_c.t[l])->t + (np-k)*nt*nr ;
172
173 // copy :
174 for (int j=0 ; j<nt ; j++) {
175 for (int i=0 ; i<nr ; i++) {
176 *ptx = *ptx_symy ;
177 ptx++ ;
178 ptx_symy++ ;
179 }
180 }
181 }
182
183 } // End of the loop on the domains where the evaluation had to be performed
184
185 // In the remaining domains, uu is set to zero:
186 // -------------------------------------------
187
188 uu.annule(nzet, nz - 1) ;
189
190
191}
192
193void Map_radial::reevaluate_symy(const Map* mp_prev0, int nzet, Scalar& uu) const {
194
195 const Map_radial* mp_prev = dynamic_cast<const Map_radial*>(mp_prev0) ;
196
197 if (mp_prev == 0x0) {
198 cout <<
199 "Map_radial::reevaluate_symy : the mapping mp_prev does not belong"
200 << endl ;
201 cout << " to the class Map_radial !" << endl ;
202 abort() ;
203 }
204
205 int nz = mg->get_nzone() ;
206
207 // Protections
208 // -----------
209 assert(uu.get_mp() == *this) ;
210 assert(uu.get_dzpuis() == 0) ;
211 assert(uu.get_etat() != ETATNONDEF) ;
212 assert(mp_prev->mg == mg) ;
213 assert(nzet <= nz) ;
214 assert(mg->get_type_p() == NONSYM) ;
215
216
217 // Maybe nothing to do ?
218 if ( uu.get_etat() == ETATZERO ) {
219 return ;
220 }
221
222 assert(uu.get_etat() == ETATQCQ) ;
223 uu.set_spectral_va().coef() ; // the coefficients of uu are required
224
225 // Copy of the coefficients
226 Mtbl_cf va_cf = *(uu.set_spectral_va().c_cf) ;
227
228 // Preparation of uu.va for the storage of the new values.
229 // ------------------------------------------------------
230
231 uu.set_spectral_va().set_etat_c_qcq() ; // Allocates the memory for the Mtbl uu.va.c
232 // if it does not exist already
233 // Destroys the coefficients
234
235 Mtbl& va_c = *(uu.set_spectral_va().c) ; // Reference to the Mtbl uu.va.c
236
237 va_c.set_etat_qcq() ; // Allocates the memory for the Tbl's in each
238 // domain if they do not exist already
239
240
241 // Values of the Coord r
242 // ---------------------
243
244 if ( r.c == 0x0 ) r.fait() ;
245 const Mtbl& rc = *(r.c) ;
246
247 // Precision for val_lx_jk :
248 // ------------------------
249 int nitermax = 100 ; // Maximum number of iterations in the secant method
250 int niter ; // Number of iterations effectively used
251 double precis = 1.e-15 ; // Absolute precision in the secant method
252 Param par ;
253 par.add_int(nitermax) ;
254 par.add_int_mod(niter) ;
255 par.add_double(precis) ;
256
257 // Loop on the domains where the evaluation is to be performed
258 // -----------------------------------------------------------
259
260 for (int l=0; l<nzet; l++) {
261 int nr = mg->get_nr(l) ;
262 int nt = mg->get_nt(l) ;
263 int np = mg->get_np(l) ;
264
265 va_c.t[l]->set_etat_qcq() ; // Allocates the array of double to
266 // store the result
267
268 double* ptx = (va_c.t[l])->t ; // Pointer on the allocated array
269
270 double* pr = (rc.t[l])->t ; // Pointer on the values of r
271
272 // Loop on half the grid points in the considered arrival domain
273 // (the other half will be obtained by symmetry with respect to
274 // the y=0 plane).
275
276 for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
277 for (int j=0; j<nt; j++) {
278 for (int i=0; i<nr; i++) {
279
280 // domain l0 and value of the coordinate xi0 of the
281 // considered point in the previous mapping
282
283 int l0 ;
284 double xi0 ;
285 mp_prev->val_lx_jk(*pr, j, k, par, l0, xi0) ;
286
287 // Value of uu at this point
288
289 *ptx = va_cf.val_point_jk_symy(l0, xi0, j, k) ;
290
291 // next point
292 pr++ ;
293 ptx++ ;
294 }
295 }
296 }
297
298 // The remaining points are obtained by symmetry with rspect to the
299 // y=0 plane
300
301 for (int k=np/2+1 ; k<np ; k++) {
302
303 // pointer on the value (already computed) at the point symmetric
304 // with respect to the plane y=0
305 double* ptx_symy = (va_c.t[l])->t + (np-k)*nt*nr ;
306
307 // copy :
308 for (int j=0 ; j<nt ; j++) {
309 for (int i=0 ; i<nr ; i++) {
310 *ptx = *ptx_symy ;
311 ptx++ ;
312 ptx_symy++ ;
313 }
314 }
315 }
316
317 } // End of the loop on the domains where the evaluation had to be performed
318
319 // In the remaining domains, uu is set to zero:
320 // -------------------------------------------
321
322 uu.annule(nzet, nz - 1) ;
323
324
325}
326}
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_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.
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_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
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_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...
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