LORENE
tenseur_sym.C
1 /*
2  * Methods of class Tenseur_sym
3  *
4  * (see file tenseur.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 1999-2001 Philippe Grandclement
10  * Copyright (c) 2000-2001 Eric Gourgoulhon
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: tenseur_sym.C,v 1.11 2021/04/14 16:20:53 j_novak Exp $
35  * $Log: tenseur_sym.C,v $
36  * Revision 1.11 2021/04/14 16:20:53 j_novak
37  * Corrected uncorrect recent change
38  *
39  * Revision 1.10 2021/04/13 11:24:53 j_novak
40  * Corrected a bug in Etoile_rot::fait_shift() which was missing the np=1 case.
41  *
42  * Revision 1.9 2016/12/05 16:18:17 j_novak
43  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
44  *
45  * Revision 1.8 2014/10/13 08:53:42 j_novak
46  * Lorene classes and functions now belong to the namespace Lorene.
47  *
48  * Revision 1.7 2014/10/06 15:13:19 j_novak
49  * Modified #include directives to use c++ syntax.
50  *
51  * Revision 1.6 2003/03/03 19:39:58 f_limousin
52  * Modification of an assert to have a check on a triad and not only on a pointer.
53  *
54  * Revision 1.5 2002/10/16 14:37:14 j_novak
55  * Reorganization of #include instructions of standard C++, in order to
56  * use experimental version 3 of gcc.
57  *
58  * Revision 1.4 2002/09/06 14:49:25 j_novak
59  * Added method lie_derive for Tenseur and Tenseur_sym.
60  * Corrected various errors for derive_cov and arithmetic.
61  *
62  * Revision 1.3 2002/08/14 13:46:15 j_novak
63  * Derived quantities of a Tenseur can now depend on several Metrique's
64  *
65  * Revision 1.2 2002/08/07 16:14:11 j_novak
66  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
67  *
68  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
69  * LORENE
70  *
71  * Revision 2.3 2001/10/10 13:55:23 eric
72  * Modif Joachim: pow(3, *) --> pow(3., *)
73  *
74  * Revision 2.2 2000/02/09 19:33:38 eric
75  * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
76  * argument des constructeurs.
77  *
78  * Revision 2.1 2000/01/11 11:14:41 eric
79  * Gestion de la base vectorielle (triad).
80  *
81  * Revision 2.0 1999/12/02 17:18:37 phil
82  * *** empty log message ***
83  *
84  *
85  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym.C,v 1.11 2021/04/14 16:20:53 j_novak Exp $
86  *
87  */
88 
89 // Headers C
90 #include <cstdlib>
91 #include <cassert>
92 #include <cmath>
93 
94 // Headers Lorene
95 #include "tenseur.h"
96 #include "metrique.h"
97 
98  //--------------//
99  // Constructors //
100  //--------------//
101 
102 // Standard constructor
103 // --------------------
104 namespace Lorene {
105 Tenseur_sym::Tenseur_sym(const Map& map, int val, const Itbl& tipe,
106  const Base_vect& triad_i, const Metrique* met,
107  double weight)
108  : Tenseur(map, val, tipe, int(pow(3., val-2)) * 6, triad_i,
109  met, weight) {
110 
111  assert (val >= 2) ;
112 }
113 
114 // Standard constructor when all the indices are of the same type
115 // --------------------------------------------------------------
116 Tenseur_sym::Tenseur_sym(const Map& map, int val, int tipe,
117  const Base_vect& triad_i, const Metrique* met,
118  double weight)
119  : Tenseur(map, val, tipe, int(pow(3., val-2)) * 6, triad_i,
120  met, weight) {
121 
122  assert (val >= 2) ;
123 }
124 
125 // Copy constructor
126 // ----------------
128  Tenseur (*source.mp, source.valence, source.type_indice,
129  int(pow(3., source.valence-2)*6), *(source.triad), source.metric,
130  source.poids) {
131 
132  assert (valence >= 2) ;
133  for (int i=0 ; i<n_comp ; i++) {
134  int place_source = source.donne_place(donne_indices(i)) ;
135  if (source.c[place_source] == 0x0)
136  c[i] = 0x0 ;
137  else
138  c[i] = new Cmp (*source.c[place_source]) ;
139  }
140  etat = source.etat ;
141  assert(source.met_depend != 0x0) ;
142  assert(source.p_derive_cov != 0x0) ;
143  assert(source.p_derive_con != 0x0) ;
144  assert(source.p_carre_scal != 0x0) ;
145 
146  if (source.p_gradient != 0x0)
147  p_gradient = new Tenseur_sym (*source.p_gradient) ;
148 
149  for (int i=0; i<N_MET_MAX; i++) {
150  met_depend[i] = source.met_depend[i] ;
151  if (met_depend[i] != 0x0) {
152 
153  set_dependance (*met_depend[i]) ;
154 
155  if (source.p_derive_cov[i] != 0x0)
156  p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
157  if (source.p_derive_con[i] != 0x0)
158  p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
159  if (source.p_carre_scal[i] != 0x0)
160  p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
161  }
162  }
163 }
164 
165 
166 // Constructor from a Tenseur
167 // --------------------------
169  Tenseur (*source.mp, source.valence, source.type_indice,
170  int(pow(3., source.valence-2)*6), *(source.triad), source.metric,
171  source.poids) {
172 
173  assert (valence >= 2) ;
174 
175  for (int i=0 ; i<n_comp ; i++) {
176  int place_source = source.donne_place(donne_indices(i)) ;
177  if (source.c[place_source] == 0x0)
178  c[i] = 0x0 ;
179  else
180  c[i] = new Cmp (*source.c[place_source]) ;
181  }
182 
183  etat = source.etat ;
184 
185  assert(source.met_depend != 0x0) ;
186  assert(source.p_derive_cov != 0x0) ;
187  assert(source.p_derive_con != 0x0) ;
188  assert(source.p_carre_scal != 0x0) ;
189 
190  if (source.p_gradient != 0x0)
191  p_gradient = new Tenseur (*source.p_gradient) ;
192 
193 }
194 
195 
196 // Constructor from a file
197 // -----------------------
198 Tenseur_sym::Tenseur_sym(const Map& map, const Base_vect& triad_i, FILE* fd,
199  const Metrique* met)
200  : Tenseur(map, triad_i, fd, met) {
201 
202  assert (valence >= 2) ;
203  assert (n_comp == int(pow(3., valence-2))*6) ;
204 }
205 
206  //--------------//
207  // Destructor //
208  //--------------//
209 
211 
212 
213 
214 
215 
216 int Tenseur_sym::donne_place (const Itbl& idx) const {
217 
218  assert (idx.get_ndim() == 1) ;
219  assert (idx.get_dim(0) == valence) ;
220  for (int i=0 ; i<valence ; i++)
221  assert ((idx(i) >= 0) && (idx(i) < 3)) ;
222 
223 
224  // Gestion des deux derniers indices :
225  int last = idx(valence-1) ;
226  int lastm1 = idx(valence-2) ;
227  if (last < lastm1) {
228  int auxi = last ;
229  last = lastm1 ;
230  lastm1 = auxi ;
231  }
232 
233  int place_fin ;
234  switch (lastm1) {
235  case 0 : {
236  place_fin = last ;
237  break ;
238  }
239  case 1 : {
240  place_fin = 2+last ;
241  break ;
242  }
243  case 2 : {
244  place_fin = 5 ;
245  break ;
246  }
247  default : {
248  abort() ;
249  }
250  }
251 
252  int res = 0 ;
253  for (int i=0 ; i<valence-2 ; i++)
254  res = 3*res+idx(i) ;
255 
256  res = 6*res + place_fin ;
257 
258  return res ;
259 }
260 
261 Itbl Tenseur_sym::donne_indices (int place) const {
262  Itbl res(valence) ;
263  res.set_etat_qcq() ;
264  assert ((place>=0) && (place<n_comp)) ;
265 
266  int reste = div(place, 6).rem ;
267  place = int((place-reste)/6) ;
268 
269  for (int i=valence-3 ; i>=0 ; i--) {
270  res.set(i) = div(place, 3).rem ;
271  place = int((place-res(i))/3) ;
272  }
273 
274  if (reste<3) {
275  res.set(valence-2) = 0 ;
276  res.set(valence-1) = reste ;
277  }
278 
279  if ((reste>2) && (reste<5)) {
280  res.set(valence-2) = 1 ;
281  res.set(valence-1) = reste - 2 ;
282  }
283 
284  if (reste == 5) {
285  res.set(valence-2) = 2 ;
286  res.set(valence-1) = 2 ;
287  }
288 
289  return res ;
290 }
291 
293 
294  assert (valence == t.get_valence()) ;
295 
296  triad = t.triad ;
297  poids = t.poids ;
298  metric = t.metric ;
299 
300  for (int i=0 ; i<valence ; i++)
301  assert (type_indice(i) == t.type_indice(i)) ;
302 
303  switch (t.get_etat()) {
304  case ETATNONDEF: {
305  set_etat_nondef() ;
306  break ;
307  }
308 
309  case ETATZERO: {
310  set_etat_zero() ;
311  break ;
312  }
313 
314  case ETATQCQ: {
315  set_etat_qcq() ;
316  for (int i=0 ; i<n_comp ; i++) {
317  int place_t = t.donne_place(donne_indices(i)) ;
318  if (t.c[place_t] == 0x0)
319  c[i] = 0x0 ;
320  else
321  *c[i] = *t.c[place_t] ;
322  }
323  break ;
324  }
325 
326  default: {
327  cout << "Unknown state in Tenseur_sym::operator= " << endl ;
328  abort() ;
329  break ;
330  }
331  }
332 }
333 
334 void Tenseur_sym::operator= (const Tenseur_sym& t) {
335  Tenseur::operator=(t) ;
336 }
337 
339 
340  assert (etat != ETATNONDEF) ;
341 
342  if (p_gradient != 0x0)
343  return ;
344  else {
345 
346  // Construction du resultat :
347  Itbl tipe (valence+1) ;
348  tipe.set_etat_qcq() ;
349  tipe.set(0) = COV ;
350  for (int i=0 ; i<valence ; i++)
351  tipe.set(i+1) = type_indice(i) ;
352 
353  // Vectorial basis
354  // ---------------
355  assert(*triad == mp->get_bvect_cart()) ;
356 
357  p_gradient = new Tenseur_sym(*mp, valence+1, tipe,
358  mp->get_bvect_cart(), metric, poids) ;
359 
360 
361  if (etat == ETATZERO)
363  else {
365  // Boucle sur les indices :
366 
367  Itbl indices_source(valence) ;
368  indices_source.set_etat_qcq() ;
369 
370  for (int j=0 ; j<p_gradient->n_comp ; j++) {
371 
372  Itbl indices_res(p_gradient->donne_indices(j)) ;
373 
374  for (int m=0 ; m<valence ; m++)
375  indices_source.set(m) = indices_res(m+1) ;
376 
377  p_gradient->set(indices_res) =
378  (*this)(indices_source).deriv(indices_res(0)) ;
379  }
380  }
381  }
382 }
383 
384 
385 void Tenseur_sym::fait_derive_cov (const Metrique& metre, int ind) const {
386 
387  assert (etat != ETATNONDEF) ;
388  assert (valence != 0) ;
389 
390  if (p_derive_cov[ind] != 0x0)
391  return ;
392  else {
393  p_derive_cov[ind] = new Tenseur_sym (gradient()) ;
394 
395  if ((valence != 0) && (etat != ETATZERO)) {
396 
397  assert( *(metre.gamma().get_triad()) == *triad ) ;
398 
399  Tenseur* auxi ;
400  for (int i=0 ; i<valence ; i++) {
401 
402  if (type_indice(i) == COV) {
403  auxi = new Tenseur(contract(metre.gamma(), 0,(*this), i)) ;
404 
405  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
406  indices_gamma.set_etat_qcq() ;
407  //On range comme il faut :
408  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
409 
410  Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
411  indices_gamma.set(0) = indices(0) ;
412  indices_gamma.set(1) = indices(i+1) ;
413  for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
414  if (idx<=i+1)
415  indices_gamma.set(idx) = indices(idx-1) ;
416  else
417  indices_gamma.set(idx) = indices(idx) ;
418 
419  p_derive_cov[ind]->set(indices) -= (*auxi)(indices_gamma) ;
420  }
421  }
422  else {
423  auxi = new Tenseur(contract(metre.gamma(), 1, (*this), i)) ;
424 
425  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
426  indices_gamma.set_etat_qcq() ;
427 
428  //On range comme il faut :
429  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
430 
431  Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
432  indices_gamma.set(0) = indices(i+1) ;
433  indices_gamma.set(1) = indices(0) ;
434  for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
435  if (idx<=i+1)
436  indices_gamma.set(idx) = indices(idx-1) ;
437  else
438  indices_gamma.set(idx) = indices(idx) ;
439  p_derive_cov[ind]->set(indices) += (*auxi)(indices_gamma) ;
440  }
441  }
442  delete auxi ;
443  }
444  }
445  if ((poids != 0.)&&(etat != ETATZERO))
446  *p_derive_cov[ind] = *p_derive_cov[ind] -
447  poids*contract(metre.gamma(), 0, 2) * (*this) ;
448 
449  }
450 }
451 
452 
453 
454 void Tenseur_sym::fait_derive_con (const Metrique& metre, int ind) const {
455 
456  if (p_derive_con[ind] != 0x0)
457  return ;
458  else {
459  // On calcul la derivee covariante :
460  if (valence != 0)
461  p_derive_con[ind] = new Tenseur_sym
462  (contract(metre.con(), 1, derive_cov(metre), 0)) ;
463 
464  else
465  p_derive_con[ind] = new Tenseur_sym
466  (contract(metre.con(), 1, gradient(), 0)) ;
467  }
468 }
469 }
Itbl type_indice
Array of size valence contening the type of each index, COV for a covariant one and CON for a contrav...
Definition: tenseur.h:321
const Map *const mp
Reference mapping.
Definition: tenseur.h:309
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tenseur.h:315
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Cmp ** c
The components.
Definition: tenseur.h:325
int n_comp
Number of components, depending on the symmetry.
Definition: tenseur.h:323
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition: tenseur.h:1256
int get_ndim() const
Gives the number of dimensions (ie dim.ndim )
Definition: itbl.h:323
Lorene prototypes.
Definition: app_hor.h:67
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.
Definition: tenseur.C:706
Base class for coordinate mappings.
Definition: map.h:688
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:324
virtual Itbl donne_indices(int place) const
Returns the indices of a component given by its position in the Cmp 1-D array c . ...
Definition: tenseur.C:720
Basic integer array class.
Definition: itbl.h:122
virtual ~Tenseur_sym()
Destructor.
Definition: tenseur_sym.C:210
virtual Itbl donne_indices(int place) const
Returns the indices of a component given by its position in the Cmp 1-D array c . ...
Definition: tenseur_sym.C:261
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
virtual void fait_derive_cov(const Metrique &met, int i) const
Calculates, if needed, the covariant derivative of *this , with respect to met .
Definition: tenseur_sym.C:385
int get_valence() const
Returns the valence.
Definition: tenseur.h:713
virtual void operator=(const Tenseur &tens)
Assignment to another Tenseur.
Definition: tenseur.C:734
virtual void operator=(const Tenseur &)
Assignment from a Tenseur .
Definition: tenseur_sym.C:292
virtual void fait_gradient() const
Calculates, if needed, the gradient of *this .
Definition: tenseur_sym.C:338
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Tenseur * p_gradient
Pointer on the gradient of *this .
Definition: tenseur.h:346
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.
Definition: tenseur_sym.C:216
int valence
Valence.
Definition: tenseur.h:310
double poids
For tensor densities: the weight.
Definition: tenseur.h:326
Tenseur ** p_derive_con
Array of pointers on the contravariant derivatives of *this with respect to the corresponding metric...
Definition: tenseur.h:368
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: itbl.C:264
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Metrique ** met_depend
Array of pointers on the Metrique &#39;s used to calculate derivatives members.
Definition: tenseur.h:340
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition: tenseur.C:225
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
void set_dependance(const Metrique &met) const
To be used to describe the fact that the derivatives members have been calculated with met ...
Definition: tenseur.C:624
Tenseur_sym(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i, const Metrique *met=0x0, double weight=0)
Standard constructor.
Definition: tenseur_sym.C:105
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:328
int get_dim(int i) const
Gives the i th dimension (ie {tt dim.dim[i] )
Definition: itbl.h:326
Tenseur ** p_derive_cov
Array of pointers on the covariant derivatives of *this with respect to the corresponding metric in ...
Definition: tenseur.h:361
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
const Tenseur & derive_cov(const Metrique &met) const
Returns the covariant derivative of *this , with respect to met .
Definition: tenseur.C:1570
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
virtual void fait_derive_con(const Metrique &, int i) const
Calculates, if needed, the contravariant derivative of *this , with respect to met ...
Definition: tenseur_sym.C:454
Tenseur ** p_carre_scal
Array of pointers on the scalar squares of *this with respect to the corresponding metric in *met_de...
Definition: tenseur.h:375
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Definition: tenseur.C:666
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558