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 
42 namespace Lorene {
43 Mtbl* 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 
113 Mtbl* 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 
151 Mtbl* 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 
189 Mtbl* 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 
207 Mtbl* 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 
225 Mtbl* 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 
243 Mtbl* 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 
273 Mtbl* 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 
303 Mtbl* 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 
323 Mtbl* 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 
357 Mtbl* 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 
391 Mtbl* 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 
425 Mtbl* 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 
465 Mtbl* 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 
543 Mtbl* 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 
607 Mtbl* 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 
695 Mtbl* 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 
785 Mtbl* 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 
880 Mtbl* 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 
976 Mtbl* 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 
1064 Mtbl* 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 
1153 Mtbl* 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 
1171 Mtbl* 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 
1247 Mtbl* map_eps_fait_d2rdtdx(const Map* cvi) {
1248 cerr << "Map_eps::map_eps_fait_d2rdtdx pas fait" << endl;
1249 abort() ;
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 
1318 Mtbl* map_eps_fait_sstd2rdpdx(const Map* cvi) {
1319 cerr << "map_eps_fait_sstd2rdpdx pas fait" << endl;
1320 abort() ;
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 
1408 Mtbl* map_eps_fait_sr2d2rdt2(const Map* cvi) {
1409 cerr << "map_eps_fait_sr2d2rdt2 pas fait" << endl;
1410 abort() ;
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::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::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:688
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Lorene::Grille3d::x
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:215
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::Map_eps::get_alpha
const Valeur & get_alpha() const
Returns the reference on the Tbl alpha.
Definition: map_eps.C:243
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:176
Lorene::Grille3d::tet
double * tet
Array of values of at the nt collocation points.
Definition: grilles.h:217
Lorene::Map_eps
Definition: map_star.h:322
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
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
Multi-domain grid.
Definition: grilles.h:279
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
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_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::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Lorene::Mtbl::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the Mtbl is defined.
Definition: mtbl.h:274
Lorene::Grille3d
3D grid class in one domain.
Definition: grilles.h:200