LORENE
map_et_fait_der.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_der.C,v 1.8 2016/12/05 16:17:57 j_novak Exp $
27  * $Log: map_et_fait_der.C,v $
28  * Revision 1.8 2016/12/05 16:17:57 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.7 2014/10/13 08:53:04 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.6 2014/10/06 15:13:13 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.5 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.4 2008/08/27 08:49:01 jl_cornou
42  * Added R_JACO02 case
43  *
44  * Revision 1.3 2003/10/15 10:40:26 e_gourgoulhon
45  * Added new Coord's drdt and stdrdp.
46  * Changed cast (const Map_et *) into static_cast<const Map_et*>.
47  *
48  * Revision 1.2 2002/10/16 14:36:41 j_novak
49  * Reorganization of #include instructions of standard C++, in order to
50  * use experimental version 3 of gcc.
51  *
52  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
53  * LORENE
54  *
55  * Revision 1.2 1999/12/20 17:27:37 eric
56  * Modif lapr_tp.
57  *
58  * Revision 1.1 1999/11/24 11:23:06 eric
59  * Initial revision
60  *
61  *
62  * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_fait_der.C,v 1.8 2016/12/05 16:17:57 j_novak Exp $
63  *
64  */
65 
66 // Includes
67 #include <cassert>
68 #include <cstdlib>
69 #include <cmath>
70 
71 #include "map.h"
72 
73 
75 // Derivees d'ordre 1 du changement de variables //
77 
78 /*
79  ************************************************************************
80  * 1/(dR/dx) ( -1/(dU/dx) ds la ZEC )
81  ************************************************************************
82  */
83 
84 namespace Lorene {
85 Mtbl* map_et_fait_dxdr(const Map* cvi) {
86 
87  // recup du changement de variable
88  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
89  const Mg3d* mg = cv->get_mg() ;
90  int nz = mg->get_nzone() ;
91 
92  // Le resultat
93  Mtbl* mti = new Mtbl(mg) ;
94  mti->set_etat_qcq() ;
95 
96  // Pour le confort
97  const double* alpha = cv->alpha ;
98  const Valeur& ff = cv->ff ;
99  const Valeur& gg = cv->gg ;
100 
101  for (int l=0 ; l<nz ; l++) {
102 
103  int nr = mg->get_nr(l);
104  int nt = mg->get_nt(l) ;
105  int np = mg->get_np(l) ;
106 
107  const Tbl& da = *((cv->daa)[l]) ;
108  const Tbl& db = *((cv->dbb)[l]) ;
109 
110  Tbl* tb = (mti->t)[l] ;
111  tb->set_etat_qcq() ;
112  double* p_r = tb->t ;
113 
114  switch(mg->get_type_r(l)) {
115 
116  case FIN: case RARE: {
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 = 1. /
121  ( alpha[l] * ( 1. + da(i) * ff(l, k, j, 0)
122  + db(i) * gg(l, k, j, 0) )
123  ) ;
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. /
136  ( alpha[l] * ( 1. + da(i) * ff(l, k, j, 0) )
137  ) ;
138  p_r++ ;
139  } // Fin de boucle sur r
140  } // Fin de boucle sur theta
141  } // Fin de boucle sur phi
142  break ;
143  }
144 
145  default: {
146  cout << "map_et_fait_dxdr: unknown type_r !" << endl ;
147  abort() ;
148  }
149  } // Fin du switch
150  } // Fin de boucle sur zone
151 
152  // Termine
153  return mti ;
154 }
155 
156 /*
157  *******************************************************************************
158  * dR/dtheta ( dans la ZEC: - dU/dth )
159  *******************************************************************************
160  */
161 
162 Mtbl* map_et_fait_drdt(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  // Pour le confort
174  const double* alpha = cv->alpha ;
175  const Valeur& ff = cv->ff ;
176  const Valeur& gg = cv->gg ;
177  const Valeur& dffdt = ff.dsdt() ;
178  const Valeur& dggdt = gg.dsdt() ;
179 
180 
181  for (int l=0 ; l<nz ; l++) {
182 
183  int nr = mg->get_nr(l);
184  int nt = mg->get_nt(l) ;
185  int np = mg->get_np(l) ;
186 
187  const Tbl& aa = *((cv->aa)[l]) ;
188  const Tbl& bb = *((cv->bb)[l]) ;
189 
190  Tbl* tb = (mti->t)[l] ;
191  tb->set_etat_qcq() ;
192  double* p_r = tb->t ;
193 
194  switch(mg->get_type_r(l)) {
195 
196 
197  case RARE : case FIN: {
198  for (int k=0 ; k<np ; k++) {
199  for (int j=0 ; j<nt ; j++) {
200  for (int i=0 ; i<nr ; i++) {
201  *p_r = alpha[l] * ( aa(i) * dffdt(l, k, j, 0)
202  + bb(i) * dggdt(l, k, j, 0) ) ;
203  p_r++ ;
204  } // Fin de boucle sur r
205  } // Fin de boucle sur theta
206  } // Fin de boucle sur phi
207  break ;
208  }
209 
210  case UNSURR: {
211 
212  for (int k=0 ; k<np ; k++) {
213  for (int j=0 ; j<nt ; j++) {
214  for (int i=0 ; i<nr ; i++) {
215  *p_r = - aa(i) * dffdt(l, k, j, 0) ;
216  p_r++ ;
217  } // Fin de boucle sur r
218  } // Fin de boucle sur theta
219  } // Fin de boucle sur phi
220  break ;
221  }
222 
223  default: {
224  cout << "map_et_fait_drdt: unknown type_r !" << endl ;
225  abort() ;
226  }
227  } // Fin du switch
228  } // Fin de boucle sur zone
229 
230  // Termine
231  return mti ;
232 }
233 
234 
235 /*
236  *******************************************************************************
237  * 1/sin(theta) dR/dphi ( dans la ZEC: - 1/sin(th) dU/dth )
238  *******************************************************************************
239  */
240 
241 Mtbl* map_et_fait_stdrdp(const Map* cvi) {
242 
243  // recup du changement de variable
244  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
245  const Mg3d* mg = cv->get_mg() ;
246  int nz = mg->get_nzone() ;
247 
248  // Le resultat
249  Mtbl* mti = new Mtbl(mg) ;
250  mti->set_etat_qcq() ;
251 
252  // Pour le confort
253  const double* alpha = cv->alpha ;
254  const Valeur& ff = cv->ff ;
255  const Valeur& gg = cv->gg ;
256  const Valeur& stdffdp = ff.stdsdp() ;
257  const Valeur& stdggdp = gg.stdsdp() ;
258 
259 
260  for (int l=0 ; l<nz ; l++) {
261 
262  int nr = mg->get_nr(l);
263  int nt = mg->get_nt(l) ;
264  int np = mg->get_np(l) ;
265 
266  const Tbl& aa = *((cv->aa)[l]) ;
267  const Tbl& bb = *((cv->bb)[l]) ;
268 
269  Tbl* tb = (mti->t)[l] ;
270  tb->set_etat_qcq() ;
271  double* p_r = tb->t ;
272 
273  switch(mg->get_type_r(l)) {
274 
275 
276  case RARE : case FIN: {
277  for (int k=0 ; k<np ; k++) {
278  for (int j=0 ; j<nt ; j++) {
279  for (int i=0 ; i<nr ; i++) {
280  *p_r = alpha[l] * ( aa(i) * stdffdp(l, k, j, 0)
281  + bb(i) * stdggdp(l, k, j, 0) ) ;
282  p_r++ ;
283  } // Fin de boucle sur r
284  } // Fin de boucle sur theta
285  } // Fin de boucle sur phi
286  break ;
287  }
288 
289  case UNSURR: {
290 
291  for (int k=0 ; k<np ; k++) {
292  for (int j=0 ; j<nt ; j++) {
293  for (int i=0 ; i<nr ; i++) {
294  *p_r = - aa(i) * stdffdp(l, k, j, 0) ;
295  p_r++ ;
296  } // Fin de boucle sur r
297  } // Fin de boucle sur theta
298  } // Fin de boucle sur phi
299  break ;
300  }
301 
302  default: {
303  cout << "map_et_fait_stdrdp: unknown type_r !" << endl ;
304  abort() ;
305  }
306  } // Fin du switch
307  } // Fin de boucle sur zone
308 
309  // Termine
310  return mti ;
311 }
312 
313 
314 /*
315  *******************************************************************************
316  * 1/R dR/dtheta ( dans la ZEC: - 1/U dU/dth )
317  *******************************************************************************
318  */
319 
320 Mtbl* map_et_fait_srdrdt(const Map* cvi) {
321 
322  // recup du changement de variable
323  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
324  const Mg3d* mg = cv->get_mg() ;
325  int nz = mg->get_nzone() ;
326 
327  // Le resultat
328  Mtbl* mti = new Mtbl(mg) ;
329  mti->set_etat_qcq() ;
330 
331  // Pour le confort
332  const double* alpha = cv->alpha ;
333  const double* beta = cv->beta ;
334  const Valeur& ff = cv->ff ;
335  const Valeur& gg = cv->gg ;
336  const Valeur& dffdt = ff.dsdt() ;
337  const Valeur& dggdt = gg.dsdt() ;
338 
339 
340  for (int l=0 ; l<nz ; l++) {
341 
342  const Grille3d* g = mg->get_grille3d(l) ;
343 
344  int nr = mg->get_nr(l);
345  int nt = mg->get_nt(l) ;
346  int np = mg->get_np(l) ;
347 
348  const Tbl& aa = *((cv->aa)[l]) ;
349  const Tbl& bb = *((cv->bb)[l]) ;
350 
351  Tbl* tb = (mti->t)[l] ;
352  tb->set_etat_qcq() ;
353  double* p_r = tb->t ;
354 
355  switch(mg->get_type_r(l)) {
356 
357  case RARE: {
358 
359  const Tbl& asx = cv->aasx ;
360  const Tbl& bsx = cv->bbsx ;
361 
362  for (int k=0 ; k<np ; k++) {
363  for (int j=0 ; j<nt ; j++) {
364  for (int i=0 ; i<nr ; i++) {
365  *p_r = ( asx(i) * dffdt(l, k, j, 0)
366  + bsx(i) * dggdt(l, k, j, 0) ) /
367  ( 1. + asx(i) * ff(l, k, j, 0)
368  + bsx(i) * gg(l, k, j, 0) ) ;
369  p_r++ ;
370  } // Fin de boucle sur r
371  } // Fin de boucle sur theta
372  } // Fin de boucle sur phi
373  break ;
374  }
375 
376  case FIN: {
377  for (int k=0 ; k<np ; k++) {
378  for (int j=0 ; j<nt ; j++) {
379  for (int i=0 ; i<nr ; i++) {
380  *p_r = alpha[l] * ( aa(i) * dffdt(l, k, j, 0)
381  + bb(i) * dggdt(l, k, j, 0) )
382  / ( alpha[l] * (
383  (g->x)[i] + aa(i) * ff(l, k, j, 0)
384  + bb(i) * gg(l, k, j, 0) )
385  + beta[l] ) ;
386  p_r++ ;
387  } // Fin de boucle sur r
388  } // Fin de boucle sur theta
389  } // Fin de boucle sur phi
390  break ;
391  }
392 
393  case UNSURR: {
394 
395  const Tbl& asxm1 = cv->zaasx ;
396 
397  for (int k=0 ; k<np ; k++) {
398  for (int j=0 ; j<nt ; j++) {
399  for (int i=0 ; i<nr ; i++) {
400  *p_r = - asxm1(i) * dffdt(l, k, j, 0)
401  / ( 1. + asxm1(i) * ff(l, k, j, 0) ) ;
402  p_r++ ;
403  } // Fin de boucle sur r
404  } // Fin de boucle sur theta
405  } // Fin de boucle sur phi
406  break ;
407  }
408 
409  default: {
410  cout << "map_et_fait_srdrdt: unknown type_r !" << endl ;
411  abort() ;
412  }
413  } // Fin du switch
414  } // Fin de boucle sur zone
415 
416  // Termine
417  return mti ;
418 }
419 
420 /*
421  *****************************************************************************
422  * 1/(R sin(theta)) dR/dphi ( ds la ZEC: - 1/(U sin(th)) dU/dphi )
423  *****************************************************************************
424  */
425 
426 Mtbl* map_et_fait_srstdrdp(const Map* cvi) {
427 
428  // recup du changement de variable
429  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
430  const Mg3d* mg = cv->get_mg() ;
431  int nz = mg->get_nzone() ;
432 
433  // Le resultat
434  Mtbl* mti = new Mtbl(mg) ;
435  mti->set_etat_qcq() ;
436 
437  // Pour le confort
438  const double* alpha = cv->alpha ;
439  const double* beta = cv->beta ;
440  const Valeur& ff = cv->ff ;
441  const Valeur& gg = cv->gg ;
442  const Valeur& stdffdp = ff.stdsdp() ;
443  const Valeur& stdggdp = gg.stdsdp() ;
444 
445 
446  for (int l=0 ; l<nz ; l++) {
447 
448  const Grille3d* g = mg->get_grille3d(l) ;
449 
450  int nr = mg->get_nr(l);
451  int nt = mg->get_nt(l) ;
452  int np = mg->get_np(l) ;
453 
454  const Tbl& aa = *((cv->aa)[l]) ;
455  const Tbl& bb = *((cv->bb)[l]) ;
456 
457  Tbl* tb = (mti->t)[l] ;
458  tb->set_etat_qcq() ;
459  double* p_r = tb->t ;
460 
461  switch(mg->get_type_r(l)) {
462 
463  case RARE: {
464 
465  const Tbl& asx = cv->aasx ;
466  const Tbl& bsx = cv->bbsx ;
467 
468  for (int k=0 ; k<np ; k++) {
469  for (int j=0 ; j<nt ; j++) {
470  for (int i=0 ; i<nr ; i++) {
471  *p_r = ( asx(i) * stdffdp(l, k, j, 0)
472  + bsx(i) * stdggdp(l, k, j, 0) ) /
473  ( 1. + asx(i) * ff(l, k, j, 0)
474  + bsx(i) * gg(l, k, j, 0) ) ;
475  p_r++ ;
476  } // Fin de boucle sur r
477  } // Fin de boucle sur theta
478  } // Fin de boucle sur phi
479  break ;
480  }
481 
482  case FIN: {
483  for (int k=0 ; k<np ; k++) {
484  for (int j=0 ; j<nt ; j++) {
485  for (int i=0 ; i<nr ; i++) {
486  *p_r = alpha[l] * ( aa(i) * stdffdp(l, k, j, 0)
487  + bb(i) * stdggdp(l, k, j, 0) )
488  / ( alpha[l] * (
489  (g->x)[i] + aa(i) * ff(l, k, j, 0)
490  + bb(i) * gg(l, k, j, 0) )
491  + beta[l] ) ;
492  p_r++ ;
493  } // Fin de boucle sur r
494  } // Fin de boucle sur theta
495  } // Fin de boucle sur phi
496  break ;
497  }
498 
499  case UNSURR: {
500 
501  const Tbl& asxm1 = cv->zaasx ;
502 
503  for (int k=0 ; k<np ; k++) {
504  for (int j=0 ; j<nt ; j++) {
505  for (int i=0 ; i<nr ; i++) {
506  *p_r = - asxm1(i) * stdffdp(l, k, j, 0)
507  / ( 1. + asxm1(i) * ff(l, k, j, 0) ) ;
508  p_r++ ;
509  } // Fin de boucle sur r
510  } // Fin de boucle sur theta
511  } // Fin de boucle sur phi
512  break ;
513  }
514 
515  default: {
516  cout << "map_et_fait_srstdrdp: unknown type_r !" << endl ;
517  abort() ;
518  }
519  } // Fin du switch
520  } // Fin de boucle sur zone
521 
522  return mti ;
523 }
524 
525 /*
526  *******************************************************************************
527  * 1/R^2 dR/dtheta ( dans la ZEC: - 1/U^2 dU/dth )
528  *******************************************************************************
529  */
530 
531 Mtbl* map_et_fait_sr2drdt(const Map* cvi) {
532 
533  // recup du changement de variable
534  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
535  const Mg3d* mg = cv->get_mg() ;
536  int nz = mg->get_nzone() ;
537 
538  // Le resultat
539  Mtbl* mti = new Mtbl(mg) ;
540  mti->set_etat_qcq() ;
541 
542  // Pour le confort
543  const double* alpha = cv->alpha ;
544  const double* beta = cv->beta ;
545  const Valeur& ff = cv->ff ;
546  const Valeur& gg = cv->gg ;
547  const Valeur& dffdt = ff.dsdt() ;
548  const Valeur& dggdt = gg.dsdt() ;
549 
550 
551  for (int l=0 ; l<nz ; l++) {
552 
553  const Grille3d* g = mg->get_grille3d(l) ;
554 
555  int nr = mg->get_nr(l);
556  int nt = mg->get_nt(l) ;
557  int np = mg->get_np(l) ;
558 
559  const Tbl& aa = *((cv->aa)[l]) ;
560  const Tbl& bb = *((cv->bb)[l]) ;
561 
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: {
569 
570  const Tbl& asx = cv->aasx ;
571  const Tbl& bsx = cv->bbsx ;
572  const Tbl& asx2 = cv->aasx2 ;
573  const Tbl& bsx2 = cv->bbsx2 ;
574 
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 
579  double ww = 1. + asx(i) * ff(l, k, j, 0)
580  + bsx(i) * gg(l, k, j, 0) ;
581 
582  *p_r = ( asx2(i) * dffdt(l, k, j, 0)
583  + bsx2(i) * dggdt(l, k, j, 0) ) /
584  (alpha[l] * ww*ww) ;
585  p_r++ ;
586  } // Fin de boucle sur r
587  } // Fin de boucle sur theta
588  } // Fin de boucle sur phi
589  break ;
590  }
591 
592 
593  case FIN: {
594  for (int k=0 ; k<np ; k++) {
595  for (int j=0 ; j<nt ; j++) {
596  for (int i=0 ; i<nr ; i++) {
597 
598  double ww = alpha[l] * (
599  (g->x)[i] + aa(i) * ff(l, k, j, 0)
600  + bb(i) * gg(l, k, j, 0) )
601  + beta[l] ;
602 
603  *p_r = alpha[l] * ( aa(i) * dffdt(l, k, j, 0)
604  + bb(i) * dggdt(l, k, j, 0) )
605  / ( ww*ww ) ;
606  p_r++ ;
607  } // Fin de boucle sur r
608  } // Fin de boucle sur theta
609  } // Fin de boucle sur phi
610  break ;
611  }
612 
613  case UNSURR: {
614 
615  const Tbl& asxm1 = cv->zaasx ;
616  const Tbl& asxm1car = cv->zaasx2 ;
617 
618  for (int k=0 ; k<np ; k++) {
619  for (int j=0 ; j<nt ; j++) {
620  for (int i=0 ; i<nr ; i++) {
621 
622  double ww = 1. + asxm1(i) * ff(l, k, j, 0) ;
623 
624  *p_r = - asxm1car(i) * dffdt(l, k, j, 0)
625  / ( alpha[l] * ww*ww ) ;
626  p_r++ ;
627  } // Fin de boucle sur r
628  } // Fin de boucle sur theta
629  } // Fin de boucle sur phi
630  break ;
631  }
632 
633  default: {
634  cout << "map_et_fait_sr2drdt: unknown type_r !" << endl ;
635  abort() ;
636  }
637  } // Fin du switch
638  } // Fin de boucle sur zone
639 
640  // Termine
641  return mti ;
642 }
643 
644 /*
645  *******************************************************************************
646  * 1/(R^2 sin(theta)) dR/dphi ( ds la ZEC: - 1/(U^2 sin(th)) dU/dphi )
647  *******************************************************************************
648  */
649 
650 Mtbl* map_et_fait_sr2stdrdp(const Map* cvi) {
651 
652  // recup du changement de variable
653  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
654  const Mg3d* mg = cv->get_mg() ;
655  int nz = mg->get_nzone() ;
656 
657  // Le resultat
658  Mtbl* mti = new Mtbl(mg) ;
659  mti->set_etat_qcq() ;
660 
661  // Pour le confort
662  const double* alpha = cv->alpha ;
663  const double* beta = cv->beta ;
664  const Valeur& ff = cv->ff ;
665  const Valeur& gg = cv->gg ;
666  const Valeur& stdffdp = ff.stdsdp() ;
667  const Valeur& stdggdp = gg.stdsdp() ;
668 
669 
670  for (int l=0 ; l<nz ; l++) {
671 
672  const Grille3d* g = mg->get_grille3d(l) ;
673 
674  int nr = mg->get_nr(l);
675  int nt = mg->get_nt(l) ;
676  int np = mg->get_np(l) ;
677 
678  const Tbl& aa = *((cv->aa)[l]) ;
679  const Tbl& bb = *((cv->bb)[l]) ;
680 
681  Tbl* tb = (mti->t)[l] ;
682  tb->set_etat_qcq() ;
683  double* p_r = tb->t ;
684 
685  switch(mg->get_type_r(l)) {
686 
687  case RARE: {
688 
689  const Tbl& asx = cv->aasx ;
690  const Tbl& bsx = cv->bbsx ;
691  const Tbl& asx2 = cv->aasx2 ;
692  const Tbl& bsx2 = cv->bbsx2 ;
693 
694  for (int k=0 ; k<np ; k++) {
695  for (int j=0 ; j<nt ; j++) {
696  for (int i=0 ; i<nr ; i++) {
697 
698  double ww = 1. + asx(i) * ff(l, k, j, 0)
699  + bsx(i) * gg(l, k, j, 0) ;
700 
701  *p_r = ( asx2(i) * stdffdp(l, k, j, 0)
702  + bsx2(i) * stdggdp(l, k, j, 0) ) /
703  (alpha[l] * ww*ww) ;
704  p_r++ ;
705  } // Fin de boucle sur r
706  } // Fin de boucle sur theta
707  } // Fin de boucle sur phi
708  break ;
709  }
710 
711  case FIN: {
712  for (int k=0 ; k<np ; k++) {
713  for (int j=0 ; j<nt ; j++) {
714  for (int i=0 ; i<nr ; i++) {
715 
716  double ww = alpha[l] * (
717  (g->x)[i] + aa(i) * ff(l, k, j, 0)
718  + bb(i) * gg(l, k, j, 0) )
719  + beta[l] ;
720 
721  *p_r = alpha[l] * ( aa(i) * stdffdp(l, k, j, 0)
722  + bb(i) * stdggdp(l, k, j, 0) )
723  / ( ww*ww ) ;
724  p_r++ ;
725  } // Fin de boucle sur r
726  } // Fin de boucle sur theta
727  } // Fin de boucle sur phi
728  break ;
729  }
730 
731  case UNSURR: {
732 
733  const Tbl& asxm1 = cv->zaasx ;
734  const Tbl& asxm1car = cv->zaasx2 ;
735 
736  for (int k=0 ; k<np ; k++) {
737  for (int j=0 ; j<nt ; j++) {
738  for (int i=0 ; i<nr ; i++) {
739 
740  double ww = 1. + asxm1(i) * ff(l, k, j, 0) ;
741 
742  *p_r = - asxm1car(i) * stdffdp(l, k, j, 0)
743  / ( alpha[l] * ww*ww ) ;
744  p_r++ ;
745  } // Fin de boucle sur r
746  } // Fin de boucle sur theta
747  } // Fin de boucle sur phi
748  break ;
749  }
750 
751  default: {
752  cout << "map_et_fait_sr2stdrdp: unknown type_r !" << endl ;
753  abort() ;
754  }
755  } // Fin du switch
756  } // Fin de boucle sur zone
757 
758  // Termine
759  return mti ;
760 }
761 
762 /*
763  ******************************************************************************
764  * 1/(dR/dx) R/x ds. le noyau
765  * 1/(dR/dx) R/(x + beta_l/alpha_l) ds. les coquilles
766  * 1/(dU/dx) U/(x-1) ds. la ZEC
767  ******************************************************************************
768  */
769 
770 Mtbl* map_et_fait_rsxdxdr(const Map* cvi) {
771 
772  // recup du changement de variable
773  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
774  const Mg3d* mg = cv->get_mg() ;
775  int nz = mg->get_nzone() ;
776 
777  // Le resultat
778  Mtbl* mti = new Mtbl(mg) ;
779  mti->set_etat_qcq() ;
780 
781  // Pour le confort
782  const double* alpha = cv->alpha ;
783  const double* beta = cv->beta ;
784  const Valeur& ff = cv->ff ;
785  const Valeur& gg = cv->gg ;
786 
787  for (int l=0 ; l<nz ; l++) {
788 
789  const Grille3d* g = mg->get_grille3d(l) ;
790 
791  int nr = mg->get_nr(l);
792  int nt = mg->get_nt(l) ;
793  int np = mg->get_np(l) ;
794 
795  const Tbl& aa = *((cv->aa)[l]) ;
796  const Tbl& bb = *((cv->bb)[l]) ;
797  const Tbl& da = *((cv->daa)[l]) ;
798  const Tbl& db = *((cv->dbb)[l]) ;
799 
800  Tbl* tb = (mti->t)[l] ;
801  tb->set_etat_qcq() ;
802  double* p_r = tb->t ;
803 
804  switch(mg->get_type_r(l)) {
805 
806  case RARE: {
807 
808  const Tbl& asx = cv->aasx ;
809  const Tbl& bsx = cv->bbsx ;
810 
811  for (int k=0 ; k<np ; k++) {
812  for (int j=0 ; j<nt ; j++) {
813  for (int i=0 ; i<nr ; i++) {
814 
815  *p_r = ( 1. + asx(i) * ff(l, k, j, 0)
816  + bsx(i) * gg(l, k, j, 0) ) /
817  ( 1. + da(i) * ff(l, k, j, 0)
818  + db(i) * gg(l, k, j, 0) ) ;
819  p_r++ ;
820  } // Fin de boucle sur r
821  } // Fin de boucle sur theta
822  } // Fin de boucle sur phi
823  break ;
824  }
825 
826 
827  case FIN: {
828  for (int k=0 ; k<np ; k++) {
829  for (int j=0 ; j<nt ; j++) {
830  for (int i=0 ; i<nr ; i++) {
831 
832  *p_r = ( (g->x)[i] + aa(i) * ff(l, k, j, 0)
833  + bb(i) * gg(l, k, j, 0)
834  + beta[l]/alpha[l] ) /
835  ( ( 1. + da(i) * ff(l, k, j, 0)
836  + db(i) * gg(l, k, j, 0) ) *
837  ( (g->x)[i] + beta[l]/alpha[l] ) ) ;
838 
839  p_r++ ;
840  } // Fin de boucle sur r
841  } // Fin de boucle sur theta
842  } // Fin de boucle sur phi
843  break ;
844  }
845 
846  case UNSURR: {
847 
848  const Tbl& asxm1 = cv->zaasx ;
849 
850  for (int k=0 ; k<np ; k++) {
851  for (int j=0 ; j<nt ; j++) {
852  for (int i=0 ; i<nr ; i++) {
853 
854  *p_r = ( 1. + asxm1(i) * ff(l, k, j, 0) ) /
855  ( 1. + da(i) * ff(l, k, j, 0) ) ;
856 
857  p_r++ ;
858  } // Fin de boucle sur r
859  } // Fin de boucle sur theta
860  } // Fin de boucle sur phi
861  break ;
862  }
863 
864  default: {
865  cout << "map_et_fait_rsxdxdr: unknown type_r !" << endl ;
866  abort() ;
867  }
868  } // Fin du switch
869  } // Fin de boucle sur zone
870 
871  // Termine
872  return mti ;
873 }
874 
875 /*
876  ******************************************************************************
877  * [ R/ (alpha_l x + beta_l) ]^2 (dR/dx) /alpha_l ds. le noyau et les coquilles
878  * dU/dx /alpha_l ds. la ZEC
879  ******************************************************************************
880  */
881 
882 Mtbl* map_et_fait_rsx2drdx(const Map* cvi) {
883 
884  // recup du changement de variable
885  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
886  const Mg3d* mg = cv->get_mg() ;
887  int nz = mg->get_nzone() ;
888 
889  // Le resultat
890  Mtbl* mti = new Mtbl(mg) ;
891  mti->set_etat_qcq() ;
892 
893  // Pour le confort
894  const double* alpha = cv->alpha ;
895  const double* beta = cv->beta ;
896  const Valeur& ff = cv->ff ;
897  const Valeur& gg = cv->gg ;
898 
899  for (int l=0 ; l<nz ; l++) {
900 
901  const Grille3d* g = mg->get_grille3d(l) ;
902 
903  int nr = mg->get_nr(l);
904  int nt = mg->get_nt(l) ;
905  int np = mg->get_np(l) ;
906 
907  const Tbl& aa = *((cv->aa)[l]) ;
908  const Tbl& bb = *((cv->bb)[l]) ;
909  const Tbl& da = *((cv->daa)[l]) ;
910  const Tbl& db = *((cv->dbb)[l]) ;
911 
912  Tbl* tb = (mti->t)[l] ;
913  tb->set_etat_qcq() ;
914  double* p_r = tb->t ;
915 
916  switch(mg->get_type_r(l)) {
917 
918  case RARE: {
919 
920  const Tbl& asx = cv->aasx ;
921  const Tbl& bsx = cv->bbsx ;
922 
923  for (int k=0 ; k<np ; k++) {
924  for (int j=0 ; j<nt ; j++) {
925  for (int i=0 ; i<nr ; i++) {
926 
927  double ww = 1. + asx(i) * ff(l, k, j, 0)
928  + bsx(i) * gg(l, k, j, 0) ;
929 
930  *p_r = ww * ww *
931  ( 1. + da(i) * ff(l, k, j, 0)
932  + db(i) * gg(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 
941  case FIN: {
942  for (int k=0 ; k<np ; k++) {
943  for (int j=0 ; j<nt ; j++) {
944  for (int i=0 ; i<nr ; i++) {
945 
946  double ww = ( (g->x)[i] + aa(i) * ff(l, k, j, 0)
947  + bb(i) * gg(l, k, j, 0)
948  + beta[l]/alpha[l] ) /
949  ( (g->x)[i] + beta[l]/alpha[l] ) ;
950 
951  *p_r = ww * ww *
952  ( 1. + da(i) * ff(l, k, j, 0)
953  + db(i) * gg(l, k, j, 0) ) ;
954 
955  p_r++ ;
956  } // Fin de boucle sur r
957  } // Fin de boucle sur theta
958  } // Fin de boucle sur phi
959  break ;
960  }
961 
962  case UNSURR: {
963 
964  for (int k=0 ; k<np ; k++) {
965  for (int j=0 ; j<nt ; j++) {
966  for (int i=0 ; i<nr ; i++) {
967 
968  *p_r = 1. + da(i) * ff(l, k, j, 0) ;
969 
970  p_r++ ;
971  } // Fin de boucle sur r
972  } // Fin de boucle sur theta
973  } // Fin de boucle sur phi
974  break ;
975  }
976 
977  default: {
978  cout << "map_et_fait_rsx2drdx: unknown type_r !" << endl ;
979  abort() ;
980  }
981  } // Fin du switch
982  } // Fin de boucle sur zone
983 
984  // Termine
985  return mti ;
986 }
987 
988 
989 
991 // Derivees d'ordre 2 du changement de variables //
993 
994 
995 /*
996  *******************************************************************************
997  * d^2R/dx^2 ( dans la ZEC: - d^2U/dx^2 )
998  *******************************************************************************
999  */
1000 
1001 Mtbl* map_et_fait_d2rdx2(const Map* cvi) {
1002 
1003  // recup du changement de variable
1004  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
1005  const Mg3d* mg = cv->get_mg() ;
1006  int nz = mg->get_nzone() ;
1007 
1008  // Le resultat
1009  Mtbl* mti = new Mtbl(mg) ;
1010  mti->set_etat_qcq() ;
1011 
1012  // Pour le confort
1013  const double* alpha = cv->alpha ;
1014  const Valeur& ff = cv->ff ;
1015  const Valeur& gg = cv->gg ;
1016 
1017  for (int l=0 ; l<nz ; l++) {
1018 
1019  int nr = mg->get_nr(l);
1020  int nt = mg->get_nt(l) ;
1021  int np = mg->get_np(l) ;
1022 
1023  const Tbl& dda = *((cv->ddaa)[l]) ;
1024  const Tbl& ddb = *((cv->ddbb)[l]) ;
1025 
1026  Tbl* tb = (mti->t)[l] ;
1027  tb->set_etat_qcq() ;
1028  double* p_r = tb->t ;
1029 
1030  switch(mg->get_type_r(l)) {
1031 
1032  case FIN: case RARE: {
1033  for (int k=0 ; k<np ; k++) {
1034  for (int j=0 ; j<nt ; j++) {
1035  for (int i=0 ; i<nr ; i++) {
1036 
1037  *p_r = alpha[l] * ( dda(i) * ff(l, k, j, 0)
1038  + ddb(i) * gg(l, k, j, 0) ) ;
1039  p_r++ ;
1040  } // Fin de boucle sur r
1041  } // Fin de boucle sur theta
1042  } // Fin de boucle sur phi
1043  break ;
1044  }
1045 
1046  case UNSURR: {
1047  for (int k=0 ; k<np ; k++) {
1048  for (int j=0 ; j<nt ; j++) {
1049  for (int i=0 ; i<nr ; i++) {
1050 
1051  *p_r = - alpha[l] * dda(i) * ff(l, k, j, 0) ;
1052  p_r++ ;
1053  } // Fin de boucle sur r
1054  } // Fin de boucle sur theta
1055  } // Fin de boucle sur phi
1056  break ;
1057  }
1058 
1059  default: {
1060  cout << "map_et_fait_d2rdx2: unknown type_r !" << endl ;
1061  abort() ;
1062  }
1063  } // Fin du switch
1064  } // Fin de boucle sur zone
1065 
1066  // Termine
1067  return mti ;
1068 }
1069 
1070 /*
1071  *****************************************************************************
1072  * 1/R^2 ( 1/sin(th) d/dth( sin(th) dR/dth ) + 1/sin(th)^2 d^2R/dphi^2 )
1073  *
1074  * [ dans la ZEC :
1075  * - 1/U^2 ( 1/sin(th) d/dth( sin(th) dU/dth ) + 1/sin(th)^2 d^2U/dphi^2 ) ]
1076  *****************************************************************************
1077  */
1078 
1079 Mtbl* map_et_fait_lapr_tp(const Map* cvi) {
1080 
1081  // recup du changement de variable
1082  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
1083  const Mg3d* mg = cv->get_mg() ;
1084  int nz = mg->get_nzone() ;
1085 
1086  // Le resultat
1087  Mtbl* mti = new Mtbl(mg) ;
1088  mti->set_etat_qcq() ;
1089 
1090  // Pour le confort
1091  const double* alpha = cv->alpha ;
1092  const double* beta = cv->beta ;
1093  const Valeur& ff = cv->ff ;
1094  const Valeur& gg = cv->gg ;
1095 
1096  Valeur ff_tmp = ff ;
1097  Valeur gg_tmp = gg ;
1098  ff_tmp.ylm() ;
1099  gg_tmp.ylm() ;
1100  const Valeur& lapff = ff_tmp.lapang() ;
1101  const Valeur& lapgg = gg_tmp.lapang() ;
1102 
1103 
1104  for (int l=0 ; l<nz ; l++) {
1105 
1106  const Grille3d* g = mg->get_grille3d(l) ;
1107 
1108  int nr = mg->get_nr(l);
1109  int nt = mg->get_nt(l) ;
1110  int np = mg->get_np(l) ;
1111 
1112  const Tbl& aa = *((cv->aa)[l]) ;
1113  const Tbl& bb = *((cv->bb)[l]) ;
1114 
1115  Tbl* tb = (mti->t)[l] ;
1116  tb->set_etat_qcq() ;
1117  double* p_r = tb->t ;
1118 
1119  switch(mg->get_type_r(l)) {
1120 
1121  case RARE: {
1122 
1123  const Tbl& asx = cv->aasx ;
1124  const Tbl& bsx = cv->bbsx ;
1125  const Tbl& asx2 = cv->aasx2 ;
1126  const Tbl& bsx2 = cv->bbsx2 ;
1127 
1128  for (int k=0 ; k<np ; k++) {
1129  for (int j=0 ; j<nt ; j++) {
1130  for (int i=0 ; i<nr ; i++) {
1131 
1132  double ww = 1. + asx(i) * ff(l, k, j, 0)
1133  + bsx(i) * gg(l, k, j, 0) ;
1134 
1135  *p_r = ( asx2(i) * lapff(l, k, j, 0)
1136  + bsx2(i) * lapgg(l, k, j, 0) ) /
1137  (alpha[l] * ww*ww) ;
1138  p_r++ ;
1139  } // Fin de boucle sur r
1140  } // Fin de boucle sur theta
1141  } // Fin de boucle sur phi
1142  break ;
1143  }
1144 
1145  case FIN: {
1146  for (int k=0 ; k<np ; k++) {
1147  for (int j=0 ; j<nt ; j++) {
1148  for (int i=0 ; i<nr ; i++) {
1149 
1150  double ww = alpha[l] * (
1151  (g->x)[i] + aa(i) * ff(l, k, j, 0)
1152  + bb(i) * gg(l, k, j, 0) )
1153  + beta[l] ;
1154 
1155  *p_r = alpha[l] * ( aa(i) * lapff(l, k, j, 0)
1156  + bb(i) * lapgg(l, k, j, 0) )
1157  / ( ww*ww ) ;
1158  p_r++ ;
1159  } // Fin de boucle sur r
1160  } // Fin de boucle sur theta
1161  } // Fin de boucle sur phi
1162  break ;
1163  }
1164 
1165  case UNSURR: {
1166 
1167  const Tbl& asxm1 = cv->zaasx ;
1168  const Tbl& asxm1car = cv->zaasx2 ;
1169 
1170  for (int k=0 ; k<np ; k++) {
1171  for (int j=0 ; j<nt ; j++) {
1172  for (int i=0 ; i<nr ; i++) {
1173 
1174  double ww = 1. + asxm1(i) * ff(l, k, j, 0) ;
1175 
1176  *p_r = - asxm1car(i) * lapff(l, k, j, 0)
1177  / ( alpha[l] * ww*ww ) ;
1178  p_r++ ;
1179  } // Fin de boucle sur r
1180  } // Fin de boucle sur theta
1181  } // Fin de boucle sur phi
1182  break ;
1183  }
1184 
1185  default: {
1186  cout << "map_et_fait_lapr_tp: unknown type_r !" << endl ;
1187  abort() ;
1188  }
1189  } // Fin du switch
1190  } // Fin de boucle sur zone
1191 
1192  // Termine
1193  return mti ;
1194 }
1195 
1196 /*
1197  *******************************************************************************
1198  * d^2R/dthdx ( dans la ZEC: - d^2U/dthdx )
1199  *******************************************************************************
1200  */
1201 
1202 Mtbl* map_et_fait_d2rdtdx(const Map* cvi) {
1203 
1204  // recup du changement de variable
1205  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
1206  const Mg3d* mg = cv->get_mg() ;
1207  int nz = mg->get_nzone() ;
1208 
1209  // Le resultat
1210  Mtbl* mti = new Mtbl(mg) ;
1211  mti->set_etat_qcq() ;
1212 
1213  // Pour le confort
1214  const double* alpha = cv->alpha ;
1215  const Valeur& ff = cv->ff ;
1216  const Valeur& gg = cv->gg ;
1217  const Valeur& dffdt = ff.dsdt() ;
1218  const Valeur& dggdt = gg.dsdt() ;
1219 
1220  for (int l=0 ; l<nz ; l++) {
1221 
1222  int nr = mg->get_nr(l);
1223  int nt = mg->get_nt(l) ;
1224  int np = mg->get_np(l) ;
1225 
1226  const Tbl& da = *((cv->daa)[l]) ;
1227  const Tbl& db = *((cv->dbb)[l]) ;
1228 
1229  Tbl* tb = (mti->t)[l] ;
1230  tb->set_etat_qcq() ;
1231  double* p_r = tb->t ;
1232 
1233  switch(mg->get_type_r(l)) {
1234 
1235  case FIN: case RARE: {
1236  for (int k=0 ; k<np ; k++) {
1237  for (int j=0 ; j<nt ; j++) {
1238  for (int i=0 ; i<nr ; i++) {
1239 
1240  *p_r = alpha[l] * ( da(i) * dffdt(l, k, j, 0)
1241  + db(i) * dggdt(l, k, j, 0) ) ;
1242  p_r++ ;
1243  } // Fin de boucle sur r
1244  } // Fin de boucle sur theta
1245  } // Fin de boucle sur phi
1246  break ;
1247  }
1248 
1249  case UNSURR: {
1250  for (int k=0 ; k<np ; k++) {
1251  for (int j=0 ; j<nt ; j++) {
1252  for (int i=0 ; i<nr ; i++) {
1253 
1254  *p_r = - alpha[l] * da(i) * dffdt(l, k, j, 0) ;
1255  p_r++ ;
1256  } // Fin de boucle sur r
1257  } // Fin de boucle sur theta
1258  } // Fin de boucle sur phi
1259  break ;
1260  }
1261 
1262  default: {
1263  cout << "map_et_fait_d2rdtdx: unknown type_r !" << endl ;
1264  abort() ;
1265  }
1266  } // Fin du switch
1267  } // Fin de boucle sur zone
1268 
1269  // Termine
1270  return mti ;
1271 }
1272 
1273 /*
1274  *******************************************************************************
1275  * 1/sin(th) d^2R/dphidx ( dans la ZEC: - 1/sin(th) d^2U/dphidx )
1276  *******************************************************************************
1277  */
1278 
1279 Mtbl* map_et_fait_sstd2rdpdx(const Map* cvi) {
1280 
1281  // recup du changement de variable
1282  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
1283  const Mg3d* mg = cv->get_mg() ;
1284  int nz = mg->get_nzone() ;
1285 
1286  // Le resultat
1287  Mtbl* mti = new Mtbl(mg) ;
1288  mti->set_etat_qcq() ;
1289 
1290  // Pour le confort
1291  const double* alpha = cv->alpha ;
1292  const Valeur& ff = cv->ff ;
1293  const Valeur& gg = cv->gg ;
1294  const Valeur& stdffdp = ff.stdsdp() ;
1295  const Valeur& stdggdp = gg.stdsdp() ;
1296 
1297  for (int l=0 ; l<nz ; l++) {
1298 
1299  int nr = mg->get_nr(l);
1300  int nt = mg->get_nt(l) ;
1301  int np = mg->get_np(l) ;
1302 
1303  const Tbl& da = *((cv->daa)[l]) ;
1304  const Tbl& db = *((cv->dbb)[l]) ;
1305 
1306  Tbl* tb = (mti->t)[l] ;
1307  tb->set_etat_qcq() ;
1308  double* p_r = tb->t ;
1309 
1310  switch(mg->get_type_r(l)) {
1311 
1312  case FIN: case RARE: {
1313  for (int k=0 ; k<np ; k++) {
1314  for (int j=0 ; j<nt ; j++) {
1315  for (int i=0 ; i<nr ; i++) {
1316 
1317  *p_r = alpha[l] * ( da(i) * stdffdp(l, k, j, 0)
1318  + db(i) * stdggdp(l, k, j, 0) ) ;
1319  p_r++ ;
1320  } // Fin de boucle sur r
1321  } // Fin de boucle sur theta
1322  } // Fin de boucle sur phi
1323  break ;
1324  }
1325 
1326  case UNSURR: {
1327  for (int k=0 ; k<np ; k++) {
1328  for (int j=0 ; j<nt ; j++) {
1329  for (int i=0 ; i<nr ; i++) {
1330 
1331  *p_r = - alpha[l] * da(i) * stdffdp(l, k, j, 0) ;
1332  p_r++ ;
1333  } // Fin de boucle sur r
1334  } // Fin de boucle sur theta
1335  } // Fin de boucle sur phi
1336  break ;
1337  }
1338 
1339  default: {
1340  cout << "map_et_fait_sstd2rdpdx: unknown type_r !" << endl ;
1341  abort() ;
1342  }
1343  } // Fin du switch
1344  } // Fin de boucle sur zone
1345 
1346  // Termine
1347  return mti ;
1348 }
1349 
1350 /*
1351  *******************************************************************************
1352  * 1/R^2 d^2R/dtheta^2 ( dans la ZEC: - 1/U^2 d^2U/dth^2 )
1353  *******************************************************************************
1354  */
1355 
1356 Mtbl* map_et_fait_sr2d2rdt2(const Map* cvi) {
1357 
1358  // recup du changement de variable
1359  const Map_et* cv = static_cast<const Map_et*>(cvi) ;
1360  const Mg3d* mg = cv->get_mg() ;
1361  int nz = mg->get_nzone() ;
1362 
1363  // Le resultat
1364  Mtbl* mti = new Mtbl(mg) ;
1365  mti->set_etat_qcq() ;
1366 
1367  // Pour le confort
1368  const double* alpha = cv->alpha ;
1369  const double* beta = cv->beta ;
1370  const Valeur& ff = cv->ff ;
1371  const Valeur& gg = cv->gg ;
1372  const Valeur& d2ffdt2 = ff.d2sdt2() ;
1373  const Valeur& d2ggdt2 = gg.d2sdt2() ;
1374 
1375 
1376  for (int l=0 ; l<nz ; l++) {
1377 
1378  const Grille3d* g = mg->get_grille3d(l) ;
1379 
1380  int nr = mg->get_nr(l);
1381  int nt = mg->get_nt(l) ;
1382  int np = mg->get_np(l) ;
1383 
1384  const Tbl& aa = *((cv->aa)[l]) ;
1385  const Tbl& bb = *((cv->bb)[l]) ;
1386 
1387  Tbl* tb = (mti->t)[l] ;
1388  tb->set_etat_qcq() ;
1389  double* p_r = tb->t ;
1390 
1391  switch(mg->get_type_r(l)) {
1392 
1393  case RARE: {
1394 
1395  const Tbl& asx = cv->aasx ;
1396  const Tbl& bsx = cv->bbsx ;
1397  const Tbl& asx2 = cv->aasx2 ;
1398  const Tbl& bsx2 = cv->bbsx2 ;
1399 
1400  for (int k=0 ; k<np ; k++) {
1401  for (int j=0 ; j<nt ; j++) {
1402  for (int i=0 ; i<nr ; i++) {
1403 
1404  double ww = 1. + asx(i) * ff(l, k, j, 0)
1405  + bsx(i) * gg(l, k, j, 0) ;
1406 
1407  *p_r = ( asx2(i) * d2ffdt2(l, k, j, 0)
1408  + bsx2(i) * d2ggdt2(l, k, j, 0) ) /
1409  (alpha[l] * ww*ww) ;
1410  p_r++ ;
1411  } // Fin de boucle sur r
1412  } // Fin de boucle sur theta
1413  } // Fin de boucle sur phi
1414  break ;
1415  }
1416 
1417  case FIN: {
1418  for (int k=0 ; k<np ; k++) {
1419  for (int j=0 ; j<nt ; j++) {
1420  for (int i=0 ; i<nr ; i++) {
1421 
1422  double ww = alpha[l] * (
1423  (g->x)[i] + aa(i) * ff(l, k, j, 0)
1424  + bb(i) * gg(l, k, j, 0) )
1425  + beta[l] ;
1426 
1427  *p_r = alpha[l] * ( aa(i) * d2ffdt2(l, k, j, 0)
1428  + bb(i) * d2ggdt2(l, k, j, 0) )
1429  / ( ww*ww ) ;
1430  p_r++ ;
1431  } // Fin de boucle sur r
1432  } // Fin de boucle sur theta
1433  } // Fin de boucle sur phi
1434  break ;
1435  }
1436 
1437  case UNSURR: {
1438 
1439  const Tbl& asxm1 = cv->zaasx ;
1440  const Tbl& asxm1car = cv->zaasx2 ;
1441 
1442  for (int k=0 ; k<np ; k++) {
1443  for (int j=0 ; j<nt ; j++) {
1444  for (int i=0 ; i<nr ; i++) {
1445 
1446  double ww = 1. + asxm1(i) * ff(l, k, j, 0) ;
1447 
1448  *p_r = - asxm1car(i) * d2ffdt2(l, k, j, 0)
1449  / ( alpha[l] * ww*ww ) ;
1450  p_r++ ;
1451  } // Fin de boucle sur r
1452  } // Fin de boucle sur theta
1453  } // Fin de boucle sur phi
1454  break ;
1455  }
1456 
1457  default: {
1458  cout << "map_et_fait_sr2d2rdt2: unknown type_r !" << endl ;
1459  abort() ;
1460  }
1461  } // Fin du switch
1462  } // Fin de boucle sur zone
1463 
1464  // Termine
1465  return mti ;
1466 }
1467 
1468 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::Mtbl::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the Mtbl is defined.
Definition: mtbl.h:274