LORENE
map_af_deriv.C
1/*
2 * Computations of Cmp partial derivatives for a Map_af mapping
3 */
4
5/*
6 * Copyright (c) 1999-2004 Eric Gourgoulhon
7 * Copyright (c) 1999-2001 Philippe Grandclement
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28char map_af_deriv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_deriv.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $" ;
29
30/*
31 * $Id: map_af_deriv.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
32 * $Log: map_af_deriv.C,v $
33 * Revision 1.13 2014/10/13 08:53:02 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.12 2012/01/17 10:32:09 j_penner
37 * added a derivative with respect to the computational coordinate xi
38 *
39 * Revision 1.11 2008/09/21 13:57:21 j_novak
40 * Changed the test on the CED in the derivative.
41 *
42 * Revision 1.10 2004/12/17 13:35:02 m_forot
43 * Add the case T_LEG
44 *
45 * Revision 1.9 2004/06/22 08:49:58 p_grandclement
46 * Addition of everything needed for using the logarithmic mapping
47 *
48 * Revision 1.8 2004/01/27 09:33:48 j_novak
49 * New method Map_radial::div_r_zec
50 *
51 * Revision 1.7 2004/01/26 16:16:17 j_novak
52 * Methods of gradient for Scalar s. The input can have any dzpuis.
53 *
54 * Revision 1.6 2004/01/22 16:13:00 e_gourgoulhon
55 * Case dzpuis=2 treated in dsdr, srdsdt and srstdsdp (output: dzpuis =
56 * 3).
57 * Reorganization cases dzpuis = 0 and 4.
58 *
59 * Revision 1.5 2003/11/11 15:31:43 j_novak
60 * Added a #ifnedef... to prevent warnings.
61 *
62 * Revision 1.4 2003/10/22 13:08:05 j_novak
63 * Better handling of dzpuis flags
64 *
65 * Revision 1.3 2003/10/20 19:45:27 e_gourgoulhon
66 * Treatment of dzpuis in dsdt and stdsdp.
67 *
68 * Revision 1.2 2003/10/15 10:34:07 e_gourgoulhon
69 * Added new methods dsdt and stdsdp.
70 *
71 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
72 * LORENE
73 *
74 * Revision 2.12 2000/02/25 08:59:51 eric
75 * Remplacement de ci.get_dzpuis() == 0 par ci.check_dzpuis(0).
76 * Suppression de l'affectation des dzpuis Mtbl/Mtnl_cf a la fin car
77 * c'est fait par Cmp::set_dzpuis.
78 *
79 * Revision 2.11 2000/01/26 13:09:18 eric
80 * Reprototypage complet des routines de derivation:
81 * le resultat est desormais suppose alloue a l'exterieur de la routine
82 * et est passe en argument (Cmp& resu), si bien que le prototypage
83 * complet devient:
84 * void DERIV(const Cmp& ci, Cmp& resu) const
85 *
86 * Revision 2.10 1999/11/30 12:51:32 eric
87 * Valeur::base est desormais du type Base_val et non plus Base_val*.
88 *
89 * Revision 2.9 1999/11/26 14:23:55 eric
90 * Traitement dzpuis des Cmp.
91 *
92 * Revision 2.8 1999/11/26 10:58:02 eric
93 * Traitement dzpuis.
94 *
95 * Revision 2.7 1999/11/25 16:29:29 eric
96 * Reorganisation complete du calcul des derivees partielles.
97 *
98 * Revision 2.6 1999/10/27 15:45:23 eric
99 * Suppression du membre Cmp::c.
100 *
101 * Revision 2.5 1999/10/27 08:47:03 eric
102 * Introduction de Cmp::va a la place de *(Cmp::c).
103 *
104 * Revision 2.4 1999/10/22 08:16:21 eric
105 * const Map*.
106 *
107 * Revision 2.3 1999/10/14 14:27:17 eric
108 * Methodes const.
109 *
110 * Revision 2.2 1999/10/13 15:54:40 eric
111 * Mg3d* -> const Mg3d*
112 *
113 * Revision 2.1 1999/09/17 10:01:09 phil
114 * correction pour deriv_x et deriv_y
115 *
116 * Revision 2.0 1999/09/14 16:37:06 phil
117 * *** empty log message ***
118 *
119 *
120 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_deriv.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
121 *
122 */
123
124// Header Lorene
125#include "map.h"
126#include "cmp.h"
127#include "tensor.h"
128
129
130 //-----------------------//
131 // d/d\xi //
132 //-----------------------//
133
134
135namespace Lorene {
136void Map_af::dsdxi(const Cmp& ci, Cmp& resu) const {
137
138 assert (ci.get_etat() != ETATNONDEF) ;
139 assert (ci.get_mp()->get_mg() == mg) ;
140
141
142 if (ci.get_etat() == ETATZERO) {
143 resu.set_etat_zero() ;
144 }
145 else {
146 assert( ci.get_etat() == ETATQCQ ) ;
147
148
149 (ci.va).coef() ; // (ci.va).c_cf is up to date
150
151 int nz = mg->get_nzone() ;
152 int nzm1 = nz - 1 ;
153
154 switch( ci.get_dzpuis() ) {
155
156 case 0 : {
157 resu = (ci.va).dsdx() ; // dsdx == d/d\xi
158
159 if (mg->get_type_r(nzm1) == UNSURR) {
160 resu.set_dzpuis(2) ; // r^2 d/dr has been computed in the // SOMETHING IS WRONG HERE
161 // external domain
162 }
163 break ;
164 }
165
166 case 2 : {
167 assert(mg->get_type_r(nzm1) == UNSURR) ;
168
169 Valeur tmp((ci.va).dsdx() ) ;
170 tmp.annule(nzm1) ; // zero in the CED
171
172 // Special treatment in the CED
173 Valeur tmp_ced = - (ci.va).dsdx() ;
174 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
175 tmp_ced.mult_xm1_zec() ;
176 tmp_ced.set(nzm1) -= 2. * ci.va(nzm1) ;
177
178 // Recombination shells + CED :
179 resu = tmp + tmp_ced ;
180
181 resu.set_dzpuis(3) ;
182 break ;
183 }
184
185 case 4 : {
186 assert(mg->get_type_r(nzm1) == UNSURR) ;
187
188 Valeur tmp(ci.va.dsdx()) ;
189 Valeur tmp2 = tmp ;
190 tmp2.base = (ci.va).dsdx().base ;
191 tmp.annule(nzm1) ; // not in the CED
192 tmp2.annule(0, nz-2) ; // special treatment of the CED
193 tmp2.mult_xm1_zec() ;
194 tmp2 = tmp2 / xsr ;
195 tmp2.set(nzm1) -= 4*ci.va(nzm1) ;
196 tmp2.base = ci.va.base ; //Just for the CED
197 tmp2.mult_xm1_zec() ;
198
199 resu = tmp + tmp2 / xsr ; // do not know what this is, but for now I can get away with it, no CED
200 resu.set_dzpuis(4) ;
201 break ;
202 }
203
204 default : {
205 cerr << "Map_af::dsdxi: unexpected value of input dzpuis !\n"
206 << " ci.get_dzpuis() = " << ci.get_dzpuis() << endl ;
207 abort() ;
208 break ;
209 }
210
211 }
212
213
214 (resu.va).set_base( (ci.va).dsdx().get_base() ) ; // same basis as d/dxi
215
216 }
217
218}
219
220void Map_af::dsdxi(const Scalar& uu, Scalar& resu) const {
221
222 assert (uu.get_etat() != ETATNONDEF) ;
223 assert (uu.get_mp().get_mg() == mg) ;
224
225 Mtbl unity = r/r;
226
227 if (uu.get_etat() == ETATZERO) {
228 resu.set_etat_zero() ;
229 }
230 else {
231 assert( uu.get_etat() == ETATQCQ ) ;
232
233 const Valeur& uuva = uu.get_spectral_va() ;
234
235 uuva.coef() ; // (uu.va).c_cf is up to date
236
237 int nz = mg->get_nzone() ;
238 int nzm1 = nz - 1 ;
239
240 if ( uu.get_dzpuis() == 0 ) {
241 resu = uuva.dsdx() * unity ; // dxds == d/d\xi, unity is used to set the correct formated output
242
243 if (mg->get_type_r(nzm1) == UNSURR) {
244 resu.set_dzpuis(2) ; // r^2 d/dr has been computed in the
245 // external domain
246 }
247 }
248 else {
249 int dzp = uu.get_dzpuis() ;
250
251 resu = uuva.dsdx() ;
252 if (mg->get_type_r(nzm1) == UNSURR) {
253 resu.annule_domain(nzm1) ; // zero in the CED
254
255 // Special treatment in the CED
256 Valeur tmp_ced = - uuva.dsdx() ;
257 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
258 tmp_ced.mult_xm1_zec() ;
259 tmp_ced.set(nzm1) -= dzp * uuva(nzm1) ;
260
261 // Recombination shells + CED :
262 resu.set_spectral_va() += tmp_ced ;
263 }
264 resu.set_dzpuis(dzp+1) ;
265
266 }
267
268 resu.set_spectral_base( uuva.dsdx().get_base() ) ; // same basis as d/dxi
269
270 }
271
272}
273
274 //---------------------//
275 // d/dr //
276 //---------------------//
277
278
279void Map_af::dsdr(const Cmp& ci, Cmp& resu) const {
280
281 assert (ci.get_etat() != ETATNONDEF) ;
282 assert (ci.get_mp()->get_mg() == mg) ;
283
284
285 if (ci.get_etat() == ETATZERO) {
286 resu.set_etat_zero() ;
287 }
288 else {
289 assert( ci.get_etat() == ETATQCQ ) ;
290
291
292 (ci.va).coef() ; // (ci.va).c_cf is up to date
293
294 int nz = mg->get_nzone() ;
295 int nzm1 = nz - 1 ;
296
297 switch( ci.get_dzpuis() ) {
298
299 case 0 : {
300 resu = (ci.va).dsdx() * dxdr ; // dxdr = dxi/dR, - dxi/dU (ZEC)
301
302 if (mg->get_type_r(nzm1) == UNSURR) {
303 resu.set_dzpuis(2) ; // r^2 d/dr has been computed in the
304 // external domain
305 }
306 break ;
307 }
308
309 case 2 : {
310 assert(mg->get_type_r(nzm1) == UNSURR) ;
311
312 Valeur tmp((ci.va).dsdx() * dxdr) ;
313 tmp.annule(nzm1) ; // zero in the CED
314
315 // Special treatment in the CED
316 Valeur tmp_ced = - (ci.va).dsdx() ;
317 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
318 tmp_ced.mult_xm1_zec() ;
319 tmp_ced.set(nzm1) -= 2. * ci.va(nzm1) ;
320
321 // Recombination shells + CED :
322 resu = tmp + tmp_ced ;
323
324 resu.set_dzpuis(3) ;
325 break ;
326 }
327
328 case 4 : {
329 assert(mg->get_type_r(nzm1) == UNSURR) ;
330
331 Valeur tmp(ci.va.dsdx() * dxdr) ;
332 Valeur tmp2 = tmp ;
333 tmp2.base = (ci.va).dsdx().base ;
334 tmp.annule(nzm1) ; // not in the CED
335 tmp2.annule(0, nz-2) ; // special treatment of the CED
336 tmp2.mult_xm1_zec() ;
337 tmp2 = tmp2 / xsr ;
338 tmp2.set(nzm1) -= 4*ci.va(nzm1) ;
339 tmp2.base = ci.va.base ; //Just for the CED
340 tmp2.mult_xm1_zec() ;
341
342 resu = tmp + tmp2 / xsr ;
343 resu.set_dzpuis(4) ;
344 break ;
345 }
346
347 default : {
348 cerr << "Map_af::dsdr: unexpected value of input dzpuis !\n"
349 << " ci.get_dzpuis() = " << ci.get_dzpuis() << endl ;
350 abort() ;
351 break ;
352 }
353
354 }
355
356
357 (resu.va).set_base( (ci.va).dsdx().get_base() ) ; // same basis as d/dxi
358
359 }
360
361}
362
363void Map_af::dsdr(const Scalar& uu, Scalar& resu) const {
364
365 assert (uu.get_etat() != ETATNONDEF) ;
366 assert (uu.get_mp().get_mg() == mg) ;
367
368
369 if (uu.get_etat() == ETATZERO) {
370 resu.set_etat_zero() ;
371 }
372 else {
373 assert( uu.get_etat() == ETATQCQ ) ;
374
375 const Valeur& uuva = uu.get_spectral_va() ;
376
377 uuva.coef() ; // (uu.va).c_cf is up to date
378
379 int nz = mg->get_nzone() ;
380 int nzm1 = nz - 1 ;
381
382 if ( uu.get_dzpuis() == 0 ) {
383 resu = uuva.dsdx() * dxdr ; // dxdr = dxi/dR, - dxi/dU (ZEC)
384
385 if (mg->get_type_r(nzm1) == UNSURR) {
386 resu.set_dzpuis(2) ; // r^2 d/dr has been computed in the
387 // external domain
388 }
389 }
390 else {
391 int dzp = uu.get_dzpuis() ;
392
393 resu = uuva.dsdx() * dxdr ;
394 if (mg->get_type_r(nzm1) == UNSURR) {
395 resu.annule_domain(nzm1) ; // zero in the CED
396
397 // Special treatment in the CED
398 Valeur tmp_ced = - uuva.dsdx() ;
399 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
400 tmp_ced.mult_xm1_zec() ;
401 tmp_ced.set(nzm1) -= dzp * uuva(nzm1) ;
402
403 // Recombination shells + CED :
404 resu.set_spectral_va() += tmp_ced ;
405 }
406 resu.set_dzpuis(dzp+1) ;
407
408 }
409
410 resu.set_spectral_base( uuva.dsdx().get_base() ) ; // same basis as d/dxi
411
412 }
413
414}
415// VARIABLE RADIALE
416
417void Map_af::dsdradial(const Scalar& uu, Scalar& resu) const {
418
419 assert (uu.get_etat() != ETATNONDEF) ;
420 assert (uu.get_mp().get_mg() == mg) ;
421
422
423 if (uu.get_etat() == ETATZERO) {
424 resu.set_etat_zero() ;
425 }
426 else {
427 assert( uu.get_etat() == ETATQCQ ) ;
428
429 const Valeur& uuva = uu.get_spectral_va() ;
430
431 uuva.coef() ; // (uu.va).c_cf is up to date
432
433 int nz = mg->get_nzone() ;
434 int nzm1 = nz - 1 ;
435
436 if ( uu.get_dzpuis() == 0 ) {
437 resu = uuva.dsdx() * dxdr ; // dxdr = dxi/dR, - dxi/dU (ZEC)
438
439 if (mg->get_type_r(nzm1) == UNSURR) {
440 resu.set_dzpuis(2) ; // r^2 d/dr has been computed in the
441 // external domain
442 }
443 }
444 else {
445 assert(mg->get_type_r(nzm1) == UNSURR) ;
446
447 int dzp = uu.get_dzpuis() ;
448
449 resu = uuva.dsdx() * dxdr ;
450 resu.annule_domain(nzm1) ; // zero in the CED
451
452 // Special treatment in the CED
453 Valeur tmp_ced = - uuva.dsdx() ;
454 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
455 tmp_ced.mult_xm1_zec() ;
456 tmp_ced.set(nzm1) -= dzp * uuva(nzm1) ;
457
458 // Recombination shells + CED :
459 resu.set_spectral_va() += tmp_ced ;
460
461 resu.set_dzpuis(dzp+1) ;
462
463 }
464
465 resu.set_spectral_base( uuva.dsdx().get_base() ) ; // same basis as d/dxi
466
467 }
468
469}
470
471 //------------------------//
472 // 1/r d/dtheta //
473 //------------------------//
474
475void Map_af::srdsdt(const Cmp& ci, Cmp& resu) const {
476
477 assert (ci.get_etat() != ETATNONDEF) ;
478 assert (ci.get_mp()->get_mg() == mg) ;
479
480 if (ci.get_etat() == ETATZERO) {
481 resu.set_etat_zero() ;
482 }
483 else {
484
485 assert( ci.get_etat() == ETATQCQ ) ;
486
487 (ci.va).coef() ; // (ci.va).c_cf is up to date
488
489 Valeur tmp = ci.va ;
490
491 tmp = tmp.dsdt() ; // d/dtheta
492
493 int nz = mg->get_nzone() ;
494 int nzm1 = nz - 1 ;
495
496 switch( ci.get_dzpuis() ) {
497
498 case 0 : {
499 tmp = tmp.sx() ; // 1/xi, Id, 1/(xi-1)
500
501 Base_val sauve_base( tmp.get_base() ) ;
502
503 tmp = tmp * xsr ; // xi/R, 1/R, (xi-1)/U
504
505 tmp.set_base(sauve_base) ; // The above operation does not the basis
506 resu = tmp ;
507
508 if (mg->get_type_r(nz-1) == UNSURR) {
509 resu.set_dzpuis(2) ; // r d/dtheta has been computed in
510 // the external domain
511 }
512 break ;
513 }
514
515 case 2 : {
516 assert(mg->get_type_r(nzm1) == UNSURR) ;
517
518 // Special treatment in the CED
519 Valeur tmp_ced = tmp ; // d/dtheta
520 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
521
522 tmp.annule(nzm1) ;
523 tmp = tmp.sx() ; // 1/xi, Id
524 Base_val sauve_base( tmp.get_base() ) ;
525 tmp = tmp * xsr ; // xi/R, 1/R
526 tmp.set_base(sauve_base) ; // The above operation does not the basis
527
528 // Recombination shells + CED :
529 resu = tmp + tmp_ced ;
530
531 resu.set_dzpuis(3) ;
532 break ;
533 }
534
535 case 4 : {
536 assert(mg->get_type_r(nzm1) == UNSURR) ;
537
538 // Special treatment in the CED
539 Valeur tmp_ced = tmp ; // d/dtheta
540 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
541 tmp_ced.mult_xm1_zec() ;
542
543 tmp.annule(nzm1) ;
544 tmp = tmp.sx() ; // 1/xi, Id
545 Base_val sauve_base( tmp.get_base() ) ;
546 tmp = tmp * xsr ; // xi/R, 1/R
547
548 // Recombination shells + CED :
549 resu = tmp + tmp_ced / xsr ;
550
551 resu.va.set_base( sauve_base ) ;
552 resu.set_dzpuis(4) ;
553 break ;
554 }
555
556 default : {
557 cerr << "Map_af::srdsdt: unexpected value of input dzpuis !\n"
558 << " ci.get_dzpuis() = " << ci.get_dzpuis() << endl ;
559 abort() ;
560 break ;
561 }
562
563 }
564
565 }
566
567}
568
569void Map_af::srdsdt(const Scalar& uu, Scalar& resu) const {
570 assert (uu.get_etat() != ETATNONDEF) ;
571 assert (uu.get_mp().get_mg() == mg) ;
572
573 if (uu.get_etat() == ETATZERO) {
574 resu.set_etat_zero() ;
575 }
576 else {
577
578 assert( uu.get_etat() == ETATQCQ ) ;
579
580 const Valeur& uuva = uu.get_spectral_va() ;
581 uuva.coef() ; // (uu.va).c_cf is up to date
582
583 Valeur tmp = uuva.dsdt() ; // d/dtheta
584
585 int nz = mg->get_nzone() ;
586 int nzm1 = nz - 1 ;
587
588 if (uu.get_dzpuis() == 0) {
589 tmp = tmp.sx() ; // 1/xi, Id, 1/(xi-1)
590
591 Base_val sauve_base( tmp.get_base() ) ;
592
593 tmp = tmp * xsr ; // xi/R, 1/R, (xi-1)/U
594
595 tmp.set_base(sauve_base) ; // The above operation does not change the basis
596 resu = tmp ;
597
598 if (mg->get_type_r(nz-1) == UNSURR) {
599 resu.set_dzpuis(2) ; // r d/dtheta has been computed in
600 // the external domain
601 }
602 }
603
604 else {
605 assert(mg->get_type_r(nzm1) == UNSURR) ;
606
607 int dzp = uu.get_dzpuis() ;
608 // Special treatment in the CED
609 Valeur tmp_ced = tmp ; // d/dtheta
610 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
611
612 tmp.annule(nzm1) ;
613 tmp = tmp.sx() ; // 1/xi, Id
614 Base_val sauve_base( tmp.get_base() ) ;
615 tmp = tmp * xsr ; // xi/R, 1/R
616 tmp.set_base(sauve_base) ; // The above operation does not change the basis
617
618 // Recombination shells + CED :
619 resu = tmp + tmp_ced ;
620
621 resu.set_dzpuis(dzp+1) ;
622 }
623
624 }
625
626}
627
628
629 //------------------------------------//
630 // 1/(r sin(theta)) d/dphi //
631 //------------------------------------//
632
633void Map_af::srstdsdp(const Cmp& ci, Cmp& resu) const {
634
635 assert (ci.get_etat() != ETATNONDEF) ;
636 assert (ci.get_mp()->get_mg() == mg) ;
637
638 if (ci.get_etat() == ETATZERO) {
639 resu.set_etat_zero() ;
640 }
641 else {
642
643 assert( ci.get_etat() == ETATQCQ ) ;
644
645 (ci.va).coef() ; // (ci.va).c_cf is up to date
646
647 Valeur tmp = ci.va ;
648
649 tmp = tmp.dsdp() ; // d/dphi
650 tmp = tmp.ssint() ; // 1/sin(theta)
651
652 int nz = mg->get_nzone() ;
653 int nzm1 = nz - 1 ;
654
655 switch( ci.get_dzpuis() ) {
656
657 case 0 : {
658 tmp = tmp.sx() ; // 1/xi, Id, 1/(xi-1)
659
660 Base_val sauve_base( tmp.get_base() ) ;
661
662 tmp = tmp * xsr ; // xi/R, 1/R, (xi-1)/U
663
664 tmp.set_base(sauve_base) ; // The above operation does not the basis
665 resu = tmp ;
666
667 if (mg->get_type_r(nz-1) == UNSURR) {
668 resu.set_dzpuis(2) ; // r d/dtheta has been computed in
669 // the external domain
670 }
671 break ;
672 }
673
674 case 2 : {
675 assert(mg->get_type_r(nzm1) == UNSURR) ;
676
677 // Special treatment in the CED
678 Valeur tmp_ced = tmp ; // 1/sin(theta) d/dphi
679 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
680
681 tmp.annule(nzm1) ;
682 tmp = tmp.sx() ; // 1/xi, Id
683 Base_val sauve_base( tmp.get_base() ) ;
684 tmp = tmp * xsr ; // xi/R, 1/R
685 tmp.set_base(sauve_base) ; // The above operation does not the basis
686
687 // Recombination shells + CED :
688 resu = tmp + tmp_ced ;
689
690 resu.set_dzpuis(3) ;
691 break ;
692 }
693
694 case 4 : {
695 assert(mg->get_type_r(nzm1) == UNSURR) ;
696
697 // Special treatment in the CED
698 Valeur tmp_ced = tmp ; // 1/sin(theta) d/dphi
699 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
700 tmp_ced.mult_xm1_zec() ;
701
702 tmp.annule(nzm1) ;
703 tmp = tmp.sx() ; // 1/xi, Id
704 Base_val sauve_base( tmp.get_base() ) ;
705 tmp = tmp * xsr ; // xi/R, 1/R
706
707 // Recombination shells + CED :
708 resu = tmp + tmp_ced / xsr ;
709
710 resu.va.set_base( sauve_base ) ;
711 resu.set_dzpuis(4) ;
712 break ;
713 }
714
715 default : {
716 cerr << "Map_af::srstdsdp: unexpected value of input dzpuis !\n"
717 << " ci.get_dzpuis() = " << ci.get_dzpuis() << endl ;
718 abort() ;
719 break ;
720 }
721
722 }
723
724 }
725
726
727}
728
729void Map_af::srstdsdp(const Scalar& uu, Scalar& resu) const {
730
731 assert (uu.get_etat() != ETATNONDEF) ;
732 assert (uu.get_mp().get_mg() == mg) ;
733
734 if (uu.get_etat() == ETATZERO) {
735 resu.set_etat_zero() ;
736 }
737 else {
738
739 assert( uu.get_etat() == ETATQCQ ) ;
740
741 const Valeur& uuva = uu.get_spectral_va() ;
742 uuva.coef() ; // (uu.va).c_cf is up to date
743
744 Valeur tmp = uuva.dsdp() ;
745
746 tmp = tmp.ssint() ; // 1/sin(theta)
747
748 int nz = mg->get_nzone() ;
749 int nzm1 = nz - 1 ;
750
751 if (uu.get_dzpuis() == 0) {
752
753 tmp = tmp.sx() ; // 1/xi, Id, 1/(xi-1)
754
755 Base_val sauve_base( tmp.get_base() ) ;
756
757 tmp = tmp * xsr ; // xi/R, 1/R, (xi-1)/U
758
759 tmp.set_base(sauve_base) ; // The above operation does not change the basis
760 resu = tmp ;
761
762 if (mg->get_type_r(nz-1) == UNSURR) {
763 resu.set_dzpuis(2) ; // r d/dtheta has been computed in
764 // the external domain
765 }
766 }
767
768 else {
769 assert(mg->get_type_r(nzm1) == UNSURR) ;
770
771 int dzp = uu.get_dzpuis() ;
772
773 // Special treatment in the CED
774 Valeur tmp_ced = tmp ; // 1/sin(theta) d/dphi
775 tmp_ced.annule(0, nz-2) ; // only non zero in the CED
776
777 tmp.annule(nzm1) ;
778 tmp = tmp.sx() ; // 1/xi, Id
779 Base_val sauve_base( tmp.get_base() ) ;
780 tmp = tmp * xsr ; // xi/R, 1/R
781 tmp.set_base(sauve_base) ; // The above operation does not change the basis
782
783 // Recombination shells + CED :
784 resu = tmp + tmp_ced ;
785
786 resu.set_dzpuis(dzp+1) ;
787 }
788 }
789}
790
791
792 //------------------------//
793 // d/dtheta //
794 //------------------------//
795
796
797void Map_af::dsdt(const Scalar& ci, Scalar& resu) const {
798
799 assert (ci.get_etat() != ETATNONDEF) ;
800 assert (ci.get_mp().get_mg() == mg) ;
801
802 if (ci.get_etat() == ETATZERO) {
803 resu.set_etat_zero() ;
804 }
805 else {
806
807 assert( ci.get_etat() == ETATQCQ ) ;
808
809 resu = ci.get_spectral_va().dsdt() ; // d/dtheta
810
811 }
812
813 resu.set_dzpuis( ci.get_dzpuis() ) ; // dzpuis unchanged
814
815}
816
817
818 //-----------------------------------//
819 // 1/sin(theta) d/dphi //
820 //-----------------------------------//
821
822
823void Map_af::stdsdp(const Scalar& ci, Scalar& resu) const {
824
825 assert (ci.get_etat() != ETATNONDEF) ;
826 assert (ci.get_mp().get_mg() == mg) ;
827
828 if (ci.get_etat() == ETATZERO) {
829 resu.set_etat_zero() ;
830 }
831 else {
832
833 assert( ci.get_etat() == ETATQCQ ) ;
834
835 resu = ci.get_spectral_va().stdsdp() ; // 1/sin(theta) d/dphi
836
837 }
838
839 resu.set_dzpuis( ci.get_dzpuis() ) ; // dzpuis unchanged
840
841}
842
843
844
845
846
847
848
849
850}
Bases of the spectral expansions.
Definition base_val.h:322
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
virtual void srstdsdp(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
virtual void dsdradial(const Scalar &, Scalar &) const
Computes of a Scalar.
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
virtual void dsdt(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
virtual void stdsdp(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
virtual void srdsdt(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
virtual void dsdxi(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1560
Coord r
r coordinate centered on the grid
Definition map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Multi-domain array.
Definition mtbl.h:118
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
Values and coefficients of a (real-value) function.
Definition valeur.h:287
const Valeur & dsdp() const
Returns of *this.
Definition valeur_dsdp.C:98
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition valeur_sx.C:110
void mult_xm1_zec()
Applies the following operator to *this : \ Id (r sampling = RARE, FIN ) \ (r -sampling = UNSURR )
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition valeur.h:480
const Valeur & stdsdp() const
Returns of *this.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
const Valeur & dsdt() const
Returns of *this.
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
const Valeur & dsdx() const
Returns of *this.
const Valeur & ssint() const
Returns of *this.
void coef() const
Computes the coeffcients of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene prototypes.
Definition app_hor.h:64