LORENE
cmp_import_symy.C
1 /*
2  * Member function of the Cmp class for initiating a Cmp from a Cmp defined
3  * on another mapping.
4  * Case where both Cmp's are symmetric with respect to their y=0 plane.
5  */
6 
7 /*
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 
32 /*
33  * $Id: cmp_import_symy.C,v 1.4 2016/12/05 16:17:48 j_novak Exp $
34  * $Log: cmp_import_symy.C,v $
35  * Revision 1.4 2016/12/05 16:17:48 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.3 2014/10/13 08:52:47 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.2 2014/10/06 15:13:03 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
45  * LORENE
46  *
47  * Revision 2.0 2000/03/06 10:56:16 eric
48  * *** empty log message ***
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_import_symy.C,v 1.4 2016/12/05 16:17:48 j_novak Exp $
52  *
53  */
54 
55 
56 
57 // Headers C
58 #include <cmath>
59 
60 // Headers Lorene
61 #include "cmp.h"
62 #include "param.h"
63 #include "nbr_spx.h"
64 
65  //-------------------------------//
66  // Importation in all domains //
67  //-------------------------------//
68 
69 namespace Lorene {
70 void Cmp::import_symy(const Cmp& ci) {
71 
72  int nz = mp->get_mg()->get_nzone() ;
73 
74  import_symy(nz, ci) ;
75 
76 }
77 
78  //--------------------------------------//
79  // Importation in inner domains only //
80  //--------------------------------------//
81 
82 void Cmp::import_symy(int nzet, const Cmp& cm_d) {
83 
84  const Map* mp_d = cm_d.mp ; // Departure mapping
85 
86  // Trivial case : mappings identical !
87  // -----------------------------------
88 
89  if (mp_d == mp) {
90  *this = cm_d ;
91  return ;
92  }
93 
94  // Relative orientation of the two mappings
95  // ----------------------------------------
96 
97  int align_rel = (mp->get_bvect_cart()).get_align()
98  * (mp_d->get_bvect_cart()).get_align() ;
99 
100  switch (align_rel) {
101 
102  case 1 : { // the two mappings have aligned Cartesian axis
103  import_align_symy(nzet, cm_d) ;
104  break ;
105  }
106 
107  case -1 : { // the two mappings have anti-aligned Cartesian axis
108  import_anti_symy(nzet, cm_d) ;
109  break ;
110  }
111 
112  default : {
113  cout << "Cmp::import_symy : unexpected value of align_rel : "
114  << align_rel << endl ;
115  abort() ;
116  break ;
117  }
118 
119  }
120 
121 }
122 
123 
124  //-----------------------------------------//
125  // Case of Cartesian axis anti-aligned //
126  //-----------------------------------------//
127 
128 
129 void Cmp::import_anti_symy(int nzet, const Cmp& cm_d) {
130 
131  // Trivial case : null Cmp
132  // ------------------------
133 
134  if (cm_d.etat == ETATZERO) {
135  set_etat_zero() ;
136  return ;
137  }
138 
139  const Map* mp_d = cm_d.mp ; // Departure mapping
140 
141  // Protections
142  // -----------
143  int align = (mp->get_bvect_cart()).get_align() ;
144 
145  assert( align * (mp_d->get_bvect_cart()).get_align() == -1 ) ;
146 
147  assert(cm_d.etat == ETATQCQ) ;
148 
149  if (cm_d.dzpuis != 0) {
150  cout <<
151  "Cmp::import_anti_symy : the dzpuis of the Cmp to be imported"
152  << " must be zero !" << endl ;
153  abort() ;
154  }
155 
156 
157  const Mg3d* mg_a = mp->get_mg() ;
158  assert(mg_a->get_type_p() == NONSYM) ;
159 
160  int nz_a = mg_a->get_nzone() ;
161  assert(nzet <= nz_a) ;
162 
163  const Valeur& va_d = cm_d.va ;
164  va_d.coef() ; // The coefficients are required
165 
166 
167  // Preparations for storing the result in *this
168  // --------------------------------------------
169  del_t() ; // delete all previously computed derived quantities
170 
171  set_etat_qcq() ; // Set the state to ETATQCQ
172 
173  va.set_etat_c_qcq() ; // Allocates the memory for the Mtbl va.c
174  // if it does not exist already
175  va.c->set_etat_qcq() ; // Allocates the memory for the Tbl's in each
176  // domain if they do not exist already
177 
178 
179  // Departure (x,y,z) coordinates of the origin of the Arrival mapping :
180 
181  double xx_a, yy_a, zz_a ;
182  if (align == 1) {
183  xx_a = mp_d->get_ori_x() - mp->get_ori_x() ;
184  yy_a = mp_d->get_ori_y() - mp->get_ori_y() ;
185  }
186  else {
187  xx_a = mp->get_ori_x() - mp_d->get_ori_x() ;
188  yy_a = mp->get_ori_y() - mp_d->get_ori_y() ;
189  }
190  zz_a = mp->get_ori_z() - mp_d->get_ori_z() ;
191 
192 
193  // r, theta, phi, x, y and z on the Arrival mapping
194  // update of the corresponding Coord's if necessary
195 
196  if ( (mp->r).c == 0x0 ) (mp->r).fait() ;
197  if ( (mp->tet).c == 0x0 ) (mp->tet).fait() ;
198  if ( (mp->phi).c == 0x0 ) (mp->phi).fait() ;
199  if ( (mp->x).c == 0x0 ) (mp->x).fait() ;
200  if ( (mp->y).c == 0x0 ) (mp->y).fait() ;
201  if ( (mp->z).c == 0x0 ) (mp->z).fait() ;
202 
203  const Mtbl* mr_a = (mp->r).c ;
204  const Mtbl* mtet_a = (mp->tet).c ;
205  const Mtbl* mphi_a = (mp->phi).c ;
206  const Mtbl* mx_a = (mp->x).c ;
207  const Mtbl* my_a = (mp->y).c ;
208  const Mtbl* mz_a = (mp->z).c ;
209 
210  Param par_precis ; // Required precision in the method Map::val_lx
211  int nitermax = 100 ; // Maximum number of iteration in the secant method
212  int niter ;
213  double precis = 1e-15 ; // Absolute precision in the secant method
214  par_precis.add_int(nitermax) ;
215  par_precis.add_int_mod(niter) ;
216  par_precis.add_double(precis) ;
217 
218 
219  // Loop of the Arrival domains where the computation is to be performed
220  // --------------------------------------------------------------------
221 
222  for (int l=0; l < nzet; l++) {
223 
224  int nr = mg_a->get_nr(l) ;
225  int nt = mg_a->get_nt(l) ;
226  int np = mg_a->get_np(l) ;
227 
228 
229  const double* pr_a = mr_a->t[l]->t ; // Pointer on the values of r
230  const double* ptet_a = mtet_a->t[l]->t ; // Pointer on the values of theta
231  const double* pphi_a = mphi_a->t[l]->t ; // Pointer on the values of phi
232  const double* px_a = mx_a->t[l]->t ; // Pointer on the values of X
233  const double* py_a = my_a->t[l]->t ; // Pointer on the values of Y
234  const double* pz_a = mz_a->t[l]->t ; // Pointer on the values of Z
235 
236  (va.c->t[l])->set_etat_qcq() ; // Allocates the array of double to
237  // store the result
238  double* ptx = (va.c->t[l])->t ; // Pointer on the allocated array
239 
240 
241  // Loop on half the grid points in the considered arrival domain
242  // (the other half will be obtained by symmetry with respect to
243  // the y=0 plane).
244 
245  for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
246  for (int j=0 ; j<nt ; j++) {
247  for (int i=0 ; i<nr ; i++) {
248 
249  double r = *pr_a ;
250  double rd, tetd, phid ;
251  if (r == __infinity) {
252  rd = r ;
253  tetd = *ptet_a ;
254  phid = *pphi_a + M_PI ;
255  if (phid < 0) phid += 2*M_PI ;
256  }
257  else {
258 
259  // Cartesian coordinates on the Departure mapping
260  double xd = - *px_a + xx_a ;
261  double yd = - *py_a + yy_a ;
262  double zd = *pz_a + zz_a ;
263 
264  // Spherical coordinates on the Departure mapping
265  double rhod2 = xd*xd + yd*yd ;
266  double rhod = sqrt( rhod2 ) ;
267  rd = sqrt(rhod2 + zd*zd) ;
268  tetd = atan2(rhod, zd) ;
269  phid = atan2(yd, xd) ;
270  if (phid < 0) phid += 2*M_PI ;
271  }
272 
273 
274  // NB: to increase the efficiency, the method Cmp::val_point
275  // is not invoked; the method Mtbl_cf::val_point is
276  // called directly instead.
277 
278  // Value of the grid coordinates (l,xi) corresponding to
279  // (rd,tetd,phid) :
280 
281  int ld ; // domain index
282  double xxd ; // radial coordinate xi in [0,1] or [-1,1]
283  mp_d->val_lx(rd, tetd, phid, par_precis, ld, xxd) ;
284 
285  // Value of the Departure Cmp at the obtained point:
286  *ptx = va_d.c_cf->val_point_symy(ld, xxd, tetd, phid) ;
287 
288  // Next point :
289  ptx++ ;
290  pr_a++ ;
291  ptet_a++ ;
292  pphi_a++ ;
293  px_a++ ;
294  py_a++ ;
295  pz_a++ ;
296 
297  }
298  }
299  }
300 
301  // The remaining points are obtained by symmetry with rspect to the
302  // y=0 plane
303 
304  for (int k=np/2+1 ; k<np ; k++) {
305 
306  // pointer on the value (already computed) at the point symmetric
307  // with respect to the plane y=0
308  double* ptx_symy = (va.c->t[l])->t + (np-k)*nt*nr ;
309 
310  // copy :
311  for (int j=0 ; j<nt ; j++) {
312  for (int i=0 ; i<nr ; i++) {
313  *ptx = *ptx_symy ;
314  ptx++ ;
315  ptx_symy++ ;
316  }
317  }
318  }
319 
320 
321  } // End of the loop on the Arrival domains
322 
323  // In the remaining domains, *this is set to zero:
324  // ----------------------------------------------
325 
326  if (nzet < nz_a) {
327  annule(nzet, nz_a - 1) ;
328  }
329 
330  // Treatment of dzpuis
331  // -------------------
332 
333  set_dzpuis(0) ;
334 
335 }
336 
337 
338  //-------------------------------------//
339  // Case of aligned Cartesian axis //
340  //-------------------------------------//
341 
342 
343 void Cmp::import_align_symy(int nzet, const Cmp& cm_d) {
344 
345  // Trivial case : null Cmp
346  // ------------------------
347 
348  if (cm_d.etat == ETATZERO) {
349  set_etat_zero() ;
350  return ;
351  }
352 
353  const Map* mp_d = cm_d.mp ; // Departure mapping
354 
355  // Protections
356  // -----------
357  int align = (mp->get_bvect_cart()).get_align() ;
358 
359  assert( align * (mp_d->get_bvect_cart()).get_align() == 1 ) ;
360 
361  assert(cm_d.etat == ETATQCQ) ;
362 
363  if (cm_d.dzpuis != 0) {
364  cout <<
365  "Cmp::import_align_symy : the dzpuis of the Cmp to be imported"
366  << " must be zero !" << endl ;
367  abort() ;
368  }
369 
370 
371  const Mg3d* mg_a = mp->get_mg() ;
372  assert(mg_a->get_type_p() == NONSYM) ;
373 
374  int nz_a = mg_a->get_nzone() ;
375  assert(nzet <= nz_a) ;
376 
377  const Valeur& va_d = cm_d.va ;
378  va_d.coef() ; // The coefficients are required
379 
380 
381  // Preparations for storing the result in *this
382  // --------------------------------------------
383  del_t() ; // delete all previously computed derived quantities
384 
385  set_etat_qcq() ; // Set the state to ETATQCQ
386 
387  va.set_etat_c_qcq() ; // Allocates the memory for the Mtbl va.c
388  // if it does not exist already
389  va.c->set_etat_qcq() ; // Allocates the memory for the Tbl's in each
390  // domain if they do not exist already
391 
392 
393  // Departure (x,y,z) coordinates of the origin of the Arrival mapping :
394 
395  double xx_a, yy_a, zz_a ;
396  if (align == 1) {
397  xx_a = mp->get_ori_x() - mp_d->get_ori_x() ;
398  yy_a = mp->get_ori_y() - mp_d->get_ori_y() ;
399  }
400  else {
401  xx_a = mp_d->get_ori_x() - mp->get_ori_x() ;
402  yy_a = mp_d->get_ori_y() - mp->get_ori_y() ;
403  }
404  zz_a = mp->get_ori_z() - mp_d->get_ori_z() ;
405 
406 
407  // r, theta, phi, x, y and z on the Arrival mapping
408  // update of the corresponding Coord's if necessary
409 
410  if ( (mp->r).c == 0x0 ) (mp->r).fait() ;
411  if ( (mp->tet).c == 0x0 ) (mp->tet).fait() ;
412  if ( (mp->phi).c == 0x0 ) (mp->phi).fait() ;
413  if ( (mp->x).c == 0x0 ) (mp->x).fait() ;
414  if ( (mp->y).c == 0x0 ) (mp->y).fait() ;
415  if ( (mp->z).c == 0x0 ) (mp->z).fait() ;
416 
417  const Mtbl* mr_a = (mp->r).c ;
418  const Mtbl* mtet_a = (mp->tet).c ;
419  const Mtbl* mphi_a = (mp->phi).c ;
420  const Mtbl* mx_a = (mp->x).c ;
421  const Mtbl* my_a = (mp->y).c ;
422  const Mtbl* mz_a = (mp->z).c ;
423 
424  Param par_precis ; // Required precision in the method Map::val_lx
425  int nitermax = 100 ; // Maximum number of iteration in the secant method
426  int niter ;
427  double precis = 1e-15 ; // Absolute precision in the secant method
428  par_precis.add_int(nitermax) ;
429  par_precis.add_int_mod(niter) ;
430  par_precis.add_double(precis) ;
431 
432 
433  // Loop of the Arrival domains where the computation is to be performed
434  // --------------------------------------------------------------------
435 
436  for (int l=0; l < nzet; l++) {
437 
438  int nr = mg_a->get_nr(l) ;
439  int nt = mg_a->get_nt(l) ;
440  int np = mg_a->get_np(l) ;
441 
442 
443  const double* pr_a = mr_a->t[l]->t ; // Pointer on the values of r
444  const double* ptet_a = mtet_a->t[l]->t ; // Pointer on the values of theta
445  const double* pphi_a = mphi_a->t[l]->t ; // Pointer on the values of phi
446  const double* px_a = mx_a->t[l]->t ; // Pointer on the values of X
447  const double* py_a = my_a->t[l]->t ; // Pointer on the values of Y
448  const double* pz_a = mz_a->t[l]->t ; // Pointer on the values of Z
449 
450  (va.c->t[l])->set_etat_qcq() ; // Allocates the array of double to
451  // store the result
452  double* ptx = (va.c->t[l])->t ; // Pointer on the allocated array
453 
454 
455 
456  // Loop on half the grid points in the considered arrival domain
457  // (the other half will be obtained by symmetry with respect to
458  // the y=0 plane).
459 
460  for (int k=0 ; k<np/2+1 ; k++) { // np/2+1 : half the grid
461  for (int j=0 ; j<nt ; j++) {
462  for (int i=0 ; i<nr ; i++) {
463 
464  double r = *pr_a ;
465  double rd, tetd, phid ;
466  if (r == __infinity) {
467  rd = r ;
468  tetd = *ptet_a ;
469  phid = *pphi_a ;
470  }
471  else {
472 
473  // Cartesian coordinates on the Departure mapping
474  double xd = *px_a + xx_a ;
475  double yd = *py_a + yy_a ;
476  double zd = *pz_a + zz_a ;
477 
478  // Spherical coordinates on the Departure mapping
479  double rhod2 = xd*xd + yd*yd ;
480  double rhod = sqrt( rhod2 ) ;
481  rd = sqrt(rhod2 + zd*zd) ;
482  tetd = atan2(rhod, zd) ;
483  phid = atan2(yd, xd) ;
484  if (phid < 0) phid += 2*M_PI ;
485  }
486 
487 
488  // NB: to increase the efficiency, the method Cmp::val_point
489  // is not invoked; the method Mtbl_cf::val_point is
490  // called directly instead.
491 
492  // Value of the grid coordinates (l,xi) corresponding to
493  // (rd,tetd,phid) :
494 
495  int ld ; // domain index
496  double xxd ; // radial coordinate xi in [0,1] or [-1,1]
497  mp_d->val_lx(rd, tetd, phid, par_precis, ld, xxd) ;
498 
499  // Value of the Departure Cmp at the obtained point:
500  *ptx = va_d.c_cf->val_point_symy(ld, xxd, tetd, phid) ;
501 
502  // Next point :
503  ptx++ ;
504  pr_a++ ;
505  ptet_a++ ;
506  pphi_a++ ;
507  px_a++ ;
508  py_a++ ;
509  pz_a++ ;
510 
511  }
512  }
513  }
514 
515 
516  // The remaining points are obtained by symmetry with rspect to the
517  // y=0 plane
518 
519  for (int k=np/2+1 ; k<np ; k++) {
520 
521  // pointer on the value (already computed) at the point symmetric
522  // with respect to the plane y=0
523  double* ptx_symy = (va.c->t[l])->t + (np-k)*nt*nr ;
524 
525  // copy :
526  for (int j=0 ; j<nt ; j++) {
527  for (int i=0 ; i<nr ; i++) {
528  *ptx = *ptx_symy ;
529  ptx++ ;
530  ptx_symy++ ;
531  }
532  }
533  }
534 
535  } // End of the loop on the Arrival domains
536 
537  // In the remaining domains, *this is set to zero:
538  // ----------------------------------------------
539 
540  if (nzet < nz_a) {
541  annule(nzet, nz_a - 1) ;
542  }
543 
544  // Treatment of dzpuis
545  // -------------------
546 
547  set_dzpuis(0) ;
548 
549 }
550 }
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:512
void del_t()
Logical destructor.
Definition: cmp.C:262
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
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
Multi-domain array.
Definition: mtbl.h:118
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:788
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
Base class for coordinate mappings.
Definition: map.h:688
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Coord tet
coordinate centered on the grid
Definition: map.h:737
Coord phi
coordinate centered on the grid
Definition: map.h:738
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:292
void import_align_symy(int nzet, const Cmp &ci)
Assignment to another Cmp defined on a different mapping, when the two mappings have aligned Cartesia...
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified zo...
Definition: cmp.h:461
double * t
The array of double.
Definition: tbl.h:176
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Parameter storage.
Definition: param.h:125
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: cmp.C:307
void import_symy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
double val_point_symy(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
const Map * mp
Reference mapping.
Definition: cmp.h:451
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Coord y
y coordinate centered on the grid
Definition: map.h:745
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
void import_anti_symy(int nzet, const Cmp &ci)
Assignment to another Cmp defined on a different mapping, when the two mappings have anti-aligned Car...
Coord x
x coordinate centered on the grid
Definition: map.h:744
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
double get_ori_z() const
Returns the z coordinate of the origin.
Definition: map.h:790
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:657
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Coord z
z coordinate centered on the grid
Definition: map.h:746
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
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point given by its physical coordina...