LORENE
tensor_sym.C
1 /*
2  * Methods of class Tensor_sym
3  *
4  * (see file tensor.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Eric Gourgoulhon & Jerome Novak
10  *
11  * Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
12  *
13  * This file is part of LORENE.
14  *
15  * LORENE is free software; you can redistribute it and/or modify
16  * it under the terms of the GNU General Public License as published by
17  * the Free Software Foundation; either version 2 of the License, or
18  * (at your option) any later version.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30 
31 
32 
33 
34 /*
35  * $Id: tensor_sym.C,v 1.4 2016/12/05 16:18:17 j_novak Exp $
36  * $Log: tensor_sym.C,v $
37  * Revision 1.4 2016/12/05 16:18:17 j_novak
38  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
39  *
40  * Revision 1.3 2014/10/13 08:53:44 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.2 2014/10/06 15:13:20 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.1 2004/01/04 20:51:45 e_gourgoulhon
47  * New class to deal with general tensors which are symmetric with
48  * respect to two of their indices.
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_sym.C,v 1.4 2016/12/05 16:18:17 j_novak Exp $
52  *
53  */
54 
55 // Headers C
56 #include <cstdlib>
57 #include <cassert>
58 #include <cmath>
59 
60 // Headers Lorene
61 #include "tensor.h"
62 #include "utilitaires.h"
63 
64 
65  //--------------//
66  // Constructors //
67  //--------------//
68 
69 // Standard constructor
70 // --------------------
71 namespace Lorene {
72 Tensor_sym::Tensor_sym(const Map& map, int val, const Itbl& tipe,
73  const Base_vect& triad_i, int index_sym1, int index_sym2)
74  : Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
75  id_sym1(index_sym1),
76  id_sym2(index_sym2) {
77 
78  // Des verifs :
79  assert( valence >= 2 ) ;
80  assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
81  assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
82  assert( id_sym1 != id_sym2 ) ;
83 
84  // The symmetry indices must be of same type:
85  assert( tipe(id_sym1) == tipe(id_sym2) ) ;
86 
87  // Possible re-ordering of the symmetry indices
88  if (id_sym1 > id_sym2) {
89  int tmp = id_sym1 ;
90  id_sym1 = id_sym2 ;
91  id_sym2 = tmp ;
92  }
93 
94 }
95 
96 
97 
98 // Standard constructor when all the indices are of the same type
99 // --------------------------------------------------------------
100 Tensor_sym::Tensor_sym(const Map& map, int val, int tipe,
101  const Base_vect& triad_i, int index_sym1, int index_sym2)
102  : Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
103  id_sym1(index_sym1),
104  id_sym2(index_sym2) {
105 
106  // Des verifs :
107  assert( valence >= 2 ) ;
108  assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
109  assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
110  assert( id_sym1 != id_sym2 ) ;
111 
112  // Possible re-ordering of the symmetry indices
113  if (id_sym1 > id_sym2) {
114  int tmp = id_sym1 ;
115  id_sym1 = id_sym2 ;
116  id_sym2 = tmp ;
117  }
118 
119 }
120 
121 // Constructor for a valence 3 symmetric tensor
122 // --------------------------------------------
123 Tensor_sym::Tensor_sym(const Map& map, int tipe0, int tipe1, int tipe2,
124  const Base_vect& triad_i,
125  int index_sym1, int index_sym2)
126  : Tensor(map, 3, tipe0, 18, triad_i),
127  id_sym1(index_sym1),
128  id_sym2(index_sym2) {
129 
130  assert( (tipe0==COV) || (tipe0==CON) ) ;
131  assert( (tipe1==COV) || (tipe1==CON) ) ;
132  assert( (tipe2==COV) || (tipe2==CON) ) ;
133 
134  type_indice.set(1) = tipe1 ;
135  type_indice.set(2) = tipe2 ;
136 
137  assert( (id_sym1 >=0) && (id_sym1 < 3) ) ;
138  assert( (id_sym2 >=0) && (id_sym2 < 3) ) ;
139  assert( id_sym1 != id_sym2 ) ;
140  assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
141 
142  // Possible re-ordering of the symmetry indices
143  if (id_sym1 > id_sym2) {
144  int tmp = id_sym1 ;
145  id_sym1 = id_sym2 ;
146  id_sym2 = tmp ;
147  }
148 
149 }
150 
151 
152 
153 // Copy constructor
154 // ----------------
156  : Tensor(*source.mp, source.valence, source.type_indice,
157  6*int(pow(3.,source.valence-2)) , *(source.triad)),
158  id_sym1(source.id_sym1),
159  id_sym2(source.id_sym2) {
160 
161  for (int i=0 ; i<n_comp ; i++) {
162 
163  int posi = source.position(indices(i)) ; // in case source belongs to
164  // a derived class of
165  // Tensor_sym with a different
166  // storage of components
167  *(cmp[i]) = *(source.cmp[posi]) ;
168  }
169 }
170 
171 
172 
173 
174 
175 // Constructor from a file
176 // -----------------------
177 Tensor_sym::Tensor_sym(const Map& map, const Base_vect& triad_i, FILE* fd)
178  : Tensor(map, triad_i, fd) {
179 
180  fread_be(&id_sym1, sizeof(int), 1, fd) ;
181  fread_be(&id_sym2, sizeof(int), 1, fd) ;
182 
183  assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
184 
185 }
186 
187 
188  //--------------//
189  // Destructor //
190  //--------------//
191 
193 
194 }
195 
196  //--------------//
197  // Assignment //
198  //--------------//
199 
200 
202 
203  assert (valence == tt.valence) ;
204 
205  triad = tt.triad ;
206  id_sym1 = tt.id_sym1 ;
207  id_sym2 = tt.id_sym2 ;
208 
209  for (int id=0 ; id<valence ; id++)
210  assert(tt.type_indice(id) == type_indice(id)) ;
211 
212  for (int ic=0 ; ic<n_comp ; ic++) {
213  int posi = tt.position(indices(ic)) ;
214  *cmp[ic] = *(tt.cmp[posi]) ;
215  }
216 
217  del_deriv() ;
218 
219 }
220 
221 void Tensor_sym::operator=(const Tensor& tt) {
222 
223  assert (valence == tt.get_valence()) ;
224 
225  triad = tt.get_triad() ;
226 
227  for (int id=0 ; id<valence ; id++)
228  assert(tt.type_indice(id) == type_indice(id)) ;
229 
230  // The symmetry indices must be of same type:
231  assert( tt.type_indice(id_sym1) == tt.type_indice(id_sym2) ) ;
232 
233 
234  for (int ic=0 ; ic<n_comp ; ic++) {
235  int posi = tt.position(indices(ic)) ;
236  *cmp[ic] = *(tt.cmp[posi]) ;
237  }
238 
239  del_deriv() ;
240 
241 }
242 
243 
244  //--------------//
245  // Accessor //
246  //--------------//
247 
248 int Tensor_sym::position(const Itbl& idx) const {
249 
250  // Protections:
251  assert (idx.get_ndim() == 1) ;
252  assert (idx.get_dim(0) == valence) ;
253  for (int i=0 ; i<valence ; i++) {
254  assert( (idx(i)>=1) && (idx(i)<=3) ) ;
255  }
256 
257  // The two symmetric indices are moved to the end --> new index array idx0
258  Itbl idx0(valence) ;
259  if (valence > 2) {
260  for (int id=0 ; id<id_sym1; id++) {
261  idx0.set(id) = idx(id) ;
262  }
263  for (int id=id_sym1; id<id_sym2-1; id++) {
264  idx0.set(id) = idx(id+1) ;
265  }
266  for (int id=id_sym2-1; id<valence-2; id++) {
267  idx0.set(id) = idx(id+2) ;
268  }
269  idx0.set(valence-2) = idx(id_sym1) ; //## not used
270  idx0.set(valence-1) = idx(id_sym2) ; //## in what follows
271  }
272 
273  // Values of the symmetric indices:
274  int is1 = idx(id_sym1) ;
275  int is2 = idx(id_sym2) ;
276 
277  // Reordering to ensure is1 <= is2 :
278  if (is2 < is1) {
279  int aux = is1 ;
280  is1 = is2 ;
281  is2 = aux ;
282  }
283 
284  // Position in the cmp array :
285  int pos = 0 ;
286  for (int id=0 ; id<valence-2 ; id++) {
287  pos = 3 * pos + idx0(id) - 1 ; // all the values of each non symmetric
288  // index occupy 3 "boxes"
289  }
290 
291  pos = 6 * pos ; // all the values of the two symmetric
292  // indices occupy 6 "boxes"
293  switch (is1) {
294  case 1 : {
295  pos += is2 - 1 ; // (1,1), (1,2) and (1,3) stored respectively
296  break ; // in relative position 0, 1 and 2
297  }
298  case 2 : {
299  pos += is2 + 1 ; // (2,2) and (2,3) stored respectively
300  break ; // in relative position 3 and 4
301  }
302  case 3 : {
303  pos += 5 ; // (3,3) stored in relative position 5
304  break ;
305  }
306  }
307 
308  return pos ;
309 }
310 
311 
312 
313 Itbl Tensor_sym::indices(int place) const {
314 
315  assert( (place>=0) && (place<n_comp) ) ;
316 
317  // Index set with the two symmetric indices at the end:
318 
319  Itbl idx0(valence) ;
320 
321  int reste = div(place, 6).rem ;
322  place = int((place-reste)/6) ;
323 
324  if (reste<3) {
325  idx0.set(valence-2) = 1 ;
326  idx0.set(valence-1) = reste + 1 ;
327  }
328 
329  if ( (reste>2) && (reste<5) ) {
330  idx0.set(valence-2) = 2 ;
331  idx0.set(valence-1) = reste - 1 ;
332  }
333 
334  if (reste == 5) {
335  idx0.set(valence-2) = 3 ;
336  idx0.set(valence-1) = 3 ;
337  }
338 
339  // The output is ready in the case of a valence 2 tensor:
340  if (valence == 2) return idx0 ;
341 
342  for (int id=valence-3 ; id>=0 ; id--) {
343  int ind = div(place, 3).rem ;
344  place = int((place-ind)/3) ;
345  idx0.set(id) = ind + 1 ;
346  }
347 
348  // Reorganization of the index set to put the two symmetric indices at
349  // their correct positions:
350 
351  Itbl idx(valence) ;
352 
353  for (int id=0 ; id<id_sym1; id++) {
354  idx.set(id) = idx0(id) ;
355  }
356  idx.set(id_sym1) = idx0(valence-2) ;
357 
358  for (int id=id_sym1+1; id<id_sym2; id++) {
359  idx.set(id) = idx0(id-1) ;
360  }
361  idx.set(id_sym2) = idx0(valence-1) ;
362 
363  for (int id=id_sym2+1; id<valence; id++) {
364  idx.set(id) = idx0(id-2) ;
365  }
366 
367  return idx ;
368 }
369 
370 
371  //--------------//
372  // Outputs //
373  //--------------//
374 
375 void Tensor_sym::sauve(FILE* fd) const {
376 
377  Tensor::sauve(fd) ;
378 
379  fwrite_be(&id_sym1, sizeof(int), 1, fd) ;
380  fwrite_be(&id_sym2, sizeof(int), 1, fd) ;
381 
382 }
383 
384 
385 
386 
387 
388 
389 
390 
391 
392 }
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:375
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
int n_comp
Number of stored components, depending on the symmetry.
Definition: tensor.h:318
virtual void operator=(const Tensor_sym &a)
Assignment to another Tensor_sym.
Definition: tensor_sym.C:201
int get_ndim() const
Gives the number of dimensions (ie dim.ndim )
Definition: itbl.h:323
Lorene prototypes.
Definition: app_hor.h:67
int id_sym2
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Definition: tensor.h:1062
Base class for coordinate mappings.
Definition: map.h:688
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition: tensor_sym.C:313
Basic integer array class.
Definition: itbl.h:122
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:915
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:72
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a cova...
Definition: tensor.h:316
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition: tensor.C:534
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
virtual ~Tensor_sym()
Destructor.
Definition: tensor_sym.C:192
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:73
int id_sym1
Number of the first symmetric index (0<= id_sym1 < valence )
Definition: tensor.h:1057
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tensor handling.
Definition: tensor.h:294
int get_valence() const
Returns the valence.
Definition: tensor.h:882
virtual void del_deriv() const
Deletes the derived quantities.
Definition: tensor.C:407
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:72
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:304
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:248
int get_dim(int i) const
Gives the i th dimension (ie {tt dim.dim[i] )
Definition: itbl.h:326
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050