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