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