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