LORENE
map_et_fait.C
1/*
2 * Copyright (c) 1999-2003 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: map_et_fait.C,v 1.10 2016/12/05 16:17:57 j_novak Exp $
27 * $Log: map_et_fait.C,v $
28 * Revision 1.10 2016/12/05 16:17:57 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.9 2014/10/13 08:53:04 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.8 2014/10/06 15:13:13 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.7 2013/06/05 15:10:42 j_novak
38 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
39 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
40 *
41 * Revision 1.6 2012/01/24 14:59:12 j_novak
42 * Removed functions XXX_fait_xi()
43 *
44 * Revision 1.5 2012/01/17 10:33:59 j_penner
45 * added a routine to construct the computational coordinate xi
46 *
47 * Revision 1.4 2008/08/27 08:48:42 jl_cornou
48 * Added R_JACO02 case
49 *
50 * Revision 1.3 2003/10/15 10:38:47 e_gourgoulhon
51 * Changed cast (const Map_et*) into static_cast<const Map_et*>.
52 *
53 * Revision 1.2 2002/10/16 14:36:41 j_novak
54 * Reorganization of #include instructions of standard C++, in order to
55 * use experimental version 3 of gcc.
56 *
57 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
58 * LORENE
59 *
60 * Revision 1.1 1999/11/24 11:23:00 eric
61 * Initial revision
62 *
63 *
64 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_fait.C,v 1.10 2016/12/05 16:17:57 j_novak Exp $
65 *
66 */
67
68
69// Includes
70#include <cassert>
71#include <cstdlib>
72#include <cmath>
73
74#include "map.h"
75
76 //----------------//
77 // Coord. radiale //
78 //----------------//
79
80namespace Lorene {
81Mtbl* map_et_fait_r(const Map* cvi) {
82
83 // recup du changement de variable
84 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
85 const Mg3d* mg = cv->get_mg() ;
86 int nz = mg->get_nzone() ;
87
88 // Le resultat
89 Mtbl* mti = new Mtbl(mg) ;
90 mti->set_etat_qcq() ;
91
92 // Pour le confort
93 const double* alpha = cv->alpha ;
94 const double* beta = cv->beta ;
95 const Valeur& ff = cv->ff ;
96 const Valeur& gg = cv->gg ;
97
98 for (int l=0 ; l<nz ; l++) {
99
100 const Grille3d* g = mg->get_grille3d(l) ;
101
102 const Tbl& aa = *((cv->aa)[l]) ;
103 const Tbl& bb = *((cv->bb)[l]) ;
104
105 Tbl* tb = (mti->t)[l] ;
106 tb->set_etat_qcq() ;
107 double* p_r = tb->t ;
108
109 int np = mg->get_np(l) ;
110 int nt = mg->get_nt(l) ;
111 int nr = mg->get_nr(l) ;
112
113 switch(mg->get_type_r(l)) {
114
115 case FIN: case RARE: {
116
117 for (int k=0 ; k<np ; k++) {
118 for (int j=0 ; j<nt ; j++) {
119 for (int i=0 ; i<nr ; i++) {
120 *p_r = alpha[l] * ( (g->x)[i]
121 + aa(i) * ff(l, k, j, 0)
122 + bb(i) * gg(l, k, j, 0)
123 ) + beta[l] ;
124 p_r++ ;
125 } // Fin de boucle sur r
126 } // Fin de boucle sur theta
127 } // Fin de boucle sur phi
128 break ;
129 }
130
131 case UNSURR: {
132 for (int k=0 ; k<np ; k++) {
133 for (int j=0 ; j<nt ; j++) {
134 for (int i=0 ; i<nr ; i++) {
135 *p_r = 1./( alpha[l] * (
136 (g->x)[i] + aa(i) * ff(l, k, j, 0)
137 )
138 + beta[l] ) ;
139 p_r++ ;
140 } // Fin de boucle sur r
141 } // Fin de boucle sur theta
142 } // Fin de boucle sur phi
143 break ;
144 }
145
146 default: {
147 cout << "map_et_fait_r: Unknown type_r !" << endl ;
148 abort () ;
149 }
150
151 } // Fin du switch
152 } // Fin de boucle sur zone
153
154 // Termine
155 return mti ;
156}
157
158 //--------------//
159 // Coord. Theta //
160 //--------------//
161
162Mtbl* map_et_fait_tet(const Map* cvi) {
163
164 // recup du changement de variable
165 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
166 const Mg3d* mg = cv->get_mg() ;
167 int nz = mg->get_nzone() ;
168
169 // Le resultat
170 Mtbl* mti = new Mtbl(mg) ;
171 mti->set_etat_qcq() ;
172
173 for (int l=0 ; l<nz ; l++) {
174 int nr = mg->get_nr(l);
175 int nt = mg->get_nt(l);
176 int np = mg->get_np(l);
177 const Grille3d* g = mg->get_grille3d(l) ;
178 Tbl* tb = (mti->t)[l] ;
179 tb->set_etat_qcq() ;
180 double* p_r = tb->t ;
181 for (int k=0 ; k<np ; k++) {
182 for (int j=0 ; j<nt ; j++) {
183 for (int i=0 ; i<nr ; i++) {
184 *p_r = (g->tet)[j] ;
185 p_r++ ;
186 } // Fin de boucle sur r
187 } // Fin de boucle sur theta
188 } // Fin de boucle sur phi
189 } // Fin de boucle sur zone
190
191 // Termine
192 return mti ;
193}
194
195 //------------//
196 // Coord. Phi //
197 //------------//
198
199Mtbl* map_et_fait_phi(const Map* cvi) {
200
201 // recup du changement de variable
202 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
203 const Mg3d* mg = cv->get_mg() ;
204 int nz = mg->get_nzone() ;
205
206 // Le resultat
207 Mtbl* mti = new Mtbl(mg) ;
208 mti->set_etat_qcq() ;
209
210 for (int l=0 ; l<nz ; l++) {
211 int nr = mg->get_nr(l);
212 int nt = mg->get_nt(l);
213 int np = mg->get_np(l);
214 const Grille3d* g = mg->get_grille3d(l) ;
215 Tbl* tb = (mti->t)[l] ;
216 tb->set_etat_qcq() ;
217 double* p_r = tb->t ;
218 for (int k=0 ; k<np ; k++) {
219 for (int j=0 ; j<nt ; j++) {
220 for (int i=0 ; i<nr ; i++) {
221 *p_r = (g->phi)[k] ;
222 p_r++ ;
223 } // Fin de boucle sur r
224 } // Fin de boucle sur theta
225 } // Fin de boucle sur phi
226 } // Fin de boucle sur zone
227
228 // Termine
229 return mti ;
230}
231
232 //----------//
233 // Coord. X //
234 //----------//
235
236Mtbl* map_et_fait_x(const Map* cvi) {
237
238 // recup de la grille
239 const Mg3d* mg = cvi->get_mg() ;
240
241 // Le resultat
242 Mtbl* mti = new Mtbl(mg) ;
243
244 *mti = (cvi->r) * (cvi->sint) * (cvi->cosp) ;
245
246 // Termine
247 return mti ;
248}
249
250 //----------//
251 // Coord. Y //
252 //----------//
253
254Mtbl* map_et_fait_y(const Map* cvi) {
255
256 // recup de la grille
257 const Mg3d* mg = cvi->get_mg() ;
258
259 // Le resultat
260 Mtbl* mti = new Mtbl(mg) ;
261
262 *mti = (cvi->r) * (cvi->sint) * (cvi->sinp) ;
263
264 // Termine
265 return mti ;
266}
267
268 //----------//
269 // Coord. Z //
270 //----------//
271
272Mtbl* map_et_fait_z(const Map* cvi) {
273
274 // recup de la grille
275 const Mg3d* mg = cvi->get_mg() ;
276
277 // Le resultat
278 Mtbl* mti = new Mtbl(mg) ;
279
280 *mti = (cvi->r) * (cvi->cost) ;
281
282 // Termine
283 return mti ;
284}
285
286 //--------------------//
287 // Coord. X "absolue" //
288 //--------------------//
289
290Mtbl* map_et_fait_xa(const Map* cvi) {
291
292 // recup de la grille
293 const Mg3d* mg = cvi->get_mg() ;
294
295 // Le resultat
296 Mtbl* mti = new Mtbl(mg) ;
297
298 double r_phi = cvi->get_rot_phi() ;
299 double t_x = cvi->get_ori_x() ;
300
301 if ( fabs(r_phi) < 1.e-14 ) {
302 *mti = (cvi->x) + t_x ;
303 }
304 else if ( fabs(r_phi - M_PI) < 1.e-14 ) {
305 *mti = - (cvi->x) + t_x ;
306 }
307 else {
308 Mtbl phi = cvi->phi + r_phi ;
309 *mti = (cvi->r) * (cvi->sint) * cos(phi) + t_x ;
310 }
311
312 // Termine
313 return mti ;
314}
315
316 //--------------------//
317 // Coord. Y "absolue" //
318 //--------------------//
319
320Mtbl* map_et_fait_ya(const Map* cvi) {
321
322 // recup de la grille
323 const Mg3d* mg = cvi->get_mg() ;
324
325 // Le resultat
326 Mtbl* mti = new Mtbl(mg) ;
327
328 double r_phi = cvi->get_rot_phi() ;
329 double t_y = cvi->get_ori_y() ;
330
331 if ( fabs(r_phi) < 1.e-14 ) {
332 *mti = (cvi->y) + t_y ;
333 }
334 else if ( fabs(r_phi - M_PI) < 1.e-14 ) {
335 *mti = - (cvi->y) + t_y ;
336 }
337 else {
338 Mtbl phi = cvi->phi + r_phi ;
339 *mti = (cvi->r) * (cvi->sint) * sin(phi) + t_y ;
340 }
341
342 // Termine
343 return mti ;
344}
345
346 //--------------------//
347 // Coord. Z "absolue" //
348 //--------------------//
349
350Mtbl* map_et_fait_za(const Map* cvi) {
351
352 // recup de la grille
353 const Mg3d* mg = cvi->get_mg() ;
354
355 double t_z = cvi->get_ori_z() ;
356
357 // Le resultat
358 Mtbl* mti = new Mtbl(mg) ;
359
360 *mti = (cvi->r) * (cvi->cost) + t_z ;
361
362 // Termine
363 return mti ;
364}
365
366 //---------------//
367 // Trigonometrie //
368 //---------------//
369
370Mtbl* map_et_fait_sint(const Map* cvi) {
371
372 // recup du changement de variable
373 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
374 const Mg3d* mg = cv->get_mg() ;
375 int nz = mg->get_nzone() ;
376
377 // Le resultat
378 Mtbl* mti = new Mtbl(mg) ;
379 mti->set_etat_qcq() ;
380
381 for (int l=0 ; l<nz ; l++) {
382 int nr = mg->get_nr(l);
383 int nt = mg->get_nt(l);
384 int np = mg->get_np(l);
385 const Grille3d* g = mg->get_grille3d(l) ;
386 Tbl* tb = (mti->t)[l] ;
387 tb->set_etat_qcq() ;
388 double* p_r = tb->t ;
389 for (int k=0 ; k<np ; k++) {
390 for (int j=0 ; j<nt ; j++) {
391 for (int i=0 ; i<nr ; i++) {
392 *p_r = sin(g->tet[j]) ;
393 p_r++ ;
394 } // Fin de boucle sur r
395 } // Fin de boucle sur theta
396 } // Fin de boucle sur phi
397 } // Fin de boucle sur zone
398
399 // Termine
400 return mti ;
401}
402
403Mtbl* map_et_fait_cost(const Map* cvi) {
404
405 // recup du changement de variable
406 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
407 const Mg3d* mg = cv->get_mg() ;
408 int nz = mg->get_nzone() ;
409
410 // Le resultat
411 Mtbl* mti = new Mtbl(mg) ;
412 mti->set_etat_qcq() ;
413
414 for (int l=0 ; l<nz ; l++) {
415 int nr = mg->get_nr(l);
416 int nt = mg->get_nt(l);
417 int np = mg->get_np(l);
418 const Grille3d* g = mg->get_grille3d(l) ;
419 Tbl* tb = (mti->t)[l] ;
420 tb->set_etat_qcq() ;
421 double* p_r = tb->t ;
422 for (int k=0 ; k<np ; k++) {
423 for (int j=0 ; j<nt ; j++) {
424 for (int i=0 ; i<nr ; i++) {
425 *p_r = cos(g->tet[j]) ;
426 p_r++ ;
427 } // Fin de boucle sur r
428 } // Fin de boucle sur theta
429 } // Fin de boucle sur phi
430 } // Fin de boucle sur zone
431
432 // Termine
433 return mti ;
434}
435
436Mtbl* map_et_fait_sinp(const Map* cvi) {
437
438 // recup du changement de variable
439 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
440 const Mg3d* mg = cv->get_mg() ;
441 int nz = mg->get_nzone() ;
442
443 // Le resultat
444 Mtbl* mti = new Mtbl(mg) ;
445 mti->set_etat_qcq() ;
446
447 for (int l=0 ; l<nz ; l++) {
448 int nr = mg->get_nr(l);
449 int nt = mg->get_nt(l);
450 int np = mg->get_np(l);
451 const Grille3d* g = mg->get_grille3d(l) ;
452 Tbl* tb = (mti->t)[l] ;
453 tb->set_etat_qcq() ;
454 double* p_r = tb->t ;
455 for (int k=0 ; k<np ; k++) {
456 for (int j=0 ; j<nt ; j++) {
457 for (int i=0 ; i<nr ; i++) {
458 *p_r = sin(g->phi[k]) ;
459 p_r++ ;
460 } // Fin de boucle sur r
461 } // Fin de boucle sur theta
462 } // Fin de boucle sur phi
463 } // Fin de boucle sur zone
464
465 // Termine
466 return mti ;
467}
468
469Mtbl* map_et_fait_cosp(const Map* cvi) {
470
471 // recup du changement de variable
472 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
473 const Mg3d* mg = cv->get_mg() ;
474 int nz = mg->get_nzone() ;
475
476 // Le resultat
477 Mtbl* mti = new Mtbl(mg) ;
478 mti->set_etat_qcq() ;
479
480 for (int l=0 ; l<nz ; l++) {
481 int nr = mg->get_nr(l);
482 int nt = mg->get_nt(l);
483 int np = mg->get_np(l);
484 const Grille3d* g = mg->get_grille3d(l) ;
485 Tbl* tb = (mti->t)[l] ;
486 tb->set_etat_qcq() ;
487 double* p_r = tb->t ;
488 for (int k=0 ; k<np ; k++) {
489 for (int j=0 ; j<nt ; j++) {
490 for (int i=0 ; i<nr ; i++) {
491 *p_r = cos(g->phi[k]) ;
492 p_r++ ;
493 } // Fin de boucle sur r
494 } // Fin de boucle sur theta
495 } // Fin de boucle sur phi
496 } // Fin de boucle sur zone
497
498 // Termine
499 return mti ;
500}
501
502/*
503 ************************************************************************
504 * x/R dans le noyau, 1/R dans les coquilles, (x-1)/U dans la ZEC
505 ************************************************************************
506 */
507
508Mtbl* map_et_fait_xsr(const Map* cvi) {
509
510 // recup du changement de variable
511 const Map_et* cv = static_cast<const Map_et*>(cvi) ;
512 const Mg3d* mg = cv->get_mg() ;
513 int nz = mg->get_nzone() ;
514
515 // Le resultat
516 Mtbl* mti = new Mtbl(mg) ;
517 mti->set_etat_qcq() ;
518
519 // Pour le confort
520 const double* alpha = cv->alpha ;
521 const double* beta = cv->beta ;
522 const Valeur& ff = cv->ff ;
523 const Valeur& gg = cv->gg ;
524 const Tbl& asx = cv->aasx ;
525 const Tbl& bsx = cv->bbsx ;
526 const Tbl& asxm1 = cv->zaasx ;
527
528 for (int l=0 ; l<nz ; l++) {
529 int nr = mg->get_nr(l);
530 int nt = mg->get_nt(l) ;
531 int np = mg->get_np(l) ;
532 const Grille3d* g = mg->get_grille3d(l) ;
533
534 const Tbl& aa = *((cv->aa)[l]) ;
535 const Tbl& bb = *((cv->bb)[l]) ;
536
537 Tbl* tb = (mti->t)[l] ;
538 tb->set_etat_qcq() ;
539 double* p_r = tb->t ;
540
541 switch(mg->get_type_r(l)) {
542
543 case RARE: {
544 assert(beta[l]==0) ;
545 for (int k=0 ; k<np ; k++) {
546 for (int j=0 ; j<nt ; j++) {
547 for (int i=0 ; i<nr ; i++) {
548 *p_r = 1. / ( alpha[l] * ( 1. + asx(i) * ff(l, k, j, 0)
549 + bsx(i) * gg(l, k, j, 0)
550 ) ) ;
551 p_r++ ;
552 } // Fin de boucle sur r
553 } // Fin de boucle sur theta
554 } // Fin de boucle sur phi
555 break ;
556 }
557
558 case FIN: {
559 for (int k=0 ; k<np ; k++) {
560 for (int j=0 ; j<nt ; j++) {
561 for (int i=0 ; i<nr ; i++) {
562 *p_r = 1. / ( alpha[l] * ( (g->x)[i]
563 + aa(i) * ff(l, k, j, 0)
564 + bb(i) * gg(l, k, j, 0)
565 ) + beta[l] );
566 p_r++ ;
567 } // Fin de boucle sur r
568 } // Fin de boucle sur theta
569 } // Fin de boucle sur phi
570 break ;
571 }
572
573 case UNSURR: {
574 assert(beta[l] == - alpha[l]) ;
575 for (int k=0 ; k<np ; k++) {
576 for (int j=0 ; j<nt ; j++) {
577 for (int i=0 ; i<nr ; i++) {
578 *p_r = 1. / ( alpha[l] * ( 1.
579 + asxm1(i) * ff(l, k, j, 0)
580 ) ) ;
581 p_r++ ;
582 } // Fin de boucle sur r
583 } // Fin de boucle sur theta
584 } // Fin de boucle sur phi
585 break ;
586 }
587
588 default: {
589 cout << "map_et_fait_xsr: unknown type_r !" << endl ;
590 abort() ;
591 }
592
593 } // Fin du switch
594 } // Fin de boucle sur zone
595
596 // Termine
597 return mti ;
598
599}
600
601}
Lorene::Grille3d
3D grid class in one domain.
Definition grilles.h:200
Lorene::Map_et
Radial mapping of rather general form.
Definition map.h:2770
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
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::phi
Coord phi
coordinate centered on the grid
Definition map.h:732