LORENE
cmp_import_asymy.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 antisymmetric 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_asymy.C,v 1.4 2016/12/05 16:17:48 j_novak Exp $
34 * $Log: cmp_import_asymy.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:07 eric
48 * *** empty log message ***
49 *
50 *
51 * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_import_asymy.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_asymy(const Cmp& ci) {
71
72 int nz = mp->get_mg()->get_nzone() ;
73
74 import_asymy(nz, ci) ;
75
76}
77
78 //--------------------------------------//
79 // Importation in inner domains only //
80 //--------------------------------------//
81
82void Cmp::import_asymy(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_asymy(nzet, cm_d) ;
104 break ;
105 }
106
107 case -1 : { // the two mappings have anti-aligned Cartesian axis
108 import_anti_asymy(nzet, cm_d) ;
109 break ;
110 }
111
112 default : {
113 cout << "Cmp::import_asymy : 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_asymy(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_asymy : 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
161 int nz_a = mg_a->get_nzone() ;
162 assert(nzet <= nz_a) ;
163
164 const Valeur& va_d = cm_d.va ;
165 va_d.coef() ; // The coefficients are required
166
167
168 // Preparations for storing the result in *this
169 // --------------------------------------------
170 del_t() ; // delete all previously computed derived quantities
171
172 set_etat_qcq() ; // Set the state to ETATQCQ
173
174 va.set_etat_c_qcq() ; // Allocates the memory for the Mtbl va.c
175 // if it does not exist already
176 va.c->set_etat_qcq() ; // Allocates the memory for the Tbl's in each
177 // domain if they do not exist already
178
179
180 // Departure (x,y,z) coordinates of the origin of the Arrival mapping :
181
182 double xx_a, yy_a, zz_a ;
183 if (align == 1) {
184 xx_a = mp_d->get_ori_x() - mp->get_ori_x() ;
185 yy_a = mp_d->get_ori_y() - mp->get_ori_y() ;
186 }
187 else {
188 xx_a = mp->get_ori_x() - mp_d->get_ori_x() ;
189 yy_a = mp->get_ori_y() - mp_d->get_ori_y() ;
190 }
191 zz_a = mp->get_ori_z() - mp_d->get_ori_z() ;
192
193
194 // r, theta, phi, x, y and z on the Arrival mapping
195 // update of the corresponding Coord's if necessary
196
197 if ( (mp->r).c == 0x0 ) (mp->r).fait() ;
198 if ( (mp->tet).c == 0x0 ) (mp->tet).fait() ;
199 if ( (mp->phi).c == 0x0 ) (mp->phi).fait() ;
200 if ( (mp->x).c == 0x0 ) (mp->x).fait() ;
201 if ( (mp->y).c == 0x0 ) (mp->y).fait() ;
202 if ( (mp->z).c == 0x0 ) (mp->z).fait() ;
203
204 const Mtbl* mr_a = (mp->r).c ;
205 const Mtbl* mtet_a = (mp->tet).c ;
206 const Mtbl* mphi_a = (mp->phi).c ;
207 const Mtbl* mx_a = (mp->x).c ;
208 const Mtbl* my_a = (mp->y).c ;
209 const Mtbl* mz_a = (mp->z).c ;
210
211 Param par_precis ; // Required precision in the method Map::val_lx
212 int nitermax = 100 ; // Maximum number of iteration in the secant method
213 int niter ;
214 double precis = 1e-15 ; // Absolute precision in the secant method
215 par_precis.add_int(nitermax) ;
216 par_precis.add_int_mod(niter) ;
217 par_precis.add_double(precis) ;
218
219
220 // Loop of the Arrival domains where the computation is to be performed
221 // --------------------------------------------------------------------
222
223 for (int l=0; l < nzet; l++) {
224
225 int nr = mg_a->get_nr(l) ;
226 int nt = mg_a->get_nt(l) ;
227 int np = mg_a->get_np(l) ;
228 int ntnr = nt*nr ;
229
230 const double* pr_a = mr_a->t[l]->t ; // Pointer on the values of r
231 const double* ptet_a = mtet_a->t[l]->t ; // Pointer on the values of theta
232 const double* pphi_a = mphi_a->t[l]->t ; // Pointer on the values of phi
233 const double* px_a = mx_a->t[l]->t ; // Pointer on the values of X
234 const double* py_a = my_a->t[l]->t ; // Pointer on the values of Y
235 const double* pz_a = mz_a->t[l]->t ; // Pointer on the values of Z
236
237 (va.c->t[l])->set_etat_qcq() ; // Allocates the array of double to
238 // store the result
239 double* ptx = (va.c->t[l])->t ; // Pointer on the allocated array
240
241
242 // Loop on half the grid points in the considered arrival domain
243 // (the other half will be obtained by antisymmetry with respect to
244 // the y=0 plane).
245
246 // Case k=0 (phi=0) : the function is zero (by antisymmetry)
247 for (int i=0; i<ntnr; i++) {
248 *ptx = 0 ;
249 ptx++ ; // next point
250 }
251
252 // Go to k=1 :
253 pr_a += ntnr ;
254 ptet_a += ntnr ;
255 pphi_a += ntnr ;
256 px_a += ntnr ;
257 py_a += ntnr ;
258 pz_a += ntnr ;
259
260 for (int k=1 ; k<np/2 ; k++) { // np/2 : ~ half the grid
261 for (int j=0 ; j<nt ; j++) {
262 for (int i=0 ; i<nr ; i++) {
263
264 double r = *pr_a ;
265 double rd, tetd, phid ;
266 if (r == __infinity) {
267 rd = r ;
268 tetd = *ptet_a ;
269 phid = *pphi_a + M_PI ;
270 if (phid < 0) phid += 2*M_PI ;
271 }
272 else {
273
274 // Cartesian coordinates on the Departure mapping
275 double xd = - *px_a + xx_a ;
276 double yd = - *py_a + yy_a ;
277 double zd = *pz_a + zz_a ;
278
279 // Spherical coordinates on the Departure mapping
280 double rhod2 = xd*xd + yd*yd ;
281 double rhod = sqrt( rhod2 ) ;
282 rd = sqrt(rhod2 + zd*zd) ;
283 tetd = atan2(rhod, zd) ;
284 phid = atan2(yd, xd) ;
285 if (phid < 0) phid += 2*M_PI ;
286 }
287
288
289 // NB: to increase the efficiency, the method Cmp::val_point
290 // is not invoked; the method Mtbl_cf::val_point is
291 // called directly instead.
292
293 // Value of the grid coordinates (l,xi) corresponding to
294 // (rd,tetd,phid) :
295
296 int ld ; // domain index
297 double xxd ; // radial coordinate xi in [0,1] or [-1,1]
298 mp_d->val_lx(rd, tetd, phid, par_precis, ld, xxd) ;
299
300 // Value of the Departure Cmp at the obtained point:
301 *ptx = va_d.c_cf->val_point_asymy(ld, xxd, tetd, phid) ;
302
303 // Next point :
304 ptx++ ;
305 pr_a++ ;
306 ptet_a++ ;
307 pphi_a++ ;
308 px_a++ ;
309 py_a++ ;
310 pz_a++ ;
311
312 }
313 }
314 }
315
316 // Case k=np/2 (phi=pi) : the function is zero (by antisymmetry)
317 for (int i=0; i<ntnr; i++) {
318 *ptx = 0 ;
319 ptx++ ; // next point
320 }
321
322 // Go to k=np/2+1 :
323 pr_a += ntnr ;
324 ptet_a += ntnr ;
325 pphi_a += ntnr ;
326 px_a += ntnr ;
327 py_a += ntnr ;
328 pz_a += ntnr ;
329
330 // The remaining points are obtained by antisymmetry with rspect to the
331 // y=0 plane
332
333 for (int k=np/2+1 ; k<np ; k++) {
334
335 // pointer on the value (already computed) at the point symmetric
336 // with respect to the plane y=0
337 double* ptx_symy = (va.c->t[l])->t + (np-k)*nt*nr ;
338
339 // copy :
340 for (int j=0 ; j<nt ; j++) {
341 for (int i=0 ; i<nr ; i++) {
342 *ptx = - (*ptx_symy) ;
343 ptx++ ;
344 ptx_symy++ ;
345 }
346 }
347 }
348
349
350 } // End of the loop on the Arrival domains
351
352 // In the remaining domains, *this is set to zero:
353 // ----------------------------------------------
354
355 if (nzet < nz_a) {
356 annule(nzet, nz_a - 1) ;
357 }
358
359 // Treatment of dzpuis
360 // -------------------
361
362 set_dzpuis(0) ;
363
364}
365
366
367 //-------------------------------------//
368 // Case of aligned Cartesian axis //
369 //-------------------------------------//
370
371
372void Cmp::import_align_asymy(int nzet, const Cmp& cm_d) {
373
374 // Trivial case : null Cmp
375 // ------------------------
376
377 if (cm_d.etat == ETATZERO) {
378 set_etat_zero() ;
379 return ;
380 }
381
382 const Map* mp_d = cm_d.mp ; // Departure mapping
383
384 // Protections
385 // -----------
386 int align = (mp->get_bvect_cart()).get_align() ;
387
388 assert( align * (mp_d->get_bvect_cart()).get_align() == 1 ) ;
389
390 assert(cm_d.etat == ETATQCQ) ;
391
392 if (cm_d.dzpuis != 0) {
393 cout <<
394 "Cmp::import_align_asymy : the dzpuis of the Cmp to be imported"
395 << " must be zero !" << endl ;
396 abort() ;
397 }
398
399
400 const Mg3d* mg_a = mp->get_mg() ;
401 assert(mg_a->get_type_p() == NONSYM) ;
402
403 int nz_a = mg_a->get_nzone() ;
404 assert(nzet <= nz_a) ;
405
406 const Valeur& va_d = cm_d.va ;
407 va_d.coef() ; // The coefficients are required
408
409
410 // Preparations for storing the result in *this
411 // --------------------------------------------
412 del_t() ; // delete all previously computed derived quantities
413
414 set_etat_qcq() ; // Set the state to ETATQCQ
415
416 va.set_etat_c_qcq() ; // Allocates the memory for the Mtbl va.c
417 // if it does not exist already
418 va.c->set_etat_qcq() ; // Allocates the memory for the Tbl's in each
419 // domain if they do not exist already
420
421
422 // Departure (x,y,z) coordinates of the origin of the Arrival mapping :
423
424 double xx_a, yy_a, zz_a ;
425 if (align == 1) {
426 xx_a = mp->get_ori_x() - mp_d->get_ori_x() ;
427 yy_a = mp->get_ori_y() - mp_d->get_ori_y() ;
428 }
429 else {
430 xx_a = mp_d->get_ori_x() - mp->get_ori_x() ;
431 yy_a = mp_d->get_ori_y() - mp->get_ori_y() ;
432 }
433 zz_a = mp->get_ori_z() - mp_d->get_ori_z() ;
434
435
436 // r, theta, phi, x, y and z on the Arrival mapping
437 // update of the corresponding Coord's if necessary
438
439 if ( (mp->r).c == 0x0 ) (mp->r).fait() ;
440 if ( (mp->tet).c == 0x0 ) (mp->tet).fait() ;
441 if ( (mp->phi).c == 0x0 ) (mp->phi).fait() ;
442 if ( (mp->x).c == 0x0 ) (mp->x).fait() ;
443 if ( (mp->y).c == 0x0 ) (mp->y).fait() ;
444 if ( (mp->z).c == 0x0 ) (mp->z).fait() ;
445
446 const Mtbl* mr_a = (mp->r).c ;
447 const Mtbl* mtet_a = (mp->tet).c ;
448 const Mtbl* mphi_a = (mp->phi).c ;
449 const Mtbl* mx_a = (mp->x).c ;
450 const Mtbl* my_a = (mp->y).c ;
451 const Mtbl* mz_a = (mp->z).c ;
452
453 Param par_precis ; // Required precision in the method Map::val_lx
454 int nitermax = 100 ; // Maximum number of iteration in the secant method
455 int niter ;
456 double precis = 1e-15 ; // Absolute precision in the secant method
457 par_precis.add_int(nitermax) ;
458 par_precis.add_int_mod(niter) ;
459 par_precis.add_double(precis) ;
460
461
462 // Loop of the Arrival domains where the computation is to be performed
463 // --------------------------------------------------------------------
464
465 for (int l=0; l < nzet; l++) {
466
467 int nr = mg_a->get_nr(l) ;
468 int nt = mg_a->get_nt(l) ;
469 int np = mg_a->get_np(l) ;
470 int ntnr = nt*nr ;
471
472 const double* pr_a = mr_a->t[l]->t ; // Pointer on the values of r
473 const double* ptet_a = mtet_a->t[l]->t ; // Pointer on the values of theta
474 const double* pphi_a = mphi_a->t[l]->t ; // Pointer on the values of phi
475 const double* px_a = mx_a->t[l]->t ; // Pointer on the values of X
476 const double* py_a = my_a->t[l]->t ; // Pointer on the values of Y
477 const double* pz_a = mz_a->t[l]->t ; // Pointer on the values of Z
478
479 (va.c->t[l])->set_etat_qcq() ; // Allocates the array of double to
480 // store the result
481 double* ptx = (va.c->t[l])->t ; // Pointer on the allocated array
482
483
484
485 // Loop on half the grid points in the considered arrival domain
486 // (the other half will be obtained by antisymmetry with respect to
487 // the y=0 plane).
488
489 // Case k=0 (phi=0) : the function is zero (by antisymmetry)
490 for (int i=0; i<ntnr; i++) {
491 *ptx = 0 ;
492 ptx++ ; // next point
493 }
494
495 // Go to k=1 :
496 pr_a += ntnr ;
497 ptet_a += ntnr ;
498 pphi_a += ntnr ;
499 px_a += ntnr ;
500 py_a += ntnr ;
501 pz_a += ntnr ;
502
503 for (int k=1 ; k<np/2 ; k++) { // np/2 : ~ half the grid
504 for (int j=0 ; j<nt ; j++) {
505 for (int i=0 ; i<nr ; i++) {
506
507 double r = *pr_a ;
508 double rd, tetd, phid ;
509 if (r == __infinity) {
510 rd = r ;
511 tetd = *ptet_a ;
512 phid = *pphi_a ;
513 }
514 else {
515
516 // Cartesian coordinates on the Departure mapping
517 double xd = *px_a + xx_a ;
518 double yd = *py_a + yy_a ;
519 double zd = *pz_a + zz_a ;
520
521 // Spherical coordinates on the Departure mapping
522 double rhod2 = xd*xd + yd*yd ;
523 double rhod = sqrt( rhod2 ) ;
524 rd = sqrt(rhod2 + zd*zd) ;
525 tetd = atan2(rhod, zd) ;
526 phid = atan2(yd, xd) ;
527 if (phid < 0) phid += 2*M_PI ;
528 }
529
530
531 // NB: to increase the efficiency, the method Cmp::val_point
532 // is not invoked; the method Mtbl_cf::val_point is
533 // called directly instead.
534
535 // Value of the grid coordinates (l,xi) corresponding to
536 // (rd,tetd,phid) :
537
538 int ld ; // domain index
539 double xxd ; // radial coordinate xi in [0,1] or [-1,1]
540 mp_d->val_lx(rd, tetd, phid, par_precis, ld, xxd) ;
541
542 // Value of the Departure Cmp at the obtained point:
543 *ptx = va_d.c_cf->val_point_asymy(ld, xxd, tetd, phid) ;
544
545 // Next point :
546 ptx++ ;
547 pr_a++ ;
548 ptet_a++ ;
549 pphi_a++ ;
550 px_a++ ;
551 py_a++ ;
552 pz_a++ ;
553
554 }
555 }
556 }
557
558
559 // Case k=np/2 (phi=pi) : the function is zero (by antisymmetry)
560 for (int i=0; i<ntnr; i++) {
561 *ptx = 0 ;
562 ptx++ ; // next point
563 }
564
565 // Go to k=np/2+1 :
566 pr_a += ntnr ;
567 ptet_a += ntnr ;
568 pphi_a += ntnr ;
569 px_a += ntnr ;
570 py_a += ntnr ;
571 pz_a += ntnr ;
572
573 // The remaining points are obtained by antisymmetry with respect to the
574 // y=0 plane
575
576 for (int k=np/2+1 ; k<np ; k++) {
577
578 // pointer on the value (already computed) at the point symmetric
579 // with respect to the plane y=0
580 double* ptx_symy = (va.c->t[l])->t + (np-k)*nt*nr ;
581
582 // copy :
583 for (int j=0 ; j<nt ; j++) {
584 for (int i=0 ; i<nr ; i++) {
585 *ptx = - (*ptx_symy) ;
586 ptx++ ;
587 ptx_symy++ ;
588 }
589 }
590 }
591
592 } // End of the loop on the Arrival domains
593
594 // In the remaining domains, *this is set to zero:
595 // ----------------------------------------------
596
597 if (nzet < nz_a) {
598 annule(nzet, nz_a - 1) ;
599 }
600
601 // Treatment of dzpuis
602 // -------------------
603
604 set_dzpuis(0) ;
605
606}
607}
const Map * mp
Reference mapping.
Definition cmp.h:451
Cmp(const Map &map)
Constructor from mapping.
Definition cmp.C:211
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified z...
Definition cmp.h:461
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition cmp.h:454
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:351
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:307
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:292
void import_anti_asymy(int nzet, const Cmp &ci)
Assignment to another Cmp defined on a different mapping, when the two mappings have anti-aligned Car...
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:657
void import_align_asymy(int nzet, const Cmp &ci)
Assignment to another Cmp defined on a different mapping, when the two mappings have aligned Cartesia...
void del_t()
Logical destructor.
Definition cmp.C:262
void import_asymy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
Multi-domain grid.
Definition grilles.h:279
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:479
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:512
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:474
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:469
double val_point_asymy(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 ...
Multi-domain array.
Definition mtbl.h:118
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
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:318
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:388
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:249
double * t
The array of double.
Definition tbl.h:173
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
void coef() const
Computes the coeffcients of *this.
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene prototypes.
Definition app_hor.h:67
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Coord r
r coordinate centered on the grid
Definition map.h:730