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
69namespace Lorene {
70void 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
82void 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
129void 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
343void 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}
Lorene::Cmp::mp
const Map * mp
Reference mapping.
Definition cmp.h:451
Lorene::Cmp::Cmp
Cmp(const Map &map)
Constructor from mapping.
Definition cmp.C:211
Lorene::Cmp::import_align_symy
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...
Definition cmp_import_symy.C:343
Lorene::Cmp::dzpuis
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified z...
Definition cmp.h:461
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition cmp.h:454
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:351
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:307
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:292
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:657
Lorene::Cmp::import_symy
void import_symy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
Definition cmp_import_symy.C:70
Lorene::Cmp::del_t
void del_t()
Logical destructor.
Definition cmp.C:262
Lorene::Cmp::import_anti_symy
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...
Definition cmp_import_symy.C:129
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:479
Lorene::Mg3d::get_type_p
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:512
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:474
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:469
Lorene::Mtbl_cf::val_point_symy
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 ...
Definition mtbl_cf_vp_symy.C:72
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:318
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:388
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:249
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:151
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:730