LORENE
param_elliptic.C
1 /*
2  * Copyright (c) 2004 Philippe Grandclement
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 version 2
8  * as published by the Free Software Foundation.
9  *
10  * LORENE is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with LORENE; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18  *
19  */
20 
21 
22 
23 /*
24  * $Id: param_elliptic.C,v 1.23 2018/11/16 14:34:36 j_novak Exp $
25  * $Log: param_elliptic.C,v $
26  * Revision 1.23 2018/11/16 14:34:36 j_novak
27  * Changed minor points to avoid some compilation warnings.
28  *
29  * Revision 1.22 2016/12/05 16:18:14 j_novak
30  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31  *
32  * Revision 1.21 2014/10/13 08:53:37 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.20 2014/10/06 15:13:15 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.19 2007/05/06 10:48:14 p_grandclement
39  * Modification of a few operators for the vorton project
40  *
41  * Revision 1.18 2007/04/24 09:04:13 p_grandclement
42  * Addition of an operator for the vortons
43  *
44  * Revision 1.17 2005/05/12 09:49:44 j_novak
45  * Temptative treatment of the case where the source is null in the CED (putting
46  * dzpuis to 4). May be a bad idea...
47  *
48  * Revision 1.16 2005/02/15 15:43:17 j_novak
49  * First version of the block inversion for the vector Poisson equation (method 6).
50  *
51  * Revision 1.15 2004/12/23 16:30:16 j_novak
52  * New files and class for the solution of the rr component of the tensor Poisson
53  * equation.
54  *
55  * Revision 1.14 2004/08/24 09:14:49 p_grandclement
56  * Addition of some new operators, like Poisson in 2d... It now requieres the
57  * GSL library to work.
58  *
59  * Also, the way a variable change is stored by a Param_elliptic is changed and
60  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
61  * will requiere some modification. (It should concern only the ones about monopoles)
62  *
63  * Revision 1.13 2004/06/22 13:46:52 j_novak
64  * Deplacement du conte++ dans set_pois_vect_r
65  *
66  * Revision 1.12 2004/06/22 08:49:59 p_grandclement
67  * Addition of everything needed for using the logarithmic mapping
68  *
69  * Revision 1.11 2004/06/14 15:07:12 j_novak
70  * New methods for the construction of the elliptic operator appearing in
71  * the vector Poisson equation (acting on eta).
72  *
73  * Revision 1.10 2004/05/17 15:50:54 j_novak
74  * Removed unused nr variables
75  *
76  * Revision 1.9 2004/05/17 15:42:22 j_novak
77  * The method 1 of vector Poisson eq. solves directly for F^r.
78  * Some bugs were corrected in the operator poisson_vect.
79  *
80  * Revision 1.8 2004/05/14 08:51:02 p_grandclement
81  * *** empty log message ***
82  *
83  * Revision 1.7 2004/05/10 15:28:22 j_novak
84  * First version of functions for the solution of the r-component of the
85  * vector Poisson equation.
86  *
87  * Revision 1.6 2004/03/05 09:18:49 p_grandclement
88  * Addition of operator sec_order_r2
89  *
90  * Revision 1.5 2004/01/15 09:15:39 p_grandclement
91  * Modification and addition of the Helmholtz operators
92  *
93  * Revision 1.4 2004/01/07 14:36:38 p_grandclement
94  * Modif mineure in Param_elliptic.set_variable
95  *
96  * Revision 1.3 2003/12/11 16:11:38 e_gourgoulhon
97  * Changed #include <iostream.h> to #include "headcpp.h".
98  *
99  * Revision 1.2 2003/12/11 15:57:27 p_grandclement
100  * include stdlib.h encore ...
101  *
102  * Revision 1.1 2003/12/11 14:48:51 p_grandclement
103  * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
104  *
105  *
106  * $Header: /cvsroot/Lorene/C++/Source/Param_elliptic/param_elliptic.C,v 1.23 2018/11/16 14:34:36 j_novak Exp $
107  *
108  */
109 
110 #include "headcpp.h"
111 
112 #include <cmath>
113 #include <cstdlib>
114 
115 #include "map.h"
116 #include "ope_elementary.h"
117 #include "param_elliptic.h"
118 #include "base_val.h"
119 #include "scalar.h"
120 #include "proto.h"
121 
122 
123 namespace Lorene {
125 
126  switch (type_map) {
127  case MAP_AFF:
128  return *mp_af ;
129  break ;
130  case MAP_LOG:
131  return *mp_log ;
132  break ;
133  default:
134  cout << "Unknown mapping in Param_elliptic" << endl ;
135  abort() ;
136  return *mp_af ;
137  }
138 }
139 
140 double Param_elliptic::get_alpha(int l) const {
141 
142  switch (type_map) {
143  case MAP_AFF:
144  return mp_af->get_alpha()[l] ;
145  break ;
146  case MAP_LOG:
147  return mp_log->get_alpha(l) ;
148  break ;
149  default:
150  cout << "Unknown mapping in Param_elliptic" << endl ;
151  abort() ;
152  return 1 ;
153  }
154 }
155 
156 double Param_elliptic::get_beta(int l) const {
157 
158  switch (type_map) {
159  case MAP_AFF:
160  return mp_af->get_beta()[l] ;
161  break ;
162  case MAP_LOG:
163  return mp_log->get_beta(l) ;
164  break ;
165  default:
166  cout << "Unknown mapping in Param_elliptic" << endl ;
167  abort() ;
168  return 1 ;
169  }
170 }
171 
172 int Param_elliptic::get_type(int l) const {
173 
174  switch (type_map) {
175  case MAP_AFF:
176  return AFFINE ;
177  break ;
178  case MAP_LOG:
179  return mp_log->get_type(l) ;
180  break ;
181  default:
182  cout << "Unknown mapping in Param_elliptic" << endl ;
183  abort() ;
184  return 1 ;
185  }
186 }
187 
188 // Construction (By default it is set to Poisson with appropriate dzpuis...)
189 Param_elliptic::Param_elliptic(const Scalar& so) : var_F(so.get_mp()), var_G(so.get_mp()),
190  done_F (so.get_mp().get_mg()->get_nzone(),
191  so.get_mp().get_mg()->get_np(0) + 1,
192  so.get_mp().get_mg()->get_nt(0)),
193  done_G (so.get_mp().get_mg()->get_nzone()),
194  val_F_plus (so.get_mp().get_mg()->get_nzone(),
195  so.get_mp().get_mg()->get_np(0) + 1,
196  so.get_mp().get_mg()->get_nt(0)),
197  val_F_minus (so.get_mp().get_mg()->get_nzone(),
198  so.get_mp().get_mg()->get_np(0) + 1,
199  so.get_mp().get_mg()->get_nt(0)),
200  val_dF_plus (so.get_mp().get_mg()->get_nzone(),
201  so.get_mp().get_mg()->get_np(0) + 1,
202  so.get_mp().get_mg()->get_nt(0)),
203  val_dF_minus (so.get_mp().get_mg()->get_nzone(),
204  so.get_mp().get_mg()->get_np(0) + 1,
205  so.get_mp().get_mg()->get_nt(0)),
206  val_G_plus (so.get_mp().get_mg()->get_nzone()),
207  val_G_minus (so.get_mp().get_mg()->get_nzone()),
208  val_dG_plus (so.get_mp().get_mg()->get_nzone()),
209  val_dG_minus (so.get_mp().get_mg()->get_nzone())
210 
211 
212 {
213 
214  // On passe en Ylm
215  Scalar auxi(so) ;
216  auxi.set_spectral_va().coef() ;
217  auxi.set_spectral_va().ylm() ;
218 
219  Base_val base (auxi.get_spectral_va().base) ;
220  int dzpuis = (auxi.dz_nonzero() ? auxi.get_dzpuis() : 4) ; //## to be modified??
221 
222  // Right now, only applicable with affine mapping
223  const Map_af* map_affine = dynamic_cast <const Map_af*> (&so.get_mp()) ;
224  const Map_log* map_log = dynamic_cast <const Map_log*> (&so.get_mp()) ;
225 
226 
227  if ((map_affine == 0x0) && (map_log == 0x0)) {
228  cout << "Param_elliptic not yet defined on this type of mapping" << endl ;
229  abort() ;
230  }
231  else {
232 
233  if (map_affine != 0x0) {
234  type_map = MAP_AFF ;
235  mp_af = map_affine ;
236  mp_log = 0x0 ;
237  }
238  if (map_log != 0x0) {
239  type_map = MAP_LOG ;
240  mp_af = 0x0 ;
241  mp_log = map_log ;
242  }
243  int nz = get_mp().get_mg()->get_nzone() ;
244  int nbr = 0 ;
245  for (int l=0 ; l<nz ; l++)
246  nbr += (get_mp().get_mg()->get_np(l)+1)*
247  get_mp().get_mg()->get_nt(l) ;
248 
249  operateurs = new Ope_elementary* [nbr] ;
250 
251  for (int l=0 ; l<nbr ; l++)
252  operateurs[l] = 0x0 ;
253 
254  int nr ;
255  int base_r, m_quant, l_quant ;
256 
257  int conte = 0 ;
258  for (int l=0 ; l<nz ; l++) {
259 
260  nr = get_mp().get_mg()->get_nr(l) ;
261 
262  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
263  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
264 
266  (l, k, j, m_quant, l_quant, base_r) ;
267 
268  switch (type_map) {
269  case MAP_AFF:
270  operateurs[conte] = new
271  Ope_poisson(nr, base_r, get_alpha(l),
272  get_beta(l), l_quant, dzpuis) ;
273  break ;
274  case MAP_LOG:
275  if (mp_log->get_type(l) == AFFINE)
276  operateurs[conte] = new
277  Ope_poisson(nr, base_r, get_alpha(l),
278  get_beta(l), l_quant, dzpuis) ;
279  else
280  operateurs[conte] = new
281  Ope_sec_order (nr, base_r, get_alpha(l), get_beta(l),
282  1., 2. , l_quant) ;
283  break ;
284  default :
285  cout << "Unknown mapping in Param_elliptic::Param_elliptic"
286  << endl ;
287 
288  }
289  conte ++ ;
290  }
291  }
292  }
293 
294  // STD VARIABLE CHANGE :
295  var_F.annule_hard() ;
297 
298  var_G.set_etat_qcq() ;
299  var_G = 1 ;
301 
302  done_F.annule_hard() ;
303  done_G.annule_hard() ;
304 
309 
314 }
315 
317 
318  int nbr = 0 ;
319  for (int l=0 ; l<get_mp().get_mg()->get_nzone() ; l++)
320  nbr += (get_mp().get_mg()->get_np(l)+1)*get_mp().get_mg()->get_nt(l) ;
321 
322  for (int i=0 ; i<nbr ; i++)
323  if (operateurs[i] != 0x0)
324  delete operateurs[i] ;
325 
326  delete [] operateurs ;
327 }
328 
330 
331  int np, nt ;
332 
333  int conte = 0 ;
334  for (int l=0 ; l<get_mp().get_mg()->get_nzone() ; l++) {
335 
336  np = get_mp().get_mg()->get_np(l) ;
337  nt = get_mp().get_mg()->get_nt(l) ;
338 
339  for (int k=0 ; k<np+1 ; k++)
340  for (int j=0 ; j<nt ; j++) {
341  if ((operateurs[conte] != 0x0) && (l==zone))
342  operateurs[conte]->inc_l_quant() ;
343  conte ++ ;
344  }
345  }
346 }
347 
348 
349 void Param_elliptic::set_helmholtz_minus (int zone, double masse, Scalar& source) {
350 
351  source.set_spectral_va().coef() ;
352  source.set_spectral_va().ylm() ;
353 
354  int nz = get_mp().get_mg()->get_nzone() ;
355  int nr ;
356  int conte = 0 ;
357  int m_quant, base_r_1d, l_quant ;
358 
359  switch (type_map) {
360 
361  case MAP_AFF:
362 
363  for (int l=0 ; l<nz ; l++) {
364 
365  nr = get_mp().get_mg()->get_nr(l) ;
366 
367  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
368  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
369  if (l==zone) {
370  if (operateurs[conte] != 0x0) {
371  delete operateurs[conte] ;
373  (l, k, j, m_quant, l_quant, base_r_1d) ;
374  operateurs[conte] = new Ope_helmholtz_minus (nr, base_r_1d, l_quant, get_alpha(l),
375  get_beta(l) , masse) ;
376  }
377  }
378  conte ++ ;
379  }
380  }
381  break ;
382 
383  case MAP_LOG :
384  if (mp_log->get_type(zone) != AFFINE) {
385  cout << "Operator not define with LOG mapping..." << endl ;
386  abort() ;
387  }
388  else {
389  for (int l=0 ; l<nz ; l++) {
390 
391  nr = get_mp().get_mg()->get_nr(l) ;
392 
393  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
394  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
395  if (l==zone) {
396  if (operateurs[conte] != 0x0) {
397  delete operateurs[conte] ;
399  (l, k, j, m_quant, l_quant, base_r_1d) ;
400  operateurs[conte] = new Ope_helmholtz_minus (nr, base_r_1d, l_quant,
401  get_alpha(l), get_beta(l), masse) ;
402  }
403  }
404  conte ++ ;
405  }
406  }
407  }
408  break ;
409 
410  default :
411  cout << "Unkown mapping in set_helmhotz_minus" << endl ;
412  abort() ;
413  break ;
414  }
415 }
416 
417 void Param_elliptic::set_helmholtz_plus (int zone, double masse, Scalar& source) {
418 
419  source.set_spectral_va().coef() ;
420  source.set_spectral_va().ylm() ;
421 
422  int nz = get_mp().get_mg()->get_nzone() ;
423  int nr ;
424  int conte = 0 ;
425  int m_quant, base_r_1d, l_quant ;
426 
427  switch (type_map) {
428 
429  case MAP_AFF:
430 
431  for (int l=0 ; l<nz ; l++) {
432 
433  nr = get_mp().get_mg()->get_nr(l) ;
434 
435  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
436  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
437  if (l==zone) {
438  if (operateurs[conte] != 0x0) {
439  int old_base = operateurs[conte]->get_base_r() ;
440  // PROVISOIRE, DANS LE NOYAU SEUL LE CAS SPHERIQUE EST IMPLEMENTE
441  if ((old_base != R_CHEBI)) {
442  delete operateurs[conte] ;
444  (l, k, j, m_quant, l_quant, base_r_1d) ;
445  operateurs[conte] = new Ope_helmholtz_plus (nr, base_r_1d, l_quant, get_alpha(l),
446  get_beta(l) , masse) ;
447  }
448  }
449  }
450  conte ++ ;
451  }
452  }
453  break ;
454 
455  case MAP_LOG :
456  if (mp_log->get_type(zone) != AFFINE) {
457  cout << "Operator not define with LOG mapping..." << endl ;
458  abort() ;
459  }
460  else {
461  for (int l=0 ; l<nz ; l++) {
462 
463  nr = get_mp().get_mg()->get_nr(l) ;
464 
465  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
466  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
467  if (l==zone) {
468  if (operateurs[conte] != 0x0) {
469  int old_base = operateurs[conte]->get_base_r() ;
470  // PROVISOIRE, DANS LE NOYAU SEUL LE CAS SPHERIQUE EST IMPLEMENTE
471  if ((old_base != R_CHEBI)) {
472  delete operateurs[conte] ;
474  (l, k, j, m_quant, l_quant, base_r_1d) ;
475  operateurs[conte] = new Ope_helmholtz_plus (nr, base_r_1d, l_quant,
476  get_alpha(l), get_beta(l), masse) ;
477  }
478  }
479  }
480  conte ++ ;
481  }
482  }
483  }
484  break ;
485 
486  default :
487  cout << "Unkown mapping in set_helmhotz_minus" << endl ;
488  abort() ;
489  break ;
490  }
491 }
492 
493 void Param_elliptic::set_poisson_vect_r(int zone, bool only_l_zero) {
494 
495  if (type_map != MAP_AFF) {
496  cout << "set_poisson_vect_r only defined for an affine mapping..." << endl ;
497  abort() ;
498  }
499  else {
500 
501  int nz = get_mp().get_mg()->get_nzone() ;
502 
503  int nr ;
504 
505  int conte = 0 ;
506  for (int l=0 ; l<nz ; l++) {
507 
508  nr = get_mp().get_mg()->get_nr(l) ;
509 
510  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
511  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
512  if ((operateurs[conte] != 0x0) && (l==zone)) {
513  int old_base = operateurs[conte]->get_base_r() ;
514  Ope_poisson* op_pois =
515  dynamic_cast<Ope_poisson*>(operateurs[conte]) ;
516  assert (op_pois !=0x0) ;
517  int lq_old = op_pois->get_lquant() ;
518  int dzp = op_pois->get_dzpuis() ;
519 
520  delete operateurs[conte] ;
521  if (type_map == MAP_AFF) {
522  if ((!only_l_zero)||(lq_old == 0)) {
523  operateurs[conte] = new Ope_pois_vect_r(nr, old_base,get_alpha(l),
524  get_beta(l), lq_old, dzp) ;
525  }
526  else {
527  operateurs[conte] = 0x0 ;
528  }
529  }
530  else
531  operateurs[conte] = 0x0 ;
532  }
533  conte ++ ;
534  }
535  }
536  }
537 }
538 
540 
541  int nz = get_mp().get_mg()->get_nzone() ;
542 
543  int conte = 0 ;
544  for (int l=0 ; l<nz ; l++) {
545 
546  bool ced = (get_mp().get_mg()->get_type_r(l) == UNSURR ) ;
547  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
548  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
549  if ((operateurs[conte] != 0x0) && (l==zone)) {
550  Ope_poisson* op_pois =
551  dynamic_cast<Ope_poisson*>(operateurs[conte]) ;
552  assert (op_pois !=0x0) ;
553  int lq_old = op_pois->get_lquant() ;
554  if (lq_old == 0) {
555  delete operateurs[conte] ;
556  operateurs[conte] = 0x0 ;
557  }
558  else
559  ced ? op_pois->inc_l_quant() : op_pois->dec_l_quant() ;
560  }
561  conte ++ ;
562  }
563  }
564 }
565 
567 
568  if (type_map != MAP_AFF) {
569  cout << "set_poisson_tens_rr only defined for an affine mapping..."
570  << endl ;
571  abort() ;
572  }
573  else {
574 
575  int nz = get_mp().get_mg()->get_nzone() ;
576 
577  int nr ;
578 
579  int conte = 0 ;
580  for (int l=0 ; l<nz ; l++) {
581 
582  nr = get_mp().get_mg()->get_nr(l) ;
583 
584  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
585  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
586  if ((operateurs[conte] != 0x0) && (l==zone)) {
587  int old_base = operateurs[conte]->get_base_r() ;
588  Ope_poisson* op_pois =
589  dynamic_cast<Ope_poisson*>(operateurs[conte]) ;
590  assert (op_pois !=0x0) ;
591  int lq_old = op_pois->get_lquant() ;
592  int dzp = op_pois->get_dzpuis() ;
593 
594  delete operateurs[conte] ;
595  if (lq_old >= 2) {
596  operateurs[conte] = new Ope_pois_tens_rr
597  (nr, old_base,get_alpha(l), get_beta(l), lq_old, dzp) ;
598  }
599  else
600  operateurs[conte] = 0x0 ;
601  }
602  conte ++ ;
603  }
604  }
605  }
606 }
607 
608 void Param_elliptic::set_sec_order_r2 (int zone, double a, double b, double c){
609 
610 
611  int nz = get_mp().get_mg()->get_nzone() ;
612  int nr ;
613  int conte = 0 ;
614 
615  switch (type_map) {
616 
617  case MAP_AFF:
618 
619  for (int l=0 ; l<nz ; l++) {
620 
621  nr = get_mp().get_mg()->get_nr(l) ;
622 
623  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
624  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
625  if ((operateurs[conte] != 0x0) && (l==zone)) {
626  int old_base = operateurs[conte]->get_base_r() ;
627  // PROVISOIRE, DANS LE NOYAU SEUL LE CAS SPHERIQUE EST IMPLEMENTE
628  if ((old_base != R_CHEBI)) {
629  delete operateurs[conte] ;
630  operateurs[conte] = new Ope_sec_order_r2 (nr, old_base,
631  get_alpha(l),
632  get_beta(l), a, b, c) ;
633  }
634  }
635  conte ++ ;
636  }
637  }
638  break ;
639 
640  case MAP_LOG :
641  if (mp_log->get_type(zone) != AFFINE) {
642  cout << "Operator not define with LOG mapping..." << endl ;
643  abort() ;
644  }
645  else {
646  for (int l=0 ; l<nz ; l++) {
647 
648  nr = get_mp().get_mg()->get_nr(l) ;
649 
650  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
651  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
652  if ((operateurs[conte] != 0x0) && (l==zone)) {
653  int old_base = operateurs[conte]->get_base_r() ;
654  // PROVISOIRE, DANS LE NOYAU SEUL LE CAS SPHERIQUE EST IMPLEMENTE
655  if ((old_base != R_CHEBI)) {
656  delete operateurs[conte] ;
657  operateurs[conte] = new Ope_sec_order_r2 (nr, old_base,
658  get_alpha(l),
659  get_beta(l), a, b, c) ;
660  }
661  }
662  conte ++ ;
663  }
664  }
665  }
666  break ;
667 
668  default :
669  cout << "Unkown mapping in set_sec_order_r2" << endl ;
670  abort() ;
671  break ;
672  }
673 }
674 
675 void Param_elliptic::set_sec_order (int zone, double a, double b, double c){
676 
677  if ((type_map == MAP_AFF) || (mp_log->get_type(zone) == AFFINE)) {
678  cout << "set_sec_order only defined for a log mapping" << endl ;
679  abort() ;
680  }
681  else {
682 
683  int nz = get_mp().get_mg()->get_nzone() ;
684 
685  int nr ;
686 
687  int conte = 0 ;
688  for (int l=0 ; l<nz ; l++) {
689 
690  nr = get_mp().get_mg()->get_nr(l) ;
691 
692  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
693  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
694  if ((operateurs[conte] != 0x0) && (l==zone)) {
695 
696  int old_base = operateurs[conte]->get_base_r() ;
697  // PROVISOIRE, DANS LE NOYAU SEUL LE CAS SPHERIQUE EST IMPLEMENTE
698  if (old_base != R_CHEBI) {
699  delete operateurs[conte] ;
700  operateurs[conte] = new Ope_sec_order (nr, old_base,
701  get_alpha(l),
702  get_beta(l), a, b, c) ;
703  }
704  }
705  conte ++ ;
706  }
707  }
708  }
709 }
710 
711 void Param_elliptic::set_ope_vorton (int zone, Scalar& source) {
712 
713  source.set_spectral_va().coef() ;
714  source.set_spectral_va().ylm() ;
715 
716  int nz = get_mp().get_mg()->get_nzone() ;
717  int nr ;
718  int conte = 0 ;
719  int m_quant, base_r_1d, l_quant ;
720  int dzpuis = source.get_dzpuis() ;
721 
722  switch (type_map) {
723 
724  case MAP_AFF:
725 
726  for (int l=0 ; l<nz ; l++) {
727  int dz = (l==nz-1) ? dzpuis : 0 ;
728  nr = get_mp().get_mg()->get_nr(l) ;
729 
730  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
731  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
732  if (l==zone) {
733  if (operateurs[conte] != 0x0) {
734  delete operateurs[conte] ;
736  (l, k, j, m_quant, l_quant, base_r_1d) ;
737  operateurs[conte] = new Ope_vorton (nr, base_r_1d, get_alpha(l),
738  get_beta(l), l_quant, dz) ;
739  }
740  }
741  conte ++ ;
742  }
743  }
744  break ;
745 
746  case MAP_LOG :
747  if (mp_log->get_type(zone) != AFFINE) {
748  cout << "Operator not define with LOG mapping..." << endl ;
749  abort() ;
750  }
751  else {
752  for (int l=0 ; l<nz ; l++) {
753  int dz = (l==nz-1) ? dzpuis : 0 ;
754  nr = get_mp().get_mg()->get_nr(l) ;
755 
756  for (int k=0 ; k<get_mp().get_mg()->get_np(l)+1 ; k++)
757  for (int j=0 ; j<get_mp().get_mg()->get_nt(l) ; j++) {
758  if (l==zone) {
759  if (operateurs[conte] != 0x0) {
760  delete operateurs[conte] ;
762  (l, k, j, m_quant, l_quant, base_r_1d) ;
763  operateurs[conte] = new Ope_vorton (nr, base_r_1d,
764  get_alpha(l), get_beta(l), l_quant, dz) ;
765  }
766  }
767  conte ++ ;
768  }
769  }
770  }
771  break ;
772 
773  default :
774  cout << "Unkown mapping in set_ope_vorton" << endl ;
775  abort() ;
776  break ;
777  }
778 }
779 
781 
782  assert (so.get_etat() != ETATNONDEF) ;
783  assert (so.get_mp() == get_mp()) ;
784 
785  var_F = so ;
786  done_F.annule_hard() ;
787 }
788 
790 
791  assert (so.get_etat() != ETATNONDEF) ;
792  assert (so.get_mp() == get_mp()) ;
793 
794  var_G = so ;
795  done_G.annule_hard() ;
796 }
797 }
void set_poisson_vect_eta(int zone)
Sets the operator to be a regular elliptic operator to solve for the component of the vector Poisson...
void set_helmholtz_minus(int zone, double mas, Scalar &so)
Set the operator to in one domain (not in the nucleus).
virtual void dec_l_quant()
Decreases the quatum number l by one unit.
Definition: ope_poisson.C:158
void set_variable_F(const Scalar &)
Changes the variable function F.
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
double get_alpha(int l) const
Returns in the domain l.
Definition: map.h:3670
Lorene prototypes.
Definition: app_hor.h:67
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
void annule_hard()
Sets the Itbl to zero in a hard way.
Definition: itbl.C:275
void set_ope_vorton(int zone, Scalar &so)
Set the operator to in one domain (not implemented in the nucleus).
int get_type(int l) const
Returns the type of description in the domain l.
Definition: map.h:3674
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Itbl done_G
Stores what has been computed for G.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Param_elliptic(const Scalar &)
Standard constructor from a Scalar.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Class for operator of the type .
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:813
virtual void inc_l_quant()
Increases the quatum number l by one unit.
Definition: ope_poisson.C:139
Tbl val_F_minus
Values of F at the inner boundaries of the various domains.
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
const Map_af * mp_af
The mapping, if affine.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Class for the operator appearing for the vortons.
Tbl val_G_plus
Values of G at the outer boundaries of the various domains.
Class for the Helmholtz operator (m > 0).
Class for operator of the type .
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
Scalar var_F
Additive variable change function.
Logarithmic radial mapping.
Definition: map.h:3616
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
Tbl val_G_minus
Values of G at the inner boundaries of the various domains.
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
int type_map
Type of mapping either MAP_AFF or MAP_LOG.
Tbl val_F_plus
Values of F at the outer boundaries of the various domains.
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
void set_variable_G(const Scalar &)
Changes the variable function G.
Base class for pure radial mappings.
Definition: map.h:1557
const Map_radial & get_mp() const
Returns the mapping.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual void inc_l_quant()=0
Increases the quatum number l by one unit.
void set_poisson_tens_rr(int zone)
Sets the operator to in all domains, for only.
Class for the operator of the rr component of the divergence-free tensor Poisson equation...
int get_base_r() const
Returns base_r}.
Scalar var_G
Multiplicative variable change that must be sphericaly symetric !
void set_helmholtz_plus(int zone, double mas, Scalar &so)
Set the operator to in one domain (only in the shells).
Tbl val_dG_minus
Values of the derivative of G at the inner boundaries of the various domains.
Class for the operator of the Poisson equation (i.e.
~Param_elliptic()
Destructor.
Ope_elementary ** operateurs
Array on the elementary operators.
const Map_log * mp_log
The mapping if log type.
Class for the operator of the r component of the vector Poisson equation.
Tbl val_dG_plus
Values of the derivative of G at the outer boundaries of the various domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Bases of the spectral expansions.
Definition: base_val.h:325
Affine radial mapping.
Definition: map.h:2048
void inc_l_quant(int zone)
Increases the quantum number l in the domain zone .
bool dz_nonzero() const
Returns true if the last domain is compactified and *this is not zero in this domain.
Definition: scalar.C:820
Basic class for elementary elliptic operators.
Tbl val_dF_minus
Values of the derivative of F at the inner boundaries of the various domains.
double get_beta(int l) const
Returns in the domain l.
Definition: map.h:3672
Itbl done_F
Stores what has been computed for F.
Class for the Helmholtz operator ( ).
void set_sec_order_r2(int zone, double a, double b, double c)
Set the operator to in one domain (only in the shells).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
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
Tbl val_dF_plus
Values of the derivative of F at the outer boundaries of the various domains.
int get_dzpuis()
Returns the associated dzpuis, if in the compactified domain.
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
void set_poisson_vect_r(int zone, bool only_l_zero=false)
Sets the operator to in all domains, for ; and to in all domains otherwise.
void set_sec_order(int zone, double a, double b, double c)
Set the operator to in one domain (only in the shells).
int get_lquant()
Returns the quantum number l.
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:490