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