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 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 $" ;
33 
34 /*
35  * $Id: tensor_sym.C,v 1.3 2014/10/13 08:53:44 j_novak Exp $
36  * $Log: tensor_sym.C,v $
37  * Revision 1.3 2014/10/13 08:53:44 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.2 2014/10/06 15:13:20 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.1 2004/01/04 20:51:45 e_gourgoulhon
44  * New class to deal with general tensors which are symmetric with
45  * respect to two of their indices.
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_sym.C,v 1.3 2014/10/13 08:53:44 j_novak Exp $
49  *
50  */
51 
52 // Headers C
53 #include <cstdlib>
54 #include <cassert>
55 #include <cmath>
56 
57 // Headers Lorene
58 #include "tensor.h"
59 #include "utilitaires.h"
60 
61 
62  //--------------//
63  // Constructors //
64  //--------------//
65 
66 // Standard constructor
67 // --------------------
68 namespace Lorene {
69 Tensor_sym::Tensor_sym(const Map& map, int val, const Itbl& tipe,
70  const Base_vect& triad_i, int index_sym1, int index_sym2)
71  : Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
72  id_sym1(index_sym1),
73  id_sym2(index_sym2) {
74 
75  // Des verifs :
76  assert( valence >= 2 ) ;
77  assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
78  assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
79  assert( id_sym1 != id_sym2 ) ;
80 
81  // The symmetry indices must be of same type:
82  assert( tipe(id_sym1) == tipe(id_sym2) ) ;
83 
84  // Possible re-ordering of the symmetry indices
85  if (id_sym1 > id_sym2) {
86  int tmp = id_sym1 ;
87  id_sym1 = id_sym2 ;
88  id_sym2 = tmp ;
89  }
90 
91 }
92 
93 
94 
95 // Standard constructor when all the indices are of the same type
96 // --------------------------------------------------------------
97 Tensor_sym::Tensor_sym(const Map& map, int val, int tipe,
98  const Base_vect& triad_i, int index_sym1, int index_sym2)
99  : Tensor(map, val, tipe, 6*int(pow(3.,val-2)), triad_i),
100  id_sym1(index_sym1),
101  id_sym2(index_sym2) {
102 
103  // Des verifs :
104  assert( valence >= 2 ) ;
105  assert( (id_sym1 >=0) && (id_sym1 < valence) ) ;
106  assert( (id_sym2 >=0) && (id_sym2 < valence) ) ;
107  assert( id_sym1 != id_sym2 ) ;
108 
109  // Possible re-ordering of the symmetry indices
110  if (id_sym1 > id_sym2) {
111  int tmp = id_sym1 ;
112  id_sym1 = id_sym2 ;
113  id_sym2 = tmp ;
114  }
115 
116 }
117 
118 // Constructor for a valence 3 symmetric tensor
119 // --------------------------------------------
120 Tensor_sym::Tensor_sym(const Map& map, int tipe0, int tipe1, int tipe2,
121  const Base_vect& triad_i,
122  int index_sym1, int index_sym2)
123  : Tensor(map, 3, tipe0, 18, triad_i),
124  id_sym1(index_sym1),
125  id_sym2(index_sym2) {
126 
127  assert( (tipe0==COV) || (tipe0==CON) ) ;
128  assert( (tipe1==COV) || (tipe1==CON) ) ;
129  assert( (tipe2==COV) || (tipe2==CON) ) ;
130 
131  type_indice.set(1) = tipe1 ;
132  type_indice.set(2) = tipe2 ;
133 
134  assert( (id_sym1 >=0) && (id_sym1 < 3) ) ;
135  assert( (id_sym2 >=0) && (id_sym2 < 3) ) ;
136  assert( id_sym1 != id_sym2 ) ;
137  assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
138 
139  // Possible re-ordering of the symmetry indices
140  if (id_sym1 > id_sym2) {
141  int tmp = id_sym1 ;
142  id_sym1 = id_sym2 ;
143  id_sym2 = tmp ;
144  }
145 
146 }
147 
148 
149 
150 // Copy constructor
151 // ----------------
153  : Tensor(*source.mp, source.valence, source.type_indice,
154  6*int(pow(3.,source.valence-2)) , *(source.triad)),
155  id_sym1(source.id_sym1),
156  id_sym2(source.id_sym2) {
157 
158  for (int i=0 ; i<n_comp ; i++) {
159 
160  int posi = source.position(indices(i)) ; // in case source belongs to
161  // a derived class of
162  // Tensor_sym with a different
163  // storage of components
164  *(cmp[i]) = *(source.cmp[posi]) ;
165  }
166 }
167 
168 
169 
170 
171 
172 // Constructor from a file
173 // -----------------------
174 Tensor_sym::Tensor_sym(const Map& map, const Base_vect& triad_i, FILE* fd)
175  : Tensor(map, triad_i, fd) {
176 
177  fread_be(&id_sym1, sizeof(int), 1, fd) ;
178  fread_be(&id_sym2, sizeof(int), 1, fd) ;
179 
180  assert( type_indice(id_sym1) == type_indice(id_sym2) ) ;
181 
182 }
183 
184 
185  //--------------//
186  // Destructor //
187  //--------------//
188 
190 
191 }
192 
193  //--------------//
194  // Assignment //
195  //--------------//
196 
197 
199 
200  assert (valence == tt.valence) ;
201 
202  triad = tt.triad ;
203  id_sym1 = tt.id_sym1 ;
204  id_sym2 = tt.id_sym2 ;
205 
206  for (int id=0 ; id<valence ; id++)
207  assert(tt.type_indice(id) == type_indice(id)) ;
208 
209  for (int ic=0 ; ic<n_comp ; ic++) {
210  int posi = tt.position(indices(ic)) ;
211  *cmp[ic] = *(tt.cmp[posi]) ;
212  }
213 
214  del_deriv() ;
215 
216 }
217 
218 void Tensor_sym::operator=(const Tensor& tt) {
219 
220  assert (valence == tt.get_valence()) ;
221 
222  triad = tt.get_triad() ;
223 
224  for (int id=0 ; id<valence ; id++)
225  assert(tt.type_indice(id) == type_indice(id)) ;
226 
227  // The symmetry indices must be of same type:
228  assert( tt.type_indice(id_sym1) == tt.type_indice(id_sym2) ) ;
229 
230 
231  for (int ic=0 ; ic<n_comp ; ic++) {
232  int posi = tt.position(indices(ic)) ;
233  *cmp[ic] = *(tt.cmp[posi]) ;
234  }
235 
236  del_deriv() ;
237 
238 }
239 
240 
241  //--------------//
242  // Accessor //
243  //--------------//
244 
245 int Tensor_sym::position(const Itbl& idx) const {
246 
247  // Protections:
248  assert (idx.get_ndim() == 1) ;
249  assert (idx.get_dim(0) == valence) ;
250  for (int i=0 ; i<valence ; i++) {
251  assert( (idx(i)>=1) && (idx(i)<=3) ) ;
252  }
253 
254  // The two symmetric indices are moved to the end --> new index array idx0
255  Itbl idx0(valence) ;
256  if (valence > 2) {
257  for (int id=0 ; id<id_sym1; id++) {
258  idx0.set(id) = idx(id) ;
259  }
260  for (int id=id_sym1; id<id_sym2-1; id++) {
261  idx0.set(id) = idx(id+1) ;
262  }
263  for (int id=id_sym2-1; id<valence-2; id++) {
264  idx0.set(id) = idx(id+2) ;
265  }
266  idx0.set(valence-2) = idx(id_sym1) ; //## not used
267  idx0.set(valence-1) = idx(id_sym2) ; //## in what follows
268  }
269 
270  // Values of the symmetric indices:
271  int is1 = idx(id_sym1) ;
272  int is2 = idx(id_sym2) ;
273 
274  // Reordering to ensure is1 <= is2 :
275  if (is2 < is1) {
276  int aux = is1 ;
277  is1 = is2 ;
278  is2 = aux ;
279  }
280 
281  // Position in the cmp array :
282  int pos = 0 ;
283  for (int id=0 ; id<valence-2 ; id++) {
284  pos = 3 * pos + idx0(id) - 1 ; // all the values of each non symmetric
285  // index occupy 3 "boxes"
286  }
287 
288  pos = 6 * pos ; // all the values of the two symmetric
289  // indices occupy 6 "boxes"
290  switch (is1) {
291  case 1 : {
292  pos += is2 - 1 ; // (1,1), (1,2) and (1,3) stored respectively
293  break ; // in relative position 0, 1 and 2
294  }
295  case 2 : {
296  pos += is2 + 1 ; // (2,2) and (2,3) stored respectively
297  break ; // in relative position 3 and 4
298  }
299  case 3 : {
300  pos += 5 ; // (3,3) stored in relative position 5
301  break ;
302  }
303  }
304 
305  return pos ;
306 }
307 
308 
309 
310 Itbl Tensor_sym::indices(int place) const {
311 
312  assert( (place>=0) && (place<n_comp) ) ;
313 
314  // Index set with the two symmetric indices at the end:
315 
316  Itbl idx0(valence) ;
317 
318  int reste = div(place, 6).rem ;
319  place = int((place-reste)/6) ;
320 
321  if (reste<3) {
322  idx0.set(valence-2) = 1 ;
323  idx0.set(valence-1) = reste + 1 ;
324  }
325 
326  if ( (reste>2) && (reste<5) ) {
327  idx0.set(valence-2) = 2 ;
328  idx0.set(valence-1) = reste - 1 ;
329  }
330 
331  if (reste == 5) {
332  idx0.set(valence-2) = 3 ;
333  idx0.set(valence-1) = 3 ;
334  }
335 
336  // The output is ready in the case of a valence 2 tensor:
337  if (valence == 2) return idx0 ;
338 
339  for (int id=valence-3 ; id>=0 ; id--) {
340  int ind = div(place, 3).rem ;
341  place = int((place-ind)/3) ;
342  idx0.set(id) = ind + 1 ;
343  }
344 
345  // Reorganization of the index set to put the two symmetric indices at
346  // their correct positions:
347 
348  Itbl idx(valence) ;
349 
350  for (int id=0 ; id<id_sym1; id++) {
351  idx.set(id) = idx0(id) ;
352  }
353  idx.set(id_sym1) = idx0(valence-2) ;
354 
355  for (int id=id_sym1+1; id<id_sym2; id++) {
356  idx.set(id) = idx0(id-1) ;
357  }
358  idx.set(id_sym2) = idx0(valence-1) ;
359 
360  for (int id=id_sym2+1; id<valence; id++) {
361  idx.set(id) = idx0(id-2) ;
362  }
363 
364  return idx ;
365 }
366 
367 
368  //--------------//
369  // Outputs //
370  //--------------//
371 
372 void Tensor_sym::sauve(FILE* fd) const {
373 
374  Tensor::sauve(fd) ;
375 
376  fwrite_be(&id_sym1, sizeof(int), 1, fd) ;
377  fwrite_be(&id_sym2, sizeof(int), 1, fd) ;
378 
379 }
380 
381 
382 
383 
384 
385 
386 
387 
388 
389 }
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:372
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:312
virtual void operator=(const Tensor_sym &a)
Assignment to another Tensor_sym.
Definition: tensor_sym.C:198
int get_ndim() const
Gives the number of dimensions (ie dim.ndim )
Definition: itbl.h:323
Lorene prototypes.
Definition: app_hor.h:64
int id_sym2
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Definition: tensor.h:1049
Base class for coordinate mappings.
Definition: map.h:670
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
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:303
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
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
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 covar...
Definition: tensor.h:310
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
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
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
virtual ~Tensor_sym()
Destructor.
Definition: tensor_sym.C:189
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
int id_sym1
Number of the first symmetric index (0<= id_sym1 < valence )
Definition: tensor.h:1044
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Tensor handling.
Definition: tensor.h:288
int get_valence() const
Returns the valence.
Definition: tensor.h:869
virtual void del_deriv() const
Deletes the derived quantities.
Definition: tensor.C:398
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
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:298
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
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:1037
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295