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
32char 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// --------------------
68namespace Lorene {
69Tensor_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// --------------------------------------------------------------
97Tensor_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// --------------------------------------------
120Tensor_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 ) ;
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// ----------------
152Tensor_sym::Tensor_sym(const Tensor_sym& source)
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// -----------------------
174Tensor_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
181
182}
183
184
185 //--------------//
186 // Destructor //
187 //--------------//
188
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
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
245int 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
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
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}
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
int position(int j) const
Gives the position in the arrays step, the_time and val corresponding to the time step j.
Definition evolution.C:273
Basic integer array class.
Definition itbl.h:122
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition itbl.h:247
Base class for coordinate mappings.
Definition map.h:670
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Tensor handling.
Definition tensor.h:288
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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 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
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
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
virtual void operator=(const Tensor_sym &a)
Assignment to another Tensor_sym.
Definition tensor_sym.C:198
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
int id_sym1
Number of the first symmetric index (0<= id_sym1 < valence )
Definition tensor.h:1044
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition tensor.h:298
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
int id_sym2
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Definition tensor.h:1049
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
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
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
virtual ~Tensor_sym()
Destructor.
Definition tensor_sym.C:189
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64