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
84namespace Lorene {
85Mtbl* 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
162Mtbl* 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
241Mtbl* 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
320Mtbl* 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
426Mtbl* 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
531Mtbl* 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
650Mtbl* 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
770Mtbl* 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
882Mtbl* 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
1001Mtbl* 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
1079Mtbl* 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
1202Mtbl* 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
1279Mtbl* 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
1356Mtbl* 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::Grille3d
3D grid class in one domain.
Definition grilles.h:200
Lorene::Grille3d::get_nr
int get_nr() const
Returns nr.
Definition grilles.h:236
Lorene::Map_et
Radial mapping of rather general form.
Definition map.h:2770
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Lorene::Valeur::d2sdt2
const Valeur & d2sdt2() const
Returns of *this.
Definition valeur_d2sdt2.C:110
Lorene::Valeur::stdsdp
const Valeur & stdsdp() const
Returns of *this.
Definition valeur_stdsdp.C:63
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::dsdt
const Valeur & dsdt() const
Returns of *this.
Definition valeur_dsdt.C:115
Lorene::Valeur::lapang
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition valeur_lapang.C:75
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142