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
28char 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
74namespace Lorene {
75void 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}
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64