LORENE
tenseur_sym_operateur.C
1 /*
2  * Copyright (c) 1999-2001 Philippe Grandclement
3  * Copyright (c) 1999-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_sym_operateur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $" ;
26 
27 /*
28  * $Id: tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $
29  * $Log: tenseur_sym_operateur.C,v $
30  * Revision 1.8 2014/10/13 08:53:43 j_novak
31  * Lorene classes and functions now belong to the namespace Lorene.
32  *
33  * Revision 1.7 2014/10/06 15:13:19 j_novak
34  * Modified #include directives to use c++ syntax.
35  *
36  * Revision 1.6 2003/03/03 19:41:34 f_limousin
37  * Suppression of an assert on a metric associated with a tensor.
38  *
39  * Revision 1.5 2002/10/16 14:37:15 j_novak
40  * Reorganization of #include instructions of standard C++, in order to
41  * use experimental version 3 of gcc.
42  *
43  * Revision 1.4 2002/09/10 13:44:17 j_novak
44  * The method "manipule" of one indice has been removed for Tenseur_sym objects
45  * (the result cannot be a Tenseur_sym).
46  * The method "sans_trace" now computes the traceless part of a Tenseur (or
47  * Tenseur_sym) of valence 2.
48  *
49  * Revision 1.3 2002/09/06 14:49:26 j_novak
50  * Added method lie_derive for Tenseur and Tenseur_sym.
51  * Corrected various errors for derive_cov and arithmetic.
52  *
53  * Revision 1.2 2002/08/07 16:14:12 j_novak
54  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
55  *
56  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
57  * LORENE
58  *
59  * Revision 2.2 2000/02/09 19:32:22 eric
60  * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
61  * argument des constructeurs.
62  *
63  * Revision 2.1 2000/01/11 11:15:08 eric
64  * Gestion de la base vectorielle (triad).
65  *
66  * Revision 2.0 1999/12/02 17:19:02 phil
67  * *** empty log message ***
68  *
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $
71  *
72  */
73 
74 // Headers C
75 #include <cstdlib>
76 #include <cassert>
77 #include <cmath>
78 
79 // Headers Lorene
80 #include "tenseur.h"
81 #include "metrique.h"
82 
83 namespace Lorene {
84 Tenseur_sym operator*(const Tenseur& t1, const Tenseur_sym& t2) {
85 
86  assert ((t1.get_etat() != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
87  assert (t1.get_mp() == t2.mp) ;
88 
89  int val_res = t1.get_valence() + t2.valence ;
90  double poids_res = t1.get_poids() + t2.poids ;
91  poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
92  const Metrique* met_res = 0x0 ;
93  if (poids_res != 0.) {
94  // assert((t1.get_metric() != 0x0) || (t2.metric != 0x0)) ;
95  if (t1.get_metric() != 0x0) met_res = t1.get_metric() ;
96  else met_res = t2.metric ;
97  }
98 
99  Itbl tipe (val_res) ;
100  tipe.set_etat_qcq() ;
101  for (int i=0 ; i<t1.get_valence() ; i++)
102  tipe.set(i) = t1.get_type_indice(i) ;
103  for (int i=0 ; i<t2.valence ; i++)
104  tipe.set(i+t1.get_valence()) = t2.type_indice(i) ;
105 
106 
107  if ( t1.get_valence() != 0 ) {
108  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
109  }
110 
111  Tenseur_sym res(*t1.get_mp(), val_res, tipe, *(t2.get_triad()),
112  met_res, poids_res) ;
113 
114 
115  if ((t1.get_etat() == ETATZERO) || (t2.etat == ETATZERO))
116  res.set_etat_zero() ;
117  else {
118  res.set_etat_qcq() ;
119  Itbl jeux_indice_t1 (t1.get_valence()) ;
120  jeux_indice_t1.set_etat_qcq() ;
121  Itbl jeux_indice_t2 (t2.valence) ;
122  jeux_indice_t2.set_etat_qcq() ;
123 
124  for (int i=0 ; i<res.n_comp ; i++) {
125  Itbl jeux_indice_res(res.donne_indices(i)) ;
126  for (int j=0 ; j<t1.get_valence() ; j++)
127  jeux_indice_t1.set(j) = jeux_indice_res(j) ;
128  for (int j=0 ; j<t2.valence ; j++)
129  jeux_indice_t2.set(j) = jeux_indice_res(j+t1.get_valence()) ;
130 
131  res.set(jeux_indice_res) = t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
132  }
133  }
134  return res ;
135 }
136 
137 Tenseur_sym manipule (const Tenseur_sym& t1, const Metrique& met) {
138 
139  Tenseur* auxi ;
140  Tenseur* auxi_old = new Tenseur(t1) ;
141 
142  for (int i=0 ; i<t1.valence ; i++) {
143  auxi = new Tenseur(manipule(*auxi_old, met, i)) ;
144  delete auxi_old ;
145  auxi_old = new Tenseur(*auxi) ;
146  delete auxi ;
147  }
148 
149  Tenseur_sym result(*auxi_old) ;
150  delete auxi_old ;
151  return result ;
152 }
153 
155  const Metrique* met)
156 {
157  assert ( (t.get_etat() != ETATNONDEF) && (x.get_etat() != ETATNONDEF) ) ;
158  assert(x.get_valence() == 1) ;
159  assert(x.get_type_indice(0) == CON) ;
160  assert(x.get_poids() == 0.) ;
161  assert(t.get_mp() == x.get_mp()) ;
162 
163  int val = t.get_valence() ;
164  double poids = t.get_poids() ;
165  Itbl tipe (val+1) ;
166  tipe.set_etat_qcq() ;
167  tipe.set(0) = COV ;
168  Itbl tipx(2) ;
169  tipx.set_etat_qcq() ;
170  tipx.set(0) = COV ;
171  tipx.set(1) = CON ;
172  for (int i=0 ; i<val ; i++)
173  tipe.set(i+1) = t.get_type_indice(i) ;
174  Tenseur_sym dt(*t.get_mp(), val+1, tipe, *t.get_triad(), t.get_metric(),
175  poids) ;
176  Tenseur dx(*x.get_mp(), 2, tipx, x.get_triad()) ;
177  if (met == 0x0) {
178  dx = x.gradient() ;
179  dt = t.gradient() ;
180  }
181  else {
182  dx = x.derive_cov(*met) ;
183  dt = t.derive_cov(*met) ;
184  }
185  Tenseur_sym resu(contract(x,0,dt,0)) ;
186  Tenseur* auxi ;
187  if ((val!=0)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO)) {
188  assert(t.get_triad()->identify() == x.get_triad()->identify()) ;
189 
190  for (int i=0 ; i<val ; i++) {
191  if (t.get_type_indice(i) == COV) {
192  auxi = new Tenseur(contract(t,i,dx,1)) ;
193 
194  Itbl indices_aux(val) ;
195  indices_aux.set_etat_qcq() ;
196  for (int j=0 ; j<resu.get_n_comp() ; j++) {
197 
198  Itbl indices (resu.donne_indices(j)) ;
199  indices_aux.set(val-1) = indices(i) ;
200  for (int idx=0 ; idx<val-1 ; idx++)
201  if (idx<i)
202  indices_aux.set(idx) = indices(idx) ;
203  else
204  indices_aux.set(idx) = indices(idx+1) ;
205 
206  resu.set(indices) += (*auxi)(indices_aux) ;
207  }
208  }
209  else {
210  auxi = new Tenseur(contract(t,i,dx,0)) ;
211 
212  Itbl indices_aux(val) ;
213  indices_aux.set_etat_qcq() ;
214 
215  //On range comme il faut :
216  for (int j=0 ; j<resu.get_n_comp() ; j++) {
217 
218  Itbl indices (resu.donne_indices(j)) ;
219  indices_aux.set(val-1) = indices(i) ;
220  for (int idx=0 ; idx<val-1 ; idx++)
221  if (idx<i)
222  indices_aux.set(idx) = indices(idx) ;
223  else
224  indices_aux.set(idx) = indices(idx+1) ;
225  resu.set(indices) -= (*auxi)(indices_aux) ;
226  }
227  }
228  delete auxi ;
229  }
230  }
231  if ((poids != 0.)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO))
232  resu = resu + poids*contract(dx,0,1)*t ;
233 
234  return resu ;
235 }
236 
237 Tenseur_sym sans_trace(const Tenseur_sym& t, const Metrique& metre)
238 {
239  assert(t.get_etat() != ETATNONDEF) ;
240  assert(metre.get_etat() != ETATNONDEF) ;
241  assert(t.get_valence() == 2) ;
242 
243  Tenseur_sym resu(t) ;
244  if (resu.get_etat() == ETATZERO) return resu ;
245  assert(resu.get_etat() == ETATQCQ) ;
246 
247  int t0 = t.get_type_indice(0) ;
248  int t1 = t.get_type_indice(1) ;
249  Itbl mix(2) ;
250  mix.set_etat_qcq() ;
251  mix.set(0) = (t0 == t1 ? -t0 : t0) ;
252  mix.set(1) = t1 ;
253 
254  Tenseur tmp(*t.get_mp(), 2, mix, *t.get_triad(), t.get_metric(),
255  t.get_poids()) ;
256  if (t0 == t1)
257  tmp = manipule(t, metre, 0) ;
258  else
259  tmp = t ;
260 
261  Tenseur trace(contract(tmp, 0, 1)) ;
262 
263  if (t0 == t1) {
264  switch (t0) {
265  case COV : {
266  resu = resu - 1./3.*trace * metre.cov() ;
267  break ;
268  }
269  case CON : {
270  resu = resu - 1./3.*trace * metre.con() ;
271  break ;
272  }
273  default :
274  cout << "Erreur bizarre dans sans_trace!" << endl ;
275  abort() ;
276  break ;
277  }
278  }
279  else {
280  Tenseur_sym delta(*t.get_mp(), 2, mix, *t.get_triad(),
281  t.get_metric(), t.get_poids()) ;
282  delta.set_etat_qcq() ;
283  for (int i=0; i<3; i++)
284  for (int j=i; j<3; j++)
285  delta.set(i,j) = (i==j ? 1 : 0) ;
286  resu = resu - trace/3. * delta ;
287  }
288  resu.set_std_base() ;
289  return resu ;
290 }
291 }
Lorene::Tenseur::etat
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:321
Lorene::Tenseur::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tenseur::get_type_indice
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:726
Lorene::Tenseur::get_metric
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition: tenseur.h:745
Lorene::Tenseur::valence
int valence
Valence.
Definition: tenseur.h:307
Lorene::operator*
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Definition: base_val_mult.C:102
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
Lorene::Itbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: itbl.C:261
Lorene::Tenseur::poids
double poids
For tensor densities: the weight.
Definition: tenseur.h:323
Lorene::Tenseur::type_indice
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
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
Lorene::Itbl
Basic integer array class.
Definition: itbl.h:122
Lorene::Base_vect::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
Lorene::manipule
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Definition: tenseur_operateur.C:509
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
Lorene::Tenseur::get_valence
int get_valence() const
Returns the valence.
Definition: tenseur.h:710
Lorene::lie_derive
Tenseur lie_derive(const Tenseur &t, const Tenseur &x, const Metrique *=0x0)
Lie Derivative of t with respect to x .
Definition: tenseur_operateur.C:816
Lorene::Tenseur_sym
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition: tenseur.h:1253
Lorene::Tenseur::derive_cov
const Tenseur & derive_cov(const Metrique &met) const
Returns the covariant derivative of *this , with respect to met .
Definition: tenseur.C:1554
Lorene::sans_trace
Tenseur sans_trace(const Tenseur &tens, const Metrique &metre)
Computes the traceless part of a Tenseur of valence 2.
Definition: tenseur_operateur.C:897
Lorene::Tenseur::get_mp
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:699
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Lorene::Tenseur_sym::donne_indices
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_sym.C:252
Lorene::Tenseur::mp
const Map *const mp
Reference mapping.
Definition: tenseur.h:306
Lorene::Tenseur::get_poids
double get_poids() const
Returns the weight.
Definition: tenseur.h:738
Lorene::Tenseur::get_n_comp
int get_n_comp() const
Returns the number of components.
Definition: tenseur.h:713
Lorene::Tenseur::metric
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:325
Lorene::Tenseur::n_comp
int n_comp
Number of components, depending on the symmetry.
Definition: tenseur.h:320
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279