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