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
77namespace Lorene {
78void 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::Vector
Tensor field of valence 1.
Definition vector.h:188
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:521
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:101
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:69
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:301
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:321
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:309
Lorene
Lorene prototypes.
Definition app_hor.h:67