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