LORENE
map_star_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_star_fait_r(const Map* cvi) {
44 
45  // recup du changement de variable
46  const Map_star* cv = static_cast<const Map_star*>(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_star_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_star_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_star_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_star_fait_tet(const Map* cvi) {
114 
115  // recup du changement de variable
116  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_phi(const Map* cvi) {
152 
153  // recup du changement de variable
154  const Map_star* cv = static_cast<const Map_star*>(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_star_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_star_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_star_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_star_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_star_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_star_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_star_fait_sint(const Map* cvi) {
324 
325  // recup du changement de variable
326  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_cost(const Map* cvi) {
358 
359  // recup du changement de variable
360  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_sinp(const Map* cvi) {
392 
393  // recup du changement de variable
394  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_cosp(const Map* cvi) {
426 
427  // recup du changement de variable
428  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_xsr(const Map* cvi) {
466 
467  // recup du changement de variable
468  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_xsr: Compactified domain not allowed !" << endl;
522  abort() ;
523  break ;
524 
525  default:
526  cout << "map_star_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_star_fait_dxdr(const Map* cvi) {
544 
545  // recup du changement de variable
546  const Map_star* cv = static_cast<const Map_star*>(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_star_fait_dxdr: Shells not implemented yet..." << endl;
581  // abort() ;
582  // break ;
583  // }
584 
585  case UNSURR:
586  cout << "map_star_fait_dxdr: Compactified domain not allowed !" << endl;
587  abort() ;
588  break ;
589 
590 
591  default:
592  cout << "map_star_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_star_fait_drdt(const Map* cvi) {
608 
609  // recup du changement de variable
610  const Map_star* cv = static_cast<const Map_star*>(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  const Valeur& dalphadt = alpha.dsdt() ;
621  const Valeur& beta = cv->get_beta() ;
622  const Valeur& dbetadt = beta.dsdt() ;
623 
624 
625  for (int l=0 ; l<nz ; l++) {
626 
627  int nr = mg->get_nr(l);
628  int nt = mg->get_nt(l) ;
629  int np = mg->get_np(l) ;
630  const Grille3d* g = mg->get_grille3d(l) ;
631  Tbl* tb = (mti->t)[l] ;
632  tb->set_etat_qcq() ;
633  double* p_r = tb->t ;
634 
635  switch(mg->get_type_r(l)) {
636 
637 
638  case RARE: case FIN: {
639  for (int k=0 ; k<np ; k++) {
640  for (int j=0 ; j<nt ; j++) {
641  for (int i=0 ; i<nr ; i++) {
642  *p_r = dalphadt(l,k,j,0) * (g->x)[i] + dbetadt(l,k,j,0) ;
643  p_r++ ;
644  } // Fin de boucle sur r
645  } // Fin de boucle sur theta
646  } // Fin de boucle sur phi
647  break ;
648  }
649  // case FIN:{
650  // cout << "map_star_fait_drdt: Shells not implemented yet..." << endl;
651  // abort() ;
652  // break ;
653  // }
654 
655  case UNSURR: {
656  cout << "map_star_fait_drdt: Compactified domain not allowed !" << endl;
657  abort() ;
658  break ;
659  }
660 
661  default: {
662  cout << "map_star_fait_drdt: unknown type_r !" << endl ;
663  abort() ;
664  }
665  } // Fin du switch
666  } // Fin de boucle sur zone
667 
668  // Termine
669 
670  return mti ;
671 }
672 
673 /*
674  ************************************************************************
675  * 1/sin(theta) dR/dphi
676  ************************************************************************
677  */
678 
679 Mtbl* map_star_fait_stdrdp(const Map* cvi) {
680 
681  // recup du changement de variable
682  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
683  const Mg3d* mg = cv->get_mg() ;
684  int nz = mg->get_nzone() ;
685 
686  // Le resultat
687  Mtbl* mti = new Mtbl(mg) ;
688  mti->set_etat_qcq() ;
689 
690  // Pour le confort
691  const Valeur& alpha = cv->get_alpha() ;
692  const Valeur& stalphadp = alpha.stdsdp() ;
693  const Valeur& beta = cv->get_beta() ;
694  const Valeur& stbetadp = beta.stdsdp() ;
695 
696 
697  for (int l=0 ; l<nz ; l++) {
698 
699  int nr = mg->get_nr(l);
700  int nt = mg->get_nt(l) ;
701  int np = mg->get_np(l) ;
702  const Grille3d* g = mg->get_grille3d(l) ;
703  Tbl* tb = (mti->t)[l] ;
704  tb->set_etat_qcq() ;
705  double* p_r = tb->t ;
706 
707  switch(mg->get_type_r(l)) {
708 
709 
710  case RARE: case FIN: {
711  for (int k=0 ; k<np ; k++) {
712  for (int j=0 ; j<nt ; j++) {
713  for (int i=0 ; i<nr ; i++) {
714  *p_r = stalphadp(l,k,j,0) * (g->x)[i] + stbetadp(l,k,j,0) ;
715  p_r++ ;
716  } // Fin de boucle sur r
717  } // Fin de boucle sur theta
718  } // Fin de boucle sur phi
719  break ;
720  }
721 
722  // case FIN:{
723  // cout << "map_star_fait_stdrdp: Shells not implemented yet..." << endl;
724  // abort() ;
725  // break ;
726  // }
727 
728  case UNSURR: {
729  cout << "map_star_fait_stdrdp: Compactified domain not allowed !" << endl;
730  abort() ;
731  break ;
732  }
733 
734  default: {
735  cout << "map_star_fait_stdrdp: unknown type_r !" << endl ;
736  abort() ;
737  }
738  } // Fin du switch
739  } // Fin de boucle sur zone
740 
741  // Termine
742 
743  return mti ;
744 }
745 
746 /*
747  ************************************************************************
748  * 1/R dR/dtheta
749  ************************************************************************
750  */
751 
752 Mtbl* map_star_fait_srdrdt(const Map* cvi) {
753 
754  // recup du changement de variable
755  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
756  const Mg3d* mg = cv->get_mg() ;
757  int nz = mg->get_nzone() ;
758 
759  // Le resultat
760  Mtbl* mti = new Mtbl(mg) ;
761  mti->set_etat_qcq() ;
762 
763  // Pour le confort
764  const Valeur& alpha = cv->get_alpha() ;
765  const Valeur& dalphadt = alpha.dsdt() ;
766  const Valeur& beta = cv->get_beta() ;
767  const Valeur& dbetadt = beta.dsdt() ;
768 
769 
770  for (int l=0 ; l<nz ; l++) {
771 
772  int nr = mg->get_nr(l);
773  int nt = mg->get_nt(l) ;
774  int np = mg->get_np(l) ;
775  const Grille3d* g = mg->get_grille3d(l) ;
776  Tbl* tb = (mti->t)[l] ;
777  tb->set_etat_qcq() ;
778  double* p_r = tb->t ;
779 
780  switch(mg->get_type_r(l)) {
781 
782 
783  case RARE : {
784  for (int k=0 ; k<np ; k++) {
785  for (int j=0 ; j<nt ; j++) {
786  for (int i=0 ; i<nr ; i++) {
787  *p_r = dalphadt(l,k,j,0) / alpha(l,k,j,0) ;
788  p_r++ ;
789  } // Fin de boucle sur r
790  } // Fin de boucle sur theta
791  } // Fin de boucle sur phi
792  break ;
793  }
794 
795  case FIN:{
796  for (int k=0 ; k<np ; k++) {
797  for (int j=0 ; j<nt ; j++) {
798  for (int i=0 ; i<nr ; i++) {
799  *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)) ;
800  p_r++ ;
801  } // Fin de boucle sur r
802  } // Fin de boucle sur theta
803  } // Fin de boucle sur phi
804  break ;
805  }
806 
807  case UNSURR: {
808  cout << "map_star_fait_srdrdt: Compactified domain not allowed !" << endl;
809  abort() ;
810  break ;
811  }
812 
813  default: {
814  cout << "map_star_fait_srdrdt: unknown type_r !" << endl ;
815  abort() ;
816  }
817  } // Fin du switch
818  } // Fin de boucle sur zone
819 
820  // Termine
821 
822  return mti ;
823 }
824 
825 /*
826  ************************************************************************
827  * 1/(R sin(theta)) dR/dphi
828  ************************************************************************
829  */
830 
831 Mtbl* map_star_fait_srstdrdp(const Map* cvi) {
832 
833  // recup du changement de variable
834  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
835  const Mg3d* mg = cv->get_mg() ;
836  int nz = mg->get_nzone() ;
837 
838  // Le resultat
839  Mtbl* mti = new Mtbl(mg) ;
840  mti->set_etat_qcq() ;
841 
842  // Pour le confort
843  const Valeur& alpha = cv->get_alpha() ;
844  const Valeur& stalphadp = alpha.stdsdp() ;
845  const Valeur& beta = cv->get_beta() ;
846  const Valeur& stbetadp = beta.stdsdp() ;
847 
848 
849  for (int l=0 ; l<nz ; l++) {
850 
851  int nr = mg->get_nr(l);
852  int nt = mg->get_nt(l) ;
853  int np = mg->get_np(l) ;
854  const Grille3d* g = mg->get_grille3d(l) ;
855  Tbl* tb = (mti->t)[l] ;
856  tb->set_etat_qcq() ;
857  double* p_r = tb->t ;
858 
859  switch(mg->get_type_r(l)) {
860 
861 
862  case RARE : {
863  for (int k=0 ; k<np ; k++) {
864  for (int j=0 ; j<nt ; j++) {
865  for (int i=0 ; i<nr ; i++) {
866  *p_r = stalphadp(l,k,j,0) / alpha(l,k,j,0) ;
867  p_r++ ;
868  } // Fin de boucle sur r
869  } // Fin de boucle sur theta
870  } // Fin de boucle sur phi
871  break ;
872  }
873 
874  case FIN:{
875  for (int k=0 ; k<np ; k++) {
876  for (int j=0 ; j<nt ; j++) {
877  for (int i=0 ; i<nr ; i++) {
878  *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)) ;
879  p_r++ ;
880  } // Fin de boucle sur r
881  } // Fin de boucle sur theta
882  } // Fin de boucle sur phi
883  break ;
884  }
885 
886  case UNSURR: {
887  cout << "map_star_fait_srstdrdp: Compactified domain not allowed !" << endl;
888  abort() ;
889  break ;
890  }
891 
892  default: {
893  cout << "map_star_fait_srstdrdp: unknown type_r !" << endl ;
894  abort() ;
895  }
896  } // Fin du switch
897  } // Fin de boucle sur zone
898 
899  // Termine
900 
901  return mti ;
902 }
903 
904 /*
905  ************************************************************************
906  * 1/R^2 dR/dtheta
907  ************************************************************************
908  */
909 
910 Mtbl* map_star_fait_sr2drdt(const Map* cvi) {
911 
912  // recup du changement de variable
913  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
914  const Mg3d* mg = cv->get_mg() ;
915  int nz = mg->get_nzone() ;
916 
917  // Le resultat
918  Mtbl* mti = new Mtbl(mg) ;
919  mti->set_etat_qcq() ;
920 
921  // Pour le confort
922  const Valeur& alpha = cv->get_alpha() ;
923  const Valeur& dalphadt = alpha.dsdt() ;
924 
925 
926  for (int l=0 ; l<nz ; l++) {
927 
928  int nr = mg->get_nr(l);
929  int nt = mg->get_nt(l) ;
930  int np = mg->get_np(l) ;
931  const Grille3d* g = mg->get_grille3d(l) ;
932  Tbl* tb = (mti->t)[l] ;
933  tb->set_etat_qcq() ;
934  double* p_r = tb->t ;
935 
936  switch(mg->get_type_r(l)) {
937 
938 
939  case RARE : {
940  for (int k=0 ; k<np ; k++) {
941  for (int j=0 ; j<nt ; j++) {
942  for (int i=0 ; i<nr ; i++) {
943  double ww = alpha(l,k,j,0) ;
944  *p_r = dalphadt(l,k,j,0) / (ww*ww * (g->x)[i]) ;
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 FIN:{
953  cout << "map_star_fait_sr2drdt: Shells not implemented yet..." << endl;
954  abort() ;
955  break ;
956  }
957 
958  case UNSURR: {
959  cout << "map_star_fait_sr2drdt: Compactified domain not allowed !" << endl;
960  abort() ;
961  break ;
962  }
963 
964  default: {
965  cout << "map_star_fait_sr2drdt: unknown type_r !" << endl ;
966  abort() ;
967  }
968  } // Fin du switch
969  } // Fin de boucle sur zone
970 
971  // Termine
972 
973  return mti ;
974 }
975 
976 /*
977  ************************************************************************
978  * 1/(R^2 sin(theta)) dR/dphi
979  ************************************************************************
980  */
981 
982 Mtbl* map_star_fait_sr2stdrdp(const Map* cvi) {
983 
984  // recup du changement de variable
985  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
986  const Mg3d* mg = cv->get_mg() ;
987  int nz = mg->get_nzone() ;
988 
989  // Le resultat
990  Mtbl* mti = new Mtbl(mg) ;
991  mti->set_etat_qcq() ;
992 
993  // Pour le confort
994  const Valeur& alpha = cv->get_alpha() ;
995  const Valeur& stalphadp = alpha.stdsdp() ;
996 
997 
998  for (int l=0 ; l<nz ; l++) {
999 
1000  int nr = mg->get_nr(l);
1001  int nt = mg->get_nt(l) ;
1002  int np = mg->get_np(l) ;
1003  const Grille3d* g = mg->get_grille3d(l) ;
1004  Tbl* tb = (mti->t)[l] ;
1005  tb->set_etat_qcq() ;
1006  double* p_r = tb->t ;
1007 
1008  switch(mg->get_type_r(l)) {
1009 
1010 
1011  case RARE : {
1012  for (int k=0 ; k<np ; k++) {
1013  for (int j=0 ; j<nt ; j++) {
1014  for (int i=0 ; i<nr ; i++) {
1015  double ww = alpha(l,k,j,0) ;
1016  *p_r = stalphadp(l,k,j,0) / (ww*ww * (g->x)[i]) ;
1017  p_r++ ;
1018  } // Fin de boucle sur r
1019  } // Fin de boucle sur theta
1020  } // Fin de boucle sur phi
1021  break ;
1022  }
1023 
1024  case FIN:{
1025  cout << "map_star_fait_sr2stdrdp: Shells not implemented yet..." << endl;
1026  abort() ;
1027  break ;
1028  }
1029 
1030  case UNSURR: {
1031  cout << "map_star_fait_sr2stdrdp: Compactified domain not allowed !" << endl;
1032  abort() ;
1033  break ;
1034  }
1035 
1036  default: {
1037  cout << "map_star_fait_sr2stdrdp: unknown type_r !" << endl ;
1038  abort() ;
1039  }
1040  } // Fin du switch
1041  } // Fin de boucle sur zone
1042 
1043  // Termine
1044 
1045  return mti ;
1046 }
1047 
1048 /*
1049  ************************************************************************
1050  * d^2R/dx^2
1051  ************************************************************************
1052  */
1053 
1054 Mtbl* map_star_fait_d2rdx2(const Map* cvi) {
1055 
1056  // recup de la grille
1057  const Mg3d* mg = cvi->get_mg() ;
1058 
1059  // Le resultat est nul :
1060  Mtbl* mti = new Mtbl(mg) ;
1061  mti->set_etat_zero() ;
1062 
1063  return mti ;
1064 }
1065 
1066 /*
1067  *****************************************************************************
1068  * 1/R^2 ( 1/sin(th) d/dth( sin(th) dR/dth ) + 1/sin(th)^2 d^2R/dphi^2 )
1069  *****************************************************************************
1070  */
1071 
1072 Mtbl* map_star_fait_lapr_tp(const Map* cvi) {
1073 
1074  // recup du changement de variable
1075  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
1076  const Mg3d* mg = cv->get_mg() ;
1077  int nz = mg->get_nzone() ;
1078 
1079  // Le resultat
1080  Mtbl* mti = new Mtbl(mg) ;
1081  mti->set_etat_qcq() ;
1082 
1083  // Pour le confort
1084  const Valeur& alpha = cv->get_alpha() ;
1085 
1086  Valeur alpha_tmp = alpha ;
1087  alpha_tmp.ylm() ;
1088  const Valeur& lapalpha = alpha_tmp.lapang() ;
1089 
1090 
1091  for (int l=0 ; l<nz ; l++) {
1092 
1093  const Grille3d* g = mg->get_grille3d(l) ;
1094 
1095  int nr = mg->get_nr(l);
1096  int nt = mg->get_nt(l) ;
1097  int np = mg->get_np(l) ;
1098 
1099  Tbl* tb = (mti->t)[l] ;
1100  tb->set_etat_qcq() ;
1101  double* p_r = tb->t ;
1102 
1103  switch(mg->get_type_r(l)) {
1104 
1105  case RARE: {
1106 
1107  for (int k=0 ; k<np ; k++) {
1108  for (int j=0 ; j<nt ; j++) {
1109  for (int i=0 ; i<nr ; i++) {
1110  *p_r = (g->x)[i]*lapalpha(l,k,j,0) ;
1111  p_r++ ;
1112  } // Fin de boucle sur r
1113  } // Fin de boucle sur theta
1114  } // Fin de boucle sur phi
1115  break ;
1116  }
1117 
1118  case FIN:{
1119  cout << "map_star_fait_lapr_tp: Shells not implemented yet..." << endl;
1120  abort() ;
1121  break ;
1122  }
1123 
1124  case UNSURR: {
1125  cout << "map_star_fait_lapr_tp: Compactified domain not allowed !" << endl;
1126  abort() ;
1127  break ;
1128  }
1129 
1130  default: {
1131  cout << "map_star_fait_lapr_tp: unknown type_r !" << endl ;
1132  abort() ;
1133  }
1134  } // Fin du switch
1135  } // Fin de boucle sur zone
1136 
1137  // Termine
1138 
1139  return mti ;
1140 }
1141 
1142 /*
1143  ************************************************************************
1144  * d^2R/dthdx
1145  ************************************************************************
1146  */
1147 
1148 Mtbl* map_star_fait_d2rdtdx(const Map* cvi) {
1149 
1150  // recup du changement de variable
1151  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
1152  const Mg3d* mg = cv->get_mg() ;
1153  int nz = mg->get_nzone() ;
1154 
1155  // Le resultat
1156  Mtbl* mti = new Mtbl(mg) ;
1157  mti->set_etat_qcq() ;
1158 
1159  // Pour le confort
1160  const Valeur& alpha = cv->get_alpha() ;
1161  const Valeur& d2alphadxdt = alpha.dsdt().dsdx() ;
1162 
1163 
1164  for (int l=0 ; l<nz ; l++) {
1165 
1166  int nr = mg->get_nr(l);
1167  int nt = mg->get_nt(l) ;
1168  int np = mg->get_np(l) ;
1169 
1170  Tbl* tb = (mti->t)[l] ;
1171  tb->set_etat_qcq() ;
1172  double* p_r = tb->t ;
1173 
1174  switch(mg->get_type_r(l)) {
1175 
1176 
1177  case RARE : {
1178  for (int k=0 ; k<np ; k++) {
1179  for (int j=0 ; j<nt ; j++) {
1180  for (int i=0 ; i<nr ; i++) {
1181  *p_r = d2alphadxdt(l,k,j,0)/alpha(l,k,j,0) ;
1182  p_r++ ;
1183  } // Fin de boucle sur r
1184  } // Fin de boucle sur theta
1185  } // Fin de boucle sur phi
1186  break ;
1187  }
1188  case FIN:{
1189  cout << "map_star_fait_d2rdtdx: Shells not implemented yet..." << endl;
1190  abort() ;
1191  break ;
1192  }
1193 
1194  case UNSURR: {
1195  cout << "map_star_fait_d2rdtdx: Compactified domain not allowed !" << endl;
1196  abort() ;
1197  break ;
1198  }
1199 
1200  default: {
1201  cout << "map_star_fait_d2rdtdx: unknown type_r !" << endl ;
1202  abort() ;
1203  }
1204  } // Fin du switch
1205  } // Fin de boucle sur zone
1206 
1207  // Termine
1208 
1209  return mti ;
1210 }
1211 
1212 /*
1213  ************************************************************************
1214  * 1/sin(th) d^2R/dphidx
1215  ************************************************************************
1216  */
1217 
1218 Mtbl* map_star_fait_sstd2rdpdx(const Map* cvi) {
1219 
1220  // recup du changement de variable
1221  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
1222  const Mg3d* mg = cv->get_mg() ;
1223  int nz = mg->get_nzone() ;
1224 
1225  // Le resultat
1226  Mtbl* mti = new Mtbl(mg) ;
1227  mti->set_etat_qcq() ;
1228 
1229  // Pour le confort
1230  const Valeur& alpha = cv->get_alpha() ;
1231  const Valeur& stalphadp = alpha.stdsdp() ;
1232 
1233 
1234  for (int l=0 ; l<nz ; l++) {
1235 
1236  int nr = mg->get_nr(l);
1237  int nt = mg->get_nt(l) ;
1238  int np = mg->get_np(l) ;
1239  // const Grille3d* g = mg->get_grille3d(l) ;
1240  Tbl* tb = (mti->t)[l] ;
1241  tb->set_etat_qcq() ;
1242  double* p_r = tb->t ;
1243 
1244  switch(mg->get_type_r(l)) {
1245 
1246 
1247  case RARE : {
1248  for (int k=0 ; k<np ; k++) {
1249  for (int j=0 ; j<nt ; j++) {
1250  for (int i=0 ; i<nr ; i++) {
1251  *p_r = stalphadp(l,k,j,0) ;
1252  p_r++ ;
1253  } // Fin de boucle sur r
1254  } // Fin de boucle sur theta
1255  } // Fin de boucle sur phi
1256  break ;
1257  }
1258 
1259  case FIN:{
1260  cout << "map_star_fait_sstd2rdpdx: Shells not implemented yet..." << endl;
1261  abort() ;
1262  break ;
1263  }
1264 
1265  case UNSURR: {
1266  cout << "map_star_fait_sstd2rdpdx: Compactified domain not allowed !" << endl;
1267  abort() ;
1268  break ;
1269  }
1270 
1271  default: {
1272  cout << "map_star_fait_sstd2rdpdx: unknown type_r !" << endl ;
1273  abort() ;
1274  }
1275  } // Fin du switch
1276  } // Fin de boucle sur zone
1277 
1278  // Termine
1279 
1280 
1281  return mti ;
1282 }
1283 
1284 /*
1285  ************************************************************************
1286  * 1/R^2 d^2R/dtheta^2
1287  ************************************************************************
1288  */
1289 
1290 Mtbl* map_star_fait_sr2d2rdt2(const Map* cvi) {
1291 
1292  // recup du changement de variable
1293  const Map_star* cv = static_cast<const Map_star*>(cvi) ;
1294  const Mg3d* mg = cv->get_mg() ;
1295  int nz = mg->get_nzone() ;
1296 
1297  // Le resultat
1298  Mtbl* mti = new Mtbl(mg) ;
1299  mti->set_etat_qcq() ;
1300 
1301  // Pour le confort
1302  const Valeur& alpha = cv->get_alpha() ;
1303  const Valeur& d2alphadt2 = alpha.d2sdt2() ;
1304 
1305 
1306  for (int l=0 ; l<nz ; l++) {
1307 
1308  int nr = mg->get_nr(l);
1309  int nt = mg->get_nt(l) ;
1310  int np = mg->get_np(l) ;
1311  const Grille3d* g = mg->get_grille3d(l) ;
1312  Tbl* tb = (mti->t)[l] ;
1313  tb->set_etat_qcq() ;
1314  double* p_r = tb->t ;
1315 
1316  switch(mg->get_type_r(l)) {
1317 
1318 
1319  case RARE : {
1320  for (int k=0 ; k<np ; k++) {
1321  for (int j=0 ; j<nt ; j++) {
1322  for (int i=0 ; i<nr ; i++) {
1323  double ww = alpha(l,k,j,0) ;
1324  *p_r = d2alphadt2(l,k,j,0) / (ww*ww * (g->x)[i]) ;
1325  p_r++ ;
1326  } // Fin de boucle sur r
1327  } // Fin de boucle sur theta
1328  } // Fin de boucle sur phi
1329  break ;
1330  }
1331  case FIN:{
1332  cout << "map_star_fait_drdt: Shells not implemented yet..." << endl;
1333  abort() ;
1334  break ;
1335  }
1336 
1337  case UNSURR: {
1338  cout << "map_star_fait_drdt: Compactified domain not allowed !" << endl;
1339  abort() ;
1340  break ;
1341  }
1342 
1343  default: {
1344  cout << "map_star_fait_drdt: unknown type_r !" << endl ;
1345  abort() ;
1346  }
1347  } // Fin du switch
1348  } // Fin de boucle sur zone
1349 
1350  // Termine
1351 
1352  return mti ;
1353 }
1354 
1355 }
Lorene::Map_star
Definition: map_star.h:23
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::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::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::Map_star::get_alpha
const Valeur & get_alpha() const
Returns the reference on the Tbl alpha.
Definition: map_star.C:237
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