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
32char 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//---------------------------------------
191namespace Lorene {
192bool 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// ------------------------------
209Tenseur::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// ---------------------------------------------------
224Tenseur::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// --------------------
246Tenseur::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// -------------------------------------------------------
268Tenseur::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// --------------------------------------------------------------
297Tenseur::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// ----------------
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
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// -----------------------------------
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
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// -----------------------
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) {
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// ------------------------------------------
448Tenseur::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// ---------------------------------------
485Tenseur::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// -------------------------------------------------------------------
509Tenseur::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
554void 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] ;
569 }
570}
571
572
573void Tenseur::del_derive () const {
574 for (int i=0; i<N_MET_MAX; 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
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
592 for (int i=0; i<N_MET_MAX; i++)
594 p_gradient = 0x0 ;
595 p_gradient_spher = 0x0 ;
596}
597
598int 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
608void 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
675
676 triad = &bi ;
677
678}
679
681
682 poids = weight ;
683}
684
685void Tenseur::set_metric(const Metrique& met) {
686
687 metric = &met ;
688}
689
690int 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
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
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: {
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++) {
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
762
763 assert (valence == 0) ;
764 poids = 0. ;
765 metric = 0x0 ;
766
767 switch (ci.get_etat()) {
768 case ETATNONDEF: {
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
792void 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
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 :
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 :
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 :
864Cmp& 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
884Cmp& 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
901
902 annule(l, l) ;
903}
904
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) ;
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 :
943const 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
972const 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
1000const 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
1034const 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
1071const 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 :
1256ostream& 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
1325void 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
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
1376 indices_source.set_etat_qcq() ;
1377
1378 for (int j=0 ; j<p_gradient->n_comp ; j++) {
1379
1381
1382 for (int m=0 ; m<valence ; m++)
1383 indices_source.set(m) = indices_res(m+1) ;
1384
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
1429void 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
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
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))
1492 poids*contract(metre.gamma(), 0, 2) * (*this) ;
1493
1494 }
1495}
1496
1497
1498
1499void 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
1517void 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++) {
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
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
1554const Tenseur& Tenseur::derive_cov (const Metrique& metre) const {
1555
1556 if (valence == 0)
1557 return gradient() ;
1558 else {
1560 int j = get_place_met(metre) ;
1561 assert(j!=-1) ;
1562 if (p_derive_cov[j] == 0x0)
1564 return *p_derive_cov[j] ;
1565 }
1566}
1567
1568const Tenseur& Tenseur::derive_con (const Metrique& metre) const {
1570 int j = get_place_met(metre) ;
1571 assert(j!=-1) ;
1572 if (p_derive_con[j] == 0x0)
1574 return *p_derive_con[j] ;
1575}
1576
1577const Tenseur& Tenseur::carre_scal (const Metrique& metre) const {
1579 int j = get_place_met(metre) ;
1580 assert(j!=-1) ;
1581 if (p_carre_scal[j] == 0x0)
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
static Base_vect * bvect_from_file(FILE *)
Construction of a vectorial basis from a file (see sauve(FILE* ) ).
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)
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...
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...
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...
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...
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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 get_dim(int i) const
Gives the i th dimension (ie {tt dim.dim[i] )
Definition itbl.h:326
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition itbl.h:247
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
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition tenseur.h:1253
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 .
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
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