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