LORENE
vector_poisson_block.C
1 /*
2  * Method for vector Poisson equation inverting eqs. for V^r and eta as a block.
3  *
4  * (see file vector.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: vector_poisson_block.C,v 1.10 2016/12/05 16:18:18 j_novak Exp $
32  * $Log: vector_poisson_block.C,v $
33  * Revision 1.10 2016/12/05 16:18:18 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.9 2014/10/13 08:53:45 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.8 2014/10/06 15:13:21 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.7 2012/10/12 11:43:38 j_novak
43  * Removed some headers
44  *
45  * Revision 1.6 2007/12/21 10:45:06 j_novak
46  * Better treatment of spherical symmetry
47  *
48  * Revision 1.5 2007/09/05 12:35:18 j_novak
49  * Homogeneous solutions are no longer obtained through the analytic formula, but
50  * by solving (again) the operator system with almost zero as r.h.s. This is to
51  * lower the condition number of the matching system.
52  *
53  * Revision 1.4 2007/01/23 17:08:46 j_novak
54  * New function pois_vect_r0.C to solve the l=0 part of the vector Poisson
55  * equation, which involves only the r-component.
56  *
57  * Revision 1.3 2005/12/30 13:39:38 j_novak
58  * Changed the Galerkin base in the nucleus to (hopefully) stabilise the solver
59  * when used in an iteration. Similar changes in the CED too.
60  *
61  * Revision 1.2 2005/02/15 15:43:18 j_novak
62  * First version of the block inversion for the vector Poisson equation (method 6).
63  *
64  *
65  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_block.C,v 1.10 2016/12/05 16:18:18 j_novak Exp $
66  *
67  */
68 
69 // C headers
70 #include <cmath>
71 
72 // Lorene headers
73 #include "metric.h"
74 #include "diff.h"
75 #include "proto.h"
76 
77 namespace Lorene {
78 void Vector::poisson_block(double lam, Vector& resu) const {
79 
80  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
81 #ifndef NDEBUG
82  for (int i=0; i<3; i++)
83  assert(cmp[i]->check_dzpuis(4)) ;
84  // All this has a meaning only for spherical components:
85  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
86  assert(bvs != 0x0) ;
87  //## ... and affine mapping, for the moment!
88  assert( mpaff != 0x0) ;
89 #endif
90 
91  if (fabs(lam + 1.) < 1.e-8) {
92  const Metric_flat& mets = mp->flat_met_spher() ;
93  Vector_divfree sou(*mp, *triad, mets) ;
94  for (int i=1; i<=3; i++) sou.set(i) = *cmp[i-1] ;
95  resu = sou.poisson() ;
96  return ;
97  }
98 
99  // Some working objects
100  //---------------------
101  const Mg3d& mg = *mpaff->get_mg() ;
102  int nz = mg.get_nzone() ; int nzm1 = nz - 1;
103  int nt = mg.get_nt(0) ;
104  int np = mg.get_np(0) ;
105  assert(mg.get_type_r(0) == RARE) ;
106  assert(mg.get_type_r(nzm1) == UNSURR) ;
107  Scalar S_r = *cmp[0] ;
108  Scalar S_eta = eta() ;
109  Scalar het(*mpaff) ;
110  Scalar vr(*mpaff) ;
111  bool all_zero = false ;
112  if (S_r.get_etat() == ETATZERO) {
113  if (S_eta.get_etat() == ETATZERO) {
114  vr.set_etat_zero() ;
115  het.set_etat_zero() ;
116  all_zero = true ;
117  }
118  else {
119  S_r.annule_hard() ;
120  S_r.set_spectral_base(S_eta.get_spectral_base()) ;
121  }
122  }
123  if ((S_eta.get_etat() == ETATZERO)&&(!all_zero)) {
124  S_eta.annule_hard() ;
125  S_eta.set_spectral_base(S_r.get_spectral_base()) ;
126  }
127  if (!all_zero) {
128  S_r.set_spectral_va().ylm() ;
129  S_eta.set_spectral_va().ylm() ;
130  const Base_val& base = S_eta.get_spectral_va().base ;
131  Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
132  Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
133  Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
134  Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
135  Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
136  Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
137  Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
138  Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
139  Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
140  Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
141 
142  // Solution of the l=0 part for Vr
143  //---------------------------------
144  Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
145  sou_l0.set_spectral_va().ylm() ;
146  Scalar vrl0 = pois_vect_r0(sou_l0) ;
147 
148  // Build-up & inversion of the system for (eta, V^r) in each domain
149  //-----------------------------------------------------------------
150 
151  // Nucleus
152  //--------
153  int nr = mg.get_nr(0) ;
154  double alpha = mpaff->get_alpha()[0] ; double alp2 = alpha*alpha ;
155  double beta = mpaff->get_beta()[0] ;
156  int l_q = 0 ; int m_q = 0; int base_r = 0 ;
157 
158  // Loop on l and m
159  //----------------
160  for (int k=0 ; k<np+1 ; k++) {
161  for (int j=0 ; j<nt ; j++) {
162  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
163  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0) ) {
164  int aa = 0 ; int bb = 0 ; int nr0 = 0 ;
165  if (base_r == R_CHEBP) {
166  nr0 = nr - 1 ;
167  aa = 0 ; bb = 1 ;
168  }
169  else {
170  assert (base_r == R_CHEBI) ;
171  nr0 = nr - 2 ;
172  aa = 2 ; bb = 1 ;
173  }
174  int d0 = nr - nr0 ;
175  int nrtot = 2*nr0 ;
176  Matrice oper(2*nr, 2*nr) ; oper.set_etat_qcq() ;
177  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
178  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
179  Diff_sxdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
180  Diff_sx2 s2(base_r, nr) ; const Matrice& ms2 = s2.get_matrice() ;
181 
182  // Building the operator
183  //----------------------
184  for (int lin=0; lin<nr; lin++) { //eq.1
185  for (int col=0; col<nr; col++)
186  oper.set(lin,col) =
187  (md2(lin,col) + 2*mxd(lin,col)
188  -(lam+1)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
189  for (int col=0; col<nr; col++)
190  oper.set(lin,col+nr) =
191  (lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
192  }
193  for (int lin=0; lin<nr; lin++) { //eq.2
194  for (int col=0; col<nr; col++)
195  oper.set(lin+nr,col) =
196  (-lam*l_q*(l_q+1)*mxd(lin,col)
197  +(lam+2)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
198  for (int col=0; col<nr; col++)
199  oper.set(lin+nr, col+nr) =
200  ((lam+1)*(md2(lin,col) + 2*mxd(lin,col))
201  - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
202  }
203  bool pb_eta = ( ( fabs( lam*double(l_q+3) + 2 ) < 0.01) && (l_q <=2) ) ;
204  if (!pb_eta) {
205  for (int col=0; col<nr; col++)
206  oper.set(nr0-1, col) = 1 ;
207  for (int col=0; col<nr; col++)
208  oper.set(nr0-1,nr+col) = 0 ;
209  }
210  if ((l_q > 2)||pb_eta) {
211  for (int col=0; col<nr; col++)
212  oper.set(nr+nr0-1, col) = 0 ;
213  for (int col=0; col<nr; col++)
214  oper.set(nr+nr0-1,nr+col) = 1 ;
215  }
216 
217  Matrice op2(nrtot, nrtot) ;
218  op2.set_etat_qcq() ;
219  for (int i=0; i<nr0; i++) {
220  for (int col=0; col<nr0;col++)
221  op2.set(i,col) = (aa*col+bb)*oper(i,col+1)
222  + (aa*(col+1)+bb)*oper(i,col) ;
223  for (int col=0; col<nr0;col++)
224  op2.set(i,col+nr0) = (aa*col+bb)*oper(i,col+nr+1)
225  + (aa*(col+1)+bb)*oper(i,col+nr) ;
226  }
227 
228  for (int i=nr0; i<nrtot; i++) {
229  for (int col=0; col<nr0;col++)
230  op2.set(i,col) = (aa*col+bb)*oper(i+d0,col+1)
231  + (aa*(col+1)+bb)*oper(i+d0,col) ;
232  for (int col=0; col<nr0;col++)
233  op2.set(i,col+nr0) = (aa*col+bb)*oper(i+d0,col+nr+1)
234  + (aa*(col+1)+bb)*oper(i+d0,col+nr) ;
235  }
236  op2.set_lu() ;
237 
238  // Filling the r.h.s
239  //------------------
240  for (int i=0; i<nr0; i++) //eq.1
241  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
242  if (!pb_eta) sec_membre.set(nr0-1) = 0 ;
243  for (int i=0; i<nr0; i++) //eq.2
244  sec_membre.set(i+nr0)
245  = (*S_r.get_spectral_va().c_cf)(0, k, j, i) ;
246  if ((l_q > 2)||pb_eta) sec_membre.set(nrtot-1) = 0 ;
247 
248  // Inversion of the "big" operator
249  //--------------------------------
250  Tbl big_res = op2.inverse(sec_membre) ;
251 
252  // Putting coefficients of eta and Vr to individual arrays
253  //--------------------------------------------------------
254  sol_part_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
255  for (int i=1; i<nr0; i++)
256  sol_part_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
257  + (aa*(i+1)+bb)*big_res(i);
258  sol_part_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
259  sol_part_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
260  for (int i=1; i<nr0; i++)
261  sol_part_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
262  + (aa*(i+1)+bb)*big_res(nr0+i);
263  sol_part_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
264  if (base_r == R_CHEBI) {
265  sol_part_eta.set(0, k, j, nr-1) = 0 ;
266  sol_part_vr.set(0, k, j, nr-1) = 0 ;
267  }
268 
269  // Homogeneous solutions (only r^(l-1) and r^(l+1) in the nucleus)
270  if (l_q<=2) {
271  double fac_eta = lam*double(l_q+3) + 2. ;
272  double fac_vr = double(l_q+1)*(lam*l_q - 2.) ;
273  Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
274  Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
275  for (int i=0 ; i<nr ; i++) {
276  sol_hom_un_eta.set(0, k, j, i) = sol_hom1(i) ;
277  sol_hom_un_vr.set(0, k, j, i) = l_q*sol_hom1(i) ;
278  sol_hom_trois_eta.set(0, k, j, i) = fac_eta*sol_hom2(i) ;
279  sol_hom_trois_vr.set(0, k, j, i) = fac_vr*sol_hom2(i) ;
280  }
281  }
282  else {
283  for (int i=0; i<nrtot; i++) sec_membre.set(i) = 0 ;
284  // First homogeneous sol.
285  sec_membre.set(nr0-1) = 1 ;
286  big_res = op2.inverse(sec_membre) ;
287 
288  sol_hom_un_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
289  for (int i=1; i<nr0; i++)
290  sol_hom_un_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
291  + (aa*(i+1)+bb)*big_res(i);
292  sol_hom_un_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
293  sol_hom_un_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
294  for (int i=1; i<nr0; i++)
295  sol_hom_un_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
296  + (aa*(i+1)+bb)*big_res(nr0+i);
297  sol_hom_un_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
298  if (base_r == R_CHEBI) {
299  sol_hom_un_eta.set(0, k, j, nr-1) = 0 ;
300  sol_hom_un_vr.set(0, k, j, nr-1) = 0 ;
301  }
302 
303  sec_membre.set(nr0-1) = 0 ;
304  sec_membre.set(nrtot-1) = 1 ;
305  big_res = op2.inverse(sec_membre) ;
306 
307  sol_hom_trois_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
308  for (int i=1; i<nr0; i++)
309  sol_hom_trois_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
310  + (aa*(i+1)+bb)*big_res(i);
311  sol_hom_trois_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1)+bb) ;
312  sol_hom_trois_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
313  for (int i=1; i<nr0; i++)
314  sol_hom_trois_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
315  + (aa*(i+1)+bb)*big_res(nr0+i);
316  sol_hom_trois_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
317  if (base_r == R_CHEBI) {
318  sol_hom_trois_eta.set(0, k, j, nr-1) = 0 ;
319  sol_hom_trois_vr.set(0, k, j, nr-1) = 0 ;
320  }
321 
322  }
323  }
324  }
325  }
326 
327 
328  // Shells
329  //-------
330  for (int zone=1 ; zone<nzm1 ; zone++) {
331  nr = mg.get_nr(zone) ;
332  assert (nr > 5) ;
333  assert(nt == mg.get_nt(zone)) ;
334  assert(np == mg.get_np(zone)) ;
335  alpha = mpaff->get_alpha()[zone] ;
336  beta = mpaff->get_beta()[zone] ;
337  double ech = beta / alpha ;
338 
339  // Loop on l and m
340  //----------------
341  for (int k=0 ; k<np+1 ; k++) {
342  for (int j=0 ; j<nt ; j++) {
343  base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
344  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
345  int dege1 = 2 ; //degeneracy of eq.1
346  int dege2 = 2 ; //degeneracy of eq.2
347  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
348  int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
349  int nrtot = 2*nr ;
350  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
351  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
352  Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
353  Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
354  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
355  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
356  Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
357  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
358 
359  // Building the operator
360  //----------------------
361  for (int lin=0; lin<nr_eq1; lin++) {
362  for (int col=0; col<nr; col++)
363  oper.set(lin,col)
364  = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
365  + 2*(mxd(lin,col) + ech*md(lin,col))
366  - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
367  for (int col=0; col<nr; col++)
368  oper.set(lin,col+nr)
369  = lam*(mxd(lin,col) + ech*md(lin,col))
370  + 2*(1+lam)*mid(lin,col) ;
371  }
372  oper.set(nr_eq1, 0) = 1 ;
373  for (int col=1; col<2*nr; col++)
374  oper.set(nr_eq1, col) = 0 ;
375  oper.set(nr_eq1+1, 0) = 0 ;
376  oper.set(nr_eq1+1, 1) = 1 ;
377  for (int col=2; col<2*nr; col++)
378  oper.set(nr_eq1+1, col) = 0 ;
379  for (int lin=0; lin<nr_eq2; lin++) {
380  for (int col=0; col<nr; col++)
381  oper.set(lin+nr,col)
382  = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
383  - (lam+2)*mid(lin,col)) ;
384  for (int col=0; col<nr; col++)
385  oper.set(lin+nr, col+nr)
386  = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
387  + ech*ech*md2(lin,col)
388  + 2*(mxd(lin,col) + ech*md(lin,col)))
389  -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
390  }
391  for (int col=0; col<nr; col++)
392  oper.set(nr+nr_eq2, col) = 0 ;
393  oper.set(nr+nr_eq2, nr) = 1 ;
394  for (int col=nr+1; col<2*nr; col++)
395  oper.set(nr+nr_eq2, col) = 0 ;
396  for (int col=0; col<nr+1; col++)
397  oper.set(nr+nr_eq2+1, col) = 0 ;
398  oper.set(nr+nr_eq2+1, nr+1) = 1 ;
399  for (int col=nr+2; col<2*nr; col++)
400  oper.set(nr+nr_eq2+1, col) = 0 ;
401 
402  oper.set_lu() ;
403 
404  // Filling the r.h.s
405  //------------------
406  Tbl sr(nr) ; sr.set_etat_qcq() ;
407  Tbl seta(nr) ; seta.set_etat_qcq() ;
408  for (int i=0; i<nr; i++) {
409  sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
410  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
411  }
412  Tbl xsr= sr ; Tbl x2sr= sr ;
413  Tbl xseta= seta ; Tbl x2seta = seta ;
414  multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
415  multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
416 
417  for (int i=0; i<nr_eq1; i++)
418  sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
419  + beta*beta*seta(i);
420  sec_membre.set(nr_eq1) = 0 ;
421  sec_membre.set(nr_eq1+1) = 0 ;
422  for (int i=0; i<nr_eq2; i++)
423  sec_membre.set(i+nr) = beta*beta*sr(i)
424  + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
425  sec_membre.set(nr+nr_eq2) = 0 ;
426  sec_membre.set(nr+nr_eq2+1) = 0 ;
427 
428  // Inversion of the "big" operator
429  //--------------------------------
430  Tbl big_res = oper.inverse(sec_membre) ;
431  for (int i=0; i<nr; i++) {
432  sol_part_eta.set(zone, k, j, i) = big_res(i) ;
433  sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
434  }
435 
436  // Getting the homogeneous solutions
437  //----------------------------------
438  sec_membre.annule_hard() ;
439  sec_membre.set(nr_eq1) = 1 ;
440  big_res = oper.inverse(sec_membre) ;
441  for (int i=0 ; i<nr ; i++) {
442  sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
443  sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
444  }
445  sec_membre.set(nr_eq1) = 0 ;
446  sec_membre.set(nr_eq1+1) = 1 ;
447  big_res = oper.inverse(sec_membre) ;
448  for (int i=0 ; i<nr ; i++) {
449  sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
450  sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
451  }
452  sec_membre.set(nr_eq1+1) = 0 ;
453  sec_membre.set(nr+nr_eq2) = 1 ;
454  big_res = oper.inverse(sec_membre) ;
455  for (int i=0 ; i<nr ; i++) {
456  sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
457  sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
458  }
459  sec_membre.set(nr+nr_eq2) = 0 ;
460  sec_membre.set(nr+nr_eq2+1) = 1 ;
461  big_res = oper.inverse(sec_membre) ;
462  for (int i=0 ; i<nr ; i++) {
463  sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
464  sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
465  }
466 
467  }
468  }
469  }
470  }
471 
472  // Compactified external domain
473  //-----------------------------
474  nr = mg.get_nr(nzm1) ;
475  assert(nt == mg.get_nt(nzm1)) ;
476  assert(np == mg.get_np(nzm1)) ;
477  alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
478  assert (nr > 4) ;
479 
480  // Loop on l and m
481  //----------------
482  for (int k=0 ; k<np+1 ; k++) {
483  for (int j=0 ; j<nt ; j++) {
484  base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
485  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
486  int dege1 = 3; //degeneracy of eq.1
487  int dege2 = (l_q == 1 ? 2 : 3); //degeneracy of eq.2
488  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
489  int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
490  int nrtot = 2*nr ;
491  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
492  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
493  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
494  Diff_sxdsdx sxd(base_r, nr) ; const Matrice& mxd = sxd.get_matrice() ;
495  Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
496 
497  // Building the operator
498  //----------------------
499  for (int lin=0; lin<nr_eq1; lin++) {
500  for (int col=0; col<nr; col++)
501  oper.set(lin,col)
502  = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
503  for (int col=0; col<nr; col++)
504  oper.set(lin,col+nr) =
505  (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
506  }
507  for (int col=0; col<nr; col++)
508  oper.set(nr_eq1, col) = 1 ;
509  for (int col=nr; col<nrtot; col++)
510  oper.set(nr_eq1, col) = 0 ;
511  int pari = -1 ;
512  for (int col=0; col<nr; col++) {
513  oper.set(nr_eq1+1, col) = pari*col*col ;
514  pari = pari ;
515  }
516  for (int col=nr; col<nrtot; col++)
517  oper.set(nr_eq1+1, col) = 0 ;
518  oper.set(nr_eq1+2, 0) = 1 ;
519  for (int col=1; col<nrtot; col++)
520  oper.set(nr_eq1+2, col) = 0 ;
521  for (int lin=0; lin<nr_eq2; lin++) {
522  for (int col=0; col<nr; col++)
523  oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
524  + (lam+2)*ms2(lin,col))) / alp2 ;
525  for (int col=0; col<nr; col++)
526  oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
527  - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
528  }
529  for (int col=0; col<nr; col++)
530  oper.set(nr+nr_eq2, col) = 0 ;
531  for (int col=nr; col<nrtot; col++)
532  oper.set(nr+nr_eq2, col) = 1 ;
533  for (int col=0; col<nr; col++)
534  oper.set(nr+nr_eq2+1, col) = 0 ;
535  pari = -1 ;
536  for (int col=0; col<nr; col++) {
537  oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
538  pari = pari ;
539  }
540  if (l_q > 1) {
541  for (int col=0; col<nr; col++)
542  oper.set(nr+nr_eq2+2, col) = 0 ;
543  oper.set(nr+nr_eq2+2, nr) = 1 ;
544  for (int col=nr+1; col<nrtot; col++)
545  oper.set(nr+nr_eq2+2, col) = 0 ;
546  }
547  oper.set_lu() ;
548 
549  // Filling the r.h.s
550  //------------------
551  for (int i=0; i<nr_eq1; i++)
552  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
553  for (int i=nr_eq1; i<nr; i++)
554  sec_membre.set(i) = 0 ;
555  for (int i=0; i<nr_eq2; i++)
556  sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
557  for (int i=nr_eq2; i<nr; i++)
558  sec_membre.set(nr+i) = 0 ;
559  Tbl big_res = oper.inverse(sec_membre) ;
560 
561  // Putting coefficients of h and v to individual arrays
562  //-----------------------------------------------------
563  for (int i=0; i<nr; i++) {
564  sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
565  sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
566  }
567 
568  // Homogeneous solutions (only 1/r^(l+2) and 1/r^l in the CED)
569  //------------------------------------------------------------
570  if (l_q == 1) {
571  Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
572  Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
573  double fac_eta = lam + 2. ;
574  double fac_vr = 2*lam + 2. ;
575  for (int i=0 ; i<nr ; i++) {
576  sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
577  sol_hom_quatre_eta.set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
578  sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
579  sol_hom_quatre_vr.set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
580  }
581  }
582  else {
583  sec_membre.annule_hard() ;
584  sec_membre.set(nr-1) = 1 ;
585  big_res = oper.inverse(sec_membre) ;
586 
587  for (int i=0; i<nr; i++) {
588  sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
589  sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
590  }
591  sec_membre.set(nr-1) = 0 ;
592  sec_membre.set(2*nr-1) = 1 ;
593  big_res = oper.inverse(sec_membre) ;
594  for (int i=0; i<nr; i++) {
595  sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
596  sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
597  }
598  }
599  }
600  }
601  }
602 
603  // Now let's match everything ...
604  //-------------------------------
605 
606  // Resulting V^r & eta
607  vr.set_etat_qcq() ;
608  vr.set_spectral_base(base) ;
609  vr.set_spectral_va().set_etat_cf_qcq() ;
610  Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
611  cf_vr.annule_hard() ;
612  het.set_etat_qcq() ;
613  het.set_spectral_base(base) ;
614  het.set_spectral_va().set_etat_cf_qcq() ;
615  Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
616  cf_eta.annule_hard() ;
617  int taille = 4*nzm1 ;
618  Tbl sec_membre(taille) ;
619  Matrice systeme(taille, taille) ;
620  systeme.set_etat_qcq() ;
621  int ligne ; int colonne ;
622 
623  // Loop on l and m
624  //----------------
625  for (int k=0 ; k<np+1 ; k++)
626  for (int j=0 ; j<nt ; j++) {
627  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
628  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
629 
630  ligne = 0 ;
631  colonne = 0 ;
632  systeme.annule_hard() ;
633  sec_membre.annule_hard() ;
634 
635  //Nucleus
636  nr = mg.get_nr(0) ;
637  alpha = mpaff->get_alpha()[0] ;
638  // value of at x=1 of eta ...
639  systeme.set(ligne, colonne) = sol_hom_un_eta.val_out_bound_jk(0, j, k) ;
640  systeme.set(ligne, colonne+1) = sol_hom_trois_eta.val_out_bound_jk(0, j, k) ;
641  sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
642  ligne++ ;
643  // ... and of its couterpart for V^r
644  systeme.set(ligne, colonne) = sol_hom_un_vr.val_out_bound_jk(0, j, k) ;
645  systeme.set(ligne, colonne+1) = sol_hom_trois_vr.val_out_bound_jk(0, j, k) ;
646  sec_membre.set(ligne) = -sol_part_vr.val_out_bound_jk(0,j,k) ;
647  ligne++ ;
648 
649  //derivatives
650  int pari = (base_r == R_CHEBP ? 0 : 1) ;
651  for (int i=0; i<nr; i++) {
652  systeme.set(ligne, colonne)
653  += (2*i+pari)*(2*i+pari)*sol_hom_un_eta(0, k, j, i)/alpha ;
654  systeme.set(ligne, colonne+1)
655  += (2*i+pari)*(2*i+pari)*sol_hom_trois_eta(0, k, j, i)/alpha ;
656  sec_membre.set(ligne)
657  -= (2*i+pari)*(2*i+pari)* sol_part_eta(0, k, j, i)/alpha ;
658  }
659  ligne++ ;
660  // ... and of its couterpart for V^r
661  for (int i=0; i<nr; i++) {
662  systeme.set(ligne, colonne)
663  += (2*i+pari)*(2*i+pari)*sol_hom_un_vr(0, k, j, i)/alpha ;
664  systeme.set(ligne, colonne+1)
665  += (2*i+pari)*(2*i+pari)*sol_hom_trois_vr(0, k, j, i)/alpha ;
666  sec_membre.set(ligne)
667  -= (2*i+pari)*(2*i+pari)* sol_part_vr(0, k, j, i)/alpha ;
668  }
669  colonne += 2 ;
670 
671  //shells
672  for (int zone=1 ; zone<nzm1 ; zone++) {
673  nr = mg.get_nr(zone) ;
674  alpha = mpaff->get_alpha()[zone] ;
675  ligne -= 3 ;
676  systeme.set(ligne, colonne)
677  = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
678  systeme.set(ligne, colonne+1)
679  = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
680  systeme.set(ligne, colonne+2)
681  = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
682  systeme.set(ligne, colonne+3)
683  = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
684  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
685  ligne++ ;
686  // ... and their couterparts for V^r
687  systeme.set(ligne, colonne)
688  = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
689  systeme.set(ligne, colonne+1)
690  = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
691  systeme.set(ligne, colonne+2)
692  = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
693  systeme.set(ligne, colonne+3)
694  = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
695  sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
696  ligne++ ;
697 
698  //derivative of (x+echelle)^(l-1) at -1
699  pari = -1 ;
700  for (int i=0; i<nr; i++) {
701  systeme.set(ligne, colonne)
702  -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
703  systeme.set(ligne, colonne+1)
704  -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
705  systeme.set(ligne, colonne+2)
706  -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
707  systeme.set(ligne, colonne+3)
708  -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
709  sec_membre.set(ligne)
710  += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
711  pari = -pari ;
712  }
713  ligne++ ;
714  // ... and their couterparts for V^r
715  pari = -1 ;
716  for (int i=0; i<nr; i++) {
717  systeme.set(ligne, colonne)
718  -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
719  systeme.set(ligne, colonne+1)
720  -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
721  systeme.set(ligne, colonne+2)
722  -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
723  systeme.set(ligne, colonne+3)
724  -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
725  sec_membre.set(ligne)
726  += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
727  pari = -pari ;
728  }
729  ligne++ ;
730 
731  systeme.set(ligne, colonne)
732  += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
733  systeme.set(ligne, colonne+1)
734  += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
735  systeme.set(ligne, colonne+2)
736  += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
737  systeme.set(ligne, colonne+3)
738  += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
739  sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
740  ligne++ ;
741  // ... and their couterparts for V^r
742  systeme.set(ligne, colonne)
743  += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
744  systeme.set(ligne, colonne+1)
745  += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
746  systeme.set(ligne, colonne+2)
747  += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
748  systeme.set(ligne, colonne+3)
749  += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
750  sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
751  ligne++ ;
752 
753  //derivative at 1
754  pari = 1 ;
755  for (int i=0; i<nr; i++) {
756  systeme.set(ligne, colonne)
757  += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
758  systeme.set(ligne, colonne+1)
759  += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
760  systeme.set(ligne, colonne+2)
761  += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
762  systeme.set(ligne, colonne+3)
763  += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
764  sec_membre.set(ligne)
765  -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
766  pari = pari ;
767  }
768  ligne++ ;
769  // ... and their couterparts for V^r
770  pari = 1 ;
771  for (int i=0; i<nr; i++) {
772  systeme.set(ligne, colonne)
773  += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
774  systeme.set(ligne, colonne+1)
775  += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
776  systeme.set(ligne, colonne+2)
777  += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
778  systeme.set(ligne, colonne+3)
779  += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
780  sec_membre.set(ligne)
781  -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
782  pari = pari ;
783  }
784  colonne += 4 ;
785  }
786  //Compactified external domain
787  nr = mg.get_nr(nzm1) ;
788 
789  alpha = mpaff->get_alpha()[nzm1] ;
790  ligne -= 3 ;
791  //value of (x-1)^(l+2) at -1 :
792  systeme.set(ligne, colonne)
793  -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
794  systeme.set(ligne, colonne+1)
795  -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
796  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
797  //... and of its couterpart for V^r
798  systeme.set(ligne+1, colonne)
799  -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
800  systeme.set(ligne+1, colonne+1)
801  -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
802  sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
803 
804  ligne += 2 ;
805  //derivative of (x-1)^(l+2) at -1 :
806  pari = 1 ;
807  for (int i=0; i<nr; i++) {
808  systeme.set(ligne, colonne)
809  -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
810  systeme.set(ligne, colonne+1)
811  -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
812  sec_membre.set(ligne)
813  += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
814  pari = -pari ;
815  }
816  ligne++ ;
817  // ... and their couterparts for V^r
818  pari = 1 ;
819  for (int i=0; i<nr; i++) {
820  systeme.set(ligne, colonne)
821  -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
822  systeme.set(ligne, colonne+1)
823  -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
824  sec_membre.set(ligne)
825  += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
826  pari = -pari ;
827  }
828 
829  // Solution of the system giving the coefficients for the homogeneous
830  // solutions
831  //-------------------------------------------------------------------
832  systeme.set_lu() ;
833  Tbl facteurs(systeme.inverse(sec_membre)) ;
834  int conte = 0 ;
835 
836  // everything is put to the right place, the same combination of hom.
837  // solutions (with some l or -(l+1) factors) must be used for V^r
838  //-------------------------------------------------------------------
839  nr = mg.get_nr(0) ; //nucleus
840  for (int i=0 ; i<nr ; i++) {
841  cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
842  +facteurs(conte)*sol_hom_un_eta(0, k, j, i)
843  +facteurs(conte+1)*sol_hom_trois_eta(0, k, j, i) ;
844  cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
845  +facteurs(conte)*sol_hom_un_vr(0, k, j, i)
846  +facteurs(conte+1)*sol_hom_trois_vr(0, k, j, i) ;
847  }
848  conte += 2 ;
849  for (int zone=1 ; zone<nzm1 ; zone++) { //shells
850  nr = mg.get_nr(zone) ;
851  for (int i=0 ; i<nr ; i++) {
852  cf_eta.set(zone, k, j, i) =
853  sol_part_eta(zone, k, j, i)
854  +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
855  +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
856  +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
857  +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
858  cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
859  +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
860  +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
861  +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
862  +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
863  }
864  conte+=4 ;
865  }
866  nr = mg.get_nr(nz-1) ; //compactified external domain
867  for (int i=0 ; i<nr ; i++) {
868  cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
869  +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
870  +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
871  cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
872  +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
873  +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
874 
875  }
876  } // End of nullite_plm
877  } //End of loop on theta
878  vr.set_spectral_va().ylm_i() ;
879  vr += vrl0 ;
880  het.set_spectral_va().ylm_i() ;
881  }
882  resu.set_vr_eta_mu(vr, het, mu().poisson()) ;
883  if ((nt==1)&&(np==1)) {
884  resu.set(2).set_etat_zero() ;
885  resu.set(3).set_etat_zero() ;
886  }
887 
888  return ;
889 
890 }
891 }
Lorene prototypes.
Definition: app_hor.h:67
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:72
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:107
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:324