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