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 
23 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 $" ;
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 
58 namespace Lorene {
59 void 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 
193 void 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
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
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.
Definition: valeur_coef.C:148
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene prototypes.
Definition: app_hor.h:64