LORENE
cmp_pde_frontiere.C
1 /*
2  * Copyright (c) 2000-2001 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 as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char cmp_pde_frontiere_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_pde_frontiere.C,v 1.7 2014/10/13 08:52:48 j_novak Exp $" ;
24 
25 /*
26  * $Id: cmp_pde_frontiere.C,v 1.7 2014/10/13 08:52:48 j_novak Exp $
27  * $Log: cmp_pde_frontiere.C,v $
28  * Revision 1.7 2014/10/13 08:52:48 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.6 2005/02/18 13:14:08 j_novak
32  * Changing of malloc/free to new/delete + suppression of some unused variables
33  * (trying to avoid compilation warnings).
34  *
35  * Revision 1.5 2004/11/23 12:49:58 f_limousin
36  * Intoduce function poisson_dir_neu(...) to solve a scalar poisson
37  * equation with a mixed boundary condition (Dirichlet + Neumann).
38  *
39  * Revision 1.4 2004/05/20 07:04:02 k_taniguchi
40  * Recovery of "->get_angu()" in the assertion of Map_af::poisson_frontiere
41  * because "limite" is the boundary value.
42  *
43  * Revision 1.3 2004/03/31 11:18:42 f_limousin
44  * Methods Map_et::poisson_interne, Map_af::poisson_interne and
45  * Cmp::poisson_neumann_interne have been implemented to solve the
46  * continuity equation for strange stars.
47  *
48  * Revision 1.2 2003/10/03 15:58:45 j_novak
49  * Cleaning of some headers
50  *
51  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
52  * LORENE
53  *
54  * Revision 2.6 2000/05/22 16:07:03 phil
55  * *** empty log message ***
56  *
57  * Revision 2.5 2000/05/22 16:03:48 phil
58  * ajout du cas dzpuis = 3
59  *
60  * Revision 2.4 2000/04/27 15:18:27 phil
61  * ajout des procedures relatives a la resolution dans une seule zone avec deux conditions limites.
62  *
63  * Revision 2.3 2000/03/31 15:59:54 phil
64  * gestion des cas ou la source est nulle.
65  *
66  * Revision 2.2 2000/03/20 13:08:53 phil
67  * *** empty log message ***
68  *
69  * Revision 2.1 2000/03/17 17:33:05 phil
70  * *** empty log message ***
71  *
72  * Revision 2.0 2000/03/17 17:25:08 phil
73  * *** empty log message ***
74  *
75  *
76  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_pde_frontiere.C,v 1.7 2014/10/13 08:52:48 j_novak Exp $
77  *
78  */
79 
80 // Header Lorene:
81 #include "scalar.h"
82 #include "param.h"
83 #include "cmp.h"
84 
85 namespace Lorene {
86 Mtbl_cf sol_poisson_frontiere(const Map_af&, const Mtbl_cf&, const Mtbl_cf&,
87  int, int, int, double = 0.,
88  double = 0.) ;
89 
90 Mtbl_cf sol_poisson_frontiere_double (const Map_af&, const Mtbl_cf&, const Mtbl_cf&,
91  const Mtbl_cf&, int) ;
92 
93 Mtbl_cf sol_poisson_interne(const Map_af&, const Mtbl_cf&, const Mtbl_cf&) ;
94 
95 Cmp Cmp::poisson_dirichlet (const Valeur& limite, int num_front) const {
96 
97  Cmp resu(*mp) ;
98  mp->poisson_frontiere (*this, limite, 1, num_front, resu) ;
99  return resu ;
100 }
101 
102 
103 Cmp Cmp::poisson_neumann (const Valeur& limite, int num_front) const {
104 
105  Cmp resu(*mp) ;
106  mp->poisson_frontiere (*this, limite, 2, num_front, resu) ;
107  return resu ;
108 }
109 
111  Param& par, Cmp& resu) const {
112 
113  mp->poisson_interne (*this, limite, par, resu) ;
114  return resu ;
115 }
116 
117 Cmp Cmp::poisson_frontiere_double (const Valeur& lim_func, const Valeur& lim_der,
118  int num_zone) const {
119  Cmp resu(*mp) ;
120  mp->poisson_frontiere_double (*this, lim_func, lim_der, num_zone, resu) ;
121  return resu ;
122 }
123 
124 void Map_et::poisson_frontiere(const Cmp&, const Valeur&, int, int, Cmp&, double, double) const {
125  cout << "Procedure non implantee ! " << endl ;
126  abort() ;
127 }
128 
129 void Map_et::poisson_frontiere_double (const Cmp&, const Valeur&, const Valeur&,
130  int, Cmp&) const {
131  cout << "Procedure non implantee ! " << endl ;
132  abort() ;
133 }
134 
135 void Map_af::poisson_frontiere(const Cmp& source, const Valeur& limite,
136  int type_raccord, int num_front,
137  Cmp& pot, double fact_dir, double fact_neu) const {
138 
139  assert(source.get_etat() != ETATNONDEF) ;
140  assert(source.get_mp()->get_mg() == mg) ;
141  assert(pot.get_mp()->get_mg() == mg) ;
142 
143  assert( source.check_dzpuis(2) || source.check_dzpuis(4)
144  || source.check_dzpuis(3)) ;
145  assert ((type_raccord == 1) || (type_raccord==2)|| (type_raccord==3)) ;
146  int dzpuis ;
147 
148  if (source.dz_nonzero()){
149  dzpuis = source.get_dzpuis() ;
150  }
151  else{
152  dzpuis = 4 ;
153  }
154 
155  // Spherical harmonic expansion of the source
156  // ------------------------------------------
157 
158  Valeur sourva = source.va ;
159  sourva.coef() ;
160  sourva.ylm() ;
161 
162  // Pour gerer le cas ou source est dans ETAT_ZERO...
163  Mtbl_cf so_cf (sourva.get_mg(), sourva.base) ;
164  if (sourva.get_etat() == ETATZERO) {
165  so_cf.set_etat_zero() ;
166  }
167  else
168  so_cf = *sourva.c_cf ;
169 
170 
171  Valeur conditions (limite) ;
172  conditions.coef() ;
173  conditions.ylm() ; // spherical harmonic transforms
174 
175  // Pour gerer le cas ou condition est dans ETAT_ZERO...
176  Mtbl_cf auxiliaire (conditions.get_mg(), conditions.base) ;
177  if (conditions.get_etat() == ETATZERO)
178  auxiliaire.set_etat_zero() ;
179  else
180  auxiliaire = *conditions.c_cf ;
181 
182  Mtbl_cf resu (sourva.get_mg(), sourva.base) ;
183 
184  if (type_raccord == 3){
185 
186  // Call to the Mtbl_cf version
187  // ---------------------------
188  resu = sol_poisson_frontiere(*this, so_cf, auxiliaire,
189  type_raccord, num_front, dzpuis,
190  fact_dir, fact_neu) ;
191  }
192  else{
193  resu = sol_poisson_frontiere(*this, so_cf, auxiliaire,
194  type_raccord, num_front, dzpuis) ;
195  }
196 
197  // Final result returned as a Cmp
198  // ------------------------------
199 
200  pot.set_etat_zero() ; // to call Cmp::del_t().
201  pot.set_etat_qcq() ;
202  pot.va = resu ;
203  (pot.va).ylm_i() ; // On repasse en base standard.
204  pot.set_dzpuis(0) ;
205 }
206 
207 
208 void Map_af::poisson_frontiere_double (const Cmp& source, const Valeur& lim_func,
209  const Valeur& lim_der, int num_zone, Cmp& pot) const {
210 
211  assert(source.get_etat() != ETATNONDEF) ;
212  assert(source.get_mp()->get_mg() == mg) ;
213  assert(pot.get_mp()->get_mg() == mg) ;
214  assert (source.get_mp()->get_mg()->get_angu() == lim_func.get_mg()) ;
215  assert (source.get_mp()->get_mg()->get_angu() == lim_der.get_mg()) ;
216 
217  // Spherical harmonic expansion of the source
218  // ------------------------------------------
219 
220  Valeur sourva = source.va ;
221  sourva.coef() ;
222  sourva.ylm() ;
223 
224  // Pour gerer le cas ou source est dans ETAT_ZERO...
225  Mtbl_cf so_cf (sourva.get_mg(), sourva.base) ;
226  if (sourva.get_etat() == ETATZERO) {
227  so_cf.set_etat_zero() ;
228  }
229  else
230  so_cf = *sourva.c_cf ;
231 
232 
233  Valeur cond_func (lim_func) ;
234  cond_func.coef() ;
235  cond_func.ylm() ; // spherical harmonic transforms
236 
237  // Pour gerer le cas ou condition est dans ETAT_ZERO...
238  Mtbl_cf auxi_func (cond_func.get_mg(), cond_func.base) ;
239  if (cond_func.get_etat() == ETATZERO)
240  auxi_func.set_etat_zero() ;
241  else
242  auxi_func = *cond_func.c_cf ;
243 
244  Valeur cond_der (lim_der) ;
245  cond_der.coef() ;
246  cond_der.ylm() ; // spherical harmonic transforms
247 
248  // Pour gerer le cas ou condition est dans ETAT_ZERO...
249  Mtbl_cf auxi_der (cond_der.get_mg(), cond_der.base) ;
250  if (cond_der.get_etat() == ETATZERO)
251  auxi_der.set_etat_zero() ;
252  else
253  auxi_der = *cond_der.c_cf ;
254 
255 
256 
257  // Call to the Mtbl_cf version
258  // ---------------------------
259 
260  Mtbl_cf resu = sol_poisson_frontiere_double (*this, so_cf, auxi_func,
261  auxi_der, num_zone) ;
262 
263  // Final result returned as a Cmp
264  // ------------------------------
265 
266  pot.set_etat_zero() ; // to call Cmp::del_t().
267  pot.set_etat_qcq() ;
268  pot.va = resu ;
269  (pot.va).ylm_i() ; // On repasse en base standard.
270 }
271 
272 
273 
274 void Map_et::poisson_interne(const Cmp& source, const Valeur& limite,
275  Param& par, Cmp& uu) const {
276 
277  assert(source.get_etat() != ETATNONDEF) ;
278  assert(source.get_mp() == this) ;
279 
280  assert(uu.get_mp() == this) ;
281 
282  int nz = mg->get_nzone() ;
283 
284  //-------------------------------
285  // Computation of the prefactor a ---> Cmp apre
286  //-------------------------------
287 
288  Mtbl unjj = 1 + srdrdt*srdrdt + srstdrdp*srstdrdp ;
289 
290  Mtbl apre1(*mg) ;
291  apre1.set_etat_qcq() ;
292  for (int l=0; l<nz; l++) {
293  *(apre1.t[l]) = alpha[l]*alpha[l] ;
294  }
295 
296  apre1 = apre1 * dxdr * dxdr * unjj ;
297 
298  Cmp apre(*this) ;
299  apre = apre1 ;
300 
301  Tbl amax0 = max(apre1) ; // maximum values in each domain
302 
303  // The maximum values of a in each domain are put in a Mtbl
304  Mtbl amax1(*mg) ;
305  amax1.set_etat_qcq() ;
306  for (int l=0; l<nz; l++) {
307  *(amax1.t[l]) = amax0(l) ;
308  }
309 
310  Cmp amax(*this) ;
311  amax = amax1 ;
312 
313 
314  //-------------------
315  // Initializations
316  //-------------------
317 
318  int nitermax = par.get_int() ;
319  int& niter = par.get_int_mod() ;
320  double lambda = par.get_double() ;
321  double unmlambda = 1. - lambda ;
322  double precis = par.get_double(1) ;
323 
324  Cmp& ssj = par.get_cmp_mod() ;
325 
326  Cmp ssjm1 = ssj ;
327  Cmp ssjm2 = ssjm1 ;
328 
329  Valeur& vuu = uu.va ;
330 
331  Valeur vuujm1(*mg) ;
332  if (uu.get_etat() == ETATZERO) {
333  vuujm1 = 1 ; // to take relative differences
334  vuujm1.set_base( vuu.base ) ;
335  }
336  else{
337  vuujm1 = vuu ;
338  }
339 
340  // Affine mapping for the Laplacian-tilde
341 
342  Map_af mpaff(*this) ;
343  Param par_nul ;
344 
345  cout << "Map_et::poisson : relat. diff. u^J <-> u^{J-1} : " << endl ;
346 
347 //==========================================================================
348 //==========================================================================
349 // Start of iteration
350 //==========================================================================
351 //==========================================================================
352 
353  Tbl tdiff(nz) ;
354  niter = 0 ;
355 
356  do {
357 
358  //====================================================================
359  // Computation of R(u) (the result is put in uu)
360  //====================================================================
361 
362 
363  //-----------------------
364  // First operations on uu
365  //-----------------------
366 
367  Valeur duudx = (uu.va).dsdx() ; // d/dx
368 
369  const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
370 
371  const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
372 
373  //------------------
374  // Angular Laplacian
375  //------------------
376 
377  Valeur sxlapang = uu.va ;
378 
379  sxlapang.ylm() ;
380 
381  sxlapang = sxlapang.lapang() ;
382 
383  sxlapang = sxlapang.sx() ; // Lap_ang(uu) /x in the nucleus
384 
385  //---------------------------------------------------------------
386  // Computation of
387  // [ 2 /(dRdx) ( A - 1 ) duu/dx + 1/R (B - 1) Lap_ang(uu) ] / x
388  //
389  // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
390  // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
391  //
392  // The result is put in uu (via vuu)
393  //---------------------------------------------------------------
394 
395  Valeur varduudx = duudx ;
396 
397  uu.set_etat_qcq() ;
398 
399  Base_val sauve_base = varduudx.base ;
400 
401  vuu = 2. * dxdr * ( rsxdxdr * unjj - 1.) * varduudx
402  + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
403 
404  vuu.set_base(sauve_base) ;
405 
406  vuu = vuu.sx() ;
407 
408  //---------------------------------------
409  // Computation of R(u)
410  //
411  // The result is put in uu (via vuu)
412  //---------------------------------------
413 
414 
415  sauve_base = vuu.base ;
416 
417  vuu = xsr * vuu
418  + 2. * dxdr * ( sr2drdt * d2uudtdx
419  + sr2stdrdp * std2uudpdx ) ;
420 
421  vuu += dxdr * ( lapr_tp + dxdr * (
422  dxdr* unjj * d2rdx2
423  - 2. * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
424  ) * duudx ;
425 
426  vuu.set_base(sauve_base) ;
427 
428  //====================================================================
429  // Computation of the effective source s^J of the ``affine''
430  // Poisson equation
431  //====================================================================
432 
433  ssj = lambda * ssjm1 + unmlambda * ssjm2 ;
434 
435  ssj = ( source + uu + (amax - apre) * ssj ) / amax ;
436 
437  (ssj.va).set_base((source.va).base) ;
438 
439  //====================================================================
440  // Resolution of the ``affine'' Poisson equation
441  //====================================================================
442 
443  mpaff.poisson_interne(ssj, limite, par_nul, uu) ;
444 
445  tdiff = diffrel(vuu, vuujm1) ;
446 
447  cout << " step " << niter << " : " ;
448  cout << tdiff(0) << " " ;
449 
450  cout << endl ;
451 
452  //=================================
453  // Updates for the next iteration
454  //=================================
455 
456  ssjm2 = ssjm1 ;
457  ssjm1 = ssj ;
458  vuujm1 = vuu ;
459 
460  niter++ ;
461 
462  } // End of iteration
463  while ( (tdiff(0) > precis) && (niter < nitermax) ) ;
464 
465 //==========================================================================
466 //==========================================================================
467 // End of iteration
468 //==========================================================================
469 //==========================================================================
470 
471 }
472 
473 
474 void Map_af::poisson_interne(const Cmp& source, const Valeur& limite,
475  Param& , Cmp& pot) const {
476 
477  assert(source.get_etat() != ETATNONDEF) ;
478  assert(source.get_mp()->get_mg() == mg) ;
479  assert(pot.get_mp()->get_mg() == mg) ;
480  assert (source.get_mp()->get_mg()->get_angu() == limite.get_mg()) ;
481 
482  // Spherical harmonic expansion of the source
483  // ------------------------------------------
484 
485  Valeur sourva = source.va ;
486  sourva.coef() ;
487  sourva.ylm() ;
488 
489  // Pour gerer le cas ou source est dans ETAT_ZERO...
490  Mtbl_cf so_cf (sourva.get_mg(), sourva.base) ;
491  if (sourva.get_etat() == ETATZERO) {
492  so_cf.set_etat_zero() ;
493  }
494  else
495  so_cf = *sourva.c_cf ;
496 
497  Valeur conditions (limite) ;
498  conditions.coef() ;
499  conditions.ylm() ; // spherical harmonic transforms
500 
501 
502  // Pour gerer le cas ou condition est dans ETAT_ZERO...
503  Mtbl_cf auxiliaire (conditions.get_mg(), conditions.base) ;
504  if (conditions.get_etat() == ETATZERO)
505  auxiliaire.set_etat_zero() ;
506  else
507  auxiliaire = *conditions.c_cf ;
508 
509  // Call to the Mtbl_cf version
510  // ---------------------------
511 
512  Mtbl_cf resu = sol_poisson_interne(*this, so_cf, auxiliaire) ;
513 
514  // Final result returned as a Cmp
515  // ------------------------------
516 
517  pot.set_etat_zero() ; // to call Cmp::del_t().
518  pot.set_etat_qcq() ;
519  pot.va = resu ;
520  (pot.va).ylm_i() ; // On repasse en base standard.
521  pot.set_dzpuis(0) ;
522 }
523 
524 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:898
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:112
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
const Cmp & dsdx() const
Returns of *this , where .
Definition: cmp_deriv.C:148
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:72
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:64
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_etat() const
Returns the logical state.
Definition: cmp.h:896
virtual void poisson_frontiere_double(const Cmp &source, const Valeur &lim_func, const Valeur &lim_der, int num_zone, Cmp &pot) const
Solver of the Poisson equation with boundary condition for the affine mapping case, cases with boundary conditions of both Dirichlet and Neumann type (no condition imposed at infinity).
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx.C:110
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1049
virtual void poisson_interne(const Cmp &source, const Valeur &limite, Param &par, Cmp &pot) const =0
Computes the solution of a Poisson equation in the shell, imposing a boundary condition at the surfac...
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
virtual void poisson_frontiere(const Cmp &, const Valeur &, int, int, Cmp &, double=0., double=0.) const
Not yet implemented.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified zon...
Definition: cmp.h:458
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
virtual void poisson_frontiere(const Cmp &, const Valeur &, int, int, Cmp &, double=0., double=0.) const
Solver of the Poisson equation with boundary condition for the affine mapping case.
bool dz_nonzero() const
Returns true if the last domain is compactified and *this is not zero in this domain.
Definition: cmp.C:660
Parameter storage.
Definition: param.h:125
const Valeur & stdsdp() const
Returns of *this.
Definition: valeur_stdsdp.C:60
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
Cmp poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Cmp::poisson() .
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
virtual void poisson_interne(const Cmp &source, const Valeur &limite, Param &par, Cmp &pot) const
Computes the solution of a Poisson equation in the shell .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: mtbl_cf.C:288
virtual void poisson_frontiere(const Cmp &source, const Valeur &limite, int raccord, int num_front, Cmp &pot, double=0., double=0.) const =0
Computes the solution of a Poisson equation from the domain num_front+1 .
Bases of the spectral expansions.
Definition: base_val.h:322
Affine radial mapping.
Definition: map.h:2027
const Map * mp
Reference mapping.
Definition: cmp.h:448
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:900
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: cmp.C:715
virtual void poisson_interne(const Cmp &source, const Valeur &limite, Param &par, Cmp &pot) const
Computes the solution of a Poisson equation in the shell, imposing a boundary condition at the surfac...
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:361
Basic array class.
Definition: tbl.h:161
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:461
Cmp poisson_neumann_interne(const Valeur &, Param &par, Cmp &resu) const
Idem as Cmp::poisson_neumann , the boundary condition is on the radial derivative of the solution...