tensor.C
/*
* Methods of class Tensor
*
* (see file tensor.h for documentation)
*
*/
/*
* Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
*
* 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 tensor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $" ;
/*
* $Id: tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $
* $Log: tensor.C,v $
* Revision 1.43 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.42 2014/10/06 15:13:20 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.41 2013/06/05 15:43:49 j_novak
* Suppression of leg_spectral_base()
*
* Revision 1.40 2013/01/11 15:44:54 j_novak
* Addition of Legendre bases (part 2).
*
* Revision 1.39 2007/12/21 16:07:08 j_novak
* Methods to filter Tensor, Vector and Sym_tensor objects.
*
* Revision 1.38 2005/10/25 08:56:38 p_grandclement
* addition of std_spectral_base in the case of odd functions near the origin
*
* Revision 1.37 2005/05/18 11:45:44 j_novak
* Added del_deriv() calls at the end of inc/dec_dzpuis.
*
* Revision 1.36 2004/07/08 12:21:53 j_novak
* Replaced tensor::annule_extern_c2 with tensor::annule_extern_cn for a
* more general transition.
*
* Revision 1.35 2004/06/17 06:55:07 e_gourgoulhon
* Added method annule_extern_c2.
*
* Revision 1.34 2004/03/04 09:47:51 e_gourgoulhon
* Method annule_domain(int, int): added call to virtual function
* del_deriv() at the end !
*
* Revision 1.33 2004/02/27 21:14:27 e_gourgoulhon
* Modif of method derive_con to create proper type of the result.
*
* Revision 1.32 2004/02/19 22:10:14 e_gourgoulhon
* Added argument "comment" in method spectral_display.
*
* Revision 1.31 2004/02/05 15:03:47 e_gourgoulhon
* Corrected bug in method derive_con().
*
* Revision 1.30 2004/01/27 13:05:11 j_novak
* Removed the method Tensor::mult_r_ced()
*
* Revision 1.29 2004/01/19 16:32:13 e_gourgoulhon
* Added operator()(int, int, int, int) and set(int, int, int, int)
* for direct access to components of valence 4 tensors.
*
* Revision 1.28 2004/01/04 20:55:23 e_gourgoulhon
* Method spectral_display(): added printing of type of class (through typeid).
*
* Revision 1.27 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.26 2003/12/27 15:01:50 e_gourgoulhon
* Method derive_con(): the position of the derivation index has
* been changed from the first one to the last one.
*
* Revision 1.25 2003/11/03 22:34:41 e_gourgoulhon
* Method dec_dzpuis: changed the name of argument dec --> decrem
* (in order not to shadow some globally defined dec).
*
* Revision 1.24 2003/10/29 11:02:54 e_gourgoulhon
* Functions dec_dzpuis and inc_dzpuis have now an integer argument to
* specify by which amount dzpuis is to be increased.
* Accordingly methods dec2_dzpuis and inc2_dzpuis have been suppressed
*
* Revision 1.23 2003/10/27 10:49:48 e_gourgoulhon
* Slightly modified operator<<.
*
* Revision 1.22 2003/10/19 19:57:00 e_gourgoulhon
* -- Added new method spectral_display
* -- slightly rearranged the operator<<
*
* Revision 1.21 2003/10/16 15:27:46 e_gourgoulhon
* Name of method annule(int ) changed to annule_domain(int ).
*
* Revision 1.20 2003/10/16 14:21:36 j_novak
* The calculation of the divergence of a Tensor is now possible.
*
* Revision 1.19 2003/10/13 13:52:40 j_novak
* Better managment of derived quantities.
*
* Revision 1.18 2003/10/11 16:47:10 e_gourgoulhon
* Suppressed the call to Ibtl::set_etat_qcq() after the construction
* of the Itbl's, thanks to the new property of the Itbl class.
*
* Revision 1.17 2003/10/11 14:48:40 e_gourgoulhon
* Line 344: suppressed assert(resu == -1)
* and added a break command to quit the loop.
*
* Revision 1.16 2003/10/08 14:24:09 j_novak
* replaced mult_r_zec with mult_r_ced
*
* Revision 1.15 2003/10/07 15:25:38 e_gourgoulhon
* Added call to del_derive_met in del_deriv().
*
* Revision 1.14 2003/10/07 09:10:00 j_novak
* Use of ::contract instead of up()
*
* Revision 1.13 2003/10/06 20:51:43 e_gourgoulhon
* In methods set: changed name "indices" to "idx" to avoid shadowing
* of class member.
*
* Revision 1.12 2003/10/06 16:17:31 j_novak
* Calculation of contravariant derivative and Ricci scalar.
*
* Revision 1.11 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.10 2003/10/05 21:11:22 e_gourgoulhon
* - Added method std_spectral_base().
* - Removed method change_triad() from this file.
*
* Revision 1.9 2003/10/03 15:09:38 j_novak
* Improved display
*
* Revision 1.8 2003/10/03 11:21:48 j_novak
* More methods for the class Metric
*
* Revision 1.7 2003/10/01 11:56:31 e_gourgoulhon
* Corrected error: '=' replaced by '==' in two assert tests.
*
* Revision 1.6 2003/09/30 08:38:23 j_novak
* added a header typeinfo
*
* Revision 1.5 2003/09/29 12:52:57 j_novak
* Methods for changing the triad are implemented.
*
* Revision 1.4 2003/09/26 14:33:53 j_novak
* Arithmetic functions for the class Tensor
*
* Revision 1.3 2003/09/24 15:10:54 j_novak
* Suppression of the etat flag in class Tensor (still present in Scalar)
*
* Revision 1.2 2003/09/23 08:52:17 e_gourgoulhon
* new version
*
* Revision 1.1 2003/09/22 12:52:51 e_gourgoulhon
* First version: not ready yet !
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
// Headers Lorene
#include "metric.h"
#include "utilitaires.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Tensor::Tensor(const Map& map, int val, const Itbl& tipe,
const Base_vect& triad_i)
: mp(&map), valence(val), triad(&triad_i), type_indice(tipe),
n_comp(int(pow(3., val))){
// Des verifs :
assert (valence >= 0) ;
assert (tipe.get_ndim() == 1) ;
assert (valence == tipe.get_dim(0)) ;
for (int i=0 ; i<valence ; i++)
assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(map) ;
set_der_0x0() ;
}
// Standard constructor with the triad passed as a pointer
// -------------------------------------------------------
Tensor::Tensor(const Map& map, int val, const Itbl& tipe,
const Base_vect* triad_i)
: mp(&map), valence(val), triad(triad_i), type_indice(tipe),
n_comp(int(pow(3., val))){
// Des verifs :
assert (valence >= 0) ;
assert (tipe.get_ndim() == 1) ;
assert (valence == tipe.get_dim(0)) ;
for (int i=0 ; i<valence ; i++)
assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(map) ;
set_der_0x0() ;
}
// Standard constructor when all the indices are of the same type
// --------------------------------------------------------------
Tensor::Tensor(const Map& map, int val, int tipe, const Base_vect& triad_i)
: mp(&map), valence(val), triad(&triad_i), type_indice(val),
n_comp(int(pow(3., val))){
// Des verifs :
assert (valence >= 0) ;
assert ((tipe == COV) || (tipe == CON)) ;
type_indice = tipe ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(map) ;
set_der_0x0() ;
}
// Copy constructor
// ----------------
Tensor::Tensor (const Tensor& source) :
mp(source.mp), valence(source.valence), triad(source.triad),
type_indice(source.type_indice) {
n_comp = int(pow(3., valence)) ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++) {
// the following writing takes into account the case where
// source belongs to a derived class of Tensor with a different
// storage of components :
int place_source = source.position(indices(i)) ;
cmp[i] = new Scalar(*source.cmp[place_source]) ;
}
set_der_0x0() ;
}
// Constructor from a file
// -----------------------
Tensor::Tensor(const Map& mapping, const Base_vect& triad_i, FILE* fd)
: mp(&mapping), triad(&triad_i), type_indice(fd){
fread_be(&valence, sizeof(int), 1, fd) ;
if (valence != 0) {
Base_vect* triad_fich = Base_vect::bvect_from_file(fd) ;
assert( *triad_fich == *triad) ;
delete triad_fich ;
}
else{
triad = 0x0 ;
}
fread_be(&n_comp, sizeof(int), 1, fd) ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(*mp, *(mp->get_mg()), fd) ;
set_der_0x0() ;
}
// Constructor for a scalar field: to be used by the derived
// class {\tt Scalar}
//-----------------------------------------------------------
Tensor::Tensor(const Map& map) : mp(&map), valence(0), triad(0x0),
type_indice(0), n_comp(1) {
cmp = new Scalar*[n_comp] ;
cmp[0] = 0x0 ;
set_der_0x0() ;
}
// Constructor used by the derived classes
// ---------------------------------------
Tensor::Tensor (const Map& map, int val, const Itbl& tipe, int compo,
const Base_vect& triad_i) :
mp(&map), valence(val), triad(&triad_i), type_indice(tipe), n_comp(compo)
{
// Des verifs :
assert (valence >= 0) ;
assert (tipe.get_ndim() == 1) ;
assert (n_comp > 0) ;
assert (valence == tipe.get_dim(0)) ;
for (int i=0 ; i<valence ; i++)
assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(map) ;
set_der_0x0() ;
}
// Constructor used by the derived classes when all the indices are of
// the same type.
// -------------------------------------------------------------------
Tensor::Tensor (const Map& map, int val, int tipe, int compo,
const Base_vect& triad_i) :
mp(&map), valence(val), triad(&triad_i), type_indice(val), n_comp(compo)
{
// Des verifs :
assert (valence >= 0) ;
assert (n_comp >= 0) ;
assert ((tipe == COV) || (tipe == CON)) ;
type_indice = tipe ;
cmp = new Scalar*[n_comp] ;
for (int i=0 ; i<n_comp ; i++)
cmp[i] = new Scalar(map) ;
set_der_0x0() ;
}
//--------------//
// Destructor //
//--------------//
Tensor::~Tensor () {
Tensor::del_deriv() ;
for (int i=0 ; i<n_comp ; i++) {
if (cmp[i] != 0x0)
delete cmp[i] ;
}
delete [] cmp ;
}
void Tensor::del_deriv() const {
for (int i=0; i<N_MET_MAX; i++)
del_derive_met(i) ;
set_der_0x0() ;
}
void Tensor::set_der_0x0() const {
for (int i=0; i<N_MET_MAX; i++)
set_der_met_0x0(i) ;
}
void Tensor::del_derive_met(int j) const {
assert( (j>=0) && (j<N_MET_MAX) ) ;
if (met_depend[j] != 0x0) {
for (int i=0 ; i<N_TENSOR_DEPEND ; i++)
if (met_depend[j]->tensor_depend[i] == this)
met_depend[j]->tensor_depend[i] = 0x0 ;
if (p_derive_cov[j] != 0x0)
delete p_derive_cov[j] ;
if (p_derive_con[j] != 0x0)
delete p_derive_con[j] ;
if (p_divergence[j] != 0x0)
delete p_divergence[j] ;
set_der_met_0x0(j) ;
}
}
void Tensor::set_der_met_0x0(int i) const {
assert( (i>=0) && (i<N_MET_MAX) ) ;
met_depend[i] = 0x0 ;
p_derive_cov[i] = 0x0 ;
p_derive_con[i] = 0x0 ;
p_divergence[i] = 0x0 ;
}
int Tensor::get_place_met(const Metric& metre) const {
int resu = -1 ;
for (int i=0; i<N_MET_MAX; i++)
if (met_depend[i] == &metre) {
resu = i ;
break ;
}
return resu ;
}
void Tensor::set_dependance (const Metric& met) const {
int nmet = 0 ;
bool deja = false ;
for (int i=0; i<N_MET_MAX; i++) {
if (met_depend[i] == &met) deja = true ;
if ((!deja) && (met_depend[i] != 0x0)) nmet++ ;
}
if (nmet == N_MET_MAX) {
cout << "Too many metrics in Tensor::set_dependances" << endl ;
abort() ;
}
if (!deja) {
int conte = 0 ;
while ((conte < N_TENSOR_DEPEND) && (met.tensor_depend[conte] != 0x0))
conte ++ ;
if (conte == N_TENSOR_DEPEND) {
cout << "Too many dependancies in Tensor::set_dependances " << endl ;
abort() ;
}
else {
met.tensor_depend[conte] = this ;
met_depend[nmet] = &met ;
}
}
}
void Tensor::set_etat_qcq() {
del_deriv() ;
for (int i=0 ; i<n_comp ; i++) {
cmp[i]->set_etat_qcq() ;
}
}
void Tensor::set_etat_nondef() {
del_deriv() ;
for (int i=0 ; i<n_comp ; i++) {
cmp[i]->set_etat_nondef() ;
}
}
void Tensor::set_etat_zero() {
del_deriv() ;
for (int i=0 ; i<n_comp ; i++) {
cmp[i]->set_etat_zero() ;
}
}
// Allocates everything
// --------------------
void Tensor::allocate_all() {
del_deriv() ;
for (int i=0 ; i<n_comp ; i++) {
cmp[i]->allocate_all() ;
}
}
void Tensor::set_triad(const Base_vect& bi) {
triad = &bi ;
}
int Tensor::position (const Itbl& idx) const {
assert (idx.get_ndim() == 1) ;
assert (idx.get_dim(0) == valence) ;
for (int i=0 ; i<valence ; i++)
assert ((idx(i)>=1) && (idx(i)<=3)) ;
int res = 0 ;
for (int i=0 ; i<valence ; i++)
res = 3*res+(idx(i)-1) ;
return res;
}
Itbl Tensor::indices (int place) const {
assert ((place >= 0) && (place < n_comp)) ;
Itbl res(valence) ;
for (int i=valence-1 ; i>=0 ; i--) {
res.set(i) = div(place, 3).rem ;
place = int((place-res(i))/3) ;
res.set(i)++ ;
}
return res ;
}
void Tensor::operator=(const Tensor& t) {
assert (valence == t.valence) ;
triad = t.triad ;
for (int i=0 ; i<valence ; i++)
assert(t.type_indice(i) == type_indice(i)) ;
for (int i=0 ; i<n_comp ; i++) {
int place_t = t.position(indices(i)) ;
*cmp[i] = *t.cmp[place_t] ;
}
del_deriv() ;
}
void Tensor::operator+=(const Tensor& t) {
assert (valence == t.valence) ;
assert (triad == t.triad) ;
for (int i=0 ; i<valence ; i++)
assert(t.type_indice(i) == type_indice(i)) ;
for (int i=0 ; i<n_comp ; i++) {
int place_t = t.position(indices(i)) ;
*cmp[i] += *t.cmp[place_t] ;
}
del_deriv() ;
}
void Tensor::operator-=(const Tensor& t) {
assert (valence == t.valence) ;
assert (triad == t.triad) ;
for (int i=0 ; i<valence ; i++)
assert(t.type_indice(i) == type_indice(i)) ;
for (int i=0 ; i<n_comp ; i++) {
int place_t = t.position(indices(i)) ;
*cmp[i] -= *t.cmp[place_t] ;
}
del_deriv() ;
}
// Affectation d'un tenseur d'ordre 2 :
Scalar& Tensor::set(int ind1, int ind2) {
assert (valence == 2) ;
Itbl ind (valence) ;
ind.set(0) = ind1 ;
ind.set(1) = ind2 ;
int place = position(ind) ;
del_deriv() ;
return *cmp[place] ;
}
// Affectation d'un tenseur d'ordre 3 :
Scalar& Tensor::set(int ind1, int ind2, int ind3) {
assert (valence == 3) ;
Itbl idx(valence) ;
idx.set(0) = ind1 ;
idx.set(1) = ind2 ;
idx.set(2) = ind3 ;
int place = position(idx) ;
del_deriv() ;
return *cmp[place] ;
}
// Affectation d'un tenseur d'ordre 4 :
Scalar& Tensor::set(int ind1, int ind2, int ind3, int ind4) {
assert (valence == 4) ;
Itbl idx(valence) ;
idx.set(0) = ind1 ;
idx.set(1) = ind2 ;
idx.set(2) = ind3 ;
idx.set(3) = ind4 ;
int place = position(idx) ;
del_deriv() ;
return *cmp[place] ;
}
// Affectation cas general
Scalar& Tensor::set(const Itbl& idx) {
assert (idx.get_ndim() == 1) ;
assert (idx.get_dim(0) == valence) ;
int place = position(idx) ;
del_deriv() ;
return *cmp[place] ;
}
// Annulation dans des domaines
void Tensor::annule_domain(int l) {
annule(l, l) ;
}
void Tensor::annule(int l_min, int l_max) {
// Cas particulier: annulation globale :
if ( (l_min == 0) && (l_max == mp->get_mg()->get_nzone()-1) ) {
set_etat_zero() ;
return ;
}
// Annulation des composantes:
for (int i=0 ; i<n_comp ; i++) {
cmp[i]->annule(l_min, l_max) ;
}
// The derived members are no longer up to date:
del_deriv() ;
}
void Tensor::annule_extern_cn(int lrac, int deg) {
// Not applicable in the nucleus nor the CED:
assert( mp->get_mg()->get_type_r(lrac) == FIN ) ;
int nz = mp->get_mg()->get_nzone() ;
#ifndef NDEBUG
if ((2*deg+1) >= mp->get_mg()->get_nr(lrac))
cout << "Tensor::annule_extern_cn : \n"
<< "WARNING!! \n"
<< "The number of coefficients in r is too low \n"
<< "to do a clean matching..." << endl ;
#endif
// Boundary of domain lrac
double r_min = mp->val_r(lrac, -1., 0., 0.) ;
double r_max = mp->val_r(lrac, 1., 0., 0.) ;
//Definition of binomial coefficients array
Itbl binom(deg+1, deg+1) ;
binom.annule_hard() ;
binom.set(0,0) = 1 ;
for (int n=1; n<=deg; n++) {
binom.set(n,0) = 1 ;
for (int k=1; k<=n; k++)
binom.set(n,k) = binom(n-1, k) + binom(n-1, k-1) ;
}
// Coefficient of the second polynomial factor
Tbl coef(deg+1) ;
coef.set_etat_qcq() ;
coef.set(deg) = 1 ;
int sg = -1 ;
for (int i=deg-1; i>=0; i--) {
coef.set(i) = double(r_max*(i+1)*coef(i+1)
+ sg*binom(deg,i)*(2*deg+1)*pow(r_min,deg-i))
/ double(deg+i+1) ;
sg *= -1 ;
}
// Normalization to have 1 at r_min
double aa = coef(deg) ;
for (int i = deg-1; i>=0; i--)
aa = r_min*aa + coef(i) ;
aa *= pow(r_min - r_max, deg+1) ;
aa = 1/aa ;
Mtbl mr = mp->r ;
Tbl rr = mr(lrac) ;
Tbl poly(rr) ;
poly = coef(deg) ;
for (int i=deg-1; i>=0; i--)
poly = rr*poly + coef(i) ;
poly *= aa*pow((rr-r_max), deg+1) ;
Scalar rac(*mp) ;
rac.allocate_all() ;
for (int l=0; l<lrac; l++) rac.set_domain(l) = 1 ;
rac.set_domain(lrac) = poly ;
rac.annule(lrac+1,nz-1) ;
rac.std_spectral_base() ;
for (int ic=0; ic<n_comp; ic++) *(cmp[ic]) *= rac ;
del_deriv() ;
}
const Scalar& Tensor::operator()(int indice1, int indice2) const {
assert(valence == 2) ;
Itbl idx(2) ;
idx.set(0) = indice1 ;
idx.set(1) = indice2 ;
return *cmp[position(idx)] ;
}
const Scalar& Tensor::operator()(int indice1, int indice2, int indice3) const {
assert(valence == 3) ;
Itbl idx(3) ;
idx.set(0) = indice1 ;
idx.set(1) = indice2 ;
idx.set(2) = indice3 ;
return *cmp[position(idx)] ;
}
const Scalar& Tensor::operator()(int indice1, int indice2, int indice3,
int indice4) const {
assert(valence == 4) ;
Itbl idx(4) ;
idx.set(0) = indice1 ;
idx.set(1) = indice2 ;
idx.set(2) = indice3 ;
idx.set(3) = indice4 ;
return *cmp[position(idx)] ;
}
const Scalar& Tensor::operator()(const Itbl& ind) const {
assert (ind.get_ndim() == 1) ;
assert (ind.get_dim(0) == valence) ;
return *cmp[position(ind)] ;
}
// Gestion de la CED :
void Tensor::dec_dzpuis(int decrem) {
for (int i=0 ; i<n_comp ; i++)
cmp[i]->dec_dzpuis(decrem) ;
del_deriv() ;
}
void Tensor::inc_dzpuis(int inc) {
for (int i=0 ; i<n_comp ; i++)
cmp[i]->inc_dzpuis(inc) ;
del_deriv() ;
}
// Le cout :
ostream& operator<<(ostream& flux, const Tensor &source ) {
flux << '\n' ;
flux << "Lorene class : " << typeid(source).name()
<< " Valence : " << source.valence << '\n' ;
if (source.get_triad() != 0x0) {
flux << "Vectorial basis (triad) on which the components are defined :"
<< '\n' ;
flux << *(source.get_triad()) << '\n' ;
}
if (source.valence != 0) {
flux << "Type of the indices : " ;
for (int i=0 ; i<source.valence ; i++) {
flux << "index " << i << " : " ;
if (source.type_indice(i) == CON)
flux << " contravariant." << '\n' ;
else
flux << " covariant." << '\n' ;
if ( i < source.valence-1 ) flux << " " ;
}
flux << '\n' ;
}
for (int i=0 ; i<source.n_comp ; i++) {
if (source.valence == 0) {
flux <<
"===================== Scalar field ========================= \n" ;
}
else {
flux << "================ Component " ;
Itbl num_indices = source.indices(i) ;
for (int j=0 ; j<source.valence ; j++) {
flux << " " << num_indices(j) ;
}
flux << " ================ \n" ;
}
flux << '\n' ;
flux << *source.cmp[i] << '\n' ;
}
return flux ;
}
void Tensor::spectral_display(const char* comment,
double thres, int precis, ostream& ost) const {
if (comment != 0x0) {
ost << comment << " : " << endl ;
}
ost << "Lorene class : " << typeid(*this).name()
<< " Valence : " << valence << '\n' ;
for (int ic=0; ic<n_comp; ic++) {
if (valence == 0) {
ost <<
"===================== Scalar field ========================= \n" ;
}
else {
ost << "================ Component " ;
Itbl num_indices = indices(ic) ;
for (int j=0 ; j<valence ; j++) {
ost << " " << num_indices(j) ;
}
ost << " ================ \n" ;
}
ost << '\n' ;
cmp[ic]->spectral_display(0x0, thres, precis, ost) ;
ost << '\n' ;
}
}
void Tensor::sauve(FILE* fd) const {
type_indice.sauve(fd) ; // type des composantes
fwrite_be(&valence, sizeof(int), 1, fd) ; // la valence
if (valence != 0) {
triad->sauve(fd) ; // Vectorial basis
}
fwrite_be(&n_comp, sizeof(int), 1, fd) ; // nbre composantes
for (int i=0 ; i<n_comp ; i++)
cmp[i]->sauve(fd) ;
}
// Sets the standard spectal bases of decomposition for each component
void Tensor::std_spectral_base() {
switch (valence) {
case 0 : {
cmp[0]->std_spectral_base() ;
break ;
}
case 1 : {
cout <<
"Tensor::std_spectral_base: should not be called on a Tensor"
<< " of valence 1 but on a Vector !" << endl ;
abort() ;
break ;
}
case 2 : {
Base_val** bases = 0x0 ;
if( triad->identify() == (mp->get_bvect_cart()).identify() ) {
bases = mp->get_mg()->std_base_vect_cart() ;
}
else {
assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
bases = mp->get_mg()->std_base_vect_spher() ;
}
Itbl ind(2) ;
for (int i=0 ; i<n_comp ; i++) {
ind = indices(i) ;
cmp[i]->set_spectral_base( (*bases[ind(0)-1]) *
(*bases[ind(1)-1]) ) ;
}
for (int i=0 ; i<3 ; i++) {
delete bases[i] ;
}
delete [] bases ;
break ;
}
default : {
cout << "Tensor::std_spectral_base: the case valence = " << valence
<< " is not treated yet !" << endl ;
abort() ;
break ;
}
}
}
// Sets the standard spectal bases of decomposition for each component (odd in the nucleus)
void Tensor::std_spectral_base_odd() {
switch (valence) {
case 0 : {
cmp[0]->std_spectral_base_odd() ;
break ;
}
default : {
cout << "Tensor::std_spectral_base_odd: the case valence = " << valence
<< " is not treated yet !" << endl ;
abort() ;
break ;
}
}
}
const Tensor& Tensor::derive_cov(const Metric& metre) const {
set_dependance(metre) ;
int j = get_place_met(metre) ;
assert ((j>=0) && (j<N_MET_MAX)) ;
if (p_derive_cov[j] == 0x0) {
p_derive_cov[j] = metre.connect().p_derive_cov(*this) ;
}
return *p_derive_cov[j] ;
}
const Tensor& Tensor::derive_con(const Metric& metre) const {
set_dependance(metre) ;
int j = get_place_met(metre) ;
assert ((j>=0) && (j<N_MET_MAX)) ;
if (p_derive_con[j] == 0x0) {
if (valence == 0) {
p_derive_con[j] =
new Vector( contract(derive_cov(metre), 0, metre.con(), 0) ) ;
}
else {
const Tensor_sym* tsym = dynamic_cast<const Tensor_sym*>(this) ;
if (tsym != 0x0) { // symmetric case, preserved by derive_con
const Tensor& dercov = derive_cov(metre) ;
Itbl type_ind = dercov.get_index_type() ;
type_ind.set(valence) = CON ;
p_derive_con[j] = new Tensor_sym(*mp, valence+1, type_ind, *triad,
tsym->sym_index1(), tsym->sym_index2()) ;
*(p_derive_con[j]) = contract(dercov, valence, metre.con(), 0) ;
// valence is the number of the last index of derive_cov
// (the "derivation" index)
}
else { // general case, no symmetry
p_derive_con[j] = new Tensor( contract(derive_cov(metre), valence,
metre.con(), 0) ) ;
// valence is the number of the last index of derive_cov
// (the "derivation" index)
}
}
}
return *p_derive_con[j] ;
}
const Tensor& Tensor::divergence(const Metric& metre) const {
set_dependance(metre) ;
int j = get_place_met(metre) ;
assert ((j>=0) && (j<N_MET_MAX)) ;
if (p_divergence[j] == 0x0) {
p_divergence[j] = metre.connect().p_divergence(*this) ;
}
return *p_divergence[j] ;
}
void Tensor::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 {
cout << "Tensor::exponential_filter_r : " << endl ;
cout << "Only Cartesian triad is implemented!" << endl ;
cout << "Exiting..." << endl ;
abort() ;
}
}
void Tensor::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 {
cout << "Tensor::exponential_filter_ylm : " << endl ;
cout << "Only Cartesian triad is implemented!" << endl ;
cout << "Exiting..." << endl ;
abort() ;
}
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Lorene::Base_vect::bvect_from_file
static Base_vect * bvect_from_file(FILE *)
Construction of a vectorial basis from a file (see sauve(FILE* ) ).
Definition base_vect_from_file.C:84
Lorene::Connection::p_derive_cov
virtual Tensor * p_derive_cov(const Tensor &tens) const
Computes the covariant derivative of a tensor (with respect to the current connection).
Definition connection.C:307
Lorene::Connection::p_divergence
virtual Tensor * p_divergence(const Tensor &tens) const
Computes the divergence of a tensor (with respect to the current connection).
Definition connection.C:463
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::get_dim
int get_dim(int i) const
Gives the i th dimension (ie {tt dim.dim[i] ).
Definition itbl.h:326
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::Itbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim ).
Definition itbl.h:323
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Lorene::Metric::connect
virtual const Connection & connect() const
Returns the connection.
Definition metric.C:301
Lorene::Metric::tensor_depend
const Tensor * tensor_depend[N_TENSOR_DEPEND]
Pointer on the dependancies, that means the array contains pointers on all the Tensor whom derivative...
Definition metric.h:139
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
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::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
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::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
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::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::fread_be
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
Lorene::fwrite_be
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene::Tensor::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 tensor.C:1079
Lorene::Tensor::operator=
virtual void operator=(const Tensor &)
Assignment to a Tensor.
Definition tensor.C:553
Lorene::Tensor_sym::sym_index1
int sym_index1() const
Number of the first symmetric index (0<= id_sym1 < valence ).
Definition tensor.h:1149
Lorene::Tensor::sauve
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
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::get_index_type
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition tensor.h:886
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::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tensor.C:519
Lorene::Tensor::allocate_all
virtual void allocate_all()
Performs the memory allocation of all the elements, down to the double arrays of the Tbl s.
Definition tensor.C:508
Lorene::Tensor::std_spectral_base_odd
virtual void std_spectral_base_odd()
Sets the standard odd spectal bases of decomposition for each component.
Definition tensor.C:982
Lorene::Tensor::p_divergence
Tensor * p_divergence[N_MET_MAX]
Array of pointers on the divergence of this with respect to various metrics.
Definition tensor.h:351
Lorene::Tensor::annule_extern_cn
void annule_extern_cn(int l_0, int deg)
Performs a smooth (C^n) transition in a given domain to zero.
Definition tensor.C:690
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::set_der_0x0
void set_der_0x0() const
Sets the pointers on derived quantities to 0x0.
Definition tensor.C:407
Lorene::Tensor::spectral_display
virtual void spectral_display(const char *comment=0x0, double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions of each component.
Definition tensor.C:874
Lorene::Tensor_sym::sym_index2
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence ).
Definition tensor.h:1154
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::p_derive_cov
Tensor * p_derive_cov[N_MET_MAX]
Array of pointers on the covariant derivatives of this with respect to various metrics.
Definition tensor.h:335
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::indices
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition tensor.C:539
Lorene::Tensor::position
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor.C:525
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_etat_nondef
virtual void set_etat_nondef()
Sets the logical state of all components to ETATNONDEF (undefined state).
Definition tensor.C:489
Lorene::Tensor::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
Lorene::Tensor::~Tensor
virtual ~Tensor()
Destructor.
Definition tensor.C:385
Lorene::Tensor::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
Lorene::Tensor::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
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::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
Lorene::Tensor::n_comp
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
Lorene::Tensor::operator-=
void operator-=(const Tensor &)
-= Tensor
Definition tensor.C:587
Lorene::Tensor::set_der_met_0x0
void set_der_met_0x0(int) const
Sets all the i-th components of met_depend , p_derive_cov , etc... to 0x0.
Definition tensor.C:433
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::operator+=
void operator+=(const Tensor &)
+= Tensor
Definition tensor.C:571
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::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
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::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 tensor.C:1066
Lorene::Tensor::p_derive_con
Tensor * p_derive_con[N_MET_MAX]
Array of pointers on the contravariant derivatives of this with respect to various metrics.
Definition tensor.h:343
Lorene::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
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::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139