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