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