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 char 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 
135 namespace Lorene {
136 void 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 
220 void 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 
279 void 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 
363 void 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 
417 void 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 
475 void 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 
569 void 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 
633 void 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 
729 void 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 
797 void 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 
823 void 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 }
Lorene::Map_af::dsdradial
virtual void dsdradial(const Scalar &, Scalar &) const
Computes of a Scalar.
Definition: map_af_deriv.C:417
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Map_af::stdsdp
virtual void stdsdp(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
Definition: map_af_deriv.C:823
Lorene::Map_af::dsdr
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:279
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:797
Lorene::Valeur::dsdp
const Valeur & dsdp() const
Returns of *this.
Definition: valeur_dsdp.C:98
Lorene::Valeur::stdsdp
const Valeur & stdsdp() const
Returns of *this.
Definition: valeur_stdsdp.C:60
Lorene::Valeur::mult_xm1_zec
void mult_xm1_zec()
Applies the following operator to *this : \ Id (r sampling = RARE, FIN ) \ (r -sampling = UNSURR )
Definition: valeur_mult_xm1.C:76
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map_af::dsdxi
virtual void dsdxi(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:136
Lorene::Map_radial::xsr
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1549
Lorene::Valeur::ssint
const Valeur & ssint() const
Returns of *this.
Definition: valeur_ssint.C:112
Lorene::Map_af::dsdt
virtual void dsdt(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
Definition: map_af_deriv.C:797
Lorene::Map_af::srstdsdp
virtual void srstdsdp(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:633
Lorene::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
Lorene::Valeur::annule
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition: valeur.C:744
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Map_radial::dxdr
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1560
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Cmp::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Valeur::dsdt
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:112
Lorene::Valeur::dsdx
const Valeur & dsdx() const
Returns of *this.
Definition: valeur_dsdx.C:111
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Valeur::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Lorene::Valeur::sx
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition: valeur_sx.C:110
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
Lorene::Map_af::srdsdt
virtual void srdsdt(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:475
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Lorene::Map::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676