sym_tensor_trans_dirac_bound2.C
/*
* Methods to impose the Dirac gauge: divergence-free condition, with interior boundary conditions added
*
* (see file sym_tensor.h for documentation).
*
*/
/*
* Copyright (c) 2006, 2007 Jerome Novak,Nicolas Vasset
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2
* as published by the Free Software Foundation.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
// char sym_tensor_trans_dirac_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $" ;
/*
* $Id: sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $
* $Log: sym_tensor_trans_dirac_bound2.C,v $
* Revision 1.5 2015/08/10 15:32:26 j_novak
* Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
*
* Revision 1.4 2014/10/13 08:53:43 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2014/10/06 15:13:19 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.2 2008/08/20 15:08:06 n_vasset
* Cleaning up the code...
*
* Revision 1.1 2007/05/04 16:45:25 n_vasset
* First version
*
* Revision 1.2 2006/10/24 13:03:19 j_novak
* New methods for the solution of the tensor wave equation. Perhaps, first
* operational version...
*
* Revision 1.1 2006/09/05 15:38:45 j_novak
* The fuctions sol_Dirac... are in a separate file, with new parameters to
* control the boundary conditions.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $
*
*/
// C headers
#include <cassert>
#include <cmath>
// Lorene headers
#include "param.h"
#include "param_elliptic.h"
#include "metric.h"
#include "diff.h"
#include "proto.h"
#include "param.h"
//----------------------------------------------------------------------------------
//
// sol_Dirac_A
//
//----------------------------------------------------------------------------------
namespace Lorene {
void Sym_tensor_trans::sol_Dirac_Abound(const Scalar& aaa, Scalar& tilde_mu, Scalar& x_new, Scalar bound_mu, const Param* par_bc) {
const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
assert(mp_aff != 0x0) ; //Only affine mapping for the moment
// WARNING!
// This will only work if we have at least 2 shells!
const Mg3d& mgrid = *mp_aff->get_mg() ;
assert(mgrid.get_type_r(0) == RARE) ;
if (aaa.get_etat() == ETATZERO) {
tilde_mu = 0 ;
x_new = 0 ;
return ;
}
double dir = 1.;
double neum = 0.;
// On suppose que le sym_tensor_entré est bien transverse;
// assert ( maxabs(contract((*this).derive_cov(mets), 1, 2)) < 0.00000000000001 ) ;
bound_mu.set_spectral_va().ylm();
assert(aaa.get_etat() != ETATNONDEF) ;
int nz = mgrid.get_nzone() ;
int nzm1 = nz - 1 ;
bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
int n_shell = ced ? nz-2 : nzm1 ;
int nz_bc = nzm1 ;
if (par_bc != 0x0)
if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
bool cedbc = (ced && (nz_bc == nzm1)) ;
#ifndef NDEBUG
if (!cedbc) {
assert(par_bc != 0x0) ;
assert(par_bc->get_n_tbl_mod() >= 3) ;
}
#endif
int nt = mgrid.get_nt(0) ;
int np = mgrid.get_np(0) ;
Scalar source = aaa ;
Scalar source_coq = aaa ;
source_coq.annule_domain(0) ;
if (ced) source_coq.annule_domain(nzm1) ;
source_coq.mult_r() ;
source.set_spectral_va().ylm() ;
source_coq.set_spectral_va().ylm() ;
Base_val base = source.get_spectral_base() ;
base.mult_x() ;
tilde_mu.annule_hard() ;
tilde_mu.set_spectral_base(base) ;
tilde_mu.set_spectral_va().set_etat_cf_qcq() ;
tilde_mu.set_spectral_va().c_cf->annule_hard() ;
x_new.annule_hard() ;
x_new.set_spectral_base(base) ;
x_new.set_spectral_va().set_etat_cf_qcq() ;
x_new.set_spectral_va().c_cf->annule_hard() ;
Mtbl_cf sol_part_mu(mgrid, base) ; sol_part_mu.annule_hard() ;
Mtbl_cf sol_part_x(mgrid, base) ; sol_part_x.annule_hard() ;
Mtbl_cf sol_hom1_mu(mgrid, base) ; sol_hom1_mu.annule_hard() ;
Mtbl_cf sol_hom1_x(mgrid, base) ; sol_hom1_x.annule_hard() ;
Mtbl_cf sol_hom2_mu(mgrid, base) ; sol_hom2_mu.annule_hard() ;
Mtbl_cf sol_hom2_x(mgrid, base) ; sol_hom2_x.annule_hard() ;
int l_q, m_q, base_r ;
//---------------
//-- NUCLEUS ---
//---------------
// On va annuler toutes les solutions dans le noyau
{ int lz = 0 ;
int nr = mgrid.get_nr(lz) ;
// double alpha = mp_aff->get_alpha()[lz] ;
// Matrice ope(2*nr, 2*nr) ;
// ope.set_etat_qcq() ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
Tbl sol(2*nr) ;
sol.set_etat_zero();
for (int i=0; i<nr; i++) {
sol_part_mu.set(lz, k, j, i) = 0 ;
sol_part_x.set(lz, k, j, i) = 0 ;
}
for (int i=0; i<nr; i++) {
sol_hom2_mu.set(lz, k, j, i) = 0 ;
sol_hom2_x.set(lz, k, j, i) = 0 ;
}
}
}
}
}
//-------------
// -- Shells --
//-------------
for (int lz=1; lz <= n_shell; lz++) {
int nr = mgrid.get_nr(lz) ;
assert(mgrid.get_nt(lz) == nt) ;
assert(mgrid.get_np(lz) == np) ;
double alpha = mp_aff->get_alpha()[lz] ;
double ech = mp_aff->get_beta()[lz] / alpha ;
Matrice ope(2*nr, 2*nr) ;
ope.set_etat_qcq() ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
+ 3*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = (2-l_q*(l_q+1))*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) ;
int ind0 = 0 ;
int ind1 = nr ;
for (int col=0; col<2*nr; col++) {
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
}
ope.set(ind0+nr-1, ind0) = 1 ;
ope.set(ind1+nr-1, ind1) = 1 ;
ope.set_lu() ;
Tbl sec(2*nr) ;
sec.set_etat_qcq() ;
for (int lin=0; lin<nr; lin++)
sec.set(lin) = 0 ;
for (int lin=0; lin<nr; lin++)
sec.set(nr+lin) = (*source_coq.get_spectral_va().c_cf)
(lz, k, j, lin) ;
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 0 ;
Tbl sol = ope.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_mu.set(lz, k, j, i) = sol(i) ;
sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
}
sec.annule_hard() ;
sec.set(ind0+nr-1) = 1 ;
sol = ope.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
}
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 1 ;
sol = ope.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom2_mu.set(lz, k, j, i) = sol(i) ;
sol_hom2_x.set(lz, k, j, i) = sol(i+nr) ;
}
}
}
}
}
//------------------------------
// Compactified external domain
//------------------------------
if (cedbc) {int lz = nzm1 ;
int nr = mgrid.get_nr(lz) ;
assert(mgrid.get_nt(lz) == nt) ;
assert(mgrid.get_np(lz) == np) ;
double alpha = mp_aff->get_alpha()[lz] ;
Matrice ope(2*nr, 2*nr) ;
ope.set_etat_qcq() ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = (2-l_q*(l_q+1))*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = -md(lin,col) ;
ope *= 1./alpha ;
int ind0 = 0 ;
int ind1 = nr ;
for (int col=0; col<2*nr; col++) {
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-2, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
}
for (int col=0; col<nr; col++) {
ope.set(ind0+nr-1, col+ind0) = 1 ;
ope.set(ind1+nr-1, col+ind1) = 1 ;
}
ope.set(ind1+nr-2, ind1+1) = 1 ;
ope.set_lu() ;
Tbl sec(2*nr) ;
sec.set_etat_qcq() ;
for (int lin=0; lin<nr; lin++)
sec.set(lin) = 0 ;
for (int lin=0; lin<nr; lin++)
sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-2) = 0 ;
sec.set(ind1+nr-1) = 0 ;
Tbl sol = ope.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_mu.set(lz, k, j, i) = sol(i) ;
sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
}
sec.annule_hard() ;
sec.set(ind1+nr-2) = 1 ;
sol = ope.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
}
}
}
}
}
// On écrit maintenant le systčme de raccord
int taille = 2*nz_bc ;
if (cedbc) taille-- ;
Mtbl_cf& mmu = *tilde_mu.set_spectral_va().c_cf ;
Mtbl_cf& mw = *x_new.set_spectral_va().c_cf ;
Tbl sec_membre(taille) ;
Matrice systeme(taille, taille) ;
int ligne ; int colonne ;
Tbl pipo(1) ;
const Tbl& mub = (cedbc ? pipo : par_bc->get_tbl_mod(2) );
double c_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(0) ) ;
double d_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(1) ) ;
double c_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(2) ) ;
double d_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(3) ) ;
Mtbl_cf dhom1_mu = sol_hom1_mu ;
Mtbl_cf dhom2_mu = sol_hom2_mu ;
Mtbl_cf dpart_mu = sol_part_mu ;
Mtbl_cf dhom1_x = sol_hom1_x ;
Mtbl_cf dhom2_x = sol_hom2_x ;
Mtbl_cf dpart_x = sol_part_x ;
dhom1_mu.dsdx() ;
dhom2_mu.dsdx() ;
dpart_mu.dsdx() ;
dhom1_x.dsdx() ;
dhom2_x.dsdx() ;
dpart_x.dsdx() ;
// Loop on l and m
//----------------
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
ligne = 0 ;
colonne = 0 ;
systeme.annule_hard() ;
sec_membre.annule_hard() ;
//First shell
// Internal boundary condition
int nr = mgrid.get_nr(1) ;
systeme.set(ligne, colonne) = dir * sol_hom1_mu.val_in_bound_jk(1, j, k) + neum* dhom1_mu.val_in_bound_jk(1,j,k);
systeme.set(ligne, colonne+1) = dir * sol_hom2_mu.val_in_bound_jk(1, j, k) + neum* dhom2_mu.val_in_bound_jk(1,j,k);
Mtbl_cf *boundmu = bound_mu.get_spectral_va().c_cf;
double blob = (*boundmu).val_in_bound_jk(1,j,k);
sec_membre.set(ligne) = - dir* sol_part_mu.val_in_bound_jk(1, j, k) - neum * dpart_mu.val_in_bound_jk(1,j,k) + blob;
ligne++;
// Condition at x=1
systeme.set(ligne, colonne) =
sol_hom1_mu.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_mu.val_out_bound_jk(1, j, k) ;
sec_membre.set(ligne) -= sol_part_mu.val_out_bound_jk(1, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_x.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_x.val_out_bound_jk(1, j, k) ;
sec_membre.set(ligne) -= sol_part_x.val_out_bound_jk(1, j, k) ;
colonne += 2 ;
//shells with nz>2
for (int zone=2 ; zone<nz_bc ; zone++) {
nr = mgrid.get_nr(zone) ;
ligne-- ;
//Condition at x = -1
systeme.set(ligne, colonne) =
- sol_hom1_mu.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_mu.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_x.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_x.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(zone, j, k) ;
ligne++ ;
// Condition at x=1
systeme.set(ligne, colonne) =
sol_hom1_mu.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_mu.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_mu.val_out_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_x.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_x.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_x.val_out_bound_jk(zone, j, k) ;
colonne += 2 ;
}
//Last domain
nr = mgrid.get_nr(nz_bc) ;
double alpha = mp_aff->get_alpha()[nz_bc] ;
ligne-- ;
//Condition at x = -1
systeme.set(ligne, colonne) =
- sol_hom1_mu.val_in_bound_jk(nz_bc, j, k) ;
if (!cedbc) systeme.set(ligne, colonne+1) =
- sol_hom2_mu.val_in_bound_jk(nz_bc, j, k) ;
sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(nz_bc, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_x.val_in_bound_jk(nz_bc, j, k) ;
if (!cedbc) systeme.set(ligne, colonne+1) =
- sol_hom2_x.val_in_bound_jk(nz_bc, j, k) ;
sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(nz_bc, j, k) ;
ligne++ ;
if (!cedbc) {// Special condition at x=1
systeme.set(ligne, colonne) =
c_mu*sol_hom1_mu.val_out_bound_jk(nz_bc, j, k)
+ d_mu*dhom1_mu.val_out_bound_jk(nz_bc, j, k) / alpha
+ c_x*sol_hom1_x.val_out_bound_jk(nz_bc, j, k)
+ d_x*dhom1_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
systeme.set(ligne, colonne+1) =
c_mu*sol_hom2_mu.val_out_bound_jk(nz_bc, j, k)
+ d_mu*dhom2_mu.val_out_bound_jk(nz_bc, j, k) / alpha
+ c_x*sol_hom2_x.val_out_bound_jk(nz_bc, j, k)
+ d_x*dhom2_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
sec_membre.set(ligne) -= c_mu*sol_part_mu.val_out_bound_jk(nz_bc, j, k)
+ d_mu*dpart_mu.val_out_bound_jk(nz_bc, j, k)/alpha
+ c_x*sol_part_mu.val_out_bound_jk(nz_bc, j, k)
+ d_x*dpart_mu.val_out_bound_jk(nz_bc, j, k)/alpha
- mub(k, j) ;
}
// System inversion: Solution of the system giving the coefficients for the homogeneous
// solutions
//-------------------------------------------------------------------
systeme.set_lu() ;
Tbl facteur = systeme.inverse(sec_membre) ;
int conte = 0 ;
// everything is put to the right place...
//----------------------------------------
nr = mgrid.get_nr(0) ; //nucleus
for (int i=0 ; i<nr ; i++) {
mmu.set(0, k, j, i) = 0 ;
mw.set(0, k, j, i) = 0 ;
}
nr = mgrid.get_nr(1) ; //first shell
for (int i=0 ; i<nr ; i++) {
mmu.set(1, k, j, i) = sol_part_mu(1, k, j, i)
+ facteur(conte)*sol_hom1_mu(1, k, j, i)
+ facteur(conte+1)*sol_hom2_mu(1, k, j, i);
mw.set(1, k, j, i) = sol_part_x(1, k, j, i)
+ facteur(conte)*sol_hom1_x(1, k, j, i)
+ facteur(conte+1)*sol_hom2_x(1,k,j,i);
}
conte+=2 ;
for (int zone=2 ; zone<=n_shell ; zone++) { //shells
nr = mgrid.get_nr(zone) ;
for (int i=0 ; i<nr ; i++) {
mmu.set(zone, k, j, i) = sol_part_mu(zone, k, j, i)
+ facteur(conte)*sol_hom1_mu(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_mu(zone, k, j, i) ;
mw.set(zone, k, j, i) = sol_part_x(zone, k, j, i)
+ facteur(conte)*sol_hom1_x(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_x(zone, k, j, i) ;
}
conte+=2 ;
}
if (cedbc) {
nr = mgrid.get_nr(nzm1) ; //compactified external domain
for (int i=0 ; i<nr ; i++) {
mmu.set(nzm1, k, j, i) = sol_part_mu(nzm1, k, j, i)
+ facteur(conte)*sol_hom1_mu(nzm1, k, j, i) ;
mw.set(nzm1, k, j, i) = sol_part_x(nzm1, k, j, i)
+ facteur(conte)*sol_hom1_x(nzm1, k, j, i) ;
}
}
} // End of nullite_plm
} //End of loop on theta
if (tilde_mu.set_spectral_va().c != 0x0)
delete tilde_mu.set_spectral_va().c ;
tilde_mu.set_spectral_va().c = 0x0 ;
tilde_mu.set_spectral_va().ylm_i() ;
if (x_new.set_spectral_va().c != 0x0)
delete x_new.set_spectral_va().c ;
x_new.set_spectral_va().c = 0x0 ;
x_new.set_spectral_va().ylm_i();
}
//----------------------------------------------------------------------------------
//
// sol_Dirac_tilde_B
//
//----------------------------------------------------------------------------------
void Sym_tensor_trans::sol_Dirac_BC2(const Scalar& bb, const Scalar& cc, const Scalar& hh,
Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_eta,double dir, double neum, double rhor, Param* par_bc, Param* par_mat) {
const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
assert(mp_aff != 0x0) ; //Only affine mapping for the moment
const Mg3d& mgrid = *mp_aff->get_mg() ;
assert(mgrid.get_type_r(0) == RARE) ;
int nz = mgrid.get_nzone() ;
int nzm1 = nz - 1 ;
bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
int n_shell = ced ? nz-2 : nzm1 ;
int nz_bc = nzm1 ;
if (par_bc != 0x0)
if (par_bc->get_n_int() > 0)
nz_bc = par_bc->get_int() ;
n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
bool cedbc = (ced && (nz_bc == nzm1)) ;
#ifndef NDEBUG
if (!cedbc) {
assert(par_bc != 0x0) ;
assert(par_bc->get_n_tbl_mod() >= 2) ;
}
#endif
int nt = mgrid.get_nt(0) ;
int np = mgrid.get_np(0) ;
int l_q, m_q, base_r;
// ON passe en ylm les quantités B et C !!! CRUCIAL !!!
Scalar bb2 = bb; bb2.set_spectral_va().ylm();
Scalar cc2 = cc; cc2.set_spectral_va().ylm();
Scalar hh2 = hh; hh2.set_spectral_va().ylm();
Base_val base = hh2.get_spectral_base() ;
Scalar tilde_b(*mp_aff); tilde_b.annule_hard(); tilde_b.std_spectral_base();
tilde_b.set_spectral_va().ylm();
// Base_val base = hrr.get_spectral_base();
// int m_q, l_q, base_r ;
for (int lz=0; lz<nz; lz++) {
int nr = mgrid.get_nr(lz);
for (int k=0; k<np+1; k++)
for (int j=0; j<nt; j++) {
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
{
for (int i=0; i<nr; i++)
tilde_b.set_spectral_va().c_cf->set(lz, k, j, i)
+= 2*(*bb2.get_spectral_va().c_cf)(lz, k, j, i)
+ (*cc2.get_spectral_va().c_cf)(lz, k, j, i)/ (2* double(l_q +1.));
}
}
}
if (tilde_b.set_spectral_va().c != 0x0)
delete tilde_b.set_spectral_va().c ;
tilde_b.set_spectral_va().c = 0x0 ;
tilde_b.set_spectral_va().coef_i();
if ( (tilde_b.get_etat() == ETATZERO) && (hh.get_etat() == ETATZERO) ) {
hrr = 0 ;
tilde_eta = 0 ;
ww = 0 ;
return ;
}
bound_eta.set_spectral_va().ylm();
Scalar tilde_b2 = tilde_b;
tilde_b2.set_spectral_va().ylm();
assert (tilde_b.get_etat() != ETATNONDEF) ;
assert (hh.get_etat() != ETATNONDEF) ;
Scalar source = tilde_b ;
Scalar source_coq = tilde_b ;
source_coq.annule_domain(0) ;
source_coq.annule_domain(nzm1) ;
source_coq.mult_r() ;
source.set_spectral_va().ylm() ;
source_coq.set_spectral_va().ylm() ;
bool bnull = (tilde_b.get_etat() == ETATZERO) ;
assert(hh.check_dzpuis(0)) ;
Scalar hoverr = hh ;
hoverr.div_r_dzpuis(2) ;
hoverr.set_spectral_va().ylm() ;
Scalar dhdr = hh.dsdr() ;
dhdr.set_spectral_va().ylm() ;
Scalar h_coq = hh ;
h_coq.set_spectral_va().ylm() ;
Scalar dh_coq = hh.dsdr() ;
dh_coq.mult_r_dzpuis(0) ;
dh_coq.set_spectral_va().ylm() ;
bool hnull = (hh.get_etat() == ETATZERO) ;
int lmax = base.give_lmax(mgrid, 0) + 1;
// Utilisation des params, facultative, permettant uniquement une optimisation....
bool need_calculation = true ;
if (par_mat != 0x0) {
bool param_new = false ;
if ((par_mat->get_n_int_mod() >= 4)
&&(par_mat->get_n_tbl_mod()>=1)
&&(par_mat->get_n_matrice_mod()>=1)
&&(par_mat->get_n_itbl_mod()>=1)) {
if (par_mat->get_int_mod(0) < nz_bc) param_new = true ;
if (par_mat->get_int_mod(1) != lmax) param_new = true ;
if (par_mat->get_int_mod(2) != mgrid.get_type_t() ) param_new = true ;
if (par_mat->get_int_mod(3) != mgrid.get_type_p() ) param_new = true ;
if (par_mat->get_itbl_mod(0)(0) != mgrid.get_nr(0)) param_new = true ;
if (fabs(par_mat->get_tbl_mod(0)(0) - mp_aff->get_alpha()[0]) > 2.e-15)
param_new = true ;
for (int l=1; l<= n_shell; l++) {
if (par_mat->get_itbl_mod(0)(l) != mgrid.get_nr(l)) param_new = true ;
if (fabs(par_mat->get_tbl_mod(0)(l) - mp_aff->get_beta()[l] /
mp_aff->get_alpha()[l]) > 2.e-15) param_new = true ;
}
if (ced) {
if (par_mat->get_itbl_mod(0)(nzm1) != mgrid.get_nr(nzm1)) param_new = true ;
if (fabs(par_mat->get_tbl_mod(0)(nzm1) - mp_aff->get_alpha()[nzm1]) > 2.e-15)
param_new = true ;
}
}
else{
param_new = true ;
}
if (param_new) {
par_mat->clean_all() ;
int* nz_bc_new = new int(nz_bc) ;
par_mat->add_int_mod(*nz_bc_new, 0) ;
int* lmax_new = new int(lmax) ;
par_mat->add_int_mod(*lmax_new, 1) ;
int* type_t_new = new int(mgrid.get_type_t()) ;
par_mat->add_int_mod(*type_t_new, 2) ;
int* type_p_new = new int(mgrid.get_type_p()) ;
par_mat->add_int_mod(*type_p_new, 3) ;
Itbl* pnr = new Itbl(nz) ;
pnr->set_etat_qcq() ;
par_mat->add_itbl_mod(*pnr) ;
for (int l=0; l<nz; l++)
pnr->set(l) = mgrid.get_nr(l) ;
Tbl* palpha = new Tbl(nz) ;
palpha->set_etat_qcq() ;
par_mat->add_tbl_mod(*palpha) ;
palpha->set(0) = mp_aff->get_alpha()[0] ;
for (int l=1; l<nzm1; l++)
palpha->set(l) = mp_aff->get_beta()[l] / mp_aff->get_alpha()[l] ;
palpha->set(nzm1) = mp_aff->get_alpha()[nzm1] ;
}
else need_calculation = false ;
}
hrr.set_etat_qcq() ;
hrr.set_spectral_base(base) ;
hrr.set_spectral_va().set_etat_cf_qcq() ;
hrr.set_spectral_va().c_cf->annule_hard() ;
tilde_eta.annule_hard() ;
tilde_eta.set_spectral_base(base) ;
tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
tilde_eta.set_spectral_va().c_cf->annule_hard() ;
ww.annule_hard() ;
ww.set_spectral_base(base) ;
ww.set_spectral_va().set_etat_cf_qcq() ;
ww.set_spectral_va().c_cf->annule_hard() ;
sol_Dirac_l01(hh2, hrr, tilde_eta, par_mat) ;
tilde_eta.annule_l(0,0, true) ;
Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
Mtbl_cf sol_part_w(mgrid, base) ; sol_part_w.annule_hard() ;
Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
Mtbl_cf sol_hom1_w(mgrid, base) ; sol_hom1_w.annule_hard() ;
Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
Mtbl_cf sol_hom2_w(mgrid, base) ; sol_hom2_w.annule_hard() ;
Mtbl_cf sol_hom3_hrr(mgrid, base) ; sol_hom3_hrr.annule_hard() ;
Mtbl_cf sol_hom3_eta(mgrid, base) ; sol_hom3_eta.annule_hard() ;
Mtbl_cf sol_hom3_w(mgrid, base) ; sol_hom3_w.annule_hard() ;
Itbl mat_done(lmax) ;
// Calcul des solutions homogenes et particulieres...
//---------------
//-- NUCLEUS ---
//---------------
{int lz = 0 ;
int nr = mgrid.get_nr(lz) ;
double alpha = mp_aff->get_alpha()[lz] ;
Matrice ope(3*nr, 3*nr) ;
int ind2 = 2*nr ;
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+2*nr) = 0 ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col) = -0.5*md(lin,col)/double(l_q+1)
- 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+2*nr) = (l_q+2)*md(lin,col)
+ l_q*(l_q+2)*ms(lin,col) ;
ope *= 1./alpha ;
for (int col=0; col<3*nr; col++)
if (l_q>2) ope.set(ind2+nr-2, col) = 0 ;
for (int col=0; col<3*nr; col++) {
ope.set(nr-1, col) = 0 ;
ope.set(2*nr-1, col) = 0 ;
ope.set(3*nr-1, col) = 0 ;
}
int pari = 1 ;
if (base_r == R_CHEBP) {
for (int col=0; col<nr; col++) {
ope.set(nr-1, col) = pari ;
ope.set(2*nr-1, col+nr) = pari ;
ope.set(3*nr-1, col+2*nr) = pari ;
pari = - pari ;
}
}
else { //In the odd case, the last coefficient must be zero!
ope.set(nr-1, nr-1) = 1 ;
ope.set(2*nr-1, 2*nr-1) = 1 ;
ope.set(3*nr-1, 3*nr-1) = 1 ;
}
if (l_q>2)
ope.set(ind2+nr-2, ind2) = 1 ;
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(3*nr) ;
sec.set_etat_qcq() ;
if (hnull) {
for (int lin=0; lin<2*nr; lin++)
sec.set(lin) = 0 ;
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
(lz, k, j, lin) ;
if (bnull) {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
+ (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
}
}
if (l_q>2) sec.set(ind2+nr-2) = 0 ;
sec.set(3*nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.annule_hard() ;
if (l_q>2) {
sec.set(ind2+nr-2) = 1 ;
sol = oper.inverse(sec) ;
}
else { //Homogeneous solution put in by hand in the case l=2
sol.annule_hard() ;
sol.set(0) = 4 ;
sol.set(nr) = 2 ;
sol.set(2*nr) = 1 ;
}
for (int i=0; i<nr; i++) {
sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
}
}
}
}
//-------------
// -- Shells --
//-------------
for (int lz=1; lz<= n_shell; lz++) {
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
int nr = mgrid.get_nr(lz) ;
int ind0 = 0 ;
int ind1 = nr ;
int ind2 = 2*nr ;
double alpha = mp_aff->get_alpha()[lz] ;
double ech = mp_aff->get_beta()[lz] / alpha ;
Matrice ope(3*nr, 3*nr) ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
+ 3*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+2*nr) = 0 ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
+ 3*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col) =
-0.5/double(l_q+1)*(mxd(lin,col) + ech*md(lin,col)
+ double(l_q+4)*mid(lin,col)) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+nr) = -2*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+2*nr) =
double(l_q+2)*(mxd(lin,col) + ech*md(lin,col)
+ l_q*mid(lin,col)) ;
for (int col=0; col<3*nr; col++) {
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
ope.set(ind2+nr-1, col) = 0 ;
}
ope.set(ind0+nr-1, ind0) = 1 ;
ope.set(ind1+nr-1, ind1) = 1 ;
ope.set(ind2+nr-1, ind2) = 1 ;
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(3*nr) ;
sec.set_etat_qcq() ;
if (hnull) {
for (int lin=0; lin<2*nr; lin++)
sec.set(lin) = 0 ;
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = (*source_coq.get_spectral_va().c_cf)
(lz, k, j, lin) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(lin+nr) = -0.5*(*h_coq.get_spectral_va().c_cf)
(lz, k, j, lin) ;
if (bnull) {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) )
+ (*source_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
}
}
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 0 ;
sec.set(ind2+nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.annule_hard() ;
sec.set(ind0+nr-1) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.set(ind1+nr-1) = 0 ;
sec.set(ind2+nr-1) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
}
}
}
}
//------------------------------
// Compactified external domain
//------------------------------
if (cedbc) {int lz = nzm1 ;
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
int nr = mgrid.get_nr(lz) ;
int ind0 = 0 ;
int ind1 = nr ;
int ind2 = 2*nr ;
double alpha = mp_aff->get_alpha()[lz] ;
Matrice ope(3*nr, 3*nr) ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+2*nr) = 0 ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col) = 0.5*md(lin,col)/double(l_q+1)
- 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+2*nr,col+2*nr) = -(l_q+2)*md(lin,col)
+ l_q*(l_q+2)*ms(lin,col) ;
ope *= 1./alpha ;
for (int col=0; col<3*nr; col++) {
ope.set(ind0+nr-2, col) = 0 ;
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-2, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
ope.set(ind2+nr-1, col) = 0 ;
}
for (int col=0; col<nr; col++) {
ope.set(ind0+nr-1, col+ind0) = 1 ;
ope.set(ind1+nr-1, col+ind1) = 1 ;
ope.set(ind2+nr-1, col+ind2) = 1 ;
}
ope.set(ind0+nr-2, ind0+1) = 1 ;
ope.set(ind1+nr-2, ind1+2) = 1 ;
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(3*nr) ;
sec.set_etat_qcq() ;
if (hnull) {
for (int lin=0; lin<2*nr; lin++)
sec.set(lin) = 0 ;
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
(lz, k, j, lin) ;
if (bnull) {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
}
else {
for (int lin=0; lin<nr; lin++)
sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
(*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
+ (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
+ (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
}
}
sec.set(ind0+nr-2) = 0 ;
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 0 ;
sec.set(ind1+nr-2) = 0 ;
sec.set(ind2+nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.annule_hard() ;
sec.set(ind0+nr-2) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
sec.set(ind0+nr-2) = 0 ;
sec.set(ind1+nr-2) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
}
}
}
}
}
// Matching system...
int taille = 3*nz_bc ; // Un de moins que dans la version complete R3
if (cedbc) taille-- ;
Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
Mtbl_cf& mw = *ww.set_spectral_va().c_cf ;
Tbl sec_membre(taille) ;
Matrice systeme(taille, taille) ;
int ligne ; int colonne ;
Tbl pipo(1) ;
const Tbl& hrrb = (cedbc ? pipo : par_bc->get_tbl_mod(1) );
double chrr = (cedbc ? 0 : par_bc->get_tbl_mod()(4) ) ;
double dhrr = (cedbc ? 0 : par_bc->get_tbl_mod()(5) ) ;
double ceta = (cedbc ? 0 : par_bc->get_tbl_mod()(6) ) ;
double deta = (cedbc ? 0 : par_bc->get_tbl_mod()(7) ) ;
double cw = (cedbc ? 0 : par_bc->get_tbl_mod()(8) ) ;
double dw = (cedbc ? 0 : par_bc->get_tbl_mod()(9) ) ;
Mtbl_cf dhom1_hrr = sol_hom1_hrr ;
Mtbl_cf dhom2_hrr = sol_hom2_hrr ;
Mtbl_cf dhom3_hrr = sol_hom3_hrr ;
Mtbl_cf dpart_hrr = sol_part_hrr ;
Mtbl_cf dhom1_eta = sol_hom1_eta ;
Mtbl_cf dhom2_eta = sol_hom2_eta ;
Mtbl_cf dhom3_eta = sol_hom3_eta ;
Mtbl_cf dpart_eta = sol_part_eta ;
Mtbl_cf dhom1_w = sol_hom1_w ;
Mtbl_cf dhom2_w = sol_hom2_w ;
Mtbl_cf dhom3_w = sol_hom3_w ;
Mtbl_cf dpart_w = sol_part_w ;
dhom1_hrr.dsdx() ; dhom1_eta.dsdx() ; dhom1_w.dsdx() ;
dhom2_hrr.dsdx() ; dhom2_eta.dsdx() ; dhom2_w.dsdx() ;
dhom3_hrr.dsdx() ; dhom3_eta.dsdx() ; dhom3_w.dsdx() ;
dpart_hrr.dsdx() ; dpart_eta.dsdx() ; dpart_w.dsdx() ;
Mtbl_cf d2hom1_hrr = dhom1_hrr ;
Mtbl_cf d2hom2_hrr = dhom2_hrr ;
Mtbl_cf d2hom3_hrr = dhom3_hrr;
Mtbl_cf d2part_hrr = dpart_hrr ;
d2hom1_hrr.dsdx(); d2hom2_hrr.dsdx(); d2hom3_hrr.dsdx(); d2part_hrr.dsdx();
// Loop on l and m//
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
ligne = 0 ;
colonne = 0 ;
systeme.annule_hard() ;
sec_membre.annule_hard() ;
//First shell
//Condition at x=-1;
int nr = mgrid.get_nr(1);
double alpha2 = mp_aff->get_alpha()[1] ;
// A compatibility condition is set in order to be coherent with second order potentials definition
const Coord& rr = (*mp_aff).r ;
Scalar rrr(*mp_aff); rrr = rr;
systeme.set(ligne, colonne)= (dhom1_w.val_in_bound_jk(1,j,k))/alpha2 + (0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom1_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom1_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom1_hrr.val_in_bound_jk(1,j,k);
systeme.set(ligne, colonne+1)= (dhom2_w.val_in_bound_jk(1,j,k))/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom2_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom2_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom2_hrr.val_in_bound_jk(1,j,k);
systeme.set(ligne, colonne+2)= (dhom3_w.val_in_bound_jk(1,j,k))/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom3_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom3_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom3_hrr.val_in_bound_jk(1,j,k) ;
sec_membre.set(ligne)= - (dpart_w.val_in_bound_jk(1,j,k)/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_part_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_part_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_part_hrr.val_in_bound_jk(1,j,k));
Mtbl_cf *Bbb = bb2.get_spectral_va().c_cf;
sec_membre.set(ligne) += (*Bbb).val_in_bound_jk(1,j,k);
ligne++ ;
// A Robyn-type boundary condition is imposed on eta
systeme.set(ligne, colonne) =
dir*sol_hom1_hrr.val_in_bound_jk(1, j, k) + neum*dhom1_hrr.val_in_bound_jk(1,j,k)/alpha2;
systeme.set(ligne, colonne+1) =
dir*sol_hom2_hrr.val_in_bound_jk(1, j, k) +neum*dhom2_hrr.val_in_bound_jk(1,j,k)/alpha2;
systeme.set(ligne, colonne+2) =
+ dir*sol_hom3_hrr.val_in_bound_jk(1, j, k) + neum*dhom3_hrr.val_in_bound_jk(1,j,k)/alpha2;
double blob = (*(bound_eta.get_spectral_va().c_cf)).val_in_bound_jk(1, j, k);
sec_membre.set(ligne) += - sol_part_hrr.val_in_bound_jk(1, j, k) - neum*dpart_hrr.val_in_bound_jk(1,j,k)/alpha2 + blob ;
ligne++ ;
//Conditions at x=1;
systeme.set(ligne, colonne) =
sol_hom1_hrr.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_hrr.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_hrr.val_out_bound_jk(1, j, k) ;
sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(1, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_eta.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_eta.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_eta.val_out_bound_jk(1, j, k) ;
sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(1, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_w.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_w.val_out_bound_jk(1, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_w.val_out_bound_jk(1, j, k) ;
sec_membre.set(ligne) -= sol_part_w.val_out_bound_jk(1, j, k) ;
colonne += 3 ;
//shells
for (int zone=2 ; zone<nz_bc ; zone++) {
nr = mgrid.get_nr(zone) ;
ligne -= 2 ;
//Condition at x = -1
systeme.set(ligne, colonne) =
- sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
- sol_hom3_hrr.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
- sol_hom3_eta.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_w.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_w.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
- sol_hom3_w.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(zone, j, k) ;
ligne++ ;
// Condition at x=1
systeme.set(ligne, colonne) =
sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_hrr.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_eta.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_w.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_w.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+2) =
sol_hom3_w.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_w.val_out_bound_jk(zone, j, k) ;
colonne += 3 ;
}
//Last domain
nr = mgrid.get_nr(nz_bc) ;
double alpha = mp_aff->get_alpha()[nz_bc] ;
ligne -= 2 ;
//Condition at x = -1
systeme.set(ligne, colonne) =
- sol_hom1_hrr.val_in_bound_jk(nz_bc, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_hrr.val_in_bound_jk(nz_bc, j, k) ;
if (!cedbc) systeme.set(ligne, colonne+2) =
- sol_hom3_hrr.val_in_bound_jk(nz_bc, j, k) ;
sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz_bc, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_eta.val_in_bound_jk(nz_bc, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_eta.val_in_bound_jk(nz_bc, j, k) ;
if (!cedbc) systeme.set(ligne, colonne+2) =
- sol_hom3_eta.val_in_bound_jk(nz_bc, j, k) ;
sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz_bc, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_w.val_in_bound_jk(nz_bc, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_w.val_in_bound_jk(nz_bc, j, k) ;
if (!cedbc) systeme.set(ligne, colonne+2) =
- sol_hom3_w.val_in_bound_jk(nz_bc, j, k) ;
sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(nz_bc, j, k) ;
ligne++ ;
if (!cedbc) {//Special condition at x=1
systeme.set(ligne, colonne) =
chrr*sol_hom1_hrr.val_out_bound_jk(nz_bc, j, k)
+ dhrr*dhom1_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
+ ceta*sol_hom1_eta.val_out_bound_jk(nz_bc, j, k)
+ deta*dhom1_eta.val_out_bound_jk(nz_bc, j, k) / alpha
+ cw*sol_hom1_w.val_out_bound_jk(nz_bc, j, k)
+ dw*dhom1_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
systeme.set(ligne, colonne+1) =
chrr*sol_hom2_hrr.val_out_bound_jk(nz_bc, j, k)
+ dhrr*dhom2_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
+ ceta*sol_hom2_eta.val_out_bound_jk(nz_bc, j, k)
+ deta*dhom2_eta.val_out_bound_jk(nz_bc, j, k) / alpha
+ cw*sol_hom2_w.val_out_bound_jk(nz_bc, j, k)
+ dw*dhom2_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
systeme.set(ligne, colonne+2) =
chrr*sol_hom3_hrr.val_out_bound_jk(nz_bc, j, k)
+ dhrr*dhom3_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
+ ceta*sol_hom3_eta.val_out_bound_jk(nz_bc, j, k)
+ deta*dhom3_eta.val_out_bound_jk(nz_bc, j, k) / alpha
+ cw*sol_hom3_w.val_out_bound_jk(nz_bc, j, k)
+ dw*dhom3_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
sec_membre.set(ligne) -= chrr*sol_part_hrr.val_out_bound_jk(nz_bc, j, k)
+ dhrr*dpart_hrr.val_out_bound_jk(nz_bc, j, k)/alpha
+ ceta*sol_part_eta.val_out_bound_jk(nz_bc, j, k)
+ deta*dpart_eta.val_out_bound_jk(nz_bc, j, k)/alpha
+ cw*sol_part_w.val_out_bound_jk(nz_bc, j, k)
+ dw*dpart_w.val_out_bound_jk(nz_bc, j, k)/alpha
- hrrb(k, j) ;
}
//System inversion: Solution of the system giving the coefficients for the homogeneous
// solutions
//-------------------------------------------------------------------
systeme.set_lu() ;
Tbl facteur = systeme.inverse(sec_membre) ;
int conte = 0 ;
// everything is put to the right place,
//---------------------------------------
nr = mgrid.get_nr(0) ; //nucleus
for (int i=0 ; i<nr ; i++) {
mhrr.set(0, k, j, i) = 0 ;
meta.set(0, k, j, i) = 0 ;
mw.set(0, k, j, i) = 0 ;
}
for (int zone=1 ; zone<=n_shell ; zone++) { //shells
nr = mgrid.get_nr(zone) ;
for (int i=0 ; i<nr ; i++) {
mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
+ facteur(conte)*sol_hom1_hrr(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_hrr(zone, k, j, i)
+ facteur(conte+2)*sol_hom3_hrr(zone, k, j, i) ;
meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
+ facteur(conte)*sol_hom1_eta(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_eta(zone, k, j, i)
+ facteur(conte+2)*sol_hom3_eta(zone, k, j, i) ;
mw.set(zone, k, j, i) = sol_part_w(zone, k, j, i)
+ facteur(conte)*sol_hom1_w(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_w(zone, k, j, i)
+ facteur(conte+2)*sol_hom3_w(zone, k, j, i) ;
}
conte+=3 ;
}
if (cedbc) {
nr = mgrid.get_nr(nzm1) ; //compactified external domain
for (int i=0 ; i<nr ; i++) {
mhrr.set(nzm1, k, j, i) = sol_part_hrr(nzm1, k, j, i)
+ facteur(conte)*sol_hom1_hrr(nzm1, k, j, i)
+ facteur(conte+1)*sol_hom2_hrr(nzm1, k, j, i) ;
meta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
+ facteur(conte)*sol_hom1_eta(nzm1, k, j, i)
+ facteur(conte+1)*sol_hom2_eta(nzm1, k, j, i) ;
mw.set(nzm1, k, j, i) = sol_part_w(nzm1, k, j, i)
+ facteur(conte)*sol_hom1_w(nzm1, k, j, i)
+ facteur(conte+1)*sol_hom2_w(nzm1, k, j, i) ;
}
}
} // End of nullite_plm
} //End of loop on theta
if (hrr.set_spectral_va().c != 0x0)
delete hrr.set_spectral_va().c;
hrr.set_spectral_va().c = 0x0 ;
hrr.set_spectral_va().ylm_i() ;
if (tilde_eta.set_spectral_va().c != 0x0)
delete tilde_eta.set_spectral_va().c ;
tilde_eta.set_spectral_va().c = 0x0 ;
tilde_eta.set_spectral_va().ylm_i() ;
if (ww.set_spectral_va().c != 0x0)
delete ww.set_spectral_va().c ;
ww.set_spectral_va().c = 0x0 ;
ww.set_spectral_va().ylm_i() ;
// // On doit resoudre séparément les cas l=0 et l=1; notamment dansces cas le systeme electrique se resout a deux equations
// // au lieu de 3.
}
void Sym_tensor_trans::sol_Dirac_l01_2(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta,
Param* par_mat) {
// void Sym_tensor_trans::sol_Dirac_tilde_B2(const Scalar& tilde_b, const Scalar& hh,
// Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_eta,double dir, double neum, double rhor, Param* par_bc, Param* par_mat) {
const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
assert(mp_aff != 0x0) ; //Only affine mapping for the moment
const Mg3d& mgrid = *mp_aff->get_mg() ;
int nz = mgrid.get_nzone() ;
assert(mgrid.get_type_r(0) == RARE) ;
assert(mgrid.get_type_r(nz-1) == UNSURR) ;
if (hh.get_etat() == ETATZERO) {
hrr.annule_hard() ;
tilde_eta.annule_hard() ;
return ;
}
int nt = mgrid.get_nt(0) ;
int np = mgrid.get_np(0) ;
Scalar source = hh ;
source.div_r_dzpuis(2) ;
Scalar source_coq = hh ;
source.set_spectral_va().ylm() ;
source_coq.set_spectral_va().ylm() ;
Base_val base = source.get_spectral_base() ;
base.mult_x() ;
int lmax = base.give_lmax(mgrid, 0) + 1;
assert (hrr.get_spectral_base() == base) ;
assert (tilde_eta.get_spectral_base() == base) ;
assert (hrr.get_spectral_va().c_cf != 0x0) ;
assert (tilde_eta.get_spectral_va().c_cf != 0x0) ;
Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
bool need_calculation = true ;
if (par_mat != 0x0)
if (par_mat->get_n_matrice_mod() > 0)
if (&par_mat->get_matrice_mod(0) != 0x0) need_calculation = false ;
int l_q, m_q, base_r ;
Itbl mat_done(lmax) ;
//---------------
//-- NUCLEUS ---
//---------------
{int lz = 0 ;
int nr = mgrid.get_nr(lz) ;
double alpha = mp_aff->get_alpha()[lz] ;
Matrice ope(2*nr, 2*nr) ;
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin, col);
ope *= 1./alpha ;
for (int col=0; col<2*nr; col++) {
ope.set(nr-1, col) = 0 ;
ope.set(2*nr-1, col) = 0 ;
}
int pari = 1 ;
if (base_r == R_CHEBP) {
for (int col=0; col<nr; col++) {
ope.set(nr-1, col) = pari ;
ope.set(2*nr-1, col+nr) = pari ;
pari = - pari ;
}
}
else { //In the odd case, the last coefficient must be zero!
ope.set(nr-1, nr-1) = 1 ;
ope.set(2*nr-1, 2*nr-1) = 1 ;
}
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(2*nr) ;
sec.set_etat_qcq() ;
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
sec.set(nr-1) = 0 ;
if (base_r == R_CHEBP) {
double h0 = 0 ; //In the l=0 case: 3*hrr(r=0) = h(r=0)
int pari = 1 ;
for (int col=0; col<nr; col++) {
h0 += pari*
(*source_coq.get_spectral_va().c_cf)(lz, k, j, col) ;
pari = - pari ;
}
sec.set(nr-1) = h0 / 3. ;
}
sec.set(2*nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
}
sec.annule_hard() ;
}
}
}
}
//-------------
// -- Shells --
//-------------
for (int lz=1; lz<nz-1; lz++) {
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
int nr = mgrid.get_nr(lz) ;
int ind0 = 0 ;
int ind1 = nr ;
assert(mgrid.get_nt(lz) == nt) ;
assert(mgrid.get_np(lz) == np) ;
double alpha = mp_aff->get_alpha()[lz] ;
double ech = mp_aff->get_beta()[lz] / alpha ;
Matrice ope(2*nr, 2*nr) ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
+ 3*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
+ 3*mid(lin, col) ;
for (int col=0; col<2*nr; col++) {
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
}
ope.set(ind0+nr-1, ind0) = 1 ;
ope.set(ind1+nr-1, ind1) = 1 ;
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(2*nr) ;
sec.set_etat_qcq() ;
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*source_coq.get_spectral_va().c_cf)
(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(nr+lin) = -0.5*(*source_coq.get_spectral_va().c_cf)
(lz, k, j, lin) ;
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
}
sec.annule_hard() ;
sec.set(ind0+nr-1) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
}
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-1) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
}
}
}
}
}
//------------------------------
// Compactified external domain
//------------------------------
{int lz = nz-1 ;
if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
int nr = mgrid.get_nr(lz) ;
int ind0 = 0 ;
int ind1 = nr ;
assert(mgrid.get_nt(lz) == nt) ;
assert(mgrid.get_np(lz) == np) ;
double alpha = mp_aff->get_alpha()[lz] ;
Matrice ope(2*nr, 2*nr) ;
for (int k=0 ; k<np+1 ; k++) {
for (int j=0 ; j<nt ; j++) {
// quantic numbers and spectral bases
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
if (need_calculation) {
ope.set_etat_qcq() ;
Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
for (int lin=0; lin<nr; lin++)
for (int col=0; col<nr; col++)
ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin, col) ;
ope *= 1./alpha ;
for (int col=0; col<2*nr; col++) {
ope.set(ind0+nr-2, col) = 0 ;
ope.set(ind0+nr-1, col) = 0 ;
ope.set(ind1+nr-2, col) = 0 ;
ope.set(ind1+nr-1, col) = 0 ;
}
for (int col=0; col<nr; col++) {
ope.set(ind0+nr-1, col+ind0) = 1 ;
ope.set(ind1+nr-1, col+ind1) = 1 ;
}
ope.set(ind0+nr-2, ind0+1) = 1 ;
ope.set(ind1+nr-2, ind1+1) = 1 ;
ope.set_lu() ;
if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
Matrice* pope = new Matrice(ope) ;
par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
mat_done.set(l_q) = 1 ;
}
} //End of case when a calculation is needed
const Matrice& oper = (par_mat == 0x0 ? ope :
par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
Tbl sec(2*nr) ;
sec.set_etat_qcq() ;
for (int lin=0; lin<nr; lin++)
sec.set(lin) = (*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
for (int lin=0; lin<nr; lin++)
sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
(lz, k, j, lin) ;
sec.set(ind0+nr-2) = 0 ;
sec.set(ind0+nr-1) = 0 ;
sec.set(ind1+nr-2) = 0 ;
sec.set(ind1+nr-1) = 0 ;
Tbl sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_part_hrr.set(lz, k, j, i) = sol(i) ;
sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
}
sec.annule_hard() ;
sec.set(ind0+nr-2) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
}
sec.set(ind0+nr-2) = 0 ;
sec.set(ind1+nr-2) = 1 ;
sol = oper.inverse(sec) ;
for (int i=0; i<nr; i++) {
sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
}
}
}
}
}
int taille = 2*(nz-1) ;
Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
Tbl sec_membre(taille) ;
Matrice systeme(taille, taille) ;
int ligne ; int colonne ;
// Loop on l and m
//----------------
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
ligne = 0 ;
colonne = 0 ;
systeme.annule_hard() ;
sec_membre.annule_hard() ;
//Nucleus
int nr = mgrid.get_nr(0) ;
sec_membre.set(ligne) = -sol_part_hrr.val_out_bound_jk(0, j, k) ;
ligne++ ;
sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
//shells
for (int zone=1 ; zone<nz-1 ; zone++) {
nr = mgrid.get_nr(zone) ;
ligne-- ;
//Condition at x = -1
systeme.set(ligne, colonne) =
- sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
ligne++ ;
// Condition at x=1
systeme.set(ligne, colonne) =
sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
systeme.set(ligne, colonne+1) =
sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
colonne += 2 ;
}
//Compactified external domain
nr = mgrid.get_nr(nz-1) ;
ligne-- ;
systeme.set(ligne, colonne) =
- sol_hom1_hrr.val_in_bound_jk(nz-1, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_hrr.val_in_bound_jk(nz-1, j, k) ;
sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz-1, j, k) ;
ligne++ ;
systeme.set(ligne, colonne) =
- sol_hom1_eta.val_in_bound_jk(nz-1, j, k) ;
systeme.set(ligne, colonne+1) =
- sol_hom2_eta.val_in_bound_jk(nz-1, j, k) ;
sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz-1, j, k) ;
// Solution of the system giving the coefficients for the homogeneous
// solutions
//-------------------------------------------------------------------
systeme.set_lu() ;
Tbl facteur = systeme.inverse(sec_membre) ;
int conte = 0 ;
// everything is put to the right place...
//----------------------------------------
nr = mgrid.get_nr(0) ; //nucleus
for (int i=0 ; i<nr ; i++) {
mhrr.set(0, k, j, i) = sol_part_hrr(0, k, j, i) ;
meta.set(0, k, j, i) = sol_part_eta(0, k, j, i) ;
}
for (int zone=1 ; zone<nz-1 ; zone++) { //shells
nr = mgrid.get_nr(zone) ;
for (int i=0 ; i<nr ; i++) {
mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
+ facteur(conte)*sol_hom1_hrr(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_hrr(zone, k, j, i) ;
meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
+ facteur(conte)*sol_hom1_eta(zone, k, j, i)
+ facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
}
conte+=2 ;
}
nr = mgrid.get_nr(nz-1) ; //compactified external domain
for (int i=0 ; i<nr ; i++) {
mhrr.set(nz-1, k, j, i) = sol_part_hrr(nz-1, k, j, i)
+ facteur(conte)*sol_hom1_hrr(nz-1, k, j, i)
+ facteur(conte+1)*sol_hom2_hrr(nz-1, k, j, i) ;
meta.set(nz-1, k, j, i) = sol_part_eta(nz-1, k, j, i)
+ facteur(conte)*sol_hom1_eta(nz-1, k, j, i)
+ facteur(conte+1)*sol_hom2_eta(nz-1, k, j, i) ;
}
} // End of nullite_plm
} //End of loop on theta
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::mult_x
void mult_x()
The basis is transformed as with a multiplication by .
Definition base_val_manip.C:157
Lorene::Base_val::give_lmax
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
Definition base_val_quantum.C:294
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:81
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
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:94
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Lorene::Diff_id::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_id.C:91
Lorene::Diff_sx
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
Lorene::Diff_sx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sx.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:98
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition itbl.C:261
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Itbl::annule_hard
void annule_hard()
Sets the Itbl to zero in a hard way.
Definition itbl.C:272
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
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:175
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:424
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition matrice.C:193
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_type_t
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition grilles.h:485
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:462
Lorene::Mg3d::get_type_p
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
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:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
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:452
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:474
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::dsdx
void dsdx()
Definition valeur_dsdx.C:153
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
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:428
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:452
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::get_n_itbl_mod
int get_n_itbl_mod() const
Returns the number of modifiable Itbl 's addresses in the list.
Definition param.C:722
Lorene::Param::get_n_matrice_mod
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
Definition param.C:859
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
Lorene::Param::get_n_tbl_mod
int get_n_tbl_mod() const
Returns the number of modifiable Tbl 's addresses in the list.
Definition param.C:584
Lorene::Param::get_itbl_mod
Itbl & get_itbl_mod(int position=0) const
Returns the reference of a stored modifiable Itbl .
Definition param.C:774
Lorene::Param::get_n_int_mod
int get_n_int_mod() const
Returns the number of modifiable int 's addresses in the list.
Definition param.C:378
Lorene::Param::get_tbl_mod
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
Definition param.C:636
Lorene::Param::clean_all
void clean_all()
Deletes all the objects stored as modifiables, i.e.
Definition param.C:174
Lorene::Param::add_matrice_mod
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Definition param.C:866
Lorene::Param::get_matrice_mod
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
Definition param.C:911
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
Lorene::Param::add_tbl_mod
void add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
Definition param.C:591
Lorene::Param::add_itbl_mod
void add_itbl_mod(Itbl &ti, int position=0)
Adds the address of a new modifiable Itbl to the list.
Definition param.C:729
Lorene::Param::get_n_int
int get_n_int() const
Returns the number of stored int 's addresses.
Definition param.C:239
Lorene::Param::get_int_mod
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition param.C:430
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::annule_l
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
Definition scalar_manip.C:111
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:147
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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:784
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1294
Lorene::Scalar::mult_r_dzpuis
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:221
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:208
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:797
Lorene::Sym_tensor_trans::sol_Dirac_BC2
void sol_Dirac_BC2(const Scalar &bb, const Scalar &cc, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &ww, Scalar bound_eta, double dir, double neum, double rhor, Param *par_bc, Param *par_mat)
Same resolution as sol_Dirac_tilde_B, but with inner boundary conditions added.
Definition sym_tensor_trans_dirac_bound2.C:580
Lorene::Sym_tensor_trans::sol_Dirac_Abound
void sol_Dirac_Abound(const Scalar &aaa, Scalar &tilde_mu, Scalar &x_new, Scalar bound_mu, const Param *par_bc)
Same resolution as sol_Dirac_A, but with inner boundary conditions added.
Definition sym_tensor_trans_dirac_bound2.C:81
Lorene::Sym_tensor_trans::sol_Dirac_l01
void sol_Dirac_l01(const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Param *par_mat) const
Solves the same system as Sym_tensor_trans::sol_Dirac_tilde_B but only for .
Definition sym_tensor_trans_dirac.C:1438
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:372
Lorene::Tbl::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
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:712
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:299
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:140
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
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:295
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene
Lorene prototypes.
Definition app_hor.h:64