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