LORENE
tenseur.C
1 /*
2  * Methods of class Tenseur
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  * Copyright (c) 2002 Jerome Novak
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 <<<<<<< tenseur.C
33 char tenseur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.14 2012/04/02 09:57:21 j_novak Exp $" ;
34 =======
35 char tenseur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $" ;
36 >>>>>>> 1.15
37 
38 /*
39 <<<<<<< tenseur.C
40  * $Id: tenseur.C,v 1.14 2012/04/02 09:57:21 j_novak Exp $
41 =======
42  * $Id: tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $
43 >>>>>>> 1.15
44  * $Log: tenseur.C,v $
45 <<<<<<< tenseur.C
46  * Revision 1.14 2012/04/02 09:57:21 j_novak
47  * Just a test
48 =======
49  * Revision 1.15 2014/10/13 08:53:41 j_novak
50  * Lorene classes and functions now belong to the namespace Lorene.
51  *
52  * Revision 1.14 2014/10/06 15:13:18 j_novak
53  * Modified #include directives to use c++ syntax.
54 >>>>>>> 1.15
55  *
56  * Revision 1.13 2003/03/03 19:37:31 f_limousin
57  * Suppression of many assert(verif()).
58  *
59  * Revision 1.12 2002/10/16 14:37:14 j_novak
60  * Reorganization of #include instructions of standard C++, in order to
61  * use experimental version 3 of gcc.
62  *
63  * Revision 1.11 2002/09/06 14:49:25 j_novak
64  * Added method lie_derive for Tenseur and Tenseur_sym.
65  * Corrected various errors for derive_cov and arithmetic.
66  *
67  * Revision 1.10 2002/08/30 13:21:38 j_novak
68  * Corrected error in constructor
69  *
70  * Revision 1.9 2002/08/14 13:46:15 j_novak
71  * Derived quantities of a Tenseur can now depend on several Metrique's
72  *
73  * Revision 1.8 2002/08/13 08:02:45 j_novak
74  * Handling of spherical vector/tensor components added in the classes
75  * Mg3d and Tenseur. Minor corrections for the class Metconf.
76  *
77  * Revision 1.7 2002/08/08 15:10:45 j_novak
78  * The flag "plat" has been added to the class Metrique to show flat metrics.
79  *
80  * Revision 1.6 2002/08/07 16:14:11 j_novak
81  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
82  *
83  * Revision 1.5 2002/08/02 15:07:41 j_novak
84  * Member function determinant has been added to the class Metrique.
85  * A better handling of spectral bases is now implemented for the class Tenseur.
86  *
87  * Revision 1.4 2002/05/07 07:36:03 e_gourgoulhon
88  * Compatibilty with xlC compiler on IBM SP2:
89  * suppressed the parentheses around argument of instruction new:
90  * e.g. t = new (Tbl *[nzone]) --> t = new Tbl*[nzone]
91  *
92  * Revision 1.3 2002/05/02 15:16:22 j_novak
93  * Added functions for more general bi-fluid EOS
94  *
95  * Revision 1.2 2001/12/04 21:27:54 e_gourgoulhon
96  *
97  * All writing/reading to a binary file are now performed according to
98  * the big endian convention, whatever the system is big endian or
99  * small endian, thanks to the functions fwrite_be and fread_be
100  *
101  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
102  * LORENE
103  *
104  * Revision 2.21 2001/10/10 13:54:40 eric
105  * Modif Joachim: pow(3, *) --> pow(3., *)
106  *
107  * Revision 2.20 2000/12/20 09:50:08 eric
108  * Correction erreur dans operator<< : la sortie doit etre flux et non cout !
109  *
110  * Revision 2.19 2000/10/12 13:11:23 eric
111  * Methode set_std_base(): traitement du cas etat = ETATZERO (return).
112  *
113  * Revision 2.18 2000/09/13 12:11:40 eric
114  * Ajout de la fonction allocate_all().
115  *
116  * Revision 2.17 2000/05/22 14:40:09 phil
117  * ajout de inc_dzpuis et dec_dzpuis
118  *
119  * Revision 2.16 2000/03/22 09:18:57 eric
120  * Traitement du cas ETATZERO dans dec2_dzpuis, inc2_dzpuis et mult_r_zec.
121  *
122  * Revision 2.15 2000/02/12 11:35:58 eric
123  * Modif de la fonction set_std_base : appel de Valeur::set_base plutot
124  * que l'affectation directe du membre Valeur::base.
125  *
126  * Revision 2.14 2000/02/10 18:30:47 eric
127  * La fonction set_triad ne fait plus que l'affectation du membre triad.
128  *
129  * Revision 2.13 2000/02/10 16:11:07 eric
130  * Ajout de la fonction change_triad.
131  *
132  * Revision 2.12 2000/02/09 19:32:39 eric
133  * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
134  * argument des constructeurs.
135  *
136  * Revision 2.11 2000/01/24 13:02:39 eric
137  * Traitement du cas triad = 0x0 dans la sauvegarde/lecture fichier
138  * Constructeur par lecture de fichier: met_depend est desormais initialise
139  * a 0x0.
140  *
141  * Revision 2.10 2000/01/20 16:02:57 eric
142  * Ajout des operator=(double ) et operator=(int ).
143  *
144  * Revision 2.9 2000/01/20 15:34:39 phil
145  * changement traid dans fait_gradient()
146  *
147  * Revision 2.8 2000/01/14 14:07:26 eric
148  * Ajout de la fonction annule.
149  *
150  * Revision 2.7 2000/01/13 14:10:53 eric
151  * Ajout du constructeur par copie d'un Cmp (pour un scalaire)
152  * ainsi que l'affectation a un Cmp.
153  *
154  * Revision 2.6 2000/01/13 13:46:38 eric
155  * Ajout du membre p_gradient_spher et des fonctions fait_gradient_spher(),
156  * gradient_spher() pour le calcul du gradient d'un scalaire en
157  * coordonnees spheriques sur la triade spherique associee.
158  *
159  * Revision 2.5 2000/01/12 13:19:04 eric
160  * Les operator::(...) renvoient desormais une reference const sur le c[...]
161  * correspondant et non plus un Cmp copie de c[...].
162  * (ceci grace a la nouvelle fonction Map::cmp_zero()).
163  *
164  * Revision 2.4 2000/01/11 11:13:59 eric
165  * Changement de nom pour la base vectorielle : base --> triad
166  *
167  * Revision 2.3 2000/01/10 17:23:07 eric
168  * Modif affichage.
169  * Methode fait_derive_cov : ajout de
170  * assert( metre.gamma().get_base() == base )
171  * Methode set_std_base : ajout de
172  * assert( *base == mp->get_bvect_cart() )
173  *
174  * Revision 2.2 2000/01/10 15:15:26 eric
175  * Ajout du membre base (base vectorielle sur laquelle sont definies
176  * les composantes).
177  *
178  * Revision 2.1 1999/12/09 12:39:23 phil
179  * changement prototypage des derivees
180  *
181  * Revision 2.0 1999/12/02 17:18:31 phil
182  * *** empty log message ***
183  *
184  *
185 <<<<<<< tenseur.C
186  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.14 2012/04/02 09:57:21 j_novak Exp $
187 =======
188  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $
189 >>>>>>> 1.15
190  *
191  */
192 
193 // Headers C
194 #include <assert.h>
195 #include <math.h>
196 
197 // Headers Lorene
198 #include "tenseur.h"
199 #include "metrique.h"
200 #include "utilitaires.h"
201 
202  //--------------//
203  // Constructors //
204  //--------------//
205 // Consistency check for tensor densities
206 //---------------------------------------
207 namespace Lorene {
208 bool Tenseur::verif() const {
209  return ( (poids == 0.) || (metric != 0x0) ) ;
210 }
211 
213  met_depend = new const Metrique*[N_MET_MAX] ;
214  p_derive_cov = new Tenseur*[N_MET_MAX] ;
215  p_derive_con = new Tenseur*[N_MET_MAX] ;
216  p_carre_scal = new Tenseur*[N_MET_MAX] ;
217  for (int i=0; i<N_MET_MAX; i++) {
218  met_depend[i] = 0x0 ;
219  }
220  set_der_0x0() ;
221 }
222 
223 // Constructor for a scalar field
224 // ------------------------------
225 Tenseur::Tenseur (const Map& map, const Metrique* met, double weight) :
226  mp(&map), valence(0), triad(0x0),
227  type_indice(0), n_comp(1), etat(ETATNONDEF), poids(weight),
228  metric(met) {
229 
230  // assert(verif()) ;
231  c = new Cmp*[n_comp] ;
232  c[0] = 0x0 ;
233  new_der_met() ;
234 }
235 
236 
237 
238 // Constructor for a scalar field and from a {\tt Cmp}
239 // ---------------------------------------------------
240 Tenseur::Tenseur (const Cmp& ci, const Metrique* met, double weight) :
241  mp(ci.get_mp()), valence(0), triad(0x0),
242  type_indice(0), n_comp(1), etat(ci.get_etat()), poids(weight),
243  metric(met){
244 
245  assert(ci.get_etat() != ETATNONDEF) ;
246  assert(verif()) ;
247 
248  c = new Cmp*[n_comp] ;
249 
250  if ( ci.get_etat() != ETATZERO ) {
251  assert( ci.get_etat() == ETATQCQ ) ;
252  c[0] = new Cmp(ci) ;
253  }
254  else {
255  c[0] = 0x0 ;
256  }
257  new_der_met() ;
258 }
259 
260 // Standard constructor
261 // --------------------
262 Tenseur::Tenseur(const Map& map, int val, const Itbl& tipe,
263  const Base_vect& triad_i, const Metrique* met, double weight)
264  : mp(&map), valence(val), triad(&triad_i), type_indice(tipe),
265  n_comp(int(pow(3., val))), etat(ETATNONDEF), poids(weight),
266  metric(met){
267 
268  // Des verifs :
269  assert (valence >= 0) ;
270  assert (tipe.get_ndim() == 1) ;
271  assert (valence == tipe.get_dim(0)) ;
272  for (int i=0 ; i<valence ; i++)
273  assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
274  assert(verif()) ;
275 
276  c = new Cmp*[n_comp] ;
277  for (int i=0 ; i<n_comp ; i++)
278  c[i] = 0x0 ;
279  new_der_met() ;
280 }
281 
282 // Standard constructor with the triad passed as a pointer
283 // -------------------------------------------------------
284 Tenseur::Tenseur(const Map& map, int val, const Itbl& tipe,
285  const Base_vect* triad_i, const Metrique* met, double weight)
286  : mp(&map), valence(val), triad(triad_i), type_indice(tipe),
287  n_comp(int(pow(3., val))), etat(ETATNONDEF), poids(weight),
288  metric(met){
289 
290  // Des verifs :
291  assert (valence >= 0) ;
292  assert (tipe.get_ndim() == 1) ;
293  assert (valence == tipe.get_dim(0)) ;
294  for (int i=0 ; i<valence ; i++)
295  assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
296  // assert(verif()) ;
297 
298  if (valence == 0) { // particular case of a scalar
299  triad = 0x0 ;
300  }
301 
302  c = new Cmp*[n_comp] ;
303  for (int i=0 ; i<n_comp ; i++)
304  c[i] = 0x0 ;
305  new_der_met() ;
306 }
307 
308 
309 
310 
311 // Standard constructor when all the indices are of the same type
312 // --------------------------------------------------------------
313 Tenseur::Tenseur(const Map& map, int val, int tipe, const Base_vect& triad_i,
314  const Metrique* met, double weight)
315  : mp(&map), valence(val), triad(&triad_i), type_indice(val),
316  n_comp(int(pow(3., val))), etat (ETATNONDEF), poids(weight),
317  metric(met){
318 
319  // Des verifs :
320  assert (valence >= 0) ;
321  assert ((tipe == COV) || (tipe == CON)) ;
322  assert(verif()) ;
324  type_indice = tipe ;
325 
326  c = new Cmp*[n_comp] ;
327  for (int i=0 ; i<n_comp ; i++)
328  c[i] = 0x0 ;
329  new_der_met() ;
330 }
331 
332 // Copy constructor
333 // ----------------
334 Tenseur::Tenseur (const Tenseur& source) :
335  mp(source.mp), valence(source.valence), triad(source.triad),
336  type_indice(source.type_indice), etat (source.etat), poids(source.poids),
337  metric(source.metric) {
338 
339  // assert(verif()) ;
340 
341  n_comp = int(pow(3., valence)) ;
342 
343  c = new Cmp*[n_comp] ;
344  for (int i=0 ; i<n_comp ; i++) {
345  int place_source = source.donne_place(donne_indices(i)) ;
346  if (source.c[place_source] == 0x0)
347  c[i] = 0x0 ;
348  else
349  c[i] = new Cmp(*source.c[place_source]) ;
350  }
351 
352  assert(source.met_depend != 0x0) ;
353  assert(source.p_derive_cov != 0x0) ;
354  assert(source.p_derive_con != 0x0) ;
355  assert(source.p_carre_scal != 0x0) ;
356  new_der_met() ;
357 
358  if (source.p_gradient != 0x0)
359  p_gradient = new Tenseur (*source.p_gradient) ;
360 
361  if (source.p_gradient_spher != 0x0)
362  p_gradient_spher = new Tenseur (*source.p_gradient_spher) ;
363 
364  for (int i=0; i<N_MET_MAX; i++) {
365  met_depend[i] = source.met_depend[i] ;
366  if (met_depend[i] != 0x0) {
367 
368  set_dependance (*met_depend[i]) ;
369 
370  if (source.p_derive_cov[i] != 0x0)
371  p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
372  if (source.p_derive_con[i] != 0x0)
373  p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
374  if (source.p_carre_scal[i] != 0x0)
375  p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
376  }
377  }
378 }
379 
380 // Constructor from a symmetric tensor
381 // -----------------------------------
382 Tenseur::Tenseur (const Tenseur_sym& source) :
383  mp(source.mp), valence(source.valence), triad(source.triad),
384  type_indice(source.type_indice), etat(source.etat),
385  poids(source.poids), metric(source.metric) {
386 
387  assert(verif()) ;
388  n_comp = int(pow(3., valence)) ;
389 
390  c = new Cmp*[n_comp] ;
391  for (int i=0 ; i<n_comp ; i++) {
392  int place_source = source.donne_place(donne_indices(i)) ;
393  if (source.c[place_source] == 0x0)
394  c[i] = 0x0 ;
395  else
396  c[i] = new Cmp(*source.c[place_source]) ;
397  }
398 
399  assert(source.met_depend != 0x0) ;
400  assert(source.p_derive_cov != 0x0) ;
401  assert(source.p_derive_con != 0x0) ;
402  assert(source.p_carre_scal != 0x0) ;
403  new_der_met() ;
404 
405  if (source.p_gradient != 0x0)
406  p_gradient = new Tenseur_sym (*source.p_gradient) ;
407 
408  for (int i=0; i<N_MET_MAX; i++) {
409  met_depend[i] = source.met_depend[i] ;
410  if (met_depend[i] != 0x0) {
411 
412  set_dependance (*met_depend[i]) ;
413 
414  if (source.p_derive_cov[i] != 0x0)
415  p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
416  if (source.p_derive_con[i] != 0x0)
417  p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
418  if (source.p_carre_scal[i] != 0x0)
419  p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
420  }
421  }
422 
423 }
424 
425 // Constructor from a file
426 // -----------------------
427 Tenseur::Tenseur(const Map& mapping, const Base_vect& triad_i, FILE* fd,
428  const Metrique* met)
429  : mp(&mapping), triad(&triad_i), type_indice(fd),
430  metric(met) {
431 
432  fread_be(&valence, sizeof(int), 1, fd) ;
433 
434  if (valence != 0) {
435  Base_vect* triad_fich = Base_vect::bvect_from_file(fd) ;
436  assert( *triad_fich == *triad) ;
437  delete triad_fich ;
438  }
439  else{
440  triad = 0x0 ;
441  }
442 
443  fread_be(&n_comp, sizeof(int), 1, fd) ;
444  fread_be(&etat, sizeof(int), 1, fd) ;
445 
446  c = new Cmp*[n_comp] ;
447  for (int i=0 ; i<n_comp ; i++)
448  c[i] = 0x0 ;
449  if (etat == ETATQCQ)
450  for (int i=0 ; i<n_comp ; i++)
451  c[i] = new Cmp(*mp, *mp->get_mg(), fd) ;
452 
453  new_der_met() ;
454 
455  if (met == 0x0)
456  poids = 0. ;
457  else
458  fread_be(&poids, sizeof(double), 1, fd) ;
459 }
460 
461 
462 // Constructor from a file for a scalar field
463 // ------------------------------------------
464 Tenseur::Tenseur (const Map& mapping, FILE* fd, const Metrique* met)
465  : mp(&mapping), type_indice(fd), metric(met){
466 
467  fread_be(&valence, sizeof(int), 1, fd) ;
468 
469  assert(valence == 0) ;
470 
471  triad = 0x0 ;
472 
473  fread_be(&n_comp, sizeof(int), 1, fd) ;
474 
475  assert(n_comp == 1) ;
476 
477  fread_be(&etat, sizeof(int), 1, fd) ;
478 
479  c = new Cmp*[n_comp] ;
480 
481  if (etat == ETATQCQ) {
482  c[0] = new Cmp(*mp, *mp->get_mg(), fd) ;
483  }
484  else{
485  c[0] = 0x0 ;
486  }
487 
488  new_der_met() ;
489 
490  if (met == 0x0)
491  poids = 0. ;
492  else
493  fread_be(&poids, sizeof(double), 1, fd) ;
494 }
495 
496 
497 
498 
499 // Constructor used by the derived classes
500 // ---------------------------------------
501 Tenseur::Tenseur (const Map& map, int val, const Itbl& tipe, int compo,
502  const Base_vect& triad_i, const Metrique* met, double weight) :
503  mp(&map), valence(val), triad(&triad_i), type_indice(tipe), n_comp(compo),
504  etat (ETATNONDEF), poids(weight), metric(met) {
505 
506  // Des verifs :
507  assert (valence >= 0) ;
508  assert (tipe.get_ndim() == 1) ;
509  assert (n_comp > 0) ;
510  assert (valence == tipe.get_dim(0)) ;
511  for (int i=0 ; i<valence ; i++)
512  assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
513  //assert(verif()) ;
514 
515  c = new Cmp*[n_comp] ;
516  for (int i=0 ; i<n_comp ; i++)
517  c[i] = 0x0 ;
518 
519  new_der_met() ;
520 }
521 
522 // Constructor used by the derived classes when all the indices are of
523 // the same type.
524 // -------------------------------------------------------------------
525 Tenseur::Tenseur (const Map& map, int val, int tipe, int compo,
526  const Base_vect& triad_i, const Metrique* met, double weight) :
527  mp(&map), valence(val), triad(&triad_i), type_indice(val), n_comp(compo),
528  etat (ETATNONDEF), poids(weight), metric(met) {
529  // Des verifs :
530  assert (valence >= 0) ;
531  assert (n_comp >= 0) ;
532  assert ((tipe == COV) || (tipe == CON)) ;
533  //assert(verif()) ;
535  type_indice = tipe ;
536 
537  c = new Cmp*[n_comp] ;
538  for (int i=0 ; i<n_comp ; i++)
539  c[i] = 0x0 ;
540 
541  new_der_met() ;
542 }
543 
544  //--------------//
545  // Destructor //
546  //--------------//
547 
548 
550 
551  del_t() ;
552  delete [] met_depend ;
553  delete [] p_derive_cov ;
554  delete [] p_derive_con ;
555  delete [] p_carre_scal ;
556  delete [] c ;
557 }
558 
559 
560 
562  del_derive() ;
563  for (int i=0 ; i<n_comp ; i++)
564  if (c[i] != 0x0) {
565  delete c[i] ;
566  c[i] = 0x0 ;
567  }
568 }
569 
570 void Tenseur::del_derive_met(int j) const {
571 
572  assert( (j>=0) && (j<N_MET_MAX) ) ;
573  // On gere la metrique ...
574  if (met_depend[j] != 0x0) {
575  for (int i=0 ; i<N_DEPEND ; i++)
576  if (met_depend[j]->dependances[i] == this)
577  met_depend[j]->dependances[i] = 0x0 ;
578  if (p_derive_cov[j] != 0x0)
579  delete p_derive_cov[j] ;
580  if (p_derive_con[j] != 0x0)
581  delete p_derive_con[j] ;
582  if (p_carre_scal[j] != 0x0)
583  delete p_carre_scal[j] ;
584  set_der_met_0x0(j) ;
585  }
586 }
587 
588 
589 void Tenseur::del_derive () const {
590  for (int i=0; i<N_MET_MAX; i++)
591  del_derive_met(i) ;
592  if (p_gradient != 0x0)
593  delete p_gradient ;
594  if (p_gradient_spher != 0x0)
595  delete p_gradient_spher ;
596  set_der_0x0() ;
597 }
598 
599 void Tenseur::set_der_met_0x0(int i) const {
600  met_depend[i] = 0x0 ;
601  p_derive_cov[i] = 0x0 ;
602  p_derive_con[i] = 0x0 ;
603  p_carre_scal[i] = 0x0 ;
604 }
605 
606 
607 void Tenseur::set_der_0x0() const {
608  for (int i=0; i<N_MET_MAX; i++)
609  set_der_met_0x0(i) ;
610  p_gradient = 0x0 ;
611  p_gradient_spher = 0x0 ;
612 }
613 
614 int Tenseur::get_place_met(const Metrique& metre) const {
615  int resu = -1 ;
616  for (int i=0; i<N_MET_MAX; i++)
617  if (met_depend[i] == &metre) {
618  assert(resu == -1) ;
619  resu = i ;
620  }
621  return resu ;
622 }
623 
624 void Tenseur::set_dependance (const Metrique& met) const {
625 
626  int nmet = 0 ;
627  bool deja = false ;
628  for (int i=0; i<N_MET_MAX; i++) {
629  if (met_depend[i] == &met) deja = true ;
630  if ((!deja) && (met_depend[i] != 0x0)) nmet++ ;
631  }
632  if (nmet == N_MET_MAX) {
633  cout << "Too many metrics in Tenseur::set_dependances" << endl ;
634  abort() ;
635  }
636  if (!deja) {
637  int conte = 0 ;
638  while ((conte < N_DEPEND) && (met.dependances[conte] != 0x0))
639  conte ++ ;
640 
641  if (conte == N_DEPEND) {
642  cout << "Too many dependancies in Tenseur::set_dependances " << endl ;
643  abort() ;
644  }
645  else {
646  met.dependances[conte] = this ;
647  met_depend[nmet] = &met ;
648  }
649  }
650 }
651 
653 
654  del_derive() ;
655  for (int i=0 ; i<n_comp ; i++)
656  if (c[i] == 0x0)
657  c[i] = new Cmp(mp) ;
658  etat = ETATQCQ ;
659 }
660 
662  del_t() ;
663  etat = ETATZERO ;
664 }
665 
667  del_t() ;
668  etat = ETATNONDEF ;
669 }
670 
671 // Allocates everything
672 // --------------------
674 
675  set_etat_qcq() ;
676  for (int i=0 ; i<n_comp ; i++) {
677  c[i]->allocate_all() ;
678  }
679 
680 }
681 
682 
683 
685 
686  bi.change_basis(*this) ;
687 
688 }
689 
690 void Tenseur::set_triad(const Base_vect& bi) {
691 
692  triad = &bi ;
693 
694 }
695 
696 void Tenseur::set_poids(double weight) {
697 
698  poids = weight ;
699 }
700 
701 void Tenseur::set_metric(const Metrique& met) {
702 
703  metric = &met ;
704 }
705 
706 int Tenseur::donne_place (const Itbl& idx) const {
707 
708  assert (idx.get_ndim() == 1) ;
709  assert (idx.get_dim(0) == valence) ;
710 
711  for (int i=0 ; i<valence ; i++)
712  assert ((idx(i)>=0) && (idx(i)<3)) ;
713  int res = 0 ;
714  for (int i=0 ; i<valence ; i++)
715  res = 3*res+idx(i) ;
716 
717  return res;
718 }
719 
720 Itbl Tenseur::donne_indices (int place) const {
721 
722  assert ((place >= 0) && (place < n_comp)) ;
723 
724  Itbl res(valence) ;
725  res.set_etat_qcq() ;
726 
727  for (int i=valence-1 ; i>=0 ; i--) {
728  res.set(i) = div(place, 3).rem ;
729  place = int((place-res(i))/3) ;
730  }
731  return res ;
732 }
733 
734 void Tenseur::operator=(const Tenseur& t) {
735 
736  assert (valence == t.valence) ;
737 
738  triad = t.triad ;
739  poids = t.poids ;
740  metric = t.metric ;
741 
742  for (int i=0 ; i<valence ; i++)
743  assert (t.type_indice(i) == type_indice(i)) ;
744 
745  switch (t.etat) {
746  case ETATNONDEF: {
747  set_etat_nondef() ;
748  break ;
749  }
750 
751  case ETATZERO: {
752  set_etat_zero() ;
753  break ;
754  }
755 
756  case ETATQCQ: {
757  set_etat_qcq() ;
758  for (int i=0 ; i<n_comp ; i++) {
759  int place_t = t.donne_place(donne_indices(i)) ;
760  if (t.c[place_t] == 0x0)
761  c[i] = 0x0 ;
762  else
763  *c[i] = *t.c[place_t] ;
764  }
765  break ;
766  }
767 
768  default: {
769  cout << "Unknown state in Tenseur::operator= " << endl ;
770  abort() ;
771  break ;
772  }
773  }
774 }
775 
776 
777 void Tenseur::operator=(const Cmp& ci) {
778 
779  assert (valence == 0) ;
780  poids = 0. ;
781  metric = 0x0 ;
782 
783  switch (ci.get_etat()) {
784  case ETATNONDEF: {
785  set_etat_nondef() ;
786  break ;
787  }
788 
789  case ETATZERO: {
790  set_etat_zero() ;
791  break ;
792  }
793 
794  case ETATQCQ: {
795  set_etat_qcq() ;
796  *(c[0]) = ci ;
797  break ;
798  }
799 
800  default: {
801  cout << "Unknown state in Tenseur::operator= " << endl ;
802  abort() ;
803  break ;
804  }
805  }
806 }
807 
808 void Tenseur::operator=(double x) {
809 
810  poids = 0. ;
811  metric = 0x0 ;
812  if (x == double(0)) {
813  set_etat_zero() ;
814  }
815  else {
816  assert(valence == 0) ;
817  set_etat_qcq() ;
818  *(c[0]) = x ;
819  }
820 
821 }
822 
823 void Tenseur::operator=(int x) {
824 
825  poids = 0. ;
826  metric = 0x0 ;
827  if (x == 0) {
828  set_etat_zero() ;
829  }
830  else {
831  assert(valence == 0) ;
832  set_etat_qcq() ;
833  *(c[0]) = x ;
834  }
835 
836 }
837 
838 
839 // Affectation d'un scalaire ...
841 
842  del_derive() ;
843  assert(etat == ETATQCQ) ;
844  assert (valence == 0) ;
845  return *c[0] ;
846 }
847 
848 
849 // Affectation d'un vecteur :
850 Cmp& Tenseur::set (int ind) {
851 
852  del_derive() ;
853  assert(valence == 1) ;
854  assert (etat == ETATQCQ) ;
855  assert ((ind >= 0) && (ind < 3)) ;
856 
857  return *c[ind] ;
858 }
859 
860 // Affectation d'un tenseur d'ordre 2 :
861 Cmp& Tenseur::set (int ind1, int ind2) {
862 
863  del_derive() ;
864  assert (valence == 2) ;
865  assert (etat == ETATQCQ) ;
866  assert ((ind1 >= 0) && (ind1 < 3)) ;
867  assert ((ind2 >= 0) && (ind2 < 3)) ;
868 
869  Itbl ind (valence) ;
870  ind.set_etat_qcq() ;
871  ind.set(0) = ind1 ;
872  ind.set(1) = ind2 ;
873 
874  int place = donne_place(ind) ;
875 
876  return *c[place] ;
877 }
878 
879 // Affectation d'un tenseur d'ordre 3 :
880 Cmp& Tenseur::set (int ind1, int ind2, int ind3) {
881 
882  del_derive() ;
883  assert (valence == 3) ;
884  assert (etat == ETATQCQ) ;
885  assert ((ind1 >= 0) && (ind1 < 3)) ;
886  assert ((ind2 >= 0) && (ind2 < 3)) ;
887  assert ((ind3 >= 0) && (ind3 < 3)) ;
888 
889  Itbl indices(valence) ;
890  indices.set_etat_qcq() ;
891  indices.set(0) = ind1 ;
892  indices.set(1) = ind2 ;
893  indices.set(2) = ind3 ;
894  int place = donne_place(indices) ;
895 
896  return *c[place] ;
897 }
898 
899 // Affectation cas general
900 Cmp& Tenseur::set(const Itbl& indices) {
901 
902  assert (indices.get_ndim() == 1) ;
903  assert (indices.get_dim(0) == valence) ;
904 
905  del_derive() ;
906  assert (etat == ETATQCQ) ;
907  for (int i=0 ; i<valence ; i++)
908  assert ((indices(i)>=0) && (indices(i)<3)) ;
909 
910  int place = donne_place(indices) ;
911 
912  return *c[place] ;
913 }
914 
915 // Annulation dans des domaines
916 void Tenseur::annule(int l) {
917 
918  annule(l, l) ;
919 }
920 
921 void Tenseur::annule(int l_min, int l_max) {
922 
923  // Cas particulier: annulation globale :
924  if ( (l_min == 0) && (l_max == mp->get_mg()->get_nzone()-1) ) {
925  set_etat_zero() ;
926  return ;
927  }
928 
929  assert( etat != ETATNONDEF ) ;
930 
931  if ( etat == ETATZERO ) {
932  return ; // rien n'a faire si c'est deja zero
933  }
934  else {
935  assert( etat == ETATQCQ ) ; // sinon...
936 
937  // Annulation des composantes:
938  for (int i=0 ; i<n_comp ; i++) {
939  c[i]->annule(l_min, l_max) ;
940  }
941 
942  // Annulation des membres derives
943  if (p_gradient != 0x0) p_gradient->annule(l_min, l_max) ;
944  if (p_gradient_spher != 0x0) p_gradient_spher->annule(l_min, l_max) ;
945  for (int j=0; j<N_MET_MAX; j++) {
946  if (p_derive_cov[j] != 0x0) p_derive_cov[j]->annule(l_min, l_max) ;
947  if (p_derive_con[j] != 0x0) p_derive_con[j]->annule(l_min, l_max) ;
948  if (p_carre_scal[j] != 0x0) p_carre_scal[j]->annule(l_min, l_max) ;
949  }
950 
951  }
952 
953 }
954 
955 
956 
957 
958 // Exctraction :
959 const Cmp& Tenseur::operator()() const {
960 
961  assert(valence == 0) ;
962 
963  if (etat == ETATQCQ) return *c[0] ; // pour la performance,
964  // ce cas est traite en premier,
965  // en dehors du switch
966  switch (etat) {
967 
968  case ETATNONDEF : {
969  cout << "Undefined Tensor in Tenseur::operator() ..." << endl ;
970  abort() ;
971  return *c[0] ; // bidon pour satisfaire le compilateur
972  }
973 
974  case ETATZERO : {
975  return mp->cmp_zero() ;
976  }
977 
978 
979  default : {
980  cout <<"Unknown state in Tenseur::operator()" << endl ;
981  abort() ;
982  return *c[0] ; // bidon pour satisfaire le compilateur
983  }
984  }
985 }
986 
987 
988 const Cmp& Tenseur::operator() (int indice) const {
989 
990  assert ((indice>=0) && (indice<3)) ;
991  assert(valence == 1) ;
992 
993  if (etat == ETATQCQ) return *c[indice] ; // pour la performance,
994  // ce cas est traite en premier,
995  // en dehors du switch
996  switch (etat) {
997 
998  case ETATNONDEF : {
999  cout << "Undefined Tensor in Tenseur::operator(int) ..." << endl ;
1000  abort() ;
1001  return *c[0] ; // bidon pour satisfaire le compilateur
1002  }
1003 
1004  case ETATZERO : {
1005  return mp->cmp_zero() ;
1006  }
1007 
1008  default : {
1009  cout <<"Unknown state in Tenseur::operator(int)" << endl ;
1010  abort() ;
1011  return *c[0] ; // bidon pour satisfaire le compilateur
1012  }
1013  }
1014 }
1015 
1016 const Cmp& Tenseur::operator() (int indice1, int indice2) const {
1017 
1018  assert ((indice1>=0) && (indice1<3)) ;
1019  assert ((indice2>=0) && (indice2<3)) ;
1020  assert(valence == 2) ;
1021 
1022  if (etat == ETATQCQ) { // pour la performance,
1023  Itbl idx(2) ; // ce cas est traite en premier,
1024  idx.set_etat_qcq() ; // en dehors du switch
1025  idx.set(0) = indice1 ;
1026  idx.set(1) = indice2 ;
1027  return *c[donne_place(idx)] ;
1028  }
1029 
1030  switch (etat) {
1031 
1032  case ETATNONDEF : {
1033  cout << "Undefined Tensor in Tenseur::operator(int, int) ..." << endl ;
1034  abort() ;
1035  return *c[0] ; // bidon pour satisfaire le compilateur
1036  }
1037 
1038  case ETATZERO : {
1039  return mp->cmp_zero() ;
1040  }
1041 
1042  default : {
1043  cout <<"Unknown state in Tenseur::operator(int, int)" << endl ;
1044  abort() ;
1045  return *c[0] ; // bidon pour satisfaire le compilateur
1046  }
1047  }
1048 }
1049 
1050 const Cmp& Tenseur::operator() (int indice1, int indice2, int indice3) const {
1051 
1052  assert ((indice1>=0) && (indice1<3)) ;
1053  assert ((indice2>=0) && (indice2<3)) ;
1054  assert ((indice3>=0) && (indice3<3)) ;
1055  assert(valence == 3) ;
1056 
1057  if (etat == ETATQCQ) { // pour la performance,
1058  Itbl idx(3) ; // ce cas est traite en premier,
1059  idx.set_etat_qcq() ; // en dehors du switch
1060  idx.set(0) = indice1 ;
1061  idx.set(1) = indice2 ;
1062  idx.set(2) = indice3 ;
1063  return *c[donne_place(idx)] ;
1064  }
1065 
1066  switch (etat) {
1067 
1068  case ETATNONDEF : {
1069  cout << "Undefined Tensor in Tenseur::operator(int, int, int) ..." << endl ;
1070  abort() ;
1071  return *c[0] ; // bidon pour satisfaire le compilateur
1072  }
1073 
1074  case ETATZERO : {
1075  return mp->cmp_zero() ;
1076  }
1077 
1078  default : {
1079  cout <<"Unknown state in Tenseur::operator(int, int, int)" << endl ;
1080  abort() ;
1081  return *c[0] ; // bidon pour satisfaire le compilateur
1082  }
1083  }
1084 }
1085 
1086 
1087 const Cmp& Tenseur::operator() (const Itbl& indices) const {
1088 
1089  assert (indices.get_ndim() == 1) ;
1090  assert (indices.get_dim(0) == valence) ;
1091  for (int i=0 ; i<valence ; i++)
1092  assert ((indices(i)>=0) && (indices(i)<3)) ;
1093 
1094  if (etat == ETATQCQ) { // pour la performance,
1095  return *c[donne_place(indices)] ; // ce cas est traite en premier,
1096  } // en dehors du switch
1097 
1098  switch (etat) {
1099 
1100  case ETATNONDEF : {
1101  cout << "Undefined Tensor in Tenseur::operator(const Itbl&) ..." << endl ;
1102  abort() ;
1103  return *c[0] ; // bidon pour satisfaire le compilateur
1104  }
1105 
1106  case ETATZERO : {
1107  return mp->cmp_zero() ;
1108  }
1109 
1110  default : {
1111  cout <<"Unknown state in Tenseur::operator(const Itbl& )" << endl ;
1112  abort() ;
1113  return *c[0] ; // bidon pour satisfaire le compilateur
1114  }
1115  }
1116 
1117 }
1118 
1119 // Gestion de la ZEC :
1121 
1122  if (etat == ETATZERO) {
1123  return ;
1124  }
1125 
1126  assert(etat == ETATQCQ) ;
1127 
1128  for (int i=0 ; i<n_comp ; i++)
1129  if (c[i] != 0x0)
1130  c[i]->dec_dzpuis() ;
1131 }
1132 
1134 
1135  if (etat == ETATZERO) {
1136  return ;
1137  }
1138 
1139  assert(etat == ETATQCQ) ;
1140 
1141  for (int i=0 ; i<n_comp ; i++)
1142  if (c[i] != 0x0)
1143  c[i]->inc_dzpuis() ;
1144 }
1145 
1147 
1148  if (etat == ETATZERO) {
1149  return ;
1150  }
1151 
1152  assert(etat == ETATQCQ) ;
1153 
1154  for (int i=0 ; i<n_comp ; i++)
1155  if (c[i] != 0x0)
1156  c[i]->dec2_dzpuis() ;
1157 }
1158 
1160 
1161  if (etat == ETATZERO) {
1162  return ;
1163  }
1164 
1165  assert(etat == ETATQCQ) ;
1166 
1167  for (int i=0 ; i<n_comp ; i++)
1168  if (c[i] != 0x0)
1169  c[i]->inc2_dzpuis() ;
1170 }
1171 
1173 
1174  if (etat == ETATZERO) {
1175  return ;
1176  }
1177 
1178  assert(etat == ETATQCQ) ;
1179 
1180  for (int i=0 ; i<n_comp ; i++)
1181  if (c[i] != 0x0)
1182  c[i]->mult_r_zec() ;
1183 }
1184 
1185 // Gestion des bases spectrales (valence <= 2)
1187 
1188  if (etat == ETATZERO) {
1189  return ;
1190  }
1191 
1192  assert(etat == ETATQCQ) ;
1193  switch (valence) {
1194 
1195  case 0 : {
1196  c[0]->std_base_scal() ;
1197  break ;
1198  }
1199 
1200  case 1 : {
1201 
1202  if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
1203 
1204  Base_val** bases = mp->get_mg()->std_base_vect_cart() ;
1205 
1206  for (int i=0 ; i<3 ; i++)
1207  (c[i]->va).set_base( *bases[i] ) ;
1208  for (int i=0 ; i<3 ; i++)
1209  delete bases[i] ;
1210  delete [] bases ;
1211  }
1212  else {
1213  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
1214  Base_val** bases = mp->get_mg()->std_base_vect_spher() ;
1215 
1216  for (int i=0 ; i<3 ; i++)
1217  (c[i]->va).set_base( *bases[i] ) ;
1218  for (int i=0 ; i<3 ; i++)
1219  delete bases[i] ;
1220  delete [] bases ;
1221  }
1222  break ;
1223 
1224  }
1225 
1226  case 2 : {
1227 
1228  if( triad->identify() == (mp->get_bvect_cart()).identify() ) {
1229 
1230  Base_val** bases = mp->get_mg()->std_base_vect_cart() ;
1231 
1232  Itbl indices (2) ;
1233  indices.set_etat_qcq() ;
1234  for (int i=0 ; i<n_comp ; i++) {
1235  indices = donne_indices(i) ;
1236  (c[i]->va).set_base( (*bases[indices(0)]) *
1237  (*bases[indices(1)]) ) ;
1238  }
1239  for (int i=0 ; i<3 ; i++)
1240  delete bases[i] ;
1241  delete [] bases ;
1242  }
1243  else {
1244  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
1245  Base_val** bases = mp->get_mg()->std_base_vect_spher() ;
1246 
1247  Itbl indices (2) ;
1248  indices.set_etat_qcq() ;
1249  for (int i=0 ; i<n_comp ; i++) {
1250  indices = donne_indices(i) ;
1251  (c[i]->va).set_base( (*bases[indices(0)]) *
1252  (*bases[indices(1)]) ) ;
1253  }
1254  for (int i=0 ; i<3 ; i++)
1255  delete bases[i] ;
1256  delete [] bases ;
1257  }
1258  break ;
1259  }
1260 
1261  default : {
1262  cout <<
1263  "Tenseur::set_std_base() : the case valence = " << valence
1264  << " is not treated !" << endl ;
1265  abort() ;
1266  break ;
1267  }
1268  }
1269 }
1270 
1271  // Le cout :
1272 ostream& operator<<(ostream& flux, const Tenseur &source ) {
1273 
1274  flux << "Valence : " << source.valence << endl ;
1275  if (source.get_poids() != 0.)
1276  flux << "Tensor density of weight " << source.poids << endl ;
1277 
1278  if (source.get_triad() != 0x0) {
1279  flux << "Vectorial basis (triad) on which the components are defined :"
1280  << endl ;
1281  flux << *(source.get_triad()) << endl ;
1282  }
1283 
1284  if (source.valence != 0)
1285  flux << "Type of the indices : " << endl ;
1286  for (int i=0 ; i<source.valence ; i++) {
1287  flux << "Index " << i << " : " ;
1288  if (source.type_indice(i) == CON)
1289  flux << " contravariant." << endl ;
1290  else
1291  flux << " covariant." << endl ;
1292  }
1293 
1294  switch (source.etat) {
1295 
1296  case ETATZERO : {
1297  flux << "Null Tenseur. " << endl ;
1298  break ;
1299  }
1300 
1301  case ETATNONDEF : {
1302  flux << "Undefined Tenseur. " << endl ;
1303  break ;
1304  }
1305 
1306  case ETATQCQ : {
1307  for (int i=0 ; i<source.n_comp ; i++) {
1308 
1309  Itbl num_indices (source.donne_indices(i)) ;
1310  flux << "Component " ;
1311 
1312  if (source.valence != 0) {
1313  for (int j=0 ; j<source.valence ; j++)
1314  flux << " " << num_indices(j) ;
1315  }
1316  else
1317  flux << " " << 0 ;
1318  flux << " : " << endl ;
1319  flux << "-------------" << endl ;
1320 
1321 
1322  if (source.c[i] != 0x0)
1323  flux << *source.c[i] << endl ;
1324  else
1325  flux << "Unknown component ... " << endl ;
1326 
1327  }
1328  break ;
1329  }
1330  default : {
1331  cout << "Unknown case in operator<< (ostream&, const Tenseur&)" << endl ;
1332  abort() ;
1333  break ;
1334  }
1335  }
1336 
1337  flux << " -----------------------------------------------------" << endl ;
1338  return flux ;
1339 }
1340 
1341 void Tenseur::sauve(FILE* fd) const {
1342 
1343  type_indice.sauve(fd) ; // type des composantes
1344  fwrite_be(&valence, sizeof(int), 1, fd) ; // la valence
1345 
1346  if (valence != 0) {
1347  triad->sauve(fd) ; // Vectorial basis
1348  }
1349 
1350  fwrite_be(&n_comp, sizeof(int), 1, fd) ; // nbre composantes
1351  fwrite_be(&etat, sizeof(int), 1, fd) ; // etat
1352 
1353  if (etat == ETATQCQ)
1354  for (int i=0 ; i<n_comp ; i++)
1355  c[i]->sauve(fd) ;
1356  if (poids != 0.)
1357  fwrite_be(&poids, sizeof(double), 1, fd) ; //poids, si pertinent
1358 }
1359 
1360 
1361 void Tenseur::fait_gradient () const {
1362 
1363  assert (etat != ETATNONDEF) ;
1364 
1365  if (p_gradient != 0x0)
1366  return ;
1367  else {
1368 
1369  // Construction du resultat :
1370  Itbl tipe (valence+1) ;
1371  tipe.set_etat_qcq() ;
1372  tipe.set(0) = COV ;
1373  for (int i=0 ; i<valence ; i++)
1374  tipe.set(i+1) = type_indice(i) ;
1375 
1376  // Vectorial basis
1377  // ---------------
1378  if (valence != 0) {
1379  // assert (*triad == mp->get_bvect_cart()) ;
1380  }
1381 
1382  p_gradient = new Tenseur(*mp, valence+1, tipe, mp->get_bvect_cart(),
1383  metric, poids) ;
1384 
1385  if (etat == ETATZERO)
1387  else {
1389  // Boucle sur les indices :
1390 
1391  Itbl indices_source(valence) ;
1392  indices_source.set_etat_qcq() ;
1393 
1394  for (int j=0 ; j<p_gradient->n_comp ; j++) {
1395 
1396  Itbl indices_res(p_gradient->donne_indices(j)) ;
1397 
1398  for (int m=0 ; m<valence ; m++)
1399  indices_source.set(m) = indices_res(m+1) ;
1400 
1401  p_gradient->set(indices_res) =
1402  (*this)(indices_source).deriv(indices_res(0)) ;
1403  }
1404  }
1405  }
1406 }
1407 
1409 
1410  assert (etat != ETATNONDEF) ;
1411 
1412  if (p_gradient_spher != 0x0)
1413  return ;
1414  else {
1415 
1416  // Construction du resultat :
1417 
1418  if (valence != 0) {
1419  cout <<
1420  "Tenseur::fait_gradient_spher : the valence must be zero !"
1421  << endl ;
1422  abort() ;
1423  }
1424 
1425  p_gradient_spher = new Tenseur(*mp, 1, COV, mp->get_bvect_spher(),
1426  metric, poids) ;
1427 
1428  if (etat == ETATZERO) {
1430  }
1431  else {
1432  assert( etat == ETATQCQ ) ;
1434 
1435  p_gradient_spher->set(0) = c[0]->dsdr() ; // d/dr
1436  p_gradient_spher->set(1) = c[0]->srdsdt() ; // 1/r d/dtheta
1437  p_gradient_spher->set(2) = c[0]->srstdsdp() ; // 1/(r sin(theta))
1438  // d/dphi
1439 
1440  }
1441  }
1442 }
1443 
1444 
1445 void Tenseur::fait_derive_cov (const Metrique& metre, int ind) const {
1446 
1447  assert (etat != ETATNONDEF) ;
1448  assert (valence != 0) ;
1449 
1450  if (p_derive_cov[ind] != 0x0)
1451  return ;
1452  else {
1453 
1454  p_derive_cov[ind] = new Tenseur (gradient()) ;
1455 
1456  if ((valence != 0) && (etat != ETATZERO)) {
1457 
1458 
1459  // assert( *(metre.gamma().get_triad()) == *triad ) ;
1460  Tenseur* auxi ;
1461  for (int i=0 ; i<valence ; i++) {
1462 
1463  if (type_indice(i) == COV) {
1464  auxi = new Tenseur(contract(metre.gamma(), 0,(*this), i)) ;
1465 
1466  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
1467  indices_gamma.set_etat_qcq() ;
1468  //On range comme il faut :
1469  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
1470 
1471  Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
1472  indices_gamma.set(0) = indices(0) ;
1473  indices_gamma.set(1) = indices(i+1) ;
1474  for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
1475  if (idx<=i+1)
1476  indices_gamma.set(idx) = indices(idx-1) ;
1477  else
1478  indices_gamma.set(idx) = indices(idx) ;
1479 
1480  p_derive_cov[ind]->set(indices) -= (*auxi)(indices_gamma) ;
1481  }
1482  }
1483  else {
1484  auxi = new Tenseur(contract(metre.gamma(), 1, (*this), i)) ;
1485 
1486  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
1487  indices_gamma.set_etat_qcq() ;
1488 
1489  //On range comme il faut :
1490  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
1491 
1492  Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
1493  indices_gamma.set(0) = indices(i+1) ;
1494  indices_gamma.set(1) = indices(0) ;
1495  for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
1496  if (idx<=i+1)
1497  indices_gamma.set(idx) = indices(idx-1) ;
1498  else
1499  indices_gamma.set(idx) = indices(idx) ;
1500  p_derive_cov[ind]->set(indices) += (*auxi)(indices_gamma) ;
1501  }
1502  }
1503  delete auxi ;
1504  }
1505  }
1506  if ((poids != 0.)&&(etat != ETATZERO))
1507  *p_derive_cov[ind] = *p_derive_cov[ind] -
1508  poids*contract(metre.gamma(), 0, 2) * (*this) ;
1509 
1510  }
1511 }
1512 
1513 
1514 
1515 void Tenseur::fait_derive_con (const Metrique& metre, int ind) const {
1516 
1517  if (p_derive_con[ind] != 0x0)
1518  return ;
1519  else {
1520  // On calcul la derivee covariante :
1521  if (valence != 0)
1522  p_derive_con[ind] = new Tenseur
1523  (contract(metre.con(), 1, derive_cov(metre), 0)) ;
1524 
1525  else {
1526  Tenseur grad (gradient()) ;
1527  grad.change_triad( *(metre.con().get_triad()) ) ;
1528  p_derive_con[ind] = new Tenseur (contract(metre.con(), 1, grad, 0)) ;
1529  }
1530  }
1531 }
1532 
1533 void Tenseur::fait_carre_scal (const Metrique& met, int ind) const {
1534 
1535  if (p_carre_scal[ind] != 0x0)
1536  return ;
1537  else {
1538  assert (valence != 0) ; // A ne pas appeler sur un scalaire ;
1539 
1540  // On bouge tous les indices :
1541  Tenseur op_t(manipule(*this, met)) ;
1542 
1543  Tenseur* auxi = new Tenseur(contract(*this, 0, op_t, 0)) ;
1544  Tenseur* auxi_old ;
1545 
1546  // On contracte tous les indices restant :
1547  for (int indice=1 ; indice<valence ; indice++) {
1548  auxi_old = new Tenseur(contract(*auxi, 0, valence-indice)) ;
1549  delete auxi ;
1550  auxi = new Tenseur(*auxi_old) ;
1551  delete auxi_old ;
1552  }
1553  p_carre_scal[ind] = new Tenseur (*auxi) ;
1554  delete auxi ;
1555  }
1556 }
1557 
1558 const Tenseur& Tenseur::gradient () const {
1559  if (p_gradient == 0x0)
1560  fait_gradient() ;
1561  return *p_gradient ;
1562 }
1563 
1565  if (p_gradient_spher == 0x0)
1567  return *p_gradient_spher ;
1568 }
1569 
1570 const Tenseur& Tenseur::derive_cov (const Metrique& metre) const {
1571 
1572  if (valence == 0)
1573  return gradient() ;
1574  else {
1575  set_dependance(metre) ;
1576  int j = get_place_met(metre) ;
1577  assert(j!=-1) ;
1578  if (p_derive_cov[j] == 0x0)
1579  fait_derive_cov (metre,j) ;
1580  return *p_derive_cov[j] ;
1581  }
1582 }
1583 
1584 const Tenseur& Tenseur::derive_con (const Metrique& metre) const {
1585  set_dependance(metre) ;
1586  int j = get_place_met(metre) ;
1587  assert(j!=-1) ;
1588  if (p_derive_con[j] == 0x0)
1589  fait_derive_con (metre, j) ;
1590  return *p_derive_con[j] ;
1591 }
1592 
1593 const Tenseur& Tenseur::carre_scal (const Metrique& metre) const {
1594  set_dependance(metre) ;
1595  int j = get_place_met(metre) ;
1596  assert(j!=-1) ;
1597  if (p_carre_scal[j] == 0x0)
1598  fait_carre_scal (metre, j) ;
1599  return *p_carre_scal[j] ;
1600 }
1601 }
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 Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:87
const Map *const mp
Reference mapping.
Definition: tenseur.h:309
double get_poids() const
Returns the weight.
Definition: tenseur.h:741
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tenseur.h:315
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
virtual ~Tenseur()
Destructor.
Definition: tenseur.C:549
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:916
int get_place_met(const Metrique &metre) const
Returns the position of the pointer on metre in the array met_depend .
Definition: tenseur.C:614
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
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1146
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1564
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:690
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Definition: map.h:825
int n_comp
Number of components, depending on the symmetry.
Definition: tenseur.h:323
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:157
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
void set_poids(double weight)
Sets the weight for a tensor density.
Definition: tenseur.C:696
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition: tenseur.h:1256
const Cmp & srstdsdp() const
Returns of *this .
Definition: cmp_deriv.C:130
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
int get_ndim() const
Gives the number of dimensions (ie dim.ndim )
Definition: itbl.h:323
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
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
virtual void sauve(FILE *) const
Save in a file.
Definition: base_vect.C:153
Base class for coordinate mappings.
Definition: map.h:688
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:324
const Cmp & srdsdt() const
Returns of *this .
Definition: cmp_deriv.C:108
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
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Definition: cmp_r_manip.C:106
void new_der_met()
Builds the arrays met_depend , p_derive_cov , p_derive_con and p_carre_scal and fills them with null ...
Definition: tenseur.C:212
Basic integer array class.
Definition: itbl.h:122
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
Definition: tenseur.C:1584
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
virtual void change_basis(Tenseur &) const =0
Change the basis in which the components of a tensor are expressed.
virtual void fait_gradient() const
Calculates, if needed, the gradient of *this .
Definition: tenseur.C:1361
void inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1133
const Cmp & operator()() const
Read only for a scalar.
Definition: tenseur.C:959
static Base_vect * bvect_from_file(FILE *)
Construction of a vectorial basis from a file (see sauve(FILE* ) ).
void inc2_dzpuis()
dzpuis += 2 ;
Definition: tenseur.C:1159
void del_t()
Logical destructor.
Definition: tenseur.C:561
virtual void operator=(const Tenseur &tens)
Assignment to another Tenseur.
Definition: tenseur.C:734
void set_metric(const Metrique &met)
Sets the pointer on the metric for a tensor density.
Definition: tenseur.C:701
void inc_dzpuis()
Increases by the value of dzpuis and changes accordingly the values of the Cmp in the external compac...
Definition: cmp_r_manip.C:169
Tenseur * p_gradient_spher
Pointer on the gradient of *this in a spherical orthonormal basis (scalar field only).
Definition: tenseur.h:353
void del_derive_met(int i) const
Logical destructor of the derivatives depending on the i-th element of *met_depend ...
Definition: tenseur.C:570
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:684
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.C:1445
void sauve(FILE *) const
Save in a file.
Definition: tenseur.C:1341
Tenseur * p_gradient
Pointer on the gradient of *this .
Definition: tenseur.h:346
void mult_r_zec()
Multiplication by r in the external zone.
Definition: tenseur.C:1172
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
bool verif() const
Returns false for a tensor density without a defined metric.
Definition: tenseur.C:208
void sauve(FILE *) const
Save in a file.
Definition: itbl.C:229
double poids
For tensor densities: the weight.
Definition: tenseur.h:326
Base_val ** std_base_vect_spher() const
Returns the standard spectral bases for the spherical components of a vector.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
Tenseur ** p_derive_con
Array of pointers on the contravariant derivatives of *this with respect to the corresponding metric...
Definition: tenseur.h:368
virtual void fait_derive_con(const Metrique &, int i) const
Calculates, if needed, the contravariant derivative of *this , with respect to met ...
Definition: tenseur.C:1515
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 .
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:183
void fait_carre_scal(const Metrique &, int i) const
Calculates, if needed, the scalar square of *this , with respect to met .
Definition: tenseur.C:1533
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
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
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:195
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
void set_der_0x0() const
Sets the pointers of all the derivatives to zero.
Definition: tenseur.C:607
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
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition: tenseur.C:225
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: cmp.C:326
Bases of the spectral expansions.
Definition: base_val.h:325
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: tenseur.C:673
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 dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1120
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
void set_der_met_0x0(int i) const
Sets the pointers of the derivatives depending on the i-th element of *met_depend to zero (as well as...
Definition: tenseur.C:599
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:328
friend Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
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
const Tenseur & carre_scal(const Metrique &) const
Returns the scalar square of *this , with respect to met .
Definition: tenseur.C:1593
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 del_derive() const
Logical destructor of all the derivatives.
Definition: tenseur.C:589
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
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
void fait_gradient_spher() const
Calculates, if needed, the gradient of *this in a spherical orthonormal basis (scalar field only)...
Definition: tenseur.C:1408
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