LORENE
vector_poisson_boundary2.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_boundary2.C,v 1.10 2016/12/05 16:18:18 j_novak Exp $
32 * $Log: vector_poisson_boundary2.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 2008/08/20 15:07:36 n_vasset
43 * Cleaning up the code...
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_boundary2.C,v 1.10 2016/12/05 16:18:18 j_novak Exp $
63 *
64 */
65
66// C headers
67#include <cassert>
68#include <cstdlib>
69#include <cmath>
70
71// Lorene headers
72#include "metric.h"
73#include "diff.h"
74#include "param_elliptic.h"
75#include "proto.h"
76#include "graphique.h"
77#include "map.h"
78#include "utilitaires.h"
79
80
81
82namespace Lorene {
83void Vector::poisson_boundary2(double lam, Vector& resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const {
84
85 const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
86#ifndef NDEBUG
87 for (int i=0; i<3; i++)
88 assert(cmp[i]->check_dzpuis(4)) ;
89 // All this has a meaning only for spherical components:
90 const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
91 assert(bvs != 0x0) ;
92 //## ... and affine mapping, for the moment!
93 assert( mpaff != 0x0) ;
94#endif
95
96 if (fabs(lam + 1.) < 1.e-8) {
97
98 cout <<" lambda = -1 is not supported by this code!" << endl;
99 abort();
100 }
101
102
103
104 // Extraction of Mtbl_cf objects from boundary informations;
105 Scalar bound_vr2 = boundvr;
106 bound_vr2.set_spectral_va().ylm();
107 Scalar bound_eta2 = boundeta;
108 bound_eta2.set_spectral_va().ylm();
109 Scalar bound_mu2 = boundmu; bound_mu2.set_spectral_va().ylm();
110
111 const Mtbl_cf *bound_vr = bound_vr2.get_spectral_va().c_cf;
112 const Mtbl_cf *bound_mu = bound_mu2.get_spectral_va().c_cf;
113 const Mtbl_cf *bound_eta = bound_eta2.get_spectral_va().c_cf;
114
115
116
117
118 // Some working objects
119 //---------------------
120 const Mg3d& mg = *mpaff->get_mg() ;
121 int nz = mg.get_nzone() ; int nzm1 = nz - 1;
122 assert(mg.get_type_r(0) == RARE) ;
123 assert(mg.get_type_r(nzm1) == UNSURR) ;
124 Scalar S_r = *cmp[0] ;
125 Scalar S_eta = eta() ;
126 Scalar het(*mpaff) ;
127 Scalar vr(*mpaff) ;
128
129 if (S_r.get_etat() == ETATZERO) {
130 if (S_eta.get_etat() == ETATZERO) {
131
132 S_r.annule_hard() ;
133 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
134 S_eta.annule_hard() ;
135 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
136 }
137 else {
138 S_r.annule_hard() ;
139 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
140 }
141 }
142 if (S_eta.get_etat() == ETATZERO) {
143 S_eta.annule_hard() ;
144 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
145 }
146
147 S_r.set_spectral_va().ylm() ;
148 S_eta.set_spectral_va().ylm() ;
149 const Base_val& base = S_eta.get_spectral_va().base ;
150 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
151 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
152 Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
153 Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
154 Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
155 Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
156 Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
157 Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
158 Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
159 Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
160
161 // Solution of the l=0 part for Vr
162 //---------------------------------
163 Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
164 sou_l0.set_spectral_va().ylm() ;
165
166 Param_elliptic param_l0(sou_l0) ;
167 for (int l=0; l<nz; l++){
168 param_l0.set_poisson_vect_r(l, true) ;
169 }
170
171
172 Scalar vrl0 = sou_l0.sol_elliptic_boundary(param_l0, *bound_vr, dir_vr, neum_vr) ;
173
174
175 // Build-up & inversion of the system for (eta, V^r) in each domain (except for the nucleus!)
176 //-----------------------------------------------------------------
177
178
179 // Shells
180 //-------
181
182 int nr = mg.get_nr(1) ;
183 int nt = mg.get_nt(1) ;
184 int np = mg.get_np(1) ;
185 double alpha = mpaff->get_alpha()[1] ; double alp2 = alpha*alpha ;
186 double beta = mpaff->get_beta()[1] ;
187
188 int l_q = 0 ; int m_q = 0; int base_r = 0 ;
189 for (int zone=1 ; zone<nzm1 ; zone++) {
190 nr = mg.get_nr(zone) ;
191 assert (nr > 5) ;
192 assert(nt == mg.get_nt(zone)) ;
193 assert(np == mg.get_np(zone)) ;
194 alpha = mpaff->get_alpha()[zone] ;
195 beta = mpaff->get_beta()[zone] ;
196 double ech = beta / alpha ;
197
198 // Loop on l and m
199 //----------------
200 for (int k=0 ; k<np+1 ; k++) {
201 for (int j=0 ; j<nt ; j++) {
202 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
203 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
204 int dege1 = 2 ; //degeneracy of eq.1
205 int dege2 = 2 ; //degeneracy of eq.2
206 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
207 int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
208 int nrtot = 2*nr ;
209 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
210 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
211 Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
212 Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
213 Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
214 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
215 Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
216 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
217
218 // Building the operator
219 //----------------------
220 for (int lin=0; lin<nr_eq1; lin++) {
221 for (int col=0; col<nr; col++)
222 oper.set(lin,col)
223 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
224 + 2*(mxd(lin,col) + ech*md(lin,col))
225 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
226 for (int col=0; col<nr; col++)
227 oper.set(lin,col+nr)
228 = lam*(mxd(lin,col) + ech*md(lin,col))
229 + 2*(1+lam)*mid(lin,col) ;
230 }
231 oper.set(nr_eq1, 0) = 1 ;
232 for (int col=1; col<2*nr; col++)
233 oper.set(nr_eq1, col) = 0 ;
234 oper.set(nr_eq1+1, 0) = 0 ;
235 oper.set(nr_eq1+1, 1) = 1 ;
236 for (int col=2; col<2*nr; col++)
237 oper.set(nr_eq1+1, col) = 0 ;
238 for (int lin=0; lin<nr_eq2; lin++) {
239 for (int col=0; col<nr; col++)
240 oper.set(lin+nr,col)
241 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
242 - (lam+2)*mid(lin,col)) ;
243 for (int col=0; col<nr; col++)
244 oper.set(lin+nr, col+nr)
245 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
246 + ech*ech*md2(lin,col)
247 + 2*(mxd(lin,col) + ech*md(lin,col)))
248 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
249 }
250 for (int col=0; col<nr; col++)
251 oper.set(nr+nr_eq2, col) = 0 ;
252 oper.set(nr+nr_eq2, nr) = 1 ;
253 for (int col=nr+1; col<2*nr; col++)
254 oper.set(nr+nr_eq2, col) = 0 ;
255 for (int col=0; col<nr+1; col++)
256 oper.set(nr+nr_eq2+1, col) = 0 ;
257 oper.set(nr+nr_eq2+1, nr+1) = 1 ;
258 for (int col=nr+2; col<2*nr; col++)
259 oper.set(nr+nr_eq2+1, col) = 0 ;
260
261 oper.set_lu() ;
262
263 // Filling the r.h.s
264 //------------------
265 Tbl sr(nr) ; sr.set_etat_qcq() ;
266 Tbl seta(nr) ; seta.set_etat_qcq() ;
267 for (int i=0; i<nr; i++) {
268 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
269 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
270 }
271 Tbl xsr= sr ; Tbl x2sr= sr ;
272 Tbl xseta= seta ; Tbl x2seta = seta ;
273 multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
274 multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
275
276 for (int i=0; i<nr_eq1; i++)
277 sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
278 + beta*beta*seta(i);
279 sec_membre.set(nr_eq1) = 0 ;
280 sec_membre.set(nr_eq1+1) = 0 ;
281 for (int i=0; i<nr_eq2; i++)
282 sec_membre.set(i+nr) = beta*beta*sr(i)
283 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
284 sec_membre.set(nr+nr_eq2) = 0 ;
285 sec_membre.set(nr+nr_eq2+1) = 0 ;
286
287 // Inversion of the "big" operator
288 //--------------------------------
289 Tbl big_res = oper.inverse(sec_membre) ;
290 for (int i=0; i<nr; i++) {
291 sol_part_eta.set(zone, k, j, i) = big_res(i) ;
292 sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
293 }
294
295 // Getting the homogeneous solutions
296 //----------------------------------
297 sec_membre.annule_hard() ;
298 sec_membre.set(nr_eq1) = 1 ;
299 big_res = oper.inverse(sec_membre) ;
300 for (int i=0 ; i<nr ; i++) {
301 sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
302 sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
303 }
304 sec_membre.set(nr_eq1) = 0 ;
305 sec_membre.set(nr_eq1+1) = 1 ;
306 big_res = oper.inverse(sec_membre) ;
307 for (int i=0 ; i<nr ; i++) {
308 sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
309 sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
310 }
311 sec_membre.set(nr_eq1+1) = 0 ;
312 sec_membre.set(nr+nr_eq2) = 1 ;
313 big_res = oper.inverse(sec_membre) ;
314 for (int i=0 ; i<nr ; i++) {
315 sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
316 sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
317 }
318 sec_membre.set(nr+nr_eq2) = 0 ;
319 sec_membre.set(nr+nr_eq2+1) = 1 ;
320 big_res = oper.inverse(sec_membre) ;
321 for (int i=0 ; i<nr ; i++) {
322 sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
323 sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
324 }
325
326 }
327 }
328 }
329 }
330
331 // Compactified external domain
332 //-----------------------------
333 nr = mg.get_nr(nzm1) ;
334 assert(nt == mg.get_nt(nzm1)) ;
335 assert(np == mg.get_np(nzm1)) ;
336 alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
337 assert (nr > 4) ;
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(nzm1, k, j, m_q, l_q, base_r) ;
344 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
345 int dege1 = 3; //degeneracy of eq.1
346 int dege2 = (l_q == 1 ? 2 : 3); //degeneracy of eq.2
347 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
348 int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
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_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
353 Diff_sxdsdx sxd(base_r, nr) ; const Matrice& mxd = sxd.get_matrice() ;
354 Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.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 = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
362 for (int col=0; col<nr; col++)
363 oper.set(lin,col+nr) =
364 (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
365 }
366 for (int col=0; col<nr; col++)
367 oper.set(nr_eq1, col) = 1 ;
368 for (int col=nr; col<nrtot; col++)
369 oper.set(nr_eq1, col) = 0 ;
370 int pari = -1 ;
371 for (int col=0; col<nr; col++) {
372 oper.set(nr_eq1+1, col) = pari*col*col ;
373 pari = pari ;
374 }
375 for (int col=nr; col<nrtot; col++)
376 oper.set(nr_eq1+1, col) = 0 ;
377 oper.set(nr_eq1+2, 0) = 1 ;
378 for (int col=1; col<nrtot; col++)
379 oper.set(nr_eq1+2, col) = 0 ;
380 for (int lin=0; lin<nr_eq2; lin++) {
381 for (int col=0; col<nr; col++)
382 oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
383 + (lam+2)*ms2(lin,col))) / alp2 ;
384 for (int col=0; col<nr; col++)
385 oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
386 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
387 }
388 for (int col=0; col<nr; col++)
389 oper.set(nr+nr_eq2, col) = 0 ;
390 for (int col=nr; col<nrtot; col++)
391 oper.set(nr+nr_eq2, col) = 1 ;
392 for (int col=0; col<nr; col++)
393 oper.set(nr+nr_eq2+1, col) = 0 ;
394 pari = -1 ;
395 for (int col=0; col<nr; col++) {
396 oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
397 pari = pari ;
398 }
399 if (l_q > 1) {
400 for (int col=0; col<nr; col++)
401 oper.set(nr+nr_eq2+2, col) = 0 ;
402 oper.set(nr+nr_eq2+2, nr) = 1 ;
403 for (int col=nr+1; col<nrtot; col++)
404 oper.set(nr+nr_eq2+2, col) = 0 ;
405 }
406 oper.set_lu() ;
407
408 // Filling the r.h.s
409 //------------------
410 for (int i=0; i<nr_eq1; i++)
411 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
412 for (int i=nr_eq1; i<nr; i++)
413 sec_membre.set(i) = 0 ;
414 for (int i=0; i<nr_eq2; i++)
415 sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
416 for (int i=nr_eq2; i<nr; i++)
417 sec_membre.set(nr+i) = 0 ;
418 Tbl big_res = oper.inverse(sec_membre) ;
419
420 // Putting coefficients of h and v to individual arrays
421 //-----------------------------------------------------
422 for (int i=0; i<nr; i++) {
423 sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
424 sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
425 }
426
427 // Homogeneous solutions (only 1/r^(l+2) and 1/r^l in the CED)
428 //------------------------------------------------------------
429 if (l_q == 1) {
430 Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
431 Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
432 double fac_eta = lam + 2. ;
433 double fac_vr = 2*lam + 2. ;
434 for (int i=0 ; i<nr ; i++) {
435 sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
436 sol_hom_quatre_eta.set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
437 sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
438 sol_hom_quatre_vr.set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
439 }
440 }
441 else {
442 sec_membre.annule_hard() ;
443 sec_membre.set(nr-1) = 1 ;
444 big_res = oper.inverse(sec_membre) ;
445
446 for (int i=0; i<nr; i++) {
447 sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
448 sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
449 }
450 sec_membre.set(nr-1) = 0 ;
451 sec_membre.set(2*nr-1) = 1 ;
452 big_res = oper.inverse(sec_membre) ;
453 for (int i=0; i<nr; i++) {
454 sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
455 sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
456 }
457 }
458 }
459 }
460 }
461
462 // Now let's match everything ...
463 //-------------------------------
464
465 // Resulting V^r & eta
466 vr.set_etat_qcq() ;
467 vr.set_spectral_base(base) ;
468 vr.set_spectral_va().set_etat_cf_qcq() ;
469 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
470 cf_vr.annule_hard() ;
471 het.set_etat_qcq() ;
472 het.set_spectral_base(base) ;
473 het.set_spectral_va().set_etat_cf_qcq() ;
474 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
475 cf_eta.annule_hard() ;
476 int taille = 4*nzm1 -2 ; // Two less than R3 case
477 Tbl sec_membre(taille) ;
478 Matrice systeme(taille, taille) ;
479 systeme.set_etat_qcq() ;
480 int ligne ; int colonne ;
481
482 // Loop on l and m
483 //----------------
484 for (int k=0 ; k<np+1 ; k++)
485 for (int j=0 ; j<nt ; j++) {
486 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
487 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
488
489 ligne = 0 ;
490 colonne = 0 ;
491 systeme.annule_hard() ;
492 sec_membre.annule_hard() ;
493
494
495 //shell 1 (boundary condition)
496 int zone1 = 1;
497 nr = mg.get_nr(zone1) ;
498 alpha = mpaff->get_alpha()[zone1] ;
499
500 // Here we prepare for a Robyn-type boundary condition on eta.
501 // Parameters are dir_eta and neum_eta
502 systeme.set(ligne, colonne)
503 = -dir_eta*sol_hom_un_eta.val_in_bound_jk(zone1, j, k) ;
504 systeme.set(ligne, colonne+1)
505 = -dir_eta * sol_hom_deux_eta.val_in_bound_jk(zone1, j, k) ;
506 systeme.set(ligne, colonne+2)
507 = -dir_eta*sol_hom_trois_eta.val_in_bound_jk(zone1, j, k) ;
508 systeme.set(ligne, colonne+3)
509 = -dir_eta*sol_hom_quatre_eta.val_in_bound_jk(zone1, j, k) ;
510
511
512
513
514 sec_membre.set(ligne) += dir_eta* sol_part_eta.val_in_bound_jk(zone1, j, k) ;
515
516
517
518 int pari = -1 ;
519 for (int i=0; i<nr; i++) {
520 systeme.set(ligne, colonne)
521 -= neum_eta*pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
522 systeme.set(ligne, colonne+1)
523 -= neum_eta*pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
524 systeme.set(ligne, colonne+2)
525 -= neum_eta*pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
526 systeme.set(ligne, colonne+3)
527 -= neum_eta*pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
528 sec_membre.set(ligne)
529 += neum_eta*pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
530 pari = -pari ;
531 }
532
533 sec_membre.set(ligne) -= (*bound_eta).val_in_bound_jk(1,j,k);
534 ligne++ ;
535
536
537
538 // ... and their couterparts for V^r
539 systeme.set(ligne, colonne)
540 = - dir_vr*sol_hom_un_vr.val_in_bound_jk(zone1, j, k) ;
541 systeme.set(ligne, colonne+1)
542 = - dir_vr*sol_hom_deux_vr.val_in_bound_jk(zone1, j, k) ;
543 systeme.set(ligne, colonne+2)
544 = -dir_vr*sol_hom_trois_vr.val_in_bound_jk(zone1, j, k) ;
545 systeme.set(ligne, colonne+3)
546 = -dir_vr*sol_hom_quatre_vr.val_in_bound_jk(zone1, j, k) ;
547 sec_membre.set(ligne) += dir_vr*sol_part_vr.val_in_bound_jk(zone1, j, k) ;
548
549
550
551
552 pari = -1 ;
553 for (int i=0; i<nr; i++) {
554 systeme.set(ligne, colonne)
555 -= neum_vr*pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
556 systeme.set(ligne, colonne+1)
557 -= neum_vr*pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
558 systeme.set(ligne, colonne+2)
559 -= neum_vr*pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
560 systeme.set(ligne, colonne+3)
561 -= neum_vr*pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
562 sec_membre.set(ligne)
563 += neum_vr*pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
564 pari = -pari ;
565 }
566
567 sec_membre.set(ligne) -= (*bound_vr).val_in_bound_jk(1,j,k);
568 ligne++ ;
569
570
571
572
573 // Values in 1: beginning with solution in eta...
574
575 systeme.set(ligne, colonne)
576 += sol_hom_un_eta.val_out_bound_jk(zone1, j, k) ;
577 systeme.set(ligne, colonne+1)
578 += sol_hom_deux_eta.val_out_bound_jk(zone1, j, k) ;
579 systeme.set(ligne, colonne+2)
580 += sol_hom_trois_eta.val_out_bound_jk(zone1, j, k) ;
581 systeme.set(ligne, colonne+3)
582 += sol_hom_quatre_eta.val_out_bound_jk(zone1, j, k) ;
583 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone1, j, k) ;
584 ligne++ ;
585 // ... and their counterparts for V^r
586 systeme.set(ligne, colonne)
587 += sol_hom_un_vr.val_out_bound_jk(zone1, j, k) ;
588 systeme.set(ligne, colonne+1)
589 += sol_hom_deux_vr.val_out_bound_jk(zone1, j, k) ;
590 systeme.set(ligne, colonne+2)
591 += sol_hom_trois_vr.val_out_bound_jk(zone1, j, k) ;
592 systeme.set(ligne, colonne+3)
593 += sol_hom_quatre_vr.val_out_bound_jk(zone1, j, k) ;
594 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone1, j, k) ;
595 ligne++ ;
596
597 //derivative at 1
598 pari = 1 ;
599 for (int i=0; i<nr; i++) {
600 systeme.set(ligne, colonne)
601 += pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
602 systeme.set(ligne, colonne+1)
603 += pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
604 systeme.set(ligne, colonne+2)
605 += pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
606 systeme.set(ligne, colonne+3)
607 += pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
608 sec_membre.set(ligne)
609 -= pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
610 pari = pari ;
611 }
612 ligne++ ;
613 // ... and their counterparts for V^r
614 pari = 1 ;
615 for (int i=0; i<nr; i++) {
616 systeme.set(ligne, colonne)
617 += pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
618 systeme.set(ligne, colonne+1)
619 += pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
620 systeme.set(ligne, colonne+2)
621 += pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
622 systeme.set(ligne, colonne+3)
623 += pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
624 sec_membre.set(ligne)
625 -= pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
626 pari = pari ;
627 }
628 colonne += 4 ;
629
630 // Other shells
631 if ( nz <= 2){
632 cout <<"WARNING!! Mapping must contain at least 2 shells!!" << endl;}
633 for (int zone=2 ; zone<nzm1 ; zone++) {
634
635 nr = mg.get_nr(zone) ;
636 alpha = mpaff->get_alpha()[zone] ;
637 ligne -= 3 ;
638 systeme.set(ligne, colonne)
639 = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
640 systeme.set(ligne, colonne+1)
641 = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
642 systeme.set(ligne, colonne+2)
643 = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
644 systeme.set(ligne, colonne+3)
645 = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
646 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
647 ligne++ ;
648 // ... and their counterparts for V^r
649 systeme.set(ligne, colonne)
650 = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
651 systeme.set(ligne, colonne+1)
652 = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
653 systeme.set(ligne, colonne+2)
654 = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
655 systeme.set(ligne, colonne+3)
656 = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
657 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
658 ligne++ ;
659
660 //derivative of (x+echelle)^(l-1) at -1
661 pari = -1 ;
662 for (int i=0; i<nr; i++) {
663 systeme.set(ligne, colonne)
664 -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
665 systeme.set(ligne, colonne+1)
666 -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
667 systeme.set(ligne, colonne+2)
668 -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
669 systeme.set(ligne, colonne+3)
670 -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
671 sec_membre.set(ligne)
672 += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
673 pari = -pari ;
674 }
675 ligne++ ;
676 // ... and their counterparts for V^r
677 pari = -1 ;
678 for (int i=0; i<nr; i++) {
679 systeme.set(ligne, colonne)
680 -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
681 systeme.set(ligne, colonne+1)
682 -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
683 systeme.set(ligne, colonne+2)
684 -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
685 systeme.set(ligne, colonne+3)
686 -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
687 sec_membre.set(ligne)
688 += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
689 pari = -pari ;
690 }
691 ligne++ ;
692
693 systeme.set(ligne, colonne)
694 += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
695 systeme.set(ligne, colonne+1)
696 += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
697 systeme.set(ligne, colonne+2)
698 += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
699 systeme.set(ligne, colonne+3)
700 += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
701 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
702 ligne++ ;
703 // ... and their counterparts for V^r
704 systeme.set(ligne, colonne)
705 += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
706 systeme.set(ligne, colonne+1)
707 += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
708 systeme.set(ligne, colonne+2)
709 += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
710 systeme.set(ligne, colonne+3)
711 += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
712 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
713 ligne++ ;
714
715 //derivative at 1
716 pari = 1 ;
717 for (int i=0; i<nr; i++) {
718 systeme.set(ligne, colonne)
719 += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
720 systeme.set(ligne, colonne+1)
721 += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
722 systeme.set(ligne, colonne+2)
723 += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
724 systeme.set(ligne, colonne+3)
725 += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
726 sec_membre.set(ligne)
727 -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
728 pari = pari ;
729 }
730 ligne++ ;
731 // ... and their couterparts for V^r
732 pari = 1 ;
733 for (int i=0; i<nr; i++) {
734 systeme.set(ligne, colonne)
735 += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
736 systeme.set(ligne, colonne+1)
737 += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
738 systeme.set(ligne, colonne+2)
739 += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
740 systeme.set(ligne, colonne+3)
741 += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
742 sec_membre.set(ligne)
743 -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
744 pari = pari ;
745 }
746 colonne += 4 ;
747 }
748
749 //Compactified external domain
750 nr = mg.get_nr(nzm1) ;
751
752 alpha = mpaff->get_alpha()[nzm1] ;
753 ligne -= 3 ;
754 //value of (x-1)^(l+2) at -1 :
755 systeme.set(ligne, colonne)
756 -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
757 systeme.set(ligne, colonne+1)
758 -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
759 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
760 //... and of its counterpart for V^r
761 systeme.set(ligne+1, colonne)
762 -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
763 systeme.set(ligne+1, colonne+1)
764 -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
765 sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
766
767 ligne += 2 ;
768 //derivative of (x-1)^(l+2) at -1 :
769 pari = 1 ;
770 for (int i=0; i<nr; i++) {
771 systeme.set(ligne, colonne)
772 -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
773 systeme.set(ligne, colonne+1)
774 -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
775 sec_membre.set(ligne)
776 += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
777 pari = -pari ;
778 }
779 ligne++ ;
780 // ... and their counterparts for V^r
781 pari = 1 ;
782 for (int i=0; i<nr; i++) {
783 systeme.set(ligne, colonne)
784 -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
785 systeme.set(ligne, colonne+1)
786 -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
787 sec_membre.set(ligne)
788 += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
789 pari = -pari ;
790 }
791
792 // Solution of the system giving the coefficients for the homogeneous
793 // solutions
794 //-------------------------------------------------------------------
795 systeme.set_lu() ;
796 Tbl facteurs(systeme.inverse(sec_membre)) ;
797 int conte = 0 ;
798
799 // everything is put to the right place, the same combination of hom.
800 // solutions (with some l or -(l+1) factors) must be used for V^r
801 //-------------------------------------------------------------------
802 nr = mg.get_nr(0) ; //nucleus
803 for (int i=0 ; i<nr ; i++) {
804 cf_eta.set(0, k, j, i) = 0.;
805 cf_vr.set(0, k, j, i) = 0.;
806 }
807
808 for (int zone=1 ; zone<nzm1 ; zone++) { //shells
809 nr = mg.get_nr(zone) ;
810 for (int i=0 ; i<nr ; i++) {
811 cf_eta.set(zone, k, j, i) =
812 sol_part_eta(zone, k, j, i)
813 +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
814 +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
815 +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
816 +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
817 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
818 +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
819 +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
820 +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
821 +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
822 }
823 conte+=4 ;
824 }
825 nr = mg.get_nr(nz-1) ; //compactified external domain
826 for (int i=0 ; i<nr ; i++) {
827 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
828 +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
829 +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
830 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
831 +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
832 +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
833
834 }
835 } // End of nullite_plm
836 } //End of loop on theta
837 vr.set_spectral_va().ylm_i() ;
838 vr += vrl0 ;
839 het.set_spectral_va().ylm_i() ;
840
841
842
843 Scalar pipo(*mpaff);
844 pipo.annule_hard();
845 pipo.std_spectral_base();
846 pipo.set_dzpuis(3);
847 Param_elliptic param_mu(mu());
848
849
850
851
852 resu.set_vr_eta_mu(vr, het, mu().sol_elliptic_boundary(param_mu, (*bound_mu), dir_mu , neum_mu)) ;
853
854 return ;
855
856}
857}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:325
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition base_val_quantum.C:84
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Diff_dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:172
Lorene::Diff_dsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_dsdx2.C:94
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Lorene::Diff_dsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_dsdx.C:97
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Lorene::Diff_sx2
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:369
Lorene::Diff_sx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sx2.C:101
Lorene::Diff_sxdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:450
Lorene::Diff_sxdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sxdsdx.C:100
Lorene::Diff_x2dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:490
Lorene::Diff_x2dsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_x2dsdx2.C:98
Lorene::Diff_xdsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:531
Lorene::Diff_xdsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_xdsdx2.C:100
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Lorene::Diff_xdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_xdsdx.C:101
Lorene::Map_af
Affine radial mapping.
Definition map.h:2042
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:608
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:604
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:178
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:427
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition matrice.C:196
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:395
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:479
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:474
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:469
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:491
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:304
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:315
Lorene::Mtbl_cf::val_out_bound_jk
double val_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Definition mtbl_cf_val_point.C:431
Lorene::Mtbl_cf::val_in_bound_jk
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Definition mtbl_cf_val_point.C:455
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition param_elliptic.h:135
Lorene::Param_elliptic::set_poisson_vect_r
void set_poisson_vect_r(int zone, bool only_l_zero=false)
Sets the operator to in all domains, for ; and to in all domains otherwise.
Definition param_elliptic.C:493
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:359
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:790
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:610
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:607
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:386
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:560
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1328
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:814
Lorene::Scalar::sol_elliptic_boundary
Scalar sol_elliptic_boundary(Param_elliptic &params, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
Definition scalar_pde.C:259
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:803
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:715
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:134
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:315
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::Vector
Vector(const Map &map, int tipe, const Base_vect &triad_i)
Standard constructor.
Definition vector.C:159
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
Lorene::Vector::poisson_boundary2
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
Definition vector_poisson_boundary2.C:83
Lorene::Vector::set_vr_eta_mu
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through potentials and (see members p_eta and p_mu ), as well as the compon...
Definition vector_etamu.C:192
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