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