tensor_sym.C
/*
* Methods of class Tensor_sym
*
* (see file tensor.h for documentation)
*
*/
/*
* Copyright (c) 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_sym_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_sym.C,v 1.3 2014/10/13 08:53:44 j_novak Exp $" ;
/*
* $Id: tensor_sym.C,v 1.3 2014/10/13 08:53:44 j_novak Exp $
* $Log: tensor_sym.C,v $
* Revision 1.3 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.2 2014/10/06 15:13:20 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.1 2004/01/04 20:51:45 e_gourgoulhon
* New class to deal with general tensors which are symmetric with
* respect to two of their indices.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_sym.C,v 1.3 2014/10/13 08:53:44 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
// Headers Lorene
#include "tensor.h"
#include "utilitaires.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Tensor_sym::Tensor_sym(const Map& map, int val, const Itbl& tipe,
const Base_vect& triad_i, int index_sym1, int index_sym2)
: Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
id_sym1(index_sym1),
id_sym2(index_sym2) {
// Des verifs :
assert( valence >= 2 ) ;
assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
assert( id_sym1 != id_sym2 ) ;
// The symmetry indices must be of same type:
assert( tipe(id_sym1) == tipe(id_sym2) ) ;
// Possible re-ordering of the symmetry indices
if (id_sym1 > id_sym2) {
int tmp = id_sym1 ;
id_sym1 = id_sym2 ;
id_sym2 = tmp ;
}
}
// Standard constructor when all the indices are of the same type
// --------------------------------------------------------------
Tensor_sym::Tensor_sym(const Map& map, int val, int tipe,
const Base_vect& triad_i, int index_sym1, int index_sym2)
: Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
id_sym1(index_sym1),
id_sym2(index_sym2) {
// Des verifs :
assert( valence >= 2 ) ;
assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
assert( id_sym1 != id_sym2 ) ;
// Possible re-ordering of the symmetry indices
if (id_sym1 > id_sym2) {
int tmp = id_sym1 ;
id_sym1 = id_sym2 ;
id_sym2 = tmp ;
}
}
// Constructor for a valence 3 symmetric tensor
// --------------------------------------------
Tensor_sym::Tensor_sym(const Map& map, int tipe0, int tipe1, int tipe2,
const Base_vect& triad_i,
int index_sym1, int index_sym2)
: Tensor(map, 3, tipe0, 18, triad_i),
id_sym1(index_sym1),
id_sym2(index_sym2) {
assert( (tipe0==COV) || (tipe0==CON) ) ;
assert( (tipe1==COV) || (tipe1==CON) ) ;
assert( (tipe2==COV) || (tipe2==CON) ) ;
type_indice.set(1) = tipe1 ;
type_indice.set(2) = tipe2 ;
assert( (id_sym1 >=0) && (id_sym1 < 3) ) ;
assert( (id_sym2 >=0) && (id_sym2 < 3) ) ;
assert( id_sym1 != id_sym2 ) ;
assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
// Possible re-ordering of the symmetry indices
if (id_sym1 > id_sym2) {
int tmp = id_sym1 ;
id_sym1 = id_sym2 ;
id_sym2 = tmp ;
}
}
// Copy constructor
// ----------------
Tensor_sym::Tensor_sym(const Tensor_sym& source)
: Tensor(*source.mp, source.valence, source.type_indice,
6*int(pow(3.,source.valence-2)) , *(source.triad)),
id_sym1(source.id_sym1),
id_sym2(source.id_sym2) {
for (int i=0 ; i<n_comp ; i++) {
int posi = source.position(indices(i)) ; // in case source belongs to
// a derived class of
// Tensor_sym with a different
// storage of components
*(cmp[i]) = *(source.cmp[posi]) ;
}
}
// Constructor from a file
// -----------------------
Tensor_sym::Tensor_sym(const Map& map, const Base_vect& triad_i, FILE* fd)
: Tensor(map, triad_i, fd) {
fread_be(&id_sym1, sizeof(int), 1, fd) ;
fread_be(&id_sym2, sizeof(int), 1, fd) ;
assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
}
//--------------//
// Destructor //
//--------------//
Tensor_sym::~Tensor_sym() {
}
//--------------//
// Assignment //
//--------------//
void Tensor_sym::operator=(const Tensor_sym& tt) {
assert (valence == tt.valence) ;
triad = tt.triad ;
id_sym1 = tt.id_sym1 ;
id_sym2 = tt.id_sym2 ;
for (int id=0 ; id<valence ; id++)
assert(tt.type_indice(id) == type_indice(id)) ;
for (int ic=0 ; ic<n_comp ; ic++) {
int posi = tt.position(indices(ic)) ;
*cmp[ic] = *(tt.cmp[posi]) ;
}
del_deriv() ;
}
void Tensor_sym::operator=(const Tensor& tt) {
assert (valence == tt.get_valence()) ;
triad = tt.get_triad() ;
for (int id=0 ; id<valence ; id++)
assert(tt.type_indice(id) == type_indice(id)) ;
// The symmetry indices must be of same type:
assert( tt.type_indice(id_sym1) == tt.type_indice(id_sym2) ) ;
for (int ic=0 ; ic<n_comp ; ic++) {
int posi = tt.position(indices(ic)) ;
*cmp[ic] = *(tt.cmp[posi]) ;
}
del_deriv() ;
}
//--------------//
// Accessor //
//--------------//
int Tensor_sym::position(const Itbl& idx) const {
// Protections:
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) ) ;
}
// The two symmetric indices are moved to the end --> new index array idx0
Itbl idx0(valence) ;
if (valence > 2) {
for (int id=0 ; id<id_sym1; id++) {
idx0.set(id) = idx(id) ;
}
for (int id=id_sym1; id<id_sym2-1; id++) {
idx0.set(id) = idx(id+1) ;
}
for (int id=id_sym2-1; id<valence-2; id++) {
idx0.set(id) = idx(id+2) ;
}
idx0.set(valence-2) = idx(id_sym1) ; //## not used
idx0.set(valence-1) = idx(id_sym2) ; //## in what follows
}
// Values of the symmetric indices:
int is1 = idx(id_sym1) ;
int is2 = idx(id_sym2) ;
// Reordering to ensure is1 <= is2 :
if (is2 < is1) {
int aux = is1 ;
is1 = is2 ;
is2 = aux ;
}
// Position in the cmp array :
int pos = 0 ;
for (int id=0 ; id<valence-2 ; id++) {
pos = 3 * pos + idx0(id) - 1 ; // all the values of each non symmetric
// index occupy 3 "boxes"
}
pos = 6 * pos ; // all the values of the two symmetric
// indices occupy 6 "boxes"
switch (is1) {
case 1 : {
pos += is2 - 1 ; // (1,1), (1,2) and (1,3) stored respectively
break ; // in relative position 0, 1 and 2
}
case 2 : {
pos += is2 + 1 ; // (2,2) and (2,3) stored respectively
break ; // in relative position 3 and 4
}
case 3 : {
pos += 5 ; // (3,3) stored in relative position 5
break ;
}
}
return pos ;
}
Itbl Tensor_sym::indices(int place) const {
assert( (place>=0) && (place<n_comp) ) ;
// Index set with the two symmetric indices at the end:
Itbl idx0(valence) ;
int reste = div(place, 6).rem ;
place = int((place-reste)/6) ;
if (reste<3) {
idx0.set(valence-2) = 1 ;
idx0.set(valence-1) = reste + 1 ;
}
if ( (reste>2) && (reste<5) ) {
idx0.set(valence-2) = 2 ;
idx0.set(valence-1) = reste - 1 ;
}
if (reste == 5) {
idx0.set(valence-2) = 3 ;
idx0.set(valence-1) = 3 ;
}
// The output is ready in the case of a valence 2 tensor:
if (valence == 2) return idx0 ;
for (int id=valence-3 ; id>=0 ; id--) {
int ind = div(place, 3).rem ;
place = int((place-ind)/3) ;
idx0.set(id) = ind + 1 ;
}
// Reorganization of the index set to put the two symmetric indices at
// their correct positions:
Itbl idx(valence) ;
for (int id=0 ; id<id_sym1; id++) {
idx.set(id) = idx0(id) ;
}
idx.set(id_sym1) = idx0(valence-2) ;
for (int id=id_sym1+1; id<id_sym2; id++) {
idx.set(id) = idx0(id-1) ;
}
idx.set(id_sym2) = idx0(valence-1) ;
for (int id=id_sym2+1; id<valence; id++) {
idx.set(id) = idx0(id-2) ;
}
return idx ;
}
//--------------//
// Outputs //
//--------------//
void Tensor_sym::sauve(FILE* fd) const {
Tensor::sauve(fd) ;
fwrite_be(&id_sym1, sizeof(int), 1, fd) ;
fwrite_be(&id_sym2, sizeof(int), 1, fd) ;
}
}
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
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::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim ).
Definition itbl.h:323
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::sauve
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
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_sym::indices
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition tensor_sym.C:310
Lorene::Tensor_sym::operator=
virtual void operator=(const Tensor_sym &a)
Assignment to another Tensor_sym.
Definition tensor_sym.C:198
Lorene::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition tensor.h:869
Lorene::Tensor_sym::sauve
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
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_sym::id_sym1
int id_sym1
Number of the first symmetric index (0<= id_sym1 < valence ).
Definition tensor.h:1044
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_sym::Tensor_sym
Tensor_sym(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i, int index_sym1, int index_sym2)
Standard constructor.
Definition tensor_sym.C:69
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
Lorene::Tensor_sym::id_sym2
int id_sym2
Number of the second symmetric index (id_sym1 < id_sym2 < valence ).
Definition tensor.h:1049
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::del_deriv
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
Lorene::Tensor_sym::position
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor_sym.C:245
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_sym::~Tensor_sym
virtual ~Tensor_sym()
Destructor.
Definition tensor_sym.C:189
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