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 
25 char 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 
118 namespace Lorene {
119 Tenseur operator*(const Tenseur& t1, const Tenseur& t2) {
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) {
136  Tenseur scal(*t1.mp, met_res, poids_res) ;
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 
169  Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
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++)
183  jeux_indice_t1.set(j) = jeux_indice_res(j) ;
184  for (int j=0 ; j<t2.valence ; j++)
185  jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
186 
187  res.set(jeux_indice_res) = t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
188  }
189  }
190  return res ;
191  }
192 }
193 
194  //------------------------------------//
195  // Tensorial product wiht desaliasing //
196  //------------------------------------//
197 
198 
199 Tenseur operator%(const Tenseur& t1, const Tenseur& t2) {
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) {
216  Tenseur scal(*t1.mp, met_res, poids_res) ;
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 
249  Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
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++)
265  jeux_indice_t1.set(j) = jeux_indice_res(j) ;
266  for (int j=0 ; j<t2.valence ; j++)
267  jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
268 
269  res.set(jeux_indice_res) = t1(jeux_indice_t1) %
270  t2(jeux_indice_t2) ;
271  }
272  }
273  return res ;
274  }
275 }
276 
277 
278 
279 Tenseur contract(const Tenseur& source, int ind_1, int ind_2) {
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 
322  Itbl jeux_indice_source(source.valence) ;
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++)
328  jeux_indice_source.set(j) = jeux_indice_res(j) ;
329  for (int j=ind_1+1 ; j<ind_2 ; j++)
330  jeux_indice_source.set(j) = jeux_indice_res(j-1) ;
331  for (int j=ind_2+1 ; j<source.valence ; j++)
332  jeux_indice_source.set(j) = jeux_indice_res(j-2) ;
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 ;
339  work = work + source(jeux_indice_source) ;
340  }
341 
342  res.set(jeux_indice_res) = work ;
343  }
344  return res ;
345 }
346 
347 
348 Tenseur 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 
384  Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
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++)
406  jeux_indice_t1.set(i) = jeux_indice_res(i) ;
407  for (int i=ind1+1 ; i<t1.valence ; i++)
408  jeux_indice_t1.set(i) = jeux_indice_res(i-1) ;
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 ;
420  work = work + t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
421  }
422 
423  res.set(jeux_indice_res) = work ;
424  }
425  return res ;
426 }
427 
428 Tenseur 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 
509 Tenseur 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 
542  Itbl place_auxi(valen) ;
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 
562 Tenseur 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 
574  Tenseur result(*auxi_old) ;
575  delete auxi_old ;
576  return result ;
577 }
578 
579 
580 Tenseur skxk(const Tenseur& source) {
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)) ;
605  Itbl indices_so (val_res+1) ;
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 ;
612  Cmp resu(source(indices_so)) ;
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 ;
625  Cmp auxiliaire (source(indices_so)) ;
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 ;
637  auxiliaire = source(indices_so) ;
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 
653 Tenseur flat_scalar_prod(const Tenseur& t1, const Tenseur& t2) {
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++) {
708  jeux_indice_t1.set(i) = jeux_indice_res(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 ;
723  work = work + t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
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++) {
790  jeux_indice_t1.set(i) = jeux_indice_res(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 ;
805  work = work + t1(jeux_indice_t1) % t2(jeux_indice_t2) ;
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 
816 Tenseur 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 
897 Tenseur 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 }
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
const Map *const mp
Reference mapping.
Definition: tenseur.h:306
double get_poids() const
Returns the weight.
Definition: tenseur.h:738
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tenseur.h:312
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:726
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
int n_comp
Number of components, depending on the symmetry.
Definition: tenseur.h:320
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition: tenseur.h:1253
Lorene prototypes.
Definition: app_hor.h:64
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:321
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...
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
Basic integer array class.
Definition: itbl.h:122
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition: cmp_arithm.C:364
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
int get_valence() const
Returns the valence.
Definition: tenseur.h:710
void mult_r()
Multiplication by r everywhere.
Definition: cmp_r_manip.C:91
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:699
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 .
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition: tenseur.h:745
int valence
Valence.
Definition: tenseur.h:307
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
Tenseur lie_derive(const Tenseur &t, const Tenseur &x, const Metrique *=0x0)
Lie Derivative of t with respect to x .
double poids
For tensor densities: the weight.
Definition: tenseur.h:323
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: itbl.C:261
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
int get_n_comp() const
Returns the number of components.
Definition: tenseur.h:713
Tenseur sans_trace(const Tenseur &tens, const Metrique &metre)
Computes the traceless part of a Tenseur of valence 2.
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:325
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
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_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .