LORENE
tenseur_arithm.C
1 /*
2  * Arithmetics functions for the Tenseur class.
3  *
4  * These functions are not member functions of the Tenseur class.
5  *
6  * (see file tenseur.h for documentation).
7  *
8  */
9 
10 /*
11  * Copyright (c) 1999-2001 Philippe Grandclement
12  * Copyright (c) 2000-2001 Eric Gourgoulhon
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation; either version 2 of the License, or
19  * (at your option) any later version.
20  *
21  * LORENE is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with LORENE; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  *
30  */
31 
32 
33 
34 /*
35  * $Id: tenseur_arithm.C,v 1.10 2016/12/05 16:18:16 j_novak Exp $
36  * $Log: tenseur_arithm.C,v $
37  * Revision 1.10 2016/12/05 16:18:16 j_novak
38  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
39  *
40  * Revision 1.9 2014/10/13 08:53:42 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.8 2014/10/06 15:13:18 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.7 2003/10/13 10:33:43 f_limousin
47  * *** empty log message ***
48  *
49  * Revision 1.6 2003/06/20 14:52:21 f_limousin
50  * Put an assert on "poids" into comments
51  *
52  * Revision 1.5 2003/03/03 19:40:52 f_limousin
53  * Suppression of an assert on a metric associated with a tensor.
54  *
55  * Revision 1.4 2002/10/16 14:37:14 j_novak
56  * Reorganization of #include instructions of standard C++, in order to
57  * use experimental version 3 of gcc.
58  *
59  * Revision 1.3 2002/09/06 14:49:25 j_novak
60  * Added method lie_derive for Tenseur and Tenseur_sym.
61  * Corrected various errors for derive_cov and arithmetic.
62  *
63  * Revision 1.2 2002/08/07 16:14:11 j_novak
64  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
65  *
66  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
67  * LORENE
68  *
69  * Revision 2.5 2000/02/09 19:30:22 eric
70  * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
71  * argument des constructeurs.
72  *
73  * Revision 2.4 2000/02/08 19:04:58 eric
74  * Les fonctions arithmetiques ne sont plus amies.
75  * Les fonctions exp, log et sqrt se trouvent desormais dans le fichier
76  * tenseur_math.C
77  * Modif de diverses operations (notament division avec double)
78  * Ajout de nouvelles operations (par ex. Tenseur + double, etc...)
79  *
80  * Revision 2.3 2000/02/01 15:40:29 eric
81  * Ajout de la fonction sqrt
82  *
83  * Revision 2.2 2000/01/11 11:14:15 eric
84  * Changement de nom pour la base vectorielle : base --> triad
85  *
86  * Revision 2.1 2000/01/10 17:25:34 eric
87  * Gestion des bases vectorielles (triades de decomposition).
88  *
89  * Revision 2.0 1999/12/02 17:18:47 phil
90  * *** empty log message ***
91  *
92  *
93  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_arithm.C,v 1.10 2016/12/05 16:18:16 j_novak Exp $
94  *
95  */
96 
97 // Headers C
98 #include <cstdlib>
99 #include <cassert>
100 #include <cmath>
101 
102 // Headers Lorene
103 #include "tenseur.h"
104 
105  //********************//
106  // OPERATEURS UNAIRES //
107  //********************//
108 
109 namespace Lorene {
111 
112  return t ;
113 
114 }
115 
117 
118  assert (t.get_etat() != ETATNONDEF) ;
119  if (t.get_etat() == ETATZERO)
120  return t ;
121  else {
122  Tenseur res(*(t.get_mp()), t.get_valence(), t.get_type_indice(),
123  t.get_triad(), t.get_metric(), t.get_poids()) ;
124 
125 
126  res.set_etat_qcq();
127 
128  for (int i=0 ; i<res.get_n_comp() ; i++) {
129  Itbl indices (res.donne_indices(i)) ;
130  res.set(indices) = -t(indices) ;
131  }
132  return res ;
133  }
134 }
135 
136  //**********//
137  // ADDITION //
138  //**********//
139 
140 Tenseur operator+(const Tenseur & t1, const Tenseur & t2) {
141 
142  assert ((t1.get_etat() != ETATNONDEF) && (t2.get_etat() != ETATNONDEF)) ;
143  assert (t1.get_valence() == t2.get_valence()) ;
144  assert (t1.get_mp() == t2.get_mp()) ;
145  if (t1.get_valence() != 0) {
146  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
147  }
148 
149  for (int i=0 ; i<t1.get_valence() ; i++)
150  assert(t1.get_type_indice(i) == t2.get_type_indice(i)) ;
151  // assert (t1.get_metric() == t2.get_metric()) ;
152  //assert (fabs(t1.get_poids() - t2.get_poids())<1.e-10) ;
153 
154  if (t1.get_etat() == ETATZERO)
155  return t2 ;
156  else if (t2.get_etat() == ETATZERO)
157  return t1 ;
158  else {
159  Tenseur res(*(t1.get_mp()), t1.get_valence(), t1.get_type_indice(),
160  t1.get_triad(), t1.get_metric(), t1.get_poids() ) ;
161 
162  res.set_etat_qcq() ;
163  for (int i=0 ; i<res.get_n_comp() ; i++) {
164  Itbl indices (res.donne_indices(i)) ;
165  res.set(indices) = t1(indices) + t2(indices) ;
166  }
167  return res ;
168  }
169 }
170 
171 
172 Tenseur operator+(const Tenseur & t1, double x) {
173 
174  assert (t1.get_etat() != ETATNONDEF) ;
175  assert (t1.get_valence() == 0) ;
176 
177  if (x == double(0)) {
178  return t1 ;
179  }
180 
181  Tenseur res( *(t1.get_mp()), t1.get_metric(), t1.get_poids() ) ;
182 
183  res.set_etat_qcq() ;
184 
185  res.set() = t1() + x ; // Cmp + double
186 
187  return res ;
188 
189 }
190 
191 
192 Tenseur operator+(double x, const Tenseur & t2) {
193 
194  return t2 + x ;
195 
196 }
197 
198 Tenseur operator+(const Tenseur & t1, int m) {
199 
200  return t1 + double(m) ;
201 
202 }
203 
204 
205 Tenseur operator+(int m, const Tenseur & t2) {
206 
207  return t2 + double(m) ;
208 
209 }
210 
211 
212 
213  //**************//
214  // SOUSTRACTION //
215  //**************//
216 
217 Tenseur operator-(const Tenseur & t1, const Tenseur & t2) {
218 
219  return (t1 + (-t2)) ;
220 
221 }
222 
223 
224 Tenseur operator-(const Tenseur & t1, double x) {
225 
226  assert (t1.get_etat() != ETATNONDEF) ;
227  assert (t1.get_valence() == 0) ;
228 
229  if (x == double(0)) {
230  return t1 ;
231  }
232 
233  Tenseur res( *(t1.get_mp()), t1.get_metric(), t1.get_poids() ) ;
234 
235  res.set_etat_qcq() ;
236 
237  res.set() = t1() - x ; // Cmp - double
238 
239  return res ;
240 
241 }
242 
243 
244 Tenseur operator-(double x, const Tenseur & t2) {
245 
246  return - (t2 - x) ;
247 
248 }
249 
250 
251 Tenseur operator-(const Tenseur & t1, int m) {
252 
253  return t1 - double(m) ;
254 
255 }
256 
257 
258 Tenseur operator-(int m, const Tenseur & t2) {
259 
260  return - (t2 - double(m)) ;
261 
262 }
263 
264 
265 
266  //****************//
267  // MULTIPLICATION //
268  //****************//
269 
270 
271 
272 Tenseur operator*(double x, const Tenseur& t) {
273 
274  assert (t.get_etat() != ETATNONDEF) ;
275  if ( (t.get_etat() == ETATZERO) || (x == double(1)) )
276  return t ;
277  else {
278  Tenseur res(*(t.get_mp()), t.get_valence(), t.get_type_indice(),
279  t.get_triad(), t.get_metric(), t.get_poids() ) ;
280 
281  if ( x == double(0) )
282  res.set_etat_zero() ;
283  else {
284  res.set_etat_qcq() ;
285  for (int i=0 ; i<res.get_n_comp() ; i++) {
286  Itbl indices (res.donne_indices(i)) ;
287  res.set(indices) = x*t(indices) ;
288  }
289  }
290  return res ;
291  }
292 }
293 
294 
295 Tenseur operator* (const Tenseur& t, double x) {
296  return x * t ;
297 }
298 
299 Tenseur operator*(int m, const Tenseur& t) {
300  return double(m) * t ;
301 }
302 
303 
304 Tenseur operator* (const Tenseur& t, int m) {
305  return double(m) * t ;
306 }
307 
308 
309  //**********//
310  // DIVISION //
311  //**********//
312 
313 Tenseur operator/ (const Tenseur& t1, const Tenseur& t2) {
314 
315  // Protections
316  assert(t1.get_etat() != ETATNONDEF) ;
317  assert(t2.get_etat() != ETATNONDEF) ;
318  assert(t2.get_valence() == 0) ; // t2 doit etre un scalaire !
319  assert(t1.get_mp() == t2.get_mp()) ;
320 
321  double poids_res = t1.get_poids() - t2.get_poids() ;
322  poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
323  const Metrique* met_res = 0x0 ;
324  if (poids_res != 0.) {
325  // assert((t1.get_metric() != 0x0) || (t2.get_metric() != 0x0)) ;
326  if (t1.get_metric() != 0x0) met_res = t1.get_metric() ;
327  else met_res = t2.get_metric() ;
328  }
329  // Cas particuliers
330  if (t2.get_etat() == ETATZERO) {
331  cout << "Division by 0 in Tenseur / Tenseur !" << endl ;
332  abort() ;
333  }
334  if (t1.get_etat() == ETATZERO) {
335  Tenseur resu(t1) ;
336  resu.set_poids(poids_res) ;
337  resu.set_metric(*met_res) ;
338  return resu ;
339  }
340 
341  // Cas general
342 
343  assert(t1.get_etat() == ETATQCQ) ; // sinon...
344  assert(t2.get_etat() == ETATQCQ) ; // sinon...
345 
346  Tenseur res(*(t1.get_mp()), t1.get_valence(), t1.get_type_indice(),
347  t1.get_triad(), met_res, poids_res) ;
348 
349  res.set_etat_qcq() ;
350  for (int i=0 ; i<res.get_n_comp() ; i++) {
351  Itbl indices (res.donne_indices(i)) ;
352  res.set(indices) = t1(indices) / t2() ; // Cmp / Cmp
353  }
354  return res ;
355 
356 }
357 
358 
359 Tenseur operator/ (const Tenseur& t, double x) {
360 
361  assert (t.get_etat() != ETATNONDEF) ;
362 
363  if ( x == double(0) ) {
364  cout << "Division by 0 in Tenseur / double !" << endl ;
365  abort() ;
366  }
367 
368  if ( (t.get_etat() == ETATZERO) || (x == double(1)) )
369  return t ;
370  else {
371  Tenseur res(*(t.get_mp()), t.get_valence(), t.get_type_indice(),
372  t.get_triad(), t.get_metric(), t.get_poids()) ;
373 
374  res.set_etat_qcq() ;
375  for (int i=0 ; i<res.get_n_comp() ; i++) {
376  Itbl indices (res.donne_indices(i)) ;
377  res.set(indices) = t(indices) / x ; // Cmp / double
378  }
379  return res ;
380  }
381 
382 }
383 
384 
385 
386 
387 Tenseur operator/ (double x, const Tenseur& t) {
388 
389  if (t.get_etat() == ETATZERO) {
390  cout << "Division by 0 in double / Tenseur !" << endl ;
391  abort() ;
392  }
393 
394  assert (t.get_etat() == ETATQCQ) ;
395  assert(t.get_valence() == 0) ; // Utilisable que sur scalaire !
396 
397  Tenseur res( *(t.get_mp()), t.get_metric(), -t.get_poids() ) ;
398  res.set_etat_qcq() ;
399  res.set() = x / t() ; // double / Cmp
400  return res ;
401 }
402 
403 
404 Tenseur operator/ (const Tenseur& t, int m) {
405 
406  return t / double(m) ;
407 }
408 
409 
410 Tenseur operator/ (int m, const Tenseur& t) {
411 
412  return double(m) / t ;
413 }
414 
415 
416 
417 
418 
419 }
double get_poids() const
Returns the weight.
Definition: tenseur.h:741
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:729
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
void set_poids(double weight)
Sets the weight for a tensor density.
Definition: tenseur.C:696
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 .
Basic integer array class.
Definition: itbl.h:122
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition: cmp_arithm.C:460
int get_valence() const
Returns the valence.
Definition: tenseur.h:713
void set_metric(const Metrique &met)
Sets the pointer on the metric for a tensor density.
Definition: tenseur.C:701
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:702
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition: tenseur.h:748
Cmp operator+(const Cmp &)
Definition: cmp_arithm.C:107
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
Cmp operator-(const Cmp &)
- Cmp
Definition: cmp_arithm.C:111
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304