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