32 char tenseur_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $" ;
183 #include "metrique.h"
184 #include "utilitaires.h"
201 for (
int i=0; i<N_MET_MAX; i++) {
210 mp(&map), valence(0), triad(0x0),
211 type_indice(0), n_comp(1), etat(ETATNONDEF), poids(weight),
225 mp(ci.get_mp()), valence(0), triad(0x0),
226 type_indice(0), n_comp(1), etat(ci.get_etat()), poids(weight),
229 assert(ci.
get_etat() != ETATNONDEF) ;
235 assert( ci.
get_etat() == ETATQCQ ) ;
247 const Base_vect& triad_i,
const Metrique* met,
double weight)
248 : mp(&map), valence(val), triad(&triad_i), type_indice(tipe),
249 n_comp(int(
pow(3., val))), etat(ETATNONDEF), poids(weight),
256 for (
int i=0 ; i<
valence ; i++)
257 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
261 for (
int i=0 ; i<
n_comp ; i++)
269 const Base_vect* triad_i,
const Metrique* met,
double weight)
270 : mp(&map), valence(val), triad(triad_i), type_indice(tipe),
271 n_comp(int(
pow(3., val))), etat(ETATNONDEF), poids(weight),
278 for (
int i=0 ; i<
valence ; i++)
279 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
287 for (
int i=0 ; i<
n_comp ; i++)
298 const Metrique* met,
double weight)
299 : mp(&map), valence(val), triad(&triad_i), type_indice(val),
300 n_comp(int(
pow(3., val))), etat (ETATNONDEF), poids(weight),
305 assert ((tipe == COV) || (tipe == CON)) ;
311 for (
int i=0 ; i<
n_comp ; i++)
319 mp(source.mp), valence(source.valence), triad(source.triad),
320 type_indice(source.type_indice), etat (source.etat), poids(source.poids),
321 metric(source.metric) {
328 for (
int i=0 ; i<
n_comp ; i++) {
330 if (source.
c[place_source] == 0x0)
333 c[i] =
new Cmp(*source.
c[place_source]) ;
348 for (
int i=0; i<N_MET_MAX; i++) {
367 mp(source.mp), valence(source.valence), triad(source.triad),
368 type_indice(source.type_indice), etat(source.etat),
369 poids(source.poids), metric(source.metric) {
375 for (
int i=0 ; i<
n_comp ; i++) {
377 if (source.
c[place_source] == 0x0)
380 c[i] =
new Cmp(*source.
c[place_source]) ;
392 for (
int i=0; i<N_MET_MAX; i++) {
413 : mp(&mapping), triad(&triad_i), type_indice(fd),
420 assert( *triad_fich == *
triad) ;
431 for (
int i=0 ; i<
n_comp ; i++)
434 for (
int i=0 ; i<
n_comp ; i++)
449 : mp(&mapping), type_indice(fd), metric(met){
465 if (
etat == ETATQCQ) {
486 const Base_vect& triad_i,
const Metrique* met,
double weight) :
487 mp(&map), valence(val), triad(&triad_i), type_indice(tipe), n_comp(compo),
488 etat (ETATNONDEF), poids(weight), metric(met) {
495 for (
int i=0 ; i<
valence ; i++)
496 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
500 for (
int i=0 ; i<
n_comp ; i++)
510 const Base_vect& triad_i,
const Metrique* met,
double weight) :
511 mp(&map), valence(val), triad(&triad_i), type_indice(val), n_comp(compo),
512 etat (ETATNONDEF), poids(weight), metric(met) {
516 assert ((tipe == COV) || (tipe == CON)) ;
522 for (
int i=0 ; i<
n_comp ; i++)
547 for (
int i=0 ; i<
n_comp ; i++)
556 assert( (j>=0) && (j<N_MET_MAX) ) ;
559 for (
int i=0 ; i<N_DEPEND ; i++)
574 for (
int i=0; i<N_MET_MAX; i++)
592 for (
int i=0; i<N_MET_MAX; i++)
600 for (
int i=0; i<N_MET_MAX; i++)
612 for (
int i=0; i<N_MET_MAX; i++) {
614 if ((!deja) && (
met_depend[i] != 0x0)) nmet++ ;
616 if (nmet == N_MET_MAX) {
617 cout <<
"Too many metrics in Tenseur::set_dependances" << endl ;
622 while ((conte < N_DEPEND) && (met.dependances[conte] != 0x0))
625 if (conte == N_DEPEND) {
626 cout <<
"Too many dependancies in Tenseur::set_dependances " << endl ;
630 met.dependances[conte] = this ;
639 for (
int i=0 ; i<
n_comp ; i++)
660 for (
int i=0 ; i<
n_comp ; i++) {
695 for (
int i=0 ; i<
valence ; i++)
696 assert ((idx(i)>=0) && (idx(i)<3)) ;
698 for (
int i=0 ; i<
valence ; i++)
706 assert ((place >= 0) && (place <
n_comp)) ;
711 for (
int i=
valence-1 ; i>=0 ; i--) {
712 res.
set(i) = div(place, 3).rem ;
713 place = int((place-res(i))/3) ;
726 for (
int i=0 ; i<
valence ; i++)
742 for (
int i=0 ; i<
n_comp ; i++) {
744 if (t.
c[place_t] == 0x0)
747 *
c[i] = *t.
c[place_t] ;
753 cout <<
"Unknown state in Tenseur::operator= " << endl ;
785 cout <<
"Unknown state in Tenseur::operator= " << endl ;
796 if (x ==
double(0)) {
827 assert(
etat == ETATQCQ) ;
838 assert (
etat == ETATQCQ) ;
839 assert ((ind >= 0) && (ind < 3)) ;
849 assert (
etat == ETATQCQ) ;
850 assert ((ind1 >= 0) && (ind1 < 3)) ;
851 assert ((ind2 >= 0) && (ind2 < 3)) ;
868 assert (
etat == ETATQCQ) ;
869 assert ((ind1 >= 0) && (ind1 < 3)) ;
870 assert ((ind2 >= 0) && (ind2 < 3)) ;
871 assert ((ind3 >= 0) && (ind3 < 3)) ;
875 indices.
set(0) = ind1 ;
876 indices.
set(1) = ind2 ;
877 indices.
set(2) = ind3 ;
890 assert (
etat == ETATQCQ) ;
891 for (
int i=0 ; i<
valence ; i++)
892 assert ((indices(i)>=0) && (indices(i)<3)) ;
913 assert(
etat != ETATNONDEF ) ;
915 if (
etat == ETATZERO ) {
919 assert(
etat == ETATQCQ ) ;
922 for (
int i=0 ; i<
n_comp ; i++) {
929 for (
int j=0; j<N_MET_MAX; j++) {
947 if (
etat == ETATQCQ)
return *
c[0] ;
953 cout <<
"Undefined Tensor in Tenseur::operator() ..." << endl ;
964 cout <<
"Unknown state in Tenseur::operator()" << endl ;
974 assert ((indice>=0) && (indice<3)) ;
977 if (
etat == ETATQCQ)
return *
c[indice] ;
983 cout <<
"Undefined Tensor in Tenseur::operator(int) ..." << endl ;
993 cout <<
"Unknown state in Tenseur::operator(int)" << endl ;
1002 assert ((indice1>=0) && (indice1<3)) ;
1003 assert ((indice2>=0) && (indice2<3)) ;
1006 if (
etat == ETATQCQ) {
1009 idx.
set(0) = indice1 ;
1010 idx.
set(1) = indice2 ;
1017 cout <<
"Undefined Tensor in Tenseur::operator(int, int) ..." << endl ;
1027 cout <<
"Unknown state in Tenseur::operator(int, int)" << endl ;
1036 assert ((indice1>=0) && (indice1<3)) ;
1037 assert ((indice2>=0) && (indice2<3)) ;
1038 assert ((indice3>=0) && (indice3<3)) ;
1041 if (
etat == ETATQCQ) {
1044 idx.
set(0) = indice1 ;
1045 idx.
set(1) = indice2 ;
1046 idx.
set(2) = indice3 ;
1053 cout <<
"Undefined Tensor in Tenseur::operator(int, int, int) ..." << endl ;
1063 cout <<
"Unknown state in Tenseur::operator(int, int, int)" << endl ;
1075 for (
int i=0 ; i<
valence ; i++)
1076 assert ((indices(i)>=0) && (indices(i)<3)) ;
1078 if (
etat == ETATQCQ) {
1085 cout <<
"Undefined Tensor in Tenseur::operator(const Itbl&) ..." << endl ;
1095 cout <<
"Unknown state in Tenseur::operator(const Itbl& )" << endl ;
1106 if (
etat == ETATZERO) {
1110 assert(
etat == ETATQCQ) ;
1112 for (
int i=0 ; i<
n_comp ; i++)
1119 if (
etat == ETATZERO) {
1123 assert(
etat == ETATQCQ) ;
1125 for (
int i=0 ; i<
n_comp ; i++)
1132 if (
etat == ETATZERO) {
1136 assert(
etat == ETATQCQ) ;
1138 for (
int i=0 ; i<
n_comp ; i++)
1145 if (
etat == ETATZERO) {
1149 assert(
etat == ETATQCQ) ;
1151 for (
int i=0 ; i<
n_comp ; i++)
1158 if (
etat == ETATZERO) {
1162 assert(
etat == ETATQCQ) ;
1164 for (
int i=0 ; i<
n_comp ; i++)
1172 if (
etat == ETATZERO) {
1176 assert(
etat == ETATQCQ) ;
1190 for (
int i=0 ; i<3 ; i++)
1191 (
c[i]->va).set_base( *bases[i] ) ;
1192 for (
int i=0 ; i<3 ; i++)
1200 for (
int i=0 ; i<3 ; i++)
1201 (
c[i]->va).set_base( *bases[i] ) ;
1202 for (
int i=0 ; i<3 ; i++)
1218 for (
int i=0 ; i<
n_comp ; i++) {
1220 (
c[i]->
va).set_base( (*bases[indices(0)]) *
1221 (*bases[indices(1)]) ) ;
1223 for (
int i=0 ; i<3 ; i++)
1233 for (
int i=0 ; i<
n_comp ; i++) {
1235 (
c[i]->
va).set_base( (*bases[indices(0)]) *
1236 (*bases[indices(1)]) ) ;
1238 for (
int i=0 ; i<3 ; i++)
1247 "Tenseur::set_std_base() : the case valence = " <<
valence
1248 <<
" is not treated !" << endl ;
1256 ostream& operator<<(ostream& flux,
const Tenseur &source ) {
1258 flux <<
"Valence : " << source.
valence << endl ;
1260 flux <<
"Tensor density of weight " << source.
poids << endl ;
1263 flux <<
"Vectorial basis (triad) on which the components are defined :"
1269 flux <<
"Type of the indices : " << endl ;
1270 for (
int i=0 ; i<source.
valence ; i++) {
1271 flux <<
"Index " << i <<
" : " ;
1273 flux <<
" contravariant." << endl ;
1275 flux <<
" covariant." << endl ;
1278 switch (source.
etat) {
1281 flux <<
"Null Tenseur. " << endl ;
1286 flux <<
"Undefined Tenseur. " << endl ;
1291 for (
int i=0 ; i<source.
n_comp ; i++) {
1294 flux <<
"Component " ;
1297 for (
int j=0 ; j<source.
valence ; j++)
1298 flux <<
" " << num_indices(j) ;
1302 flux <<
" : " << endl ;
1303 flux <<
"-------------" << endl ;
1306 if (source.
c[i] != 0x0)
1307 flux << *source.
c[i] << endl ;
1309 flux <<
"Unknown component ... " << endl ;
1315 cout <<
"Unknown case in operator<< (ostream&, const Tenseur&)" << endl ;
1321 flux <<
" -----------------------------------------------------" << endl ;
1337 if (
etat == ETATQCQ)
1338 for (
int i=0 ; i<
n_comp ; i++)
1347 assert (
etat != ETATNONDEF) ;
1357 for (
int i=0 ; i<
valence ; i++)
1369 if (
etat == ETATZERO)
1382 for (
int m=0 ; m<
valence ; m++)
1383 indices_source.
set(m) = indices_res(m+1) ;
1386 (*this)(indices_source).deriv(indices_res(0)) ;
1394 assert (
etat != ETATNONDEF) ;
1404 "Tenseur::fait_gradient_spher : the valence must be zero !"
1412 if (
etat == ETATZERO) {
1416 assert(
etat == ETATQCQ ) ;
1431 assert (
etat != ETATNONDEF) ;
1445 for (
int i=0 ; i<
valence ; i++) {
1456 indices_gamma.
set(0) = indices(0) ;
1457 indices_gamma.
set(1) = indices(i+1) ;
1460 indices_gamma.
set(idx) = indices(idx-1) ;
1462 indices_gamma.
set(idx) = indices(idx) ;
1477 indices_gamma.
set(0) = indices(i+1) ;
1478 indices_gamma.
set(1) = indices(0) ;
1481 indices_gamma.
set(idx) = indices(idx-1) ;
1483 indices_gamma.
set(idx) = indices(idx) ;
1531 for (
int indice=1 ; indice<
valence ; indice++) {
1534 auxi =
new Tenseur(*auxi_old) ;
Bases of the spectral expansions.
Vectorial bases (triads) with respect to which the tensorial components are defined.
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
static Base_vect * bvect_from_file(FILE *)
Construction of a vectorial basis from a file (see sauve(FILE* ) ).
virtual void sauve(FILE *) const
Save in a file.
virtual void change_basis(Tenseur &) const =0
Change the basis in which the components of a tensor are expressed.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
int get_etat() const
Returns the logical state.
Valeur va
The numerical value of the Cmp
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
const Cmp & srstdsdp() const
Returns of *this .
void annule(int l)
Sets the Cmp to zero in a given domain.
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
const Cmp & srdsdt() const
Returns of *this .
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
const Cmp & dsdr() const
Returns of *this .
void inc_dzpuis()
Increases by the value of dzpuis and changes accordingly the values of the Cmp in the external compac...
Basic integer array class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int & set(int i)
Read/write of a particular element (index i ) (1D case)
int get_dim(int i) const
Gives the i th dimension (ie {tt dim.dim[i] )
void sauve(FILE *) const
Save in a file.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim )
Base class for coordinate mappings.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Base_val ** std_base_vect_spher() const
Returns the standard spectral bases for the spherical components of a vector.
int get_nzone() const
Returns the number of domains.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
virtual int donne_place(const Itbl &idx) const
Returns the position in the Cmp 1-D array c of a component given by its indices.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
virtual Itbl donne_indices(int place) const
Returns the indices of a component given by its position in the Cmp 1-D array c .
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tenseur ** p_derive_con
Array of pointers on the contravariant derivatives of *this with respect to the corresponding metri...
friend Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
void fait_carre_scal(const Metrique &, int i) const
Calculates, if needed, the scalar square of *this , with respect to met .
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
void fait_gradient_spher() const
Calculates, if needed, the gradient of *this in a spherical orthonormal basis (scalar field only).
void sauve(FILE *) const
Save in a file.
virtual void operator=(const Tenseur &tens)
Assignment to another Tenseur.
Tenseur * p_gradient_spher
Pointer on the gradient of *this in a spherical orthonormal basis (scalar field only).
const Map *const mp
Reference mapping.
void inc_dzpuis()
dzpuis += 1 ;
virtual void fait_gradient() const
Calculates, if needed, the gradient of *this .
const Metrique ** met_depend
Array of pointers on the Metrique 's used to calculate derivatives members.
virtual int donne_place(const Itbl &idx) const
Returns the position in the Cmp 1-D array c of a component given by its indices.
const Tenseur & derive_cov(const Metrique &met) const
Returns the covariant derivative of *this , with respect to met .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
void annule(int l)
Sets the Tenseur to zero in a given domain.
void dec2_dzpuis()
dzpuis -= 2 ;
Tenseur ** p_carre_scal
Array of pointers on the scalar squares of *this with respect to the corresponding metric in *met_d...
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
double get_poids() const
Returns the weight.
int get_place_met(const Metrique &metre) const
Returns the position of the pointer on metre in the array met_depend .
void dec_dzpuis()
dzpuis -= 1 ;
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
virtual ~Tenseur()
Destructor.
int n_comp
Number of components, depending on the symmetry.
void new_der_met()
Builds the arrays met_depend , p_derive_cov , p_derive_con and p_carre_scal and fills them with null ...
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
void set_poids(double weight)
Sets the weight for a tensor density.
void mult_r_zec()
Multiplication by r in the external zone.
Itbl type_indice
Array of size valence contening the type of each index, COV for a covariant one and CON for a contrav...
void del_derive() const
Logical destructor of all the derivatives.
bool verif() const
Returns false for a tensor density without a defined metric.
const Cmp & operator()() const
Read only for a scalar.
void set_metric(const Metrique &met)
Sets the pointer on the metric for a tensor density.
void inc2_dzpuis()
dzpuis += 2 ;
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
void del_derive_met(int i) const
Logical destructor of the derivatives depending on the i-th element of *met_depend .
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
double poids
For tensor densities: the weight.
void set_der_met_0x0(int i) const
Sets the pointers of the derivatives depending on the i-th element of *met_depend to zero (as well as...
Tenseur ** p_derive_cov
Array of pointers on the covariant derivatives of *this with respect to the corresponding metric in...
void set_der_0x0() const
Sets the pointers of all the derivatives to zero.
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Tenseur * p_gradient
Pointer on the gradient of *this .
void set_dependance(const Metrique &met) const
To be used to describe the fact that the derivatives members have been calculated with met .
void del_t()
Logical destructor.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
virtual void fait_derive_cov(const Metrique &met, int i) const
Calculates, if needed, the covariant derivative of *this , with respect to met .
const Tenseur & carre_scal(const Metrique &) const
Returns the scalar square of *this , with respect to met .
virtual void fait_derive_con(const Metrique &, int i) const
Calculates, if needed, the contravariant derivative of *this , with respect to met .
Cmp pow(const Cmp &, int)
Power .
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
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.