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