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