LORENE
map_eps_fait.C
1/*
2 * This file is part of LORENE.
3 *
4 * LORENE is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * LORENE is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with LORENE; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 *
18 */
19
20
21// Includes
22#include <cassert>
23#include <cstdlib>
24#include <cmath>
25
26#include "mtbl.h"
27#include "map.h"
28#include "proto.h"
29#include "valeur.h"
30
31
32
33
34
35
36
37
38 //----------------//
39 // Coord. radiale //
40 //----------------//
41
42namespace Lorene {
43Mtbl* map_eps_fait_r(const Map* cvi) {
44
45 // recup du changement de variable
46 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
47 const Mg3d* mg = cv->get_mg() ;
48 int nz = mg->get_nzone() ;
49
50 // Le resultat
51 Mtbl* mti = new Mtbl(mg) ;
52 mti->set_etat_qcq() ;
53
54 // Pour le confort
55 const Valeur& alpha = cv->get_alpha() ;
56 const Valeur& beta = cv->get_beta() ;
57
58 int i, j, k ;
59 for (int l=0 ; l<nz ; l++) {
60 int ir = mg->get_nr(l);
61 int it = mg->get_nt(l) ;
62 int ip = mg->get_np(l) ;
63 const Grille3d* g = mg->get_grille3d(l) ;
64 Tbl* tb = (mti->t)[l] ;
65 tb->set_etat_qcq() ;
66 double* p_r = tb->t ;
67
68 switch(mg->get_type_r(l)) {
69 case RARE: case FIN:
70 for (k=0 ; k<ip ; k++) {
71 for (j=0 ; j<it ; j++) {
72 for (i=0 ; i<ir ; i++) {
73 *p_r = alpha(l,k,j,0) * (g->x)[i] + beta(l,k,j,0) ;
74 p_r++ ;
75 } // Fin de boucle sur r
76 } // Fin de boucle sur theta
77 } // Fin de boucle sur phi
78 break ;
79
80 // case FIN:{
81 // cout << "map_eps_fait_r: Shells not implemented yet..." << endl;
82 // abort() ;
83 // break ;
84 // }
85
86 case UNSURR:
87 for (k=0 ; k<ip ; k++) {
88 for (j=0 ; j<it ; j++) {
89 for (i=0 ; i<ir ; i++) {
90 cout << "map_eps_fait_r: Warning ! No compactified zone allowed !" << endl;
91 abort() ;
92 } // Fin de boucle sur r
93 } // Fin de boucle sur theta
94 } // Fin de boucle sur phi
95 break ;
96
97 default:
98 cout << "map_eps_fait_r: unknown type_r !\n" ;
99 abort () ;
100 exit(-1) ;
101
102 } // Fin du switch
103 } // Fin de boucle sur zone
104
105 // Termine
106 return mti ;
107}
108
109 //--------------//
110 // Coord. Theta //
111 //--------------//
112
113Mtbl* map_eps_fait_tet(const Map* cvi) {
114
115 // recup du changement de variable
116 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
117 const Mg3d* mg = cv->get_mg() ;
118 int nz = mg->get_nzone() ;
119
120 // Le resultat
121 Mtbl* mti = new Mtbl(mg) ;
122 mti->set_etat_qcq() ;
123
124 int i, j, k ;
125 for (int l=0 ; l<nz ; l++) {
126 int ir = mg->get_nr(l);
127 int it = mg->get_nt(l);
128 int ip = mg->get_np(l);
129 const Grille3d* g = mg->get_grille3d(l) ;
130 Tbl* tb = (mti->t)[l] ;
131 tb->set_etat_qcq() ;
132 double* p_r = tb->t ;
133 for (k=0 ; k<ip ; k++) {
134 for (j=0 ; j<it ; j++) {
135 for (i=0 ; i<ir ; i++) {
136 *p_r = (g->tet)[j] ;
137 p_r++ ;
138 } // Fin de boucle sur r
139 } // Fin de boucle sur theta
140 } // Fin de boucle sur phi
141 } // Fin de boucle sur zone
142
143 // Termine
144 return mti ;
145}
146
147 //------------//
148 // Coord. Phi //
149 //------------//
150
151Mtbl* map_eps_fait_phi(const Map* cvi) {
152
153 // recup du changement de variable
154 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
155 const Mg3d* mg = cv->get_mg() ;
156 int nz = mg->get_nzone() ;
157
158 // Le resultat
159 Mtbl* mti = new Mtbl(mg) ;
160 mti->set_etat_qcq() ;
161
162 int i, j, k ;
163 for (int l=0 ; l<nz ; l++) {
164 int ir = mg->get_nr(l);
165 int it = mg->get_nt(l);
166 int ip = mg->get_np(l);
167 const Grille3d* g = mg->get_grille3d(l) ;
168 Tbl* tb = (mti->t)[l] ;
169 tb->set_etat_qcq() ;
170 double* p_r = tb->t ;
171 for (k=0 ; k<ip ; k++) {
172 for (j=0 ; j<it ; j++) {
173 for (i=0 ; i<ir ; i++) {
174 *p_r = (g->phi)[k] ;
175 p_r++ ;
176 } // Fin de boucle sur r
177 } // Fin de boucle sur theta
178 } // Fin de boucle sur phi
179 } // Fin de boucle sur zone
180
181 // Termine
182 return mti ;
183}
184
185 //----------//
186 // Coord. X //
187 //----------//
188
189Mtbl* map_eps_fait_x(const Map* cvi) {
190
191 // recup de la grille
192 const Mg3d* mg = cvi->get_mg() ;
193
194 // Le resultat
195 Mtbl* mti = new Mtbl(mg) ;
196
197 *mti = (cvi->r) * (cvi->sint) * (cvi->cosp) ;
198
199 // Termine
200 return mti ;
201}
202
203 //----------//
204 // Coord. Y //
205 //----------//
206
207Mtbl* map_eps_fait_y(const Map* cvi) {
208
209 // recup de la grille
210 const Mg3d* mg = cvi->get_mg() ;
211
212 // Le resultat
213 Mtbl* mti = new Mtbl(mg) ;
214
215 *mti = (cvi->r) * (cvi->sint) * (cvi->sinp) ;
216
217 // Termine
218 return mti ;
219}
220
221 //----------//
222 // Coord. Z //
223 //----------//
224
225Mtbl* map_eps_fait_z(const Map* cvi) {
226
227 // recup de la grille
228 const Mg3d* mg = cvi->get_mg() ;
229
230 // Le resultat
231 Mtbl* mti = new Mtbl(mg) ;
232
233 *mti = (cvi->r) * (cvi->cost) ;
234
235 // Termine
236 return mti ;
237}
238
239 //--------------------//
240 // Coord. X "absolue" //
241 //--------------------//
242
243Mtbl* map_eps_fait_xa(const Map* cvi) {
244
245 // recup de la grille
246 const Mg3d* mg = cvi->get_mg() ;
247
248 // Le resultat
249 Mtbl* mti = new Mtbl(mg) ;
250
251 double r_phi = cvi->get_rot_phi() ;
252 double t_x = cvi->get_ori_x() ;
253
254 if ( fabs(r_phi) < 1.e-14 ) {
255 *mti = (cvi->x) + t_x ;
256 }
257 else if ( fabs(r_phi - M_PI) < 1.e-14 ) {
258 *mti = - (cvi->x) + t_x ;
259 }
260 else {
261 Mtbl phi = cvi->phi + r_phi ;
262 *mti = (cvi->r) * (cvi->sint) * cos(phi) + t_x ;
263 }
264
265 // Termine
266 return mti ;
267}
268
269 //--------------------//
270 // Coord. Y "absolue" //
271 //--------------------//
272
273Mtbl* map_eps_fait_ya(const Map* cvi) {
274
275 // recup de la grille
276 const Mg3d* mg = cvi->get_mg() ;
277
278 // Le resultat
279 Mtbl* mti = new Mtbl(mg) ;
280
281 double r_phi = cvi->get_rot_phi() ;
282 double t_y = cvi->get_ori_y() ;
283
284 if ( fabs(r_phi) < 1.e-14 ) {
285 *mti = (cvi->y) + t_y ;
286 }
287 else if ( fabs(r_phi - M_PI) < 1.e-14 ) {
288 *mti = - (cvi->y) + t_y ;
289 }
290 else {
291 Mtbl phi = cvi->phi + r_phi ;
292 *mti = (cvi->r) * (cvi->sint) * sin(phi) + t_y ;
293 }
294
295 // Termine
296 return mti ;
297}
298
299 //--------------------//
300 // Coord. Z "absolue" //
301 //--------------------//
302
303Mtbl* map_eps_fait_za(const Map* cvi) {
304
305 // recup de la grille
306 const Mg3d* mg = cvi->get_mg() ;
307
308 double t_z = cvi->get_ori_z() ;
309
310 // Le resultat
311 Mtbl* mti = new Mtbl(mg) ;
312
313 *mti = (cvi->r) * (cvi->cost) + t_z ;
314
315 // Termine
316 return mti ;
317}
318
319 //---------------//
320 // Trigonometrie //
321 //---------------//
322
323Mtbl* map_eps_fait_sint(const Map* cvi) {
324
325 // recup du changement de variable
326 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
327 const Mg3d* mg = cv->get_mg() ;
328 int nz = mg->get_nzone() ;
329
330 // Le resultat
331 Mtbl* mti = new Mtbl(mg) ;
332 mti->set_etat_qcq() ;
333
334 int i, j, k ;
335 for (int l=0 ; l<nz ; l++) {
336 int ir = mg->get_nr(l);
337 int it = mg->get_nt(l);
338 int ip = mg->get_np(l);
339 const Grille3d* g = mg->get_grille3d(l) ;
340 Tbl* tb = (mti->t)[l] ;
341 tb->set_etat_qcq() ;
342 double* p_r = tb->t ;
343 for (k=0 ; k<ip ; k++) {
344 for (j=0 ; j<it ; j++) {
345 for (i=0 ; i<ir ; i++) {
346 *p_r = sin(g->tet[j]) ;
347 p_r++ ;
348 } // Fin de boucle sur r
349 } // Fin de boucle sur theta
350 } // Fin de boucle sur phi
351 } // Fin de boucle sur zone
352
353 // Termine
354 return mti ;
355}
356
357Mtbl* map_eps_fait_cost(const Map* cvi) {
358
359 // recup du changement de variable
360 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
361 const Mg3d* mg = cv->get_mg() ;
362 int nz = mg->get_nzone() ;
363
364 // Le resultat
365 Mtbl* mti = new Mtbl(mg) ;
366 mti->set_etat_qcq() ;
367
368 int i, j, k ;
369 for (int l=0 ; l<nz ; l++) {
370 int ir = mg->get_nr(l);
371 int it = mg->get_nt(l);
372 int ip = mg->get_np(l);
373 const Grille3d* g = mg->get_grille3d(l) ;
374 Tbl* tb = (mti->t)[l] ;
375 tb->set_etat_qcq() ;
376 double* p_r = tb->t ;
377 for (k=0 ; k<ip ; k++) {
378 for (j=0 ; j<it ; j++) {
379 for (i=0 ; i<ir ; i++) {
380 *p_r = cos(g->tet[j]) ;
381 p_r++ ;
382 } // Fin de boucle sur r
383 } // Fin de boucle sur theta
384 } // Fin de boucle sur phi
385 } // Fin de boucle sur zone
386
387 // Termine
388 return mti ;
389}
390
391Mtbl* map_eps_fait_sinp(const Map* cvi) {
392
393 // recup du changement de variable
394 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
395 const Mg3d* mg = cv->get_mg() ;
396 int nz = mg->get_nzone() ;
397
398 // Le resultat
399 Mtbl* mti = new Mtbl(mg) ;
400 mti->set_etat_qcq() ;
401
402 int i, j, k ;
403 for (int l=0 ; l<nz ; l++) {
404 int ir = mg->get_nr(l);
405 int it = mg->get_nt(l);
406 int ip = mg->get_np(l);
407 const Grille3d* g = mg->get_grille3d(l) ;
408 Tbl* tb = (mti->t)[l] ;
409 tb->set_etat_qcq() ;
410 double* p_r = tb->t ;
411 for (k=0 ; k<ip ; k++) {
412 for (j=0 ; j<it ; j++) {
413 for (i=0 ; i<ir ; i++) {
414 *p_r = sin(g->phi[k]) ;
415 p_r++ ;
416 } // Fin de boucle sur r
417 } // Fin de boucle sur theta
418 } // Fin de boucle sur phi
419 } // Fin de boucle sur zone
420
421 // Termine
422 return mti ;
423}
424
425Mtbl* map_eps_fait_cosp(const Map* cvi) {
426
427 // recup du changement de variable
428 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
429 const Mg3d* mg = cv->get_mg() ;
430 int nz = mg->get_nzone() ;
431
432 // Le resultat
433 Mtbl* mti = new Mtbl(mg) ;
434 mti->set_etat_qcq() ;
435
436 int i, j, k ;
437 for (int l=0 ; l<nz ; l++) {
438 int ir = mg->get_nr(l);
439 int it = mg->get_nt(l);
440 int ip = mg->get_np(l);
441 const Grille3d* g = mg->get_grille3d(l) ;
442 Tbl* tb = (mti->t)[l] ;
443 tb->set_etat_qcq() ;
444 double* p_r = tb->t ;
445 for (k=0 ; k<ip ; k++) {
446 for (j=0 ; j<it ; j++) {
447 for (i=0 ; i<ir ; i++) {
448 *p_r = cos(g->phi[k]) ;
449 p_r++ ;
450 } // Fin de boucle sur r
451 } // Fin de boucle sur theta
452 } // Fin de boucle sur phi
453 } // Fin de boucle sur zone
454
455 // Termine
456 return mti ;
457}
458
459/*
460 ************************************************************************
461 * x/R dans le noyau, 1/R dans les coquilles, (x-1)/U dans la ZEC
462 ************************************************************************
463 */
464
465Mtbl* map_eps_fait_xsr(const Map* cvi) {
466
467 // recup du changement de variable
468 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
469 const Mg3d* mg = cv->get_mg() ;
470 int nz = mg->get_nzone() ;
471
472 // Le resultat
473 Mtbl* mti = new Mtbl(mg) ;
474 mti->set_etat_qcq() ;
475
476 // Pour le confort
477 const Valeur& alpha = cv->get_alpha() ;
478 const Valeur& beta = cv->get_beta() ;
479
480 int i, j, k ;
481 for (int l=0 ; l<nz ; l++) {
482 int ir = mg->get_nr(l);
483 int it = mg->get_nt(l) ;
484 int ip = mg->get_np(l) ;
485 const Grille3d* g = mg->get_grille3d(l) ;
486 Tbl* tb = (mti->t)[l] ;
487 tb->set_etat_qcq() ;
488 double* p_r = tb->t ;
489
490 switch(mg->get_type_r(l)) {
491
492 case RARE:
493 assert(beta(l).get_etat()==ETATZERO) ;
494 for (k=0 ; k<ip ; k++) {
495 for (j=0 ; j<it ; j++) {
496 for (i=0 ; i<ir ; i++) {
497 *p_r = 1. / alpha(l,k,j,0) ;
498 p_r++ ;
499 } // Fin de boucle sur r
500 } // Fin de boucle sur theta
501 } // Fin de boucle sur phi
502 break ;
503
504 case FIN:
505 for (k=0 ; k<ip ; k++) {
506 for (j=0 ; j<it ; j++) {
507 if (ir == 1) { //Some hack for angular grid case...
508 *p_r = 1. / beta(l,k,j,0) ;
509 p_r++ ;
510 }
511 else
512 for (i=0 ; i<ir ; i++) {
513 *p_r = 1. / ( alpha(l,k,j,0) * (g->x)[i] + beta(l,k,j,0) ) ;
514 p_r++ ;
515 } // Fin de boucle sur r
516 } // Fin de boucle sur theta
517 } // Fin de boucle sur phi
518 break ;
519
520 case UNSURR:
521 cout << "map_eps_fait_xsr: Compactified domain not allowed !" << endl;
522 abort() ;
523 break ;
524
525 default:
526 cout << "map_eps_fait_xsr: unknown type_r !" << endl ;
527 abort() ;
528
529 } // Fin du switch
530 } // Fin de boucle sur zone
531
532 // Termine
533 return mti ;
534
535}
536
537/*
538 ************************************************************************
539 * 1/(dR/dx) ( -1/(dU/dx) ds la ZEC )
540 ************************************************************************
541 */
542
543Mtbl* map_eps_fait_dxdr(const Map* cvi) {
544
545 // recup du changement de variable
546 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
547 const Mg3d* mg = cv->get_mg() ;
548 int nz = mg->get_nzone() ;
549
550 // Le resultat
551 Mtbl* mti = new Mtbl(mg) ;
552 mti->set_etat_qcq() ;
553
554 // Pour le confort
555 const Valeur& alpha = cv->get_alpha() ;
556
557 int i, j, k ;
558 for (int l=0 ; l<nz ; l++) {
559 int ir = mg->get_nr(l);
560 int it = mg->get_nt(l) ;
561 int ip = mg->get_np(l) ;
562 Tbl* tb = (mti->t)[l] ;
563 tb->set_etat_qcq() ;
564 double* p_r = tb->t ;
565
566 switch(mg->get_type_r(l)) {
567
568 case RARE: case FIN:
569 for (k=0 ; k<ip ; k++) {
570 for (j=0 ; j<it ; j++) {
571 for (i=0 ; i<ir ; i++) {
572 *p_r = 1. / alpha(l,k,j,0) ;
573 p_r++ ;
574 } // Fin de boucle sur r
575 } // Fin de boucle sur theta
576 } // Fin de boucle sur phi
577 break ;
578
579 // case FIN:{
580 // cout << "map_eps_fait_dxdr: Shells not implemented yet..." << endl;
581 // abort() ;
582 // break ;
583 // }
584
585 case UNSURR:
586 cout << "map_eps_fait_dxdr: Compactified domain not allowed !" << endl;
587 abort() ;
588 break ;
589
590
591 default:
592 cout << "map_eps_fait_dxdr: unknown type_r !" << endl ;
593 abort() ;
594 } // Fin du switch
595 } // Fin de boucle sur zone
596
597 // Termine
598 return mti ;
599}
600
601/*
602 ************************************************************************
603 * dR/dtheta
604 ************************************************************************
605 */
606
607Mtbl* map_eps_fait_drdt(const Map* cvi) {
608
609 // recup du changement de variable
610 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
611 const Mg3d* mg = cv->get_mg() ;
612 int nz = mg->get_nzone() ;
613
614 // Le resultat
615 Mtbl* mti = new Mtbl(mg) ;
616 mti->set_etat_qcq() ;
617
618 // Pour le confort
619 const Valeur& alpha = cv->get_alpha() ;
620 Valeur dalphadt(alpha.get_mg()) ; dalphadt.annule_hard() ;
621 for (int k=0; k<mg->get_np(0); k++)
622 for (int j=0; j<mg->get_nt(0); j++){
623 const Grille3d* g = mg->get_grille3d(0) ;
624 double theta = (g->tet)[j] ;
625 double pphi = (g->phi)[k] ;
626 double sint = sin(theta) ;
627 double cost = cos(theta) ;
628 double sinp = sin(pphi) ;
629 double cosp = cos(pphi) ;
630 double aa = cv->get_aa() ;
631 double bb = cv->get_bb() ;
632 double cc = cv->get_cc() ;
633 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
634 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
635 dalphadt.set(0, k, j, 0) = alpha3 * sint * cost * (1./(cc*cc) - sinp*sinp / (bb*bb) - cosp*cosp / (aa*aa)) ;
636 }
637 const Valeur& beta = cv->get_beta() ;
638 const Valeur& dbetadt = beta.dsdt() ;
639
640
641 for (int l=0 ; l<nz ; l++) {
642
643 int nr = mg->get_nr(l);
644 int nt = mg->get_nt(l) ;
645 int np = mg->get_np(l) ;
646 const Grille3d* g = mg->get_grille3d(l) ;
647 Tbl* tb = (mti->t)[l] ;
648 tb->set_etat_qcq() ;
649 double* p_r = tb->t ;
650
651 switch(mg->get_type_r(l)) {
652
653
654 case RARE: case FIN: {
655 for (int k=0 ; k<np ; k++) {
656 for (int j=0 ; j<nt ; j++) {
657 for (int i=0 ; i<nr ; i++) {
658 *p_r = dalphadt(l,k,j,0) * (g->x)[i] + dbetadt(l,k,j,0) ;
659 p_r++ ;
660 } // Fin de boucle sur r
661 } // Fin de boucle sur theta
662 } // Fin de boucle sur phi
663 break ;
664 }
665 // case FIN:{
666 // cout << "map_eps_fait_drdt: Shells not implemented yet..." << endl;
667 // abort() ;
668 // break ;
669 // }
670
671 case UNSURR: {
672 cout << "map_eps_fait_drdt: Compactified domain not allowed !" << endl;
673 abort() ;
674 break ;
675 }
676
677 default: {
678 cout << "map_eps_fait_drdt: unknown type_r !" << endl ;
679 abort() ;
680 }
681 } // Fin du switch
682 } // Fin de boucle sur zone
683
684 // Termine
685
686 return mti ;
687}
688
689/*
690 ************************************************************************
691 * 1/sin(theta) dR/dphi
692 ************************************************************************
693 */
694
695Mtbl* map_eps_fait_stdrdp(const Map* cvi) {
696
697 // recup du changement de variable
698 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
699 const Mg3d* mg = cv->get_mg() ;
700 int nz = mg->get_nzone() ;
701
702 // Le resultat
703 Mtbl* mti = new Mtbl(mg) ;
704 mti->set_etat_qcq() ;
705
706 // Pour le confort
707 const Valeur& alpha = cv->get_alpha() ;
708 Valeur stalphadp(alpha.get_mg()) ; stalphadp.annule_hard() ;
709 for (int k=0; k<mg->get_np(0); k++)
710 for (int j=0; j<mg->get_nt(0); j++){
711 const Grille3d* g = mg->get_grille3d(0) ;
712 double theta = (g->tet)[j] ;
713 double pphi = (g->phi)[k] ;
714 double sint = sin(theta) ;
715 double cost = cos(theta) ;
716 double sinp = sin(pphi) ;
717 double cosp = cos(pphi) ;
718 double aa = cv->get_aa() ;
719 double bb = cv->get_bb() ;
720 double cc = cv->get_cc() ;
721 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
722 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
723 stalphadp.set(0, k, j, 0) = alpha3 * sint * cosp * sinp * (1./(aa*aa) - 1./(bb*bb)) ; // Attention à la division par sin(theta) ici !!
724
725 }
726 const Valeur& beta = cv->get_beta() ;
727 const Valeur& stbetadp = beta.stdsdp() ;
728
729
730 for (int l=0 ; l<nz ; l++) {
731
732 int nr = mg->get_nr(l);
733 int nt = mg->get_nt(l) ;
734 int np = mg->get_np(l) ;
735 const Grille3d* g = mg->get_grille3d(l) ;
736 Tbl* tb = (mti->t)[l] ;
737 tb->set_etat_qcq() ;
738 double* p_r = tb->t ;
739
740 switch(mg->get_type_r(l)) {
741
742
743 case RARE: case FIN: {
744 for (int k=0 ; k<np ; k++) {
745 for (int j=0 ; j<nt ; j++) {
746 for (int i=0 ; i<nr ; i++) {
747 *p_r = stalphadp(l,k,j,0) * (g->x)[i] + stbetadp(l,k,j,0) ;
748 p_r++ ;
749 } // Fin de boucle sur r
750 } // Fin de boucle sur theta
751 } // Fin de boucle sur phi
752 break ;
753 }
754
755 // case FIN:{
756 // cout << "map_eps_fait_stdrdp: Shells not implemented yet..." << endl;
757 // abort() ;
758 // break ;
759 // }
760
761 case UNSURR: {
762 cout << "map_eps_fait_stdrdp: Compactified domain not allowed !" << endl;
763 abort() ;
764 break ;
765 }
766
767 default: {
768 cout << "map_eps_fait_stdrdp: unknown type_r !" << endl ;
769 abort() ;
770 }
771 } // Fin du switch
772 } // Fin de boucle sur zone
773
774 // Termine
775
776 return mti ;
777}
778
779/*
780 ************************************************************************
781 * 1/R dR/dtheta
782 ************************************************************************
783 */
784
785Mtbl* map_eps_fait_srdrdt(const Map* cvi) {
786
787 // recup du changement de variable
788 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
789 const Mg3d* mg = cv->get_mg() ;
790 int nz = mg->get_nzone() ;
791
792 // Le resultat
793 Mtbl* mti = new Mtbl(mg) ;
794 mti->set_etat_qcq() ;
795
796 // Pour le confort
797 const Valeur& alpha = cv->get_alpha() ;
798 Valeur dalphadt(alpha.get_mg()) ; dalphadt.annule_hard() ;
799 for (int k=0; k<mg->get_np(0); k++)
800 for (int j=0; j<mg->get_nt(0); j++){
801 const Grille3d* g = mg->get_grille3d(0) ;
802 double theta = (g->tet)[j] ;
803 double pphi = (g->phi)[k] ;
804 double sint = sin(theta) ;
805 double cost = cos(theta) ;
806 double sinp = sin(pphi) ;
807 double cosp = cos(pphi) ;
808 double aa = cv->get_aa() ;
809 double bb = cv->get_bb() ;
810 double cc = cv->get_cc() ;
811 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
812 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
813 dalphadt.set(0, k, j, 0) = alpha3 * sint * cost * (1./(cc*cc) - sinp*sinp / (bb*bb) - cosp*cosp / (aa*aa)) ;
814 }
815 const Valeur& beta = cv->get_beta() ;
816 const Valeur& dbetadt = beta.dsdt() ;
817
818
819 for (int l=0 ; l<nz ; l++) {
820
821 int nr = mg->get_nr(l);
822 int nt = mg->get_nt(l) ;
823 int np = mg->get_np(l) ;
824 const Grille3d* g = mg->get_grille3d(l) ;
825 Tbl* tb = (mti->t)[l] ;
826 tb->set_etat_qcq() ;
827 double* p_r = tb->t ;
828
829 switch(mg->get_type_r(l)) {
830
831
832 case RARE : {
833 for (int k=0 ; k<np ; k++) {
834 for (int j=0 ; j<nt ; j++) {
835 for (int i=0 ; i<nr ; i++) {
836 *p_r = dalphadt(l,k,j,0) / alpha(l,k,j,0) ;
837 p_r++ ;
838 } // Fin de boucle sur r
839 } // Fin de boucle sur theta
840 } // Fin de boucle sur phi
841 break ;
842 }
843
844 case FIN:{
845 for (int k=0 ; k<np ; k++) {
846 for (int j=0 ; j<nt ; j++) {
847 for (int i=0 ; i<nr ; i++) {
848 *p_r = (dalphadt(l,k,j,0)*(g->x)[i] + dbetadt(l,k,j,0)) / (alpha(l,k,j,0)*(g->x)[i] + beta(l,k,j,0)) ;
849 p_r++ ;
850 } // Fin de boucle sur r
851 } // Fin de boucle sur theta
852 } // Fin de boucle sur phi
853 break ;
854 }
855
856 case UNSURR: {
857 cout << "map_eps_fait_srdrdt: Compactified domain not allowed !" << endl;
858 abort() ;
859 break ;
860 }
861
862 default: {
863 cout << "map_eps_fait_srdrdt: unknown type_r !" << endl ;
864 abort() ;
865 }
866 } // Fin du switch
867 } // Fin de boucle sur zone
868
869 // Termine
870
871 return mti ;
872}
873
874/*
875 ************************************************************************
876 * 1/(R sin(theta)) dR/dphi
877 ************************************************************************
878 */
879
880Mtbl* map_eps_fait_srstdrdp(const Map* cvi) {
881
882 // recup du changement de variable
883 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
884 const Mg3d* mg = cv->get_mg() ;
885 int nz = mg->get_nzone() ;
886
887 // Le resultat
888 Mtbl* mti = new Mtbl(mg) ;
889 mti->set_etat_qcq() ;
890
891 // Pour le confort
892 const Valeur& alpha = cv->get_alpha() ;
893 Valeur stalphadp(alpha.get_mg()) ; stalphadp.annule_hard() ;
894 for (int k=0; k<mg->get_np(0); k++)
895 for (int j=0; j<mg->get_nt(0); j++){
896 const Grille3d* g = mg->get_grille3d(0) ;
897 double theta = (g->tet)[j] ;
898 double pphi = (g->phi)[k] ;
899 double sint = sin(theta) ;
900 double cost = cos(theta) ;
901 double sinp = sin(pphi) ;
902 double cosp = cos(pphi) ;
903 double aa = cv->get_aa() ;
904 double bb = cv->get_bb() ;
905 double cc = cv->get_cc() ;
906 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
907 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
908 stalphadp.set(0, k, j, 0) = alpha3 * sint * cosp * sinp * (1./(aa*aa) - 1./(bb*bb)) ; // Attention à la division par sin(theta) ici !!
909
910 }
911 const Valeur& beta = cv->get_beta() ;
912 const Valeur& stbetadp = beta.stdsdp() ;
913
914
915 for (int l=0 ; l<nz ; l++) {
916
917 int nr = mg->get_nr(l);
918 int nt = mg->get_nt(l) ;
919 int np = mg->get_np(l) ;
920 const Grille3d* g = mg->get_grille3d(l) ;
921 Tbl* tb = (mti->t)[l] ;
922 tb->set_etat_qcq() ;
923 double* p_r = tb->t ;
924
925 switch(mg->get_type_r(l)) {
926
927
928 case RARE : {
929 for (int k=0 ; k<np ; k++) {
930 for (int j=0 ; j<nt ; j++) {
931 for (int i=0 ; i<nr ; i++) {
932 *p_r = stalphadp(l,k,j,0) / alpha(l,k,j,0) ;
933 p_r++ ;
934 } // Fin de boucle sur r
935 } // Fin de boucle sur theta
936 } // Fin de boucle sur phi
937 break ;
938 }
939
940 case FIN:{
941 for (int k=0 ; k<np ; k++) {
942 for (int j=0 ; j<nt ; j++) {
943 for (int i=0 ; i<nr ; i++) {
944 *p_r = (stalphadp(l,k,j,0)*(g->x)[i] + stbetadp(l,k,j,0)) / (alpha(l,k,j,0)*(g->x)[i] + beta(l,k,j,0)) ;
945 p_r++ ;
946 } // Fin de boucle sur r
947 } // Fin de boucle sur theta
948 } // Fin de boucle sur phi
949 break ;
950 }
951
952 case UNSURR: {
953 cout << "map_eps_fait_srstdrdp: Compactified domain not allowed !" << endl;
954 abort() ;
955 break ;
956 }
957
958 default: {
959 cout << "map_eps_fait_srstdrdp: unknown type_r !" << endl ;
960 abort() ;
961 }
962 } // Fin du switch
963 } // Fin de boucle sur zone
964
965 // Termine
966
967 return mti ;
968}
969
970/*
971 ************************************************************************
972 * 1/R^2 dR/dtheta
973 ************************************************************************
974 */
975
976Mtbl* map_eps_fait_sr2drdt(const Map* cvi) {
977
978 // recup du changement de variable
979 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
980 const Mg3d* mg = cv->get_mg() ;
981 int nz = mg->get_nzone() ;
982
983 // Le resultat
984 Mtbl* mti = new Mtbl(mg) ;
985 mti->set_etat_qcq() ;
986
987 // Pour le confort
988 const Valeur& alpha = cv->get_alpha() ;
989 Valeur dalphadt(alpha.get_mg()) ; dalphadt.annule_hard() ;
990 for (int k=0; k<mg->get_np(0); k++)
991 for (int j=0; j<mg->get_nt(0); j++){
992 const Grille3d* g = mg->get_grille3d(0) ;
993 double theta = (g->tet)[j] ;
994 double pphi = (g->phi)[k] ;
995 double sint = sin(theta) ;
996 double cost = cos(theta) ;
997 double sinp = sin(pphi) ;
998 double cosp = cos(pphi) ;
999 double aa = cv->get_aa() ;
1000 double bb = cv->get_bb() ;
1001 double cc = cv->get_cc() ;
1002 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
1003 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
1004 dalphadt.set(0, k, j, 0) = alpha3 * sint * cost * (1./(cc*cc) - sinp*sinp / (bb*bb) - cosp*cosp / (aa*aa)) ;
1005 }
1006
1007
1008 for (int l=0 ; l<nz ; l++) {
1009
1010 int nr = mg->get_nr(l);
1011 int nt = mg->get_nt(l) ;
1012 int np = mg->get_np(l) ;
1013 const Grille3d* g = mg->get_grille3d(l) ;
1014 Tbl* tb = (mti->t)[l] ;
1015 tb->set_etat_qcq() ;
1016 double* p_r = tb->t ;
1017
1018 switch(mg->get_type_r(l)) {
1019
1020
1021 case RARE : {
1022 for (int k=0 ; k<np ; k++) {
1023 for (int j=0 ; j<nt ; j++) {
1024 for (int i=0 ; i<nr ; i++) {
1025 double ww = alpha(l,k,j,0) ;
1026 *p_r = dalphadt(l,k,j,0) / (ww*ww * (g->x)[i]) ;
1027 p_r++ ;
1028 } // Fin de boucle sur r
1029 } // Fin de boucle sur theta
1030 } // Fin de boucle sur phi
1031 break ;
1032 }
1033
1034 case FIN:{
1035 cout << "map_eps_fait_sr2drdt: Shells not implemented yet..." << endl;
1036 abort() ;
1037 break ;
1038 }
1039
1040 case UNSURR: {
1041 cout << "map_eps_fait_sr2drdt: Compactified domain not allowed !" << endl;
1042 abort() ;
1043 break ;
1044 }
1045
1046 default: {
1047 cout << "map_eps_fait_sr2drdt: unknown type_r !" << endl ;
1048 abort() ;
1049 }
1050 } // Fin du switch
1051 } // Fin de boucle sur zone
1052
1053 // Termine
1054
1055 return mti ;
1056}
1057
1058/*
1059 ************************************************************************
1060 * 1/(R^2 sin(theta)) dR/dphi
1061 ************************************************************************
1062 */
1063
1064Mtbl* map_eps_fait_sr2stdrdp(const Map* cvi) {
1065
1066 // recup du changement de variable
1067 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
1068 const Mg3d* mg = cv->get_mg() ;
1069 int nz = mg->get_nzone() ;
1070
1071 // Le resultat
1072 Mtbl* mti = new Mtbl(mg) ;
1073 mti->set_etat_qcq() ;
1074
1075 // Pour le confort
1076 const Valeur& alpha = cv->get_alpha() ;
1077 Valeur stalphadp(alpha.get_mg()) ; stalphadp.annule_hard() ;
1078 for (int k=0; k<mg->get_np(0); k++)
1079 for (int j=0; j<mg->get_nt(0); j++){
1080 const Grille3d* g = mg->get_grille3d(0) ;
1081 double theta = (g->tet)[j] ;
1082 double pphi = (g->phi)[k] ;
1083 double sint = sin(theta) ;
1084 double cost = cos(theta) ;
1085 double sinp = sin(pphi) ;
1086 double cosp = cos(pphi) ;
1087 double aa = cv->get_aa() ;
1088 double bb = cv->get_bb() ;
1089 double cc = cv->get_cc() ;
1090 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
1091 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
1092 stalphadp.set(0, k, j, 0) = alpha3 * sint * cosp * sinp * (1./(aa*aa) - 1./(bb*bb)) ; // Attention à la division par sin(theta) ici !!
1093
1094 }
1095
1096
1097 for (int l=0 ; l<nz ; l++) {
1098
1099 int nr = mg->get_nr(l);
1100 int nt = mg->get_nt(l) ;
1101 int np = mg->get_np(l) ;
1102 const Grille3d* g = mg->get_grille3d(l) ;
1103 Tbl* tb = (mti->t)[l] ;
1104 tb->set_etat_qcq() ;
1105 double* p_r = tb->t ;
1106
1107 switch(mg->get_type_r(l)) {
1108
1109
1110 case RARE : {
1111 for (int k=0 ; k<np ; k++) {
1112 for (int j=0 ; j<nt ; j++) {
1113 for (int i=0 ; i<nr ; i++) {
1114 double ww = alpha(l,k,j,0) ;
1115 *p_r = stalphadp(l,k,j,0) / (ww*ww * (g->x)[i]) ;
1116 p_r++ ;
1117 } // Fin de boucle sur r
1118 } // Fin de boucle sur theta
1119 } // Fin de boucle sur phi
1120 break ;
1121 }
1122
1123 case FIN:{
1124 cout << "map_eps_fait_sr2stdrdp: Shells not implemented yet..." << endl;
1125 abort() ;
1126 break ;
1127 }
1128
1129 case UNSURR: {
1130 cout << "map_eps_fait_sr2stdrdp: Compactified domain not allowed !" << endl;
1131 abort() ;
1132 break ;
1133 }
1134
1135 default: {
1136 cout << "map_eps_fait_sr2stdrdp: unknown type_r !" << endl ;
1137 abort() ;
1138 }
1139 } // Fin du switch
1140 } // Fin de boucle sur zone
1141
1142 // Termine
1143
1144 return mti ;
1145}
1146
1147/*
1148 ************************************************************************
1149 * d^2R/dx^2
1150 ************************************************************************
1151 */
1152
1153Mtbl* map_eps_fait_d2rdx2(const Map* cvi) {
1154
1155 // recup de la grille
1156 const Mg3d* mg = cvi->get_mg() ;
1157
1158 // Le resultat est nul :
1159 Mtbl* mti = new Mtbl(mg) ;
1160 mti->set_etat_zero() ;
1161
1162 return mti ;
1163}
1164
1165/*
1166 *****************************************************************************
1167 * 1/R^2 ( 1/sin(th) d/dth( sin(th) dR/dth ) + 1/sin(th)^2 d^2R/dphi^2 )
1168 *****************************************************************************
1169 */
1170
1171Mtbl* map_eps_fait_lapr_tp(const Map* cvi) {
1172
1173 // recup du changement de variable
1174 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
1175 const Mg3d* mg = cv->get_mg() ;
1176 int nz = mg->get_nzone() ;
1177
1178 // Le resultat
1179 Mtbl* mti = new Mtbl(mg) ;
1180 mti->set_etat_qcq() ;
1181
1182 // Pour le confort
1183 const Valeur& alpha = cv->get_alpha() ;
1184
1185 Valeur alpha_tmp = alpha ;
1186 alpha_tmp.ylm() ;
1187 const Valeur& lapalpha = alpha_tmp.lapang() ;
1188
1189
1190 for (int l=0 ; l<nz ; l++) {
1191
1192 const Grille3d* g = mg->get_grille3d(l) ;
1193
1194 int nr = mg->get_nr(l);
1195 int nt = mg->get_nt(l) ;
1196 int np = mg->get_np(l) ;
1197
1198 Tbl* tb = (mti->t)[l] ;
1199 tb->set_etat_qcq() ;
1200 double* p_r = tb->t ;
1201
1202 switch(mg->get_type_r(l)) {
1203
1204 case RARE: {
1205
1206 for (int k=0 ; k<np ; k++) {
1207 for (int j=0 ; j<nt ; j++) {
1208 for (int i=0 ; i<nr ; i++) {
1209 *p_r = (g->x)[i]*lapalpha(l,k,j,0) ;
1210 p_r++ ;
1211 } // Fin de boucle sur r
1212 } // Fin de boucle sur theta
1213 } // Fin de boucle sur phi
1214 break ;
1215 }
1216
1217 case FIN:{
1218 cout << "map_eps_fait_lapr_tp: Shells not implemented yet..." << endl;
1219 abort() ;
1220 break ;
1221 }
1222
1223 case UNSURR: {
1224 cout << "map_eps_fait_lapr_tp: Compactified domain not allowed !" << endl;
1225 abort() ;
1226 break ;
1227 }
1228
1229 default: {
1230 cout << "map_eps_fait_lapr_tp: unknown type_r !" << endl ;
1231 abort() ;
1232 }
1233 } // Fin du switch
1234 } // Fin de boucle sur zone
1235
1236 // Termine
1237
1238 return mti ;
1239}
1240
1241/*
1242 ************************************************************************
1243 * d^2R/dthdx
1244 ************************************************************************
1245 */
1246
1247Mtbl* map_eps_fait_d2rdtdx(const Map* cvi) {
1248cerr << "Map_eps::map_eps_fait_d2rdtdx pas fait" << endl;
1249abort() ;
1250 // recup du changement de variable
1251 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
1252 const Mg3d* mg = cv->get_mg() ;
1253 int nz = mg->get_nzone() ;
1254
1255 // Le resultat
1256 Mtbl* mti = new Mtbl(mg) ;
1257 mti->set_etat_qcq() ;
1258
1259 // Pour le confort
1260 const Valeur& alpha = cv->get_alpha() ;
1261 const Valeur& d2alphadxdt = alpha.dsdt().dsdx() ;
1262
1263
1264 for (int l=0 ; l<nz ; l++) {
1265
1266 int nr = mg->get_nr(l);
1267 int nt = mg->get_nt(l) ;
1268 int np = mg->get_np(l) ;
1269
1270 Tbl* tb = (mti->t)[l] ;
1271 tb->set_etat_qcq() ;
1272 double* p_r = tb->t ;
1273
1274 switch(mg->get_type_r(l)) {
1275
1276
1277 case RARE : {
1278 for (int k=0 ; k<np ; k++) {
1279 for (int j=0 ; j<nt ; j++) {
1280 for (int i=0 ; i<nr ; i++) {
1281 *p_r = d2alphadxdt(l,k,j,0)/alpha(l,k,j,0) ;
1282 p_r++ ;
1283 } // Fin de boucle sur r
1284 } // Fin de boucle sur theta
1285 } // Fin de boucle sur phi
1286 break ;
1287 }
1288 case FIN:{
1289 cout << "map_eps_fait_d2rdtdx: Shells not implemented yet..." << endl;
1290 abort() ;
1291 break ;
1292 }
1293
1294 case UNSURR: {
1295 cout << "map_eps_fait_d2rdtdx: Compactified domain not allowed !" << endl;
1296 abort() ;
1297 break ;
1298 }
1299
1300 default: {
1301 cout << "map_eps_fait_d2rdtdx: unknown type_r !" << endl ;
1302 abort() ;
1303 }
1304 } // Fin du switch
1305 } // Fin de boucle sur zone
1306
1307 // Termine
1308
1309 return mti ;
1310}
1311
1312/*
1313 ************************************************************************
1314 * 1/sin(th) d^2R/dphidx
1315 ************************************************************************
1316 */
1317
1318Mtbl* map_eps_fait_sstd2rdpdx(const Map* cvi) {
1319cerr << "map_eps_fait_sstd2rdpdx pas fait" << endl;
1320abort() ;
1321 // recup du changement de variable
1322 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
1323 const Mg3d* mg = cv->get_mg() ;
1324 int nz = mg->get_nzone() ;
1325
1326 // Le resultat
1327 Mtbl* mti = new Mtbl(mg) ;
1328 mti->set_etat_qcq() ;
1329
1330 // Pour le confort
1331 const Valeur& alpha = cv->get_alpha() ;
1332 Valeur stalphadp(alpha.get_mg()) ; stalphadp.annule_hard() ;
1333 for (int k=0; k<mg->get_np(0); k++)
1334 for (int j=0; j<mg->get_nt(0); j++){
1335 const Grille3d* g = mg->get_grille3d(0) ;
1336 double theta = (g->tet)[j] ;
1337 double pphi = (g->phi)[k] ;
1338 double sint = sin(theta) ;
1339 double cost = cos(theta) ;
1340 double sinp = sin(pphi) ;
1341 double cosp = cos(pphi) ;
1342 double aa = cv->get_aa() ;
1343 double bb = cv->get_bb() ;
1344 double cc = cv->get_cc() ;
1345 double alpha2 = alpha(0, k, j, 0)*alpha(0, k, j, 0) ;
1346 double alpha3 = alpha2 * alpha(0, k, j, 0) ;
1347 stalphadp.set(0, k, j, 0) = alpha3 * sint * cosp * sinp * (1./(aa*aa) - 1./(bb*bb)) ; // Attention à la division par sin(theta) ici !!
1348
1349 }
1350
1351
1352 for (int l=0 ; l<nz ; l++) {
1353
1354 int nr = mg->get_nr(l);
1355 int nt = mg->get_nt(l) ;
1356 int np = mg->get_np(l) ;
1357 const Grille3d* g = mg->get_grille3d(l) ;
1358 Tbl* tb = (mti->t)[l] ;
1359 tb->set_etat_qcq() ;
1360 double* p_r = tb->t ;
1361
1362 switch(mg->get_type_r(l)) {
1363
1364
1365 case RARE : {
1366 for (int k=0 ; k<np ; k++) {
1367 for (int j=0 ; j<nt ; j++) {
1368 for (int i=0 ; i<nr ; i++) {
1369 *p_r = stalphadp(l,k,j,0) ;
1370 p_r++ ;
1371 } // Fin de boucle sur r
1372 } // Fin de boucle sur theta
1373 } // Fin de boucle sur phi
1374 break ;
1375 }
1376
1377 case FIN:{
1378 cout << "map_eps_fait_sstd2rdpdx: Shells not implemented yet..." << endl;
1379 abort() ;
1380 break ;
1381 }
1382
1383 case UNSURR: {
1384 cout << "map_eps_fait_sstd2rdpdx: Compactified domain not allowed !" << endl;
1385 abort() ;
1386 break ;
1387 }
1388
1389 default: {
1390 cout << "map_eps_fait_sstd2rdpdx: unknown type_r !" << endl ;
1391 abort() ;
1392 }
1393 } // Fin du switch
1394 } // Fin de boucle sur zone
1395
1396 // Termine
1397
1398
1399 return mti ;
1400}
1401
1402/*
1403 ************************************************************************
1404 * 1/R^2 d^2R/dtheta^2
1405 ************************************************************************
1406 */
1407
1408Mtbl* map_eps_fait_sr2d2rdt2(const Map* cvi) {
1409cerr << "map_eps_fait_sr2d2rdt2 pas fait" << endl;
1410abort() ;
1411 // recup du changement de variable
1412 const Map_eps* cv = static_cast<const Map_eps*>(cvi) ;
1413 const Mg3d* mg = cv->get_mg() ;
1414 int nz = mg->get_nzone() ;
1415
1416 // Le resultat
1417 Mtbl* mti = new Mtbl(mg) ;
1418 mti->set_etat_qcq() ;
1419
1420 // Pour le confort
1421 const Valeur& alpha = cv->get_alpha() ;
1422 const Valeur& d2alphadt2 = alpha.d2sdt2() ;
1423
1424
1425 for (int l=0 ; l<nz ; l++) {
1426
1427 int nr = mg->get_nr(l);
1428 int nt = mg->get_nt(l) ;
1429 int np = mg->get_np(l) ;
1430 const Grille3d* g = mg->get_grille3d(l) ;
1431 Tbl* tb = (mti->t)[l] ;
1432 tb->set_etat_qcq() ;
1433 double* p_r = tb->t ;
1434
1435 switch(mg->get_type_r(l)) {
1436
1437
1438 case RARE : {
1439 for (int k=0 ; k<np ; k++) {
1440 for (int j=0 ; j<nt ; j++) {
1441 for (int i=0 ; i<nr ; i++) {
1442 double ww = alpha(l,k,j,0) ;
1443 *p_r = d2alphadt2(l,k,j,0) / (ww*ww * (g->x)[i]) ;
1444 p_r++ ;
1445 } // Fin de boucle sur r
1446 } // Fin de boucle sur theta
1447 } // Fin de boucle sur phi
1448 break ;
1449 }
1450 case FIN:{
1451 cout << "map_eps_fait_drdt: Shells not implemented yet..." << endl;
1452 abort() ;
1453 break ;
1454 }
1455
1456 case UNSURR: {
1457 cout << "map_eps_fait_drdt: Compactified domain not allowed !" << endl;
1458 abort() ;
1459 break ;
1460 }
1461
1462 default: {
1463 cout << "map_eps_fait_drdt: unknown type_r !" << endl ;
1464 abort() ;
1465 }
1466 } // Fin du switch
1467 } // Fin de boucle sur zone
1468
1469 // Termine
1470
1471 return mti ;
1472}
1473
1474}
Lorene::Grille3d
3D grid class in one domain.
Definition grilles.h:200
Lorene::Grille3d::x
double * x
Array of values of at the nr collocation points.
Definition grilles.h:215
Lorene::Grille3d::tet
double * tet
Array of values of at the nt collocation points.
Definition grilles.h:217
Lorene::Grille3d::get_nr
int get_nr() const
Returns nr.
Definition grilles.h:236
Lorene::Map_eps
Definition map.h:4210
Lorene::Map_eps::Map_eps
Map_eps(const Mg3d &mgrille, const double *r_limits)
Standard Constructor.
Definition map_eps.C:52
Lorene::Map_eps::get_alpha
const Valeur & get_alpha() const
Returns the reference on the Tbl alpha.
Definition map_eps.C:243
Lorene::Map_eps::map_eps_fait_r
friend Mtbl * map_eps_fait_r(const Map *)
< Not implemented
Definition map_eps_fait.C:43
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_grille3d
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:517
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_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::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:491
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::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:302
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::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::d2sdt2
const Valeur & d2sdt2() const
Returns of *this.
Definition valeur_d2sdt2.C:110
Lorene::Valeur::stdsdp
const Valeur & stdsdp() const
Returns of *this.
Definition valeur_stdsdp.C:63
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::dsdt
const Valeur & dsdt() const
Returns of *this.
Definition valeur_dsdt.C:115
Lorene::Valeur::dsdx
const Valeur & dsdx() const
Returns of *this.
Definition valeur_dsdx.C:114
Lorene::Valeur::lapang
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition valeur_lapang.C:75
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::cost
Coord cost
Definition map.h:734
Lorene::cosp
Coord cosp
Definition map.h:736
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
Lorene::sinp
Coord sinp
Definition map.h:735
Lorene::sint
Coord sint
Definition map.h:733