LORENE
tenseur_operateur.C
1/*
2 * Copyright (c) 1999-2001 Philippe Grandclement
3 * Copyright (c) 2000-2001 Eric Gourgoulhon
4 * Copyright (c) 2002 Jerome Novak
5 *
6 * This file is part of LORENE.
7 *
8 * LORENE is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * LORENE is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with LORENE; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 *
22 */
23
24
25char tenseur_operateur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $" ;
26
27/*
28 * $Id: tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $
29 * $Log: tenseur_operateur.C,v $
30 * Revision 1.10 2014/10/13 08:53:42 j_novak
31 * Lorene classes and functions now belong to the namespace Lorene.
32 *
33 * Revision 1.9 2014/10/06 15:13:19 j_novak
34 * Modified #include directives to use c++ syntax.
35 *
36 * Revision 1.8 2004/05/27 07:17:19 p_grandclement
37 * Correction of some shadowed variables
38 *
39 * Revision 1.7 2003/06/20 14:53:38 f_limousin
40 * Add the function contract_desal()
41 *
42 * Revision 1.6 2003/03/03 19:38:41 f_limousin
43 * Suppression of an assert on a metric associated with a tensor.
44 *
45 * Revision 1.5 2002/10/16 14:37:14 j_novak
46 * Reorganization of #include instructions of standard C++, in order to
47 * use experimental version 3 of gcc.
48 *
49 * Revision 1.4 2002/09/10 13:44:17 j_novak
50 * The method "manipule" of one indice has been removed for Tenseur_sym objects
51 * (the result cannot be a Tenseur_sym).
52 * The method "sans_trace" now computes the traceless part of a Tenseur (or
53 * Tenseur_sym) of valence 2.
54 *
55 * Revision 1.3 2002/09/06 14:49:25 j_novak
56 * Added method lie_derive for Tenseur and Tenseur_sym.
57 * Corrected various errors for derive_cov and arithmetic.
58 *
59 * Revision 1.2 2002/08/07 16:14:11 j_novak
60 * class Tenseur can now also handle tensor densities, this should be transparent to older codes
61 *
62 * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
63 * LORENE
64 *
65 * Revision 2.11 2001/08/27 10:04:21 eric
66 * Ajout de l'operator% (produit tensoriel avec desaliasing)
67 *
68 * Revision 2.10 2001/05/26 15:43:17 eric
69 * Ajout de la fonction flat_scalar_prod_desal (desaliasage)
70 *
71 * Revision 2.9 2000/10/06 15:08:40 eric
72 * Traitement des cas ETATZERO dans contract et flat_scal_prod
73 *
74 * Revision 2.8 2000/09/13 09:43:29 eric
75 * Modif skxk : appel des nouvelles fonctions Valeur::mult_cp() et
76 * Valeur::mult_sp() pour la multiplication par cos(phi) et sin(phi).
77 *
78 * Revision 2.7 2000/02/09 19:32:11 eric
79 * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
80 * argument des constructeurs.
81 *
82 * Revision 2.6 2000/02/01 14:14:25 eric
83 * Ajout de la fonction amie flat_scalar_prod.
84 *
85 * Revision 2.5 2000/01/21 12:59:18 phil
86 * ajout de skxk
87 *
88 * Revision 2.4 2000/01/14 09:29:43 eric
89 * *** empty log message ***
90 *
91 * Revision 2.3 2000/01/13 17:22:37 phil
92 * la fonction contraction de deux tenseurs ne passe plus par produit tensoriel
93 *
94 * Revision 2.2 2000/01/11 11:14:29 eric
95 * Changement de nom pour la base vectorielle : base --> triad
96 *
97 * Revision 2.1 2000/01/10 17:25:15 eric
98 * Gestion des bases vectorielles (triades de decomposition).
99 *
100 * Revision 2.0 1999/12/02 17:19:06 phil
101 * *** empty log message ***
102 *
103 *
104 * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $
105 *
106 */
107
108// Headers C
109#include <cstdlib>
110#include <cassert>
111#include <cmath>
112
113// Headers Lorene
114#include "tenseur.h"
115#include "metrique.h"
116
117
118namespace Lorene {
120
121 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
122 assert (t1.mp == t2.mp) ;
123
124 int val_res = t1.valence + t2.valence ;
125 double poids_res = t1.poids + t2.poids ;
126 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
127 const Metrique* met_res = 0x0 ;
128 if (poids_res != 0.) {
129 // assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
130 if (t1.metric != 0x0) met_res = t1.metric ;
131 else met_res = t2.metric ;
132 }
133
134 // cas scalaire :
135 if (val_res == 0) {
137 // cas ou un des deux est nul :
138 if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
139 scal.set_etat_zero() ;
140 else {
141 scal.set_etat_qcq() ;
142 scal.set() = t1() * t2() ;
143 }
144 return scal ;
145 }
146
147 else {
148
149 Itbl tipe (val_res) ;
150 tipe.set_etat_qcq() ;
151 for (int i=0 ; i<t1.valence ; i++)
152 tipe.set(i) = t1.type_indice(i) ;
153 for (int i=0 ; i<t2.valence ; i++)
154 tipe.set(i+t1.valence) = t2.type_indice(i) ;
155
156
157 if ( (t1.valence != 0) && (t2.valence != 0) ) {
158 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
159 }
160
161 const Base_vect* triad_res ;
162 if (t1.valence != 0) {
163 triad_res = t1.get_triad() ;
164 }
165 else{
166 triad_res = t2.get_triad() ;
167 }
168
170
171 if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
172 res.set_etat_zero() ;
173 else {
174 res.set_etat_qcq() ;
175 Itbl jeux_indice_t1 (t1.valence) ;
176 jeux_indice_t1.set_etat_qcq() ;
177 Itbl jeux_indice_t2 (t2.valence) ;
178 jeux_indice_t2.set_etat_qcq() ;
179
180 for (int i=0 ; i<res.n_comp ; i++) {
181 Itbl jeux_indice_res(res.donne_indices(i)) ;
182 for (int j=0 ; j<t1.valence ; j++)
184 for (int j=0 ; j<t2.valence ; j++)
185 jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
186
188 }
189 }
190 return res ;
191 }
192}
193
194 //------------------------------------//
195 // Tensorial product wiht desaliasing //
196 //------------------------------------//
197
198
200
201 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
202 assert (t1.mp == t2.mp) ;
203
204 int val_res = t1.valence + t2.valence ;
205 double poids_res = t1.poids + t2.poids ;
206 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
207 const Metrique* met_res = 0x0 ;
208 if (poids_res != 0.) {
209 // assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
210 if (t1.metric != 0x0) met_res = t1.metric ;
211 else met_res = t2.metric ;
212 }
213
214 // cas scalaire :
215 if (val_res == 0) {
217 // cas ou un des deux est nul :
218 if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
219 scal.set_etat_zero() ;
220 else {
221 scal.set_etat_qcq() ;
222 scal.set() = t1() % t2() ;
223 }
224 return scal ;
225 }
226
227 else {
228
229 Itbl tipe (val_res) ;
230 tipe.set_etat_qcq() ;
231 for (int i=0 ; i<t1.valence ; i++)
232 tipe.set(i) = t1.type_indice(i) ;
233 for (int i=0 ; i<t2.valence ; i++)
234 tipe.set(i+t1.valence) = t2.type_indice(i) ;
235
236
237 if ( (t1.valence != 0) && (t2.valence != 0) ) {
238 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
239 }
240
241 const Base_vect* triad_res ;
242 if (t1.valence != 0) {
243 triad_res = t1.get_triad() ;
244 }
245 else{
246 triad_res = t2.get_triad() ;
247 }
248
250
251
252
253 if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
254 res.set_etat_zero() ;
255 else {
256 res.set_etat_qcq() ;
257 Itbl jeux_indice_t1 (t1.valence) ;
258 jeux_indice_t1.set_etat_qcq() ;
259 Itbl jeux_indice_t2 (t2.valence) ;
260 jeux_indice_t2.set_etat_qcq() ;
261
262 for (int i=0 ; i<res.n_comp ; i++) {
263 Itbl jeux_indice_res(res.donne_indices(i)) ;
264 for (int j=0 ; j<t1.valence ; j++)
266 for (int j=0 ; j<t2.valence ; j++)
267 jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
268
271 }
272 }
273 return res ;
274 }
275}
276
277
278
280
281
282 // Les verifications :
283 assert (source.etat != ETATNONDEF) ;
284 assert ((ind_1 >= 0) && (ind_1 < source.valence)) ;
285 assert ((ind_2 >= 0) && (ind_2 < source.valence)) ;
286 assert (source.type_indice(ind_1) != source.type_indice(ind_2)) ;
287
288 // On veut ind_1 < ind_2 :
289 if (ind_1 > ind_2) {
290 int auxi = ind_2 ;
291 ind_2 = ind_1 ;
292 ind_1 = auxi ;
293 }
294
295 // On construit le resultat :
296 int val_res = source.valence - 2 ;
297
298 Itbl tipe (val_res) ;
299 tipe.set_etat_qcq() ;
300 for (int i=0 ; i<ind_1 ; i++)
301 tipe.set(i) = source.type_indice(i) ;
302 for (int i=ind_1 ; i<ind_2-1 ; i++)
303 tipe.set(i) = source.type_indice(i+1) ;
304 for (int i = ind_2-1 ; i<val_res ; i++)
305 tipe.set(i) = source.type_indice(i+2) ;
306
307 Tenseur res(*source.mp, val_res, tipe, source.triad, source.metric,
308 source.poids) ;
309
310 // Cas particulier d'une source nulle
311 if (source.etat == ETATZERO) {
312 res.set_etat_zero() ;
313 return res ;
314 }
315
316 res.set_etat_qcq() ;
317
318 Cmp work(source.mp) ;
319
320 // Boucle sur les composantes de res :
321
323 jeux_indice_source.set_etat_qcq() ;
324
325 for (int i=0 ; i<res.n_comp ; i++) {
326 Itbl jeux_indice_res (res.donne_indices(i)) ;
327 for (int j=0 ; j<ind_1 ; j++)
329 for (int j=ind_1+1 ; j<ind_2 ; j++)
331 for (int j=ind_2+1 ; j<source.valence ; j++)
333
334
335 work.set_etat_zero() ;
336 for (int j=0 ; j<3 ; j++) {
337 jeux_indice_source.set(ind_1) = j ;
338 jeux_indice_source.set(ind_2) = j ;
340 }
341
342 res.set(jeux_indice_res) = work ;
343 }
344 return res ;
345}
346
347
348Tenseur contract (const Tenseur& t1, int ind1, const Tenseur& t2, int ind2) {
349
350 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
351 // Verifs :
352 assert ((ind1>=0) && (ind1<t1.valence)) ;
353 assert ((ind2>=0) && (ind2<t2.valence)) ;
354 assert (*(t1.mp) == *(t2.mp)) ;
355
356 // Contraction possible ?
357 if ( (t1.valence != 0) && (t2.valence != 0) ) {
358 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
359 }
360 assert (t1.type_indice(ind1) != t2.type_indice(ind2)) ;
361
362 int val_res = t1.valence + t2.valence - 2;
363 double poids_res = t1.poids + t2.poids ;
364 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
365 const Metrique* met_res = 0x0 ;
366 if (poids_res != 0.) {
367 // assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
368 if (t1.metric != 0x0) met_res = t1.metric ;
369 else met_res = t2.metric ;
370 }
371 Itbl tipe(val_res) ;
372 tipe.set_etat_qcq() ;
373 for (int i=0 ; i<ind1 ; i++)
374 tipe.set(i) = t1.type_indice(i) ;
375 for (int i=ind1 ; i<t1.valence-1 ; i++)
376 tipe.set(i) = t1.type_indice(i+1) ;
377 for (int i=t1.valence-1 ; i<t1.valence+ind2-1 ; i++)
378 tipe.set(i) = t2.type_indice(i-t1.valence+1) ;
379 for (int i = t1.valence+ind2-1 ; i<val_res ; i++)
380 tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
381
382 const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
383
385
386 // Cas particulier ou l'un des deux tenseurs est nul
387 if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
388 res.set_etat_zero() ;
389 return res ;
390 }
391
392 res.set_etat_qcq() ;
393
394 Cmp work(t1.mp) ;
395
396 // Boucle sur les composantes de res :
397
398 Itbl jeux_indice_t1(t1.valence) ;
399 Itbl jeux_indice_t2(t2.valence) ;
400 jeux_indice_t1.set_etat_qcq() ;
401 jeux_indice_t2.set_etat_qcq() ;
402
403 for (int comp=0 ; comp<res.n_comp ; comp++) {
404 Itbl jeux_indice_res (res.donne_indices(comp)) ;
405 for (int i=0 ; i<ind1 ; i++)
407 for (int i=ind1+1 ; i<t1.valence ; i++)
409 for (int i=0 ; i<ind2 ; i++)
410 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-1) ;
411 for (int i=ind2+1 ; i<t2.valence ; i++)
412 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
413
414
415
416 work.set_etat_zero() ;
417 for (int j=0 ; j<3 ; j++) {
418 jeux_indice_t1.set(ind1) = j ;
419 jeux_indice_t2.set(ind2) = j ;
421 }
422
423 res.set(jeux_indice_res) = work ;
424 }
425 return res ;
426}
427
428Tenseur contract_desal (const Tenseur& t1, int ind1, const Tenseur& t2, int ind2) {
429
430 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
431 // Verifs :
432 assert ((ind1>=0) && (ind1<t1.valence)) ;
433 assert ((ind2>=0) && (ind2<t2.valence)) ;
434 assert (t1.mp == t2.mp) ;
435
436 // Contraction possible ?
437 if ( (t1.valence != 0) && (t2.valence != 0) ) {
438 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
439 }
440 assert (t1.type_indice(ind1) != t2.type_indice(ind2)) ;
441
442 int val_res = t1.valence + t2.valence - 2;
443 double poids_res = t1.poids + t2.poids ;
444 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
445 const Metrique* met_res = 0x0 ;
446 if (poids_res != 0.) {
447 // assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
448 if (t1.metric != 0x0) met_res = t1.metric ;
449 else met_res = t2.metric ;
450 }
451 Itbl tipe(val_res) ;
452 tipe.set_etat_qcq() ;
453 for (int i=0 ; i<ind1 ; i++)
454 tipe.set(i) = t1.type_indice(i) ;
455 for (int i=ind1 ; i<t1.valence-1 ; i++)
456 tipe.set(i) = t1.type_indice(i+1) ;
457 for (int i=t1.valence-1 ; i<t1.valence+ind2-1 ; i++)
458 tipe.set(i) = t2.type_indice(i-t1.valence+1) ;
459 for (int i = t1.valence+ind2-1 ; i<val_res ; i++)
460 tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
461
462 const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
463
464 Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
465
466 // Cas particulier ou l'un des deux tenseurs est nul
467 if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
468 res.set_etat_zero() ;
469 return res ;
470 }
471
472 res.set_etat_qcq() ;
473
474 Cmp work(t1.mp) ;
475
476 // Boucle sur les composantes de res :
477
478 Itbl jeux_indice_t1(t1.valence) ;
479 Itbl jeux_indice_t2(t2.valence) ;
480 jeux_indice_t1.set_etat_qcq() ;
481 jeux_indice_t2.set_etat_qcq() ;
482
483 for (int comp=0 ; comp<res.n_comp ; comp++) {
484 Itbl jeux_indice_res (res.donne_indices(comp)) ;
485 for (int i=0 ; i<ind1 ; i++)
486 jeux_indice_t1.set(i) = jeux_indice_res(i) ;
487 for (int i=ind1+1 ; i<t1.valence ; i++)
488 jeux_indice_t1.set(i) = jeux_indice_res(i-1) ;
489 for (int i=0 ; i<ind2 ; i++)
490 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-1) ;
491 for (int i=ind2+1 ; i<t2.valence ; i++)
492 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
493
494
495
496 work.set_etat_zero() ;
497 for (int j=0 ; j<3 ; j++) {
498 jeux_indice_t1.set(ind1) = j ;
499 jeux_indice_t2.set(ind2) = j ;
500 work = work + t1(jeux_indice_t1)%t2(jeux_indice_t2) ;
501 }
502
503 res.set(jeux_indice_res) = work ;
504 }
505 return res ;
506}
507
508
509Tenseur manipule(const Tenseur& t1, const Metrique& met, int place) {
510
511 assert (t1.etat != ETATNONDEF) ;
512 assert (met.get_etat() != ETATNONDEF) ;
513
514 int valen = t1.valence ;
515 assert (valen != 0) ; // Aucun interet pour un scalaire...
516 assert ((place >=0) && (place < valen)) ;
517
518 Itbl tipe (valen) ;
519 tipe.set_etat_qcq() ;
520 tipe.set(0) = -t1.type_indice(place) ;
521 for (int i=1 ; i<place+1 ; i++)
522 tipe.set(i) = t1.type_indice(i-1) ;
523 for (int i=place+1 ; i<valen ; i++)
524 tipe.set(i) = t1.type_indice(i) ;
525
526 Tenseur auxi(*t1.mp, valen, tipe, t1.triad) ;
527
528 if (t1.type_indice(place) == COV)
529 auxi = contract (met.con(), 1, t1, place) ;
530 else
531 auxi = contract (met.cov(), 1, t1, place) ;
532
533 // On doit remettre les indices a la bonne place ...
534
535 for (int i=0 ; i<valen ; i++)
536 tipe.set(i) = t1.type_indice(i) ;
537 tipe.set(place) *= -1 ;
538
539 Tenseur res(*t1.mp, valen, tipe, t1.triad, auxi.metric, auxi.poids) ;
540 res.set_etat_qcq() ;
541
543 place_auxi.set_etat_qcq() ;
544
545 for (int i=0 ; i<res.n_comp ; i++) {
546
547 Itbl place_res (res.donne_indices(i)) ;
548
549 place_auxi.set(0) = place_res(place) ;
550 for (int j=1 ; j<place+1 ; j++)
551 place_auxi.set(j) = place_res(j-1) ;
552 place_res.set(place) = place_auxi(0) ;
553 for (int j=place+1 ; j<valen ; j++)
554 place_auxi.set(j) = place_res(j);
555
556
557 res.set(place_res) = auxi(place_auxi) ;
558 }
559 return res ;
560}
561
562Tenseur manipule (const Tenseur& t1, const Metrique& met) {
563
564 Tenseur* auxi ;
565 Tenseur* auxi_old = new Tenseur(t1) ;
566
567 for (int i=0 ; i<t1.valence ; i++) {
568 auxi = new Tenseur(manipule(*auxi_old, met, i)) ;
569 delete auxi_old ;
570 auxi_old = new Tenseur(*auxi) ;
571 delete auxi ;
572 }
573
575 delete auxi_old ;
576 return result ;
577}
578
579
581
582 // Verification
583 assert (source.valence > 0) ;
584 assert (source.etat != ETATNONDEF) ;
585 assert (*source.triad == source.mp->get_bvect_cart()) ;
586
587 // Le resultat :
588 int val_res = source.valence-1 ;
589 Itbl tipe (val_res) ;
590 tipe.set_etat_qcq() ;
591 for (int i=0 ; i<val_res ; i++)
592 tipe.set(i) = source.type_indice(i) ;
593
594
595 Tenseur res (*source.mp, val_res, tipe, source.triad, source.metric,
596 source.poids) ;
597
598 if (source.etat == ETATZERO)
599 res.set_etat_zero() ;
600 else {
601 res.set_etat_qcq() ;
602
603 for (int i=0 ; i<res.n_comp ; i++) {
604 Itbl indices_res (res.donne_indices(i)) ;
606 indices_so.set_etat_qcq() ;
607 for (int j=0 ; j<val_res ; j++)
608 indices_so.set(j) = indices_res(j) ;
609 // x S_x
610 // -----
611 indices_so.set(val_res) = 0 ;
613
614 resu.mult_r() ; // Multipl. by r
615
616 // What follows is valid only for a mapping of class Map_radial :
617 assert( dynamic_cast<const Map_radial*>(source.get_mp()) != 0x0) ;
618
619 resu.va = (resu.va).mult_st() ; // Multipl. by sin(theta)
620 resu.va = (resu.va).mult_cp() ; // Multipl. by cos(phi)
621
622 // y S_y
623 // -----
624 indices_so.set(val_res) = 1 ;
626
627 auxiliaire.mult_r() ; // Multipl. by r
628
629 auxiliaire.va = (auxiliaire.va).mult_st() ; // Multipl. by sin(theta)
630 auxiliaire.va = (auxiliaire.va).mult_sp() ; // Multipl. by sin(phi)
631
632 resu = resu + auxiliaire ;
633
634 // z S_z
635 // -----
636 indices_so.set(val_res) = 2 ;
638
639 auxiliaire.mult_r() ; // Multipl. by r
640
641 auxiliaire.va = (auxiliaire.va).mult_ct() ; // Multipl. by cos(theta)
642
643 resu = resu + auxiliaire ;
644
645 res.set(indices_res) = resu ;
646 // The End
647 // -------
648 }
649 }
650 return res ;
651}
652
654
655 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
656 // Verifs :
657 assert (t1.mp == t2.mp) ;
658
659 // Contraction possible ?
660 if ( (t1.valence != 0) && (t2.valence != 0) ) {
661 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
662 }
663
664 int val_res = t1.valence + t2.valence - 2;
665 double poids_res = t1.poids + t2.poids ;
666 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
667 const Metrique* met_res = 0x0 ;
668 if (poids_res != 0.) {
669 assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
670 if (t1.metric != 0x0) met_res = t1.metric ;
671 else met_res = t2.metric ;
672 }
673 Itbl tipe(val_res) ;
674 tipe.set_etat_qcq() ;
675
676 for (int i=0 ; i<t1.valence - 1 ; i++)
677 tipe.set(i) = t1.type_indice(i) ;
678 for (int i = t1.valence-1 ; i<val_res ; i++)
679 tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
680
681 Tenseur res(*t1.mp, val_res, tipe, t1.triad, met_res, poids_res) ;
682
683 // Cas particulier ou l'un des deux tenseurs est nul
684 if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
685 res.set_etat_zero() ;
686 return res ;
687 }
688
689 res.set_etat_qcq() ;
690
691 Cmp work(t1.mp) ;
692
693 // Boucle sur les composantes de res :
694
695 Itbl jeux_indice_t1(t1.valence) ;
696 Itbl jeux_indice_t2(t2.valence) ;
697 jeux_indice_t1.set_etat_qcq() ;
698 jeux_indice_t2.set_etat_qcq() ;
699
700 for (int ir=0 ; ir<res.n_comp ; ir++) { // Boucle sur les composantes
701 // du resultat
702
703 // Indices du resultat correspondant a la position ir :
704 Itbl jeux_indice_res = res.donne_indices(ir) ;
705
706 // Premiers indices correspondant dans t1 :
707 for (int i=0 ; i<t1.valence - 1 ; i++) {
709 }
710
711 // Derniers indices correspondant dans t2 :
712 for (int i=1 ; i<t2.valence ; i++) {
713 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
714 }
715
716 work.set_etat_zero() ;
717
718 // Sommation sur le dernier indice de t1 et le premier de t2 :
719
720 for (int j=0 ; j<3 ; j++) {
721 jeux_indice_t1.set(t1.valence - 1) = j ;
722 jeux_indice_t2.set(0) = j ;
724 }
725
726 res.set(jeux_indice_res) = work ;
727
728 } // fin de la boucle sur les composantes du resultat
729
730 return res ;
731}
732
733
734
736
737 assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
738 // Verifs :
739 assert (t1.mp == t2.mp) ;
740
741 // Contraction possible ?
742 if ( (t1.valence != 0) && (t2.valence != 0) ) {
743 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
744 }
745
746 int val_res = t1.valence + t2.valence - 2;
747 double poids_res = t1.poids + t2.poids ;
748 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
749 const Metrique* met_res = 0x0 ;
750 if (poids_res != 0.) {
751 assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
752 if (t1.metric != 0x0) met_res = t1.metric ;
753 else met_res = t2.metric ;
754 }
755 Itbl tipe(val_res) ;
756 tipe.set_etat_qcq() ;
757
758 for (int i=0 ; i<t1.valence - 1 ; i++)
759 tipe.set(i) = t1.type_indice(i) ;
760 for (int i = t1.valence-1 ; i<val_res ; i++)
761 tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
762
763 Tenseur res(*t1.mp, val_res, tipe, t1.triad, met_res, poids_res) ;
764
765 // Cas particulier ou l'un des deux tenseurs est nul
766 if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
767 res.set_etat_zero() ;
768 return res ;
769 }
770
771 res.set_etat_qcq() ;
772
773 Cmp work(t1.mp) ;
774
775 // Boucle sur les composantes de res :
776
777 Itbl jeux_indice_t1(t1.valence) ;
778 Itbl jeux_indice_t2(t2.valence) ;
779 jeux_indice_t1.set_etat_qcq() ;
780 jeux_indice_t2.set_etat_qcq() ;
781
782 for (int ir=0 ; ir<res.n_comp ; ir++) { // Boucle sur les composantes
783 // du resultat
784
785 // Indices du resultat correspondant a la position ir :
786 Itbl jeux_indice_res = res.donne_indices(ir) ;
787
788 // Premiers indices correspondant dans t1 :
789 for (int i=0 ; i<t1.valence - 1 ; i++) {
791 }
792
793 // Derniers indices correspondant dans t2 :
794 for (int i=1 ; i<t2.valence ; i++) {
795 jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
796 }
797
798 work.set_etat_zero() ;
799
800 // Sommation sur le dernier indice de t1 et le premier de t2 :
801
802 for (int j=0 ; j<3 ; j++) {
803 jeux_indice_t1.set(t1.valence - 1) = j ;
804 jeux_indice_t2.set(0) = j ;
806 }
807
808 res.set(jeux_indice_res) = work ;
809
810 } // fin de la boucle sur les composantes du resultat
811
812 return res ;
813}
814
815
816Tenseur lie_derive (const Tenseur& t, const Tenseur& x, const Metrique* met)
817{
818 assert ( (t.get_etat() != ETATNONDEF) && (x.get_etat() != ETATNONDEF) ) ;
819 assert(x.get_valence() == 1) ;
820 assert(x.get_type_indice(0) == CON) ;
821 assert(x.get_poids() == 0.) ;
822 assert(t.get_mp() == x.get_mp()) ;
823
824 int val = t.get_valence() ;
825 double poids = t.get_poids() ;
826 Itbl tipe (val+1) ;
827 tipe.set_etat_qcq() ;
828 tipe.set(0) = COV ;
829 Itbl tipx(2) ;
830 tipx.set_etat_qcq() ;
831 tipx.set(0) = COV ;
832 tipx.set(1) = CON ;
833 for (int i=0 ; i<val ; i++)
834 tipe.set(i+1) = t.get_type_indice(i) ;
835 Tenseur dt(*t.get_mp(), val+1, tipe, t.get_triad(), t.get_metric(), poids) ;
836 Tenseur dx(*x.get_mp(), 2, tipx, x.get_triad()) ;
837 if (met == 0x0) {
838 dx = x.gradient() ;
839 dt = t.gradient() ;
840 }
841 else {
842 dx = x.derive_cov(*met) ;
843 dt = t.derive_cov(*met) ;
844 }
845 Tenseur resu(contract(x,0,dt,0)) ;
846 Tenseur* auxi ;
847 if ((val!=0)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO)) {
848 assert(t.get_triad()->identify() == x.get_triad()->identify()) ;
849
850 for (int i=0 ; i<val ; i++) {
851 if (t.get_type_indice(i) == COV) {
852 auxi = new Tenseur(contract(t,i,dx,1)) ;
853
854 Itbl indices_aux(val) ;
855 indices_aux.set_etat_qcq() ;
856 for (int j=0 ; j<resu.get_n_comp() ; j++) {
857
858 Itbl indices (resu.donne_indices(j)) ;
859 indices_aux.set(val-1) = indices(i) ;
860 for (int idx=0 ; idx<val-1 ; idx++)
861 if (idx<i)
862 indices_aux.set(idx) = indices(idx) ;
863 else
864 indices_aux.set(idx) = indices(idx+1) ;
865
866 resu.set(indices) += (*auxi)(indices_aux) ;
867 }
868 }
869 else {
870 auxi = new Tenseur(contract(t,i,dx,0)) ;
871
872 Itbl indices_aux(val) ;
873 indices_aux.set_etat_qcq() ;
874
875 //On range comme il faut :
876 for (int j=0 ; j<resu.get_n_comp() ; j++) {
877
878 Itbl indices (resu.donne_indices(j)) ;
879 indices_aux.set(val-1) = indices(i) ;
880 for (int idx=0 ; idx<val-1 ; idx++)
881 if (idx<i)
882 indices_aux.set(idx) = indices(idx) ;
883 else
884 indices_aux.set(idx) = indices(idx+1) ;
885 resu.set(indices) -= (*auxi)(indices_aux) ;
886 }
887 }
888 delete auxi ;
889 }
890 }
891 if ((poids != 0.)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO))
892 resu = resu + poids*contract(dx,0,1)*t ;
893
894 return resu ;
895}
896
897Tenseur sans_trace(const Tenseur& t, const Metrique& metre)
898{
899 assert(t.get_etat() != ETATNONDEF) ;
900 assert(metre.get_etat() != ETATNONDEF) ;
901 assert(t.get_valence() == 2) ;
902
903 Tenseur resu(t) ;
904 if (resu.get_etat() == ETATZERO) return resu ;
905 assert(resu.get_etat() == ETATQCQ) ;
906
907 int t0 = t.get_type_indice(0) ;
908 int t1 = t.get_type_indice(1) ;
909 Itbl mix(2) ;
910 mix.set_etat_qcq() ;
911 mix.set(0) = (t0 == t1 ? -t0 : t0) ;
912 mix.set(1) = t1 ;
913
914 Tenseur tmp(*t.get_mp(), 2, mix, *t.get_triad(), t.get_metric(),
915 t.get_poids()) ;
916 if (t0 == t1)
917 tmp = manipule(t, metre, 0) ;
918 else
919 tmp = t ;
920
921 Tenseur trace(contract(tmp, 0, 1)) ;
922
923 if (t0 == t1) {
924 switch (t0) {
925 case COV : {
926 resu = resu - 1./3.*trace * metre.cov() ;
927 break ;
928 }
929 case CON : {
930 resu = resu - 1./3.*trace * metre.con() ;
931 break ;
932 }
933 default :
934 cout << "Erreur bizarre dans sans_trace!" << endl ;
935 abort() ;
936 break ;
937 }
938 }
939 else {
940 Tenseur_sym delta(*t.get_mp(), 2, mix, *t.get_triad(),
941 t.get_metric(), t.get_poids()) ;
942 delta.set_etat_qcq() ;
943 for (int i=0; i<3; i++)
944 for (int j=i; j<3; j++)
945 delta.set(i,j) = (i==j ? 1 : 0) ;
946 resu = resu - trace/3. * delta ;
947 }
948 resu.set_std_base() ;
949 return resu ;
950}
951
952
953
954}
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition itbl.h:247
Base class for pure radial mappings.
Definition map.h:1536
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition tenseur.h:325
const Map *const mp
Reference mapping.
Definition tenseur.h:306
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
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:699
double get_poids() const
Returns the weight.
Definition tenseur.h:738
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
Itbl type_indice
Array of size valence contening the type of each index, COV for a covariant one and CON for a contrav...
Definition tenseur.h:318
int valence
Valence.
Definition tenseur.h:307
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition tenseur.h:321
double poids
For tensor densities: the weight.
Definition tenseur.h:323
int get_valence() const
Returns the valence.
Definition tenseur.h:710
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition tenseur.h:745
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition cmp_arithm.C:364
Tenseur lie_derive(const Tenseur &t, const Tenseur &x, const Metrique *=0x0)
Lie Derivative of t with respect to x .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Tenseur sans_trace(const Tenseur &tens, const Metrique &metre)
Computes the traceless part of a Tenseur of valence 2.
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Lorene prototypes.
Definition app_hor.h:64