vector.C
/*
* Methods of class Vector
*
* (see file vector.h for documentation)
*
*/
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* 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 as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* 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 vector_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector.C,v 1.30 2014/10/13 08:53:44 j_novak Exp $" ;
/*
* $Id: vector.C,v 1.30 2014/10/13 08:53:44 j_novak Exp $
* $Log: vector.C,v $
* Revision 1.30 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.29 2014/10/06 15:13:20 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.28 2008/10/29 14:09:14 jl_cornou
* Spectral bases for pseudo vectors and curl added
*
* Revision 1.27 2008/08/27 08:52:23 jl_cornou
* Added fonctions for angular potential A
*
* Revision 1.26 2007/12/21 16:07:08 j_novak
* Methods to filter Tensor, Vector and Sym_tensor objects.
*
* Revision 1.25 2005/02/14 13:01:50 j_novak
* p_eta and p_mu are members of the class Vector. Most of associated functions
* have been moved from the class Vector_divfree to the class Vector.
*
* Revision 1.24 2005/01/25 15:37:35 j_novak
* Solved some dzpuis problem...
*
* Revision 1.23 2005/01/12 16:48:23 j_novak
* Better treatment of the case where all vector components are null in
* decompose_div .
*
* Revision 1.22 2004/10/12 09:58:25 j_novak
* Better memory management.
*
* Revision 1.21 2004/10/11 09:46:31 j_novak
* Speed improvements.
*
* Revision 1.20 2004/05/09 20:55:05 e_gourgoulhon
* Added method flux.
*
* Revision 1.19 2004/03/29 11:57:45 e_gourgoulhon
* Added methods ope_killing and ope_killing_conf.
*
* Revision 1.18 2004/02/26 22:48:50 e_gourgoulhon
* -- Method divergence: call to Tensor::divergence and cast of the
* result.
* -- Added method derive_lie.
*
* Revision 1.17 2004/02/24 09:46:20 j_novak
* Correction to cope with SGI compiler's warnings.
*
* Revision 1.16 2004/02/20 10:53:04 j_novak
* Added the matching of the potential adapted to the case where the
* vector is the source of a Poisson equation (dzpuis = 4).
*
* Revision 1.15 2004/01/30 10:30:17 j_novak
* Changed dzpuis handling in Vector::decompose_div (this may be temporary).
*
* Revision 1.14 2003/12/30 23:09:47 e_gourgoulhon
* Change in methods derive_cov() and divergence() to take into account
* the change of name: Metric::get_connect() --> Metric::connect().
*
* Revision 1.13 2003/12/19 15:18:16 j_novak
* Shadow variables hunt
*
* Revision 1.12 2003/10/29 11:04:34 e_gourgoulhon
* dec2_dpzuis() replaced by dec_dzpuis(2).
* inc2_dpzuis() replaced by inc_dzpuis(2).
*
* Revision 1.11 2003/10/22 14:24:19 j_novak
* *** empty log message ***
*
* Revision 1.9 2003/10/20 13:00:38 j_novak
* Memory error corrected
*
* Revision 1.8 2003/10/20 09:32:12 j_novak
* Members p_potential and p_div_free of the Helmholtz decomposition
* + the method decompose_div(Metric).
*
* Revision 1.7 2003/10/16 14:21:37 j_novak
* The calculation of the divergence of a Tensor is now possible.
*
* Revision 1.6 2003/10/13 13:52:40 j_novak
* Better managment of derived quantities.
*
* Revision 1.5 2003/10/06 13:58:48 j_novak
* The memory management has been improved.
* Implementation of the covariant derivative with respect to the exact Tensor
* type.
*
* Revision 1.4 2003/10/05 21:14:20 e_gourgoulhon
* Added method std_spectral_base().
*
* Revision 1.3 2003/10/03 14:10:32 e_gourgoulhon
* Added constructor from Tensor.
*
* Revision 1.2 2003/10/03 14:08:46 j_novak
* Removed old change_trid...
*
* Revision 1.1 2003/09/26 08:05:31 j_novak
* New class Vector.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/vector.C,v 1.30 2014/10/13 08:53:44 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
// Headers Lorene
#include "metric.h"
#include "proto.h"
#include "matrice.h"
#include "nbr_spx.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Vector::Vector(const Map& map, int tipe, const Base_vect& triad_i)
: Tensor(map, 1, tipe, triad_i) {
set_der_0x0() ;
}
// Standard constructor with the triad passed as a pointer
// -------------------------------------------------------
Vector::Vector(const Map& map, int tipe, const Base_vect* triad_i)
: Tensor(map, 1, tipe, *triad_i) {
set_der_0x0() ;
}
// Copy constructor
// ----------------
Vector::Vector (const Vector& source) :
Tensor(source) {
assert(valence == 1) ;
set_der_0x0() ;
}
// Constructor from a {\tt Tensor}.
//--------------------------------
Vector::Vector(const Tensor& uu) : Tensor(uu) {
assert(valence == 1) ;
set_der_0x0() ;
}
// Constructor from a file
// -----------------------
Vector::Vector(const Map& mapping, const Base_vect& triad_i, FILE* fd) :
Tensor(mapping, triad_i, fd) {
assert ( (valence == 1) && (n_comp == 3) ) ;
set_der_0x0() ;
}
//--------------//
// Destructor //
//--------------//
Vector::~Vector () {
Vector::del_deriv() ;
}
//-------------------//
// Memory managment //
//-------------------//
void Vector::del_deriv() const {
for (int i=0; i<N_MET_MAX; i++)
del_derive_met(i) ;
if (p_A != 0x0) delete p_A ;
if (p_eta != 0x0) delete p_eta ;
if (p_mu != 0x0) delete p_mu ;
set_der_0x0() ;
Tensor::del_deriv() ;
}
void Vector::set_der_0x0() const {
for (int i=0; i<N_MET_MAX; i++)
set_der_met_0x0(i) ;
p_A = 0x0 ;
p_eta = 0x0 ;
p_mu = 0x0 ;
}
void Vector::del_derive_met(int j) const {
assert( (j>=0) && (j<N_MET_MAX) ) ;
if (met_depend[j] != 0x0) {
if (p_potential[j] != 0x0)
delete p_potential[j] ;
if (p_div_free[j] != 0x0)
delete p_div_free[j] ;
set_der_met_0x0(j) ;
Tensor::del_derive_met(j) ;
}
}
void Vector::set_der_met_0x0(int i) const {
assert( (i>=0) && (i<N_MET_MAX) ) ;
p_potential[i] = 0x0 ;
p_div_free[i] = 0x0 ;
}
void Vector::operator=(const Vector& t) {
triad = t.triad ;
assert(t.type_indice(0) == type_indice(0)) ;
for (int i=0 ; i<3 ; i++) {
*cmp[i] = *t.cmp[i] ;
}
del_deriv() ;
}
void Vector::operator=(const Tensor& t) {
assert (t.valence == 1) ;
triad = t.triad ;
assert(t.type_indice(0) == type_indice(0)) ;
for (int i=0 ; i<3 ; i++) {
*cmp[i] = *t.cmp[i] ;
}
del_deriv() ;
}
// Affectation d'une composante :
Scalar& Vector::set(int index) {
assert ( (index>=1) && (index<=3) ) ;
del_deriv() ;
return *cmp[index - 1] ;
}
const Scalar& Vector::operator()(int index) const {
assert ( (index>=1) && (index<=3) ) ;
return *cmp[index - 1] ;
}
// Sets the standard spectal bases of decomposition for each component
void Vector::std_spectral_base() {
Base_val** bases = 0x0 ;
if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
// Cartesian case
bases = mp->get_mg()->std_base_vect_cart() ;
}
else {
// Spherical case
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
bases = mp->get_mg()->std_base_vect_spher() ;
}
for (int i=0 ; i<3 ; i++) {
cmp[i]->set_spectral_base( *bases[i] ) ;
}
for (int i=0 ; i<3 ; i++) {
delete bases[i] ;
}
delete [] bases ;
}
// Sets the standard spectral bases of decomposition for each component for a pseudo vector
void Vector::pseudo_spectral_base() {
Base_val** bases = 0x0 ;
if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
// Cartesian case
bases = mp->get_mg()->pseudo_base_vect_cart() ;
}
else {
// Spherical case
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
bases = mp->get_mg()->pseudo_base_vect_spher() ;
}
for (int i=0 ; i<3 ; i++) {
cmp[i]->set_spectral_base( *bases[i] ) ;
}
for (int i=0 ; i<3 ; i++) {
delete bases[i] ;
}
delete [] bases ;
}
//-------------------------------//
// Computational methods //
//-------------------------------//
const Scalar& Vector::divergence(const Metric& metre) const {
const Scalar* pscal =
dynamic_cast<const Scalar*>( &(Tensor::divergence(metre)) ) ;
assert(pscal != 0x0) ;
return *pscal ;
}
Vector Vector::derive_lie(const Vector& vv) const {
Vector resu(*mp, type_indice(0), triad) ;
compute_derive_lie(vv, resu) ;
return resu ;
}
const Vector_divfree Vector::curl() const {
if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
const Metric_flat& metc = mp->flat_met_cart() ;
Vector_divfree resu(*mp, mp->get_bvect_cart(), metc) ;
resu.set(1)= cmp[3]->dsdy() - cmp[2]->dsdz();
resu.set(2)= cmp[1]->dsdz() - cmp[3]->dsdx();
resu.set(3)= cmp[2]->dsdx() - cmp[1]->dsdy();
resu.pseudo_spectral_base();
return resu ;
}
else {
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
const Metric_flat& mets = mp->flat_met_spher() ;
Vector_divfree resu(*mp, mp->get_bvect_spher(), mets);
Scalar tmp = *cmp[3] ;
tmp.div_tant();
tmp += cmp[3]->dsdt();
tmp.div_r();
resu.set(1) = tmp - cmp[2]->srstdsdp() ;
tmp = *cmp[3] ;
tmp.mult_r();
tmp = tmp.dsdr();
tmp.div_r();
resu.set(2) = cmp[1]->srstdsdp() - tmp ;
tmp = *cmp[2];
tmp.mult_r();
resu.set(3) = tmp.dsdr() - cmp[1]->dsdt() ;
resu.set(3).div_r();
resu.pseudo_spectral_base();
return resu ;
}
}
Sym_tensor Vector::ope_killing(const Metric& gam) const {
Sym_tensor resu(*mp, type_indice(0), *triad) ;
if (type_indice(0) == CON ) {
for (int i=1; i<=3; i++) {
for (int j=i; j<=3; j++) {
resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i) ;
}
}
}
else {
for (int i=1; i<=3; i++) {
for (int j=i; j<=3; j++) {
resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i) ;
}
}
}
return resu ;
}
Sym_tensor Vector::ope_killing_conf(const Metric& gam) const {
Sym_tensor resu(*mp, type_indice(0), *triad) ;
if (type_indice(0) == CON ) {
for (int i=1; i<=3; i++) {
for (int j=i; j<=3; j++) {
resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i)
- 0.6666666666666666* divergence(gam)
* gam.con()(i,j) ;
}
}
}
else {
for (int i=1; i<=3; i++) {
for (int j=i; j<=3; j++) {
resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i)
- 0.6666666666666666* derive_con(gam).trace()
* gam.cov()(i,j) ;
}
}
}
return resu ;
}
const Scalar& Vector::potential(const Metric& metre) const {
set_dependance(metre) ;
int j = get_place_met(metre) ;
assert ((j>=0) && (j<N_MET_MAX)) ;
if (p_potential[j] == 0x0) {
decompose_div(metre) ;
}
return *p_potential[j] ;
}
const Vector_divfree& Vector::div_free(const Metric& metre) const {
set_dependance(metre) ;
int j = get_place_met(metre) ;
assert ((j>=0) && (j<N_MET_MAX)) ;
if (p_div_free[j] == 0x0) {
decompose_div(metre) ;
}
return *p_div_free[j] ;
}
void Vector::decompose_div(const Metric& metre) const {
assert( type_indice(0) == CON ) ; //Only for contravariant vectors...
set_dependance(metre) ;
int ind = get_place_met(metre) ;
assert ((ind>=0) && (ind<N_MET_MAX)) ;
if ( (p_potential[ind] != 0x0) && (p_div_free[ind] != 0x0) )
return ; // Nothing to do ...
else {
if (p_div_free[ind] != 0x0)
delete p_div_free[ind] ;
if (p_potential[ind] != 0x0)
delete p_potential[ind] ;
const Mg3d* mmg = mp->get_mg() ;
int nz = mmg->get_nzone() ;
int dzp = cmp[0]->get_dzpuis() ;
#ifndef NDEBUG
bool dz_zero = cmp[0]->check_dzpuis(0) ;
bool dz_four = cmp[0]->check_dzpuis(4) ;
#endif
assert( dz_zero || dz_four) ;
assert (cmp[1]->check_dzpuis(dzp)) ;
assert (cmp[2]->check_dzpuis(dzp)) ;
Scalar dive = divergence(metre) ;
if (dive.get_etat() == ETATZERO) {
p_potential[ind] = new Scalar(*mp) ;
p_potential[ind]->set_etat_zero() ;
p_potential[ind]->set_dzpuis(dzp) ;
}
else {
//No matching of the solution at this point
p_potential[ind] = new Scalar(dive.poisson()) ;
if (dzp == 4) {
assert (mmg->get_type_r(nz-1) == UNSURR) ;
// Let's now do the matching ... in the case dzpuis = 4
const Map_af* mapping = dynamic_cast<const Map_af*>(mp) ;
assert (mapping != 0x0) ;
Valeur& val = p_potential[ind]->set_spectral_va() ;
val.ylm() ;
Mtbl_cf& mtc = *val.c_cf ;
if (nz > 1) {
int np = mmg->get_np(0) ;
int nt = mmg->get_nt(0) ;
#ifndef NDEBUG
for (int lz=0; lz<nz; lz++) {
assert (mmg->get_np(lz) == np) ;
assert (mmg->get_nt(lz) == nt) ;
}
#endif
int lmax = 0 ;
for (int k=0; k<np+1; k++)
for (int j=0; j<nt; j++)
if ( nullite_plm(j, nt, k, np, val.base)) {
int m_quant, l_quant, base_r ;
donne_lm (nz, 0, j, k, val.base, m_quant,
l_quant, base_r) ;
lmax = (l_quant > lmax ? l_quant : lmax) ;
}
Scalar** ri = new Scalar*[lmax+1] ;
Scalar** rmi = new Scalar*[lmax+1] ;
Scalar erre(*mp) ;
erre = mp->r ;
for (int l=0; l<=lmax; l++) {
ri[l] = new Scalar(*mp) ;
rmi[l] = new Scalar(*mp) ;
if (l == 0) *(ri[l]) = 1. ;
else *(ri[l]) = pow(erre, l) ;
ri[l]->annule_domain(nz - 1) ;
ri[l]->std_spectral_base() ; //exact base in r will be set later
if (l==2) *(rmi[l]) = 1 ;
else *(rmi[l]) = pow(erre, 2-l) ;
rmi[l]->annule(0,nz-2) ;
Scalar tmp = pow(erre, -l-1) ;
tmp.annule_domain(nz-1) ;
tmp.annule_domain(0) ;
*(rmi[l]) += tmp ;
rmi[l]->set_dzpuis(3) ;
rmi[l]->std_spectral_base() ;//exact base in r will be set later
}
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
if ( nullite_plm(j, nt, k, np, val.base)) {
int m_quant, l_quant, base_r, l, c ;
donne_lm (nz, 0, j, k, val.base, m_quant, l_quant, base_r) ;
bool lzero = (l_quant == 0) || (l_quant == 1) ;
int nb_hom_ced = (l_quant < 2 ? 0 : 1) ;
int taille = 2*(nz-1) - 1 + nb_hom_ced ;
Tbl deuz(taille) ;
deuz.set_etat_qcq() ;
Matrice systeme(taille,taille) ;
systeme.set_etat_qcq() ;
for (l=0; l<taille; l++)
for (c=0; c<taille; c++) systeme.set(l,c) = 0 ;
for (l=0; l<taille; l++) deuz.set(l) = 0 ;
//---------
// Nucleus
//---------
assert(mmg->get_type_r(0) == RARE) ;
assert ((base_r == R_CHEBP)||(base_r == R_CHEBI)) ;
ri[l_quant]->set_spectral_va().set_base_r(0, base_r) ;
double alpha = mapping->get_alpha()[0] ;
int nr = mmg->get_nr(0) ;
Tbl partn(nr) ;
partn.set_etat_qcq() ;
for (int i=0; i<nr; i++)
partn.set(i) = mtc(0, k, j, i) ;
l=0 ; c=0 ;
systeme.set(l,c) += pow(alpha, double(l_quant)) ;
deuz.set(l) -= val1_dern_1d(0, partn, base_r) ;
l++ ;
if (!lzero) {
systeme.set(l,c) += l_quant * pow(alpha, double(l_quant - 1)) ;
deuz.set(l) -= val1_dern_1d(1, partn, base_r) / alpha ;
}
//----------
// Shells
//----------
for (int lz=1; lz<nz-1; lz++) {
alpha = mapping->get_alpha()[lz] ;
double beta = mapping->get_beta()[lz] ;
double xm1 = beta - alpha ;
double xp1 = alpha + beta ;
nr = mmg->get_nr(lz) ;
Tbl parts(nr) ;
parts.set_etat_qcq() ;
for (int i=0; i<nr; i++)
parts.set(i) = mtc(lz, k, j, i) ;
//Function at x = -1
l-- ;
c = 2*lz - 1 ;
systeme.set(l,c) -= pow(xm1, double(l_quant)) ;
c++;
systeme.set(l,c) -= pow(xm1, double(-l_quant-1)) ;
deuz.set(l) += valm1_dern_1d(0, parts, R_CHEB) ;
if ((lz>1) || (!lzero)) {
//First derivative at x=-1
l++ ;
c-- ;
systeme.set(l,c) -= l_quant*pow(xm1, double(l_quant - 1)) ;
c++;
systeme.set(l,c) -= (-l_quant - 1)*
pow(xm1, double(-l_quant - 2)) ;
deuz.set(l) += valm1_dern_1d(1, parts, R_CHEB) / alpha ;
}
//Function at x = 1
l++ ;
c-- ;
systeme.set(l,c) += pow(xp1, double(l_quant)) ;
c++;
systeme.set(l,c) += pow(xp1, double(-l_quant-1)) ;
deuz.set(l) -= val1_dern_1d(0, parts, R_CHEB) ;
//First derivative at x = 1
l++ ;
c-- ;
systeme.set(l,c) += l_quant*pow(xp1, double(l_quant - 1)) ;
c++;
systeme.set(l,c) += (-l_quant - 1)*
pow(xp1, double(-l_quant - 2)) ;
deuz.set(l) -= val1_dern_1d(1, parts, R_CHEB) / alpha ;
} //End of the loop on different shells
//-------------------------------
// Compactified external domain
//-------------------------------
assert(mmg->get_type_r(nz-1) == UNSURR) ;
nr = mmg->get_nr(nz-1) ;
Tbl partc(nr) ;
partc.set_etat_qcq() ;
for (int i=0; i<nr; i++)
partc.set(i) = mtc(nz-1, k, j, i) ;
alpha = mapping->get_alpha()[nz-1] ;
double beta = mapping->get_beta()[nz-1] ;
double xm1 = beta - alpha ; // 1 / r_left
double part0, part1 ;
part0 = valm1_dern_1d(0, partc, R_CHEB) ;
part1 = xm1*valm1_dern_1d(1, partc, R_CHEB) / alpha ;
assert (p_potential[ind]->get_dzpuis() == 3) ;
//Function at x = -1
l--;
if (!lzero) {
c++;
systeme.set(l,c) -= pow(xm1, double(l_quant+1)) ;
}
deuz.set(l) += part0*xm1*xm1*xm1 ;
// First derivative at x = -1
l++ ;
if (!lzero) {
systeme.set(l,c) -= (-l_quant - 1)*
pow(xm1, double(l_quant + 2)) ;
}
if ( (nz > 2) || (!lzero))
deuz.set(l) += -xm1*xm1*xm1*xm1*(3*part0 + part1) ;
//--------------------------------------
// Solution of the linear system
//--------------------------------------
int inf = 1 + nb_hom_ced;
int sup = 3 - nb_hom_ced;
systeme.set_band(sup, inf) ;
systeme.set_lu() ;
Tbl facteur(systeme.inverse(deuz)) ;
ri[l_quant]->set_spectral_va().coef() ;
rmi[l_quant]->set_spectral_va().coef() ;
//Linear combination in the nucleus
nr = mmg->get_nr(0) ;
for (int i=0; i<nr; i++)
mtc.set(0, k, j, i) +=
facteur(0)*(*(ri[l_quant]->get_spectral_va().c_cf))(0, 0, 0, i) ;
//Linear combination in the shells
for (int lz=1; lz<nz-1; lz++) {
nr = mmg->get_nr(lz) ;
for (int i=0; i<nr; i++)
mtc.set(lz, k, j, i) += facteur(2*lz - 1)*
(*(ri[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
for (int i=0; i<nr; i++)
mtc.set(lz, k, j, i) += facteur(2*lz)*
(*(rmi[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
}
//Linear combination in the CED
nr = mmg->get_nr(nz-1) ;
if (!lzero) {
for (int i=0; i<nr; i++)
mtc.set(nz-1, k, j, i) +=
facteur(taille - 1)*
(*(rmi[l_quant]->get_spectral_va().c_cf))(nz-1, 0, 0, i) ;
}
} //End of nullite_plm ...
} //End of j/theta loop
} //End of k/phi loop
for (int l=0; l<=lmax; l++) {
delete ri[l] ;
delete rmi[l] ;
}
delete [] ri ;
delete [] rmi ;
} //End of the case of more than one domain
val.ylm_i() ;
} //End of the case dzp = 4
}
p_div_free[ind] = new Vector_divfree(*mp, *triad, metre) ;
Vector gradient = p_potential[ind]->derive_con(metre) ;
if (dzp != 4) gradient.dec_dzpuis(2) ;
*p_div_free[ind] = ( *this - gradient ) ;
}
}
double Vector::flux(double radius, const Metric& met) const {
assert(type_indice(0) == CON) ;
const Map_af* mp_af = dynamic_cast<const Map_af*>(mp) ;
if (mp_af == 0x0) {
cerr <<
"Vector::flux : the case of a mapping outside the class Map_af\n"
<< " is not implemented yet !" << endl ;
abort() ;
}
const Metric_flat* ff = dynamic_cast<const Metric_flat*>(&met) ;
if (ff == 0x0) {
cerr <<
"Vector::flux : the case of a non flat metric is not implemented yet !"
<< endl ;
abort() ;
}
const Base_vect_cart* bcart = dynamic_cast<const Base_vect_cart*>(triad) ;
Vector* vspher = 0x0 ;
if (bcart != 0x0) { // switch to spherical components:
vspher = new Vector(*this) ;
vspher->change_triad(mp->get_bvect_spher()) ;
}
else assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
const Vector* vv = (bcart != 0x0) ? vspher : this ;
const Scalar& vr = vv->operator()(1) ;
double resu ;
if (radius == __infinity) {
resu = mp_af->integrale_surface_infini(vr) ;
}
else {
resu = mp_af->integrale_surface(vr, radius) ;
}
return resu ;
}
void Vector::exponential_filter_r(int lzmin, int lzmax, int p,
double alpha) {
if( triad->identify() == (mp->get_bvect_cart()).identify() )
for (int i=0; i<n_comp; i++)
cmp[i]->exponential_filter_r(lzmin, lzmax, p, alpha) ;
else {
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
Scalar vr_tmp = operator()(1) ;
vr_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
Scalar eta_tmp = eta() ;
eta_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
Scalar mu_tmp = mu() ;
mu_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
}
}
void Vector::exponential_filter_ylm(int lzmin, int lzmax, int p,
double alpha) {
if( triad->identify() == (mp->get_bvect_cart()).identify() )
for (int i=0; i<n_comp; i++)
cmp[i]->exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
else {
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
Scalar vr_tmp = operator()(1) ;
vr_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
Scalar eta_tmp = eta() ;
eta_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
Scalar mu_tmp = mu() ;
mu_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
}
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_vect_cart
Cartesian vectorial bases (triads).
Definition base_vect.h:201
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
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::Map_af::integrale_surface_infini
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Definition map_af_integ_surf.C:196
Lorene::Map_af::integrale_surface
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Definition map_af_integ_surf.C:90
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::set_band
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition matrice.C:364
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
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_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::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
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::div_r
void div_r()
Division by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:133
Lorene::Scalar::div_tant
void div_tant()
Division by .
Definition scalar_th_manip.C:111
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::exponential_filter_ylm
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
Definition scalar_exp_filter.C:142
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
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::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:208
Lorene::Scalar::exponential_filter_r
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
Definition scalar_exp_filter.C:60
Lorene::Tbl
Basic array class.
Definition tbl.h:161
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
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Valeur::set_base_r
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition valeur.C:836
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Vector_divfree
Divergence-free vectors.
Definition vector.h:724
Lorene::Vector::operator()
const Scalar & operator()(int) const
Readonly access to a component.
Definition vector.C:305
Lorene::Vector::ope_killing
Sym_tensor ope_killing(const Metric &gam) const
Computes the Killing operator associated with a given metric.
Definition vector.C:438
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::decompose_div
void decompose_div(const Metric &) const
Makes the Helmholtz decomposition (see documentation of p_potential ) of this with respect to a given...
Definition vector.C:516
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::Vector::p_potential
Scalar * p_potential[N_MET_MAX]
The potential giving the gradient part in the Helmholtz decomposition of any 3D vector .
Definition vector.h:198
Lorene::Vector::del_deriv
virtual void del_deriv() const
Deletes the derived quantities.
Definition vector.C:219
Lorene::Vector::mu
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Definition vector_etamu.C:98
Lorene::Vector::del_derive_met
virtual void del_derive_met(int) const
Logical destructor of the derivatives depending on the i-th element of met_depend in the class Vector...
Definition vector.C:242
Lorene::Vector::flux
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition vector.C:807
Lorene::Vector::Vector
Vector(const Map &map, int tipe, const Base_vect &triad_i)
Standard constructor.
Definition vector.C:156
Lorene::Vector::p_div_free
Vector_divfree * p_div_free[N_MET_MAX]
The divergence-free vector of the Helmholtz decomposition of any 3D vector .
Definition vector.h:205
Lorene::Vector::curl
const Vector_divfree curl() const
The curl of this with respect to a (flat) Metric .
Definition vector.C:402
Lorene::Vector::p_A
Scalar * p_A
Field defined by.
Definition vector.h:241
Lorene::Vector::eta
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
Definition vector_etamu.C:66
Lorene::Vector::exponential_filter_r
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies exponential filters to all components (see Scalar::exponential_filter_r ).
Definition vector.C:850
Lorene::Vector::div_free
const Vector_divfree & div_free(const Metric &) const
Returns the div-free vector in the Helmholtz decomposition.
Definition vector.C:504
Lorene::Vector::set_der_0x0
void set_der_0x0() const
Sets the pointers on derived quantities to 0x0.
Definition vector.C:232
Lorene::Vector::~Vector
virtual ~Vector()
Destructor.
Definition vector.C:208
Lorene::Vector::pseudo_spectral_base
virtual void pseudo_spectral_base()
Sets the standard spectal bases of decomposition for each component for a pseudo_vector.
Definition vector.C:346
Lorene::Vector::set_der_met_0x0
void set_der_met_0x0(int) const
Sets all the i-th components of met_depend in the class Vector (p_potential , etc....
Definition vector.C:258
Lorene::Vector::p_mu
Scalar * p_mu
Field such that the angular components of the vector are written:
Definition vector.h:233
Lorene::Vector::p_eta
Scalar * p_eta
Field such that the angular components of the vector are written:
Definition vector.h:219
Lorene::Vector::ope_killing_conf
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:75
Lorene::Vector::operator=
virtual void operator=(const Vector &a)
Assignment from a Vector.
Definition vector.C:267
Lorene::Vector::potential
const Scalar & potential(const Metric &) const
Returns the potential in the Helmholtz decomposition.
Definition vector.C:492
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
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:189
Lorene::Vector::exponential_filter_ylm
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies exponential filters to all components (see Scalar::exponential_filter_ylm ).
Definition vector.C:867
Lorene::Vector::derive_lie
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition vector.C:392
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene::Tensor::del_derive_met
virtual void del_derive_met(int) const
Logical destructor of the derivatives depending on the i-th element of met_depend .
Definition tensor.C:414
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
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::get_place_met
int get_place_met(const Metric &) const
Returns the position of the pointer on metre in the array met_depend .
Definition tensor.C:443
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
Lorene::Tensor::met_depend
const Metric * met_depend[N_MET_MAX]
Array on the Metric 's which were used to compute derived quantities, like p_derive_cov ,...
Definition tensor.h:327
Lorene::Tensor::Tensor
Tensor(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition tensor.C:202
Lorene::Tensor::valence
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...).
Definition tensor.h:298
Lorene::Tensor::set_dependance
void set_dependance(const Metric &) const
To be used to describe the fact that the derivatives members have been calculated with met .
Definition tensor.C:453
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
Lorene::Tensor::compute_derive_lie
void compute_derive_lie(const Vector &v, Tensor &resu) const
Computes the Lie derivative of this with respect to some vector field v (protected method; the public...
Definition tensor_calculus.C:339
Lorene::Tensor::n_comp
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene::Tensor::del_deriv
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
Lorene::Tensor::divergence
const Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
Definition tensor.C:1055
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Lorene::Tensor::type_indice
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a cov...
Definition tensor.h:310
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139