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
23char 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
81namespace Lorene {
82Mtbl* 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
159Mtbl* 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
238Mtbl* 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
317Mtbl* 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
423Mtbl* 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
528Mtbl* 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
647Mtbl* 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
767Mtbl* 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
879Mtbl* 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
998Mtbl* 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
1076Mtbl* 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
1199Mtbl* 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
1276Mtbl* 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
1353Mtbl* 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}
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
const Mg3d * get_mg() const
Gives the Mg3d on which the Mtbl is defined.
Definition mtbl.h:274
Lorene prototypes.
Definition app_hor.h:64