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