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
32char 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
106namespace 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
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
169Tenseur 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
189Tenseur operator+(double x, const Tenseur & t2) {
190
191 return t2 + x ;
192
193}
194
195Tenseur operator+(const Tenseur & t1, int m) {
196
197 return t1 + double(m) ;
198
199}
200
201
202Tenseur operator+(int m, const Tenseur & t2) {
203
204 return t2 + double(m) ;
205
206}
207
208
209
210 //**************//
211 // SOUSTRACTION //
212 //**************//
213
215
216 return (t1 + (-t2)) ;
217
218}
219
220
221Tenseur 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
241Tenseur operator-(double x, const Tenseur & t2) {
242
243 return - (t2 - x) ;
244
245}
246
247
248Tenseur operator-(const Tenseur & t1, int m) {
249
250 return t1 - double(m) ;
251
252}
253
254
255Tenseur operator-(int m, const Tenseur & t2) {
256
257 return - (t2 - double(m)) ;
258
259}
260
261
262
263 //****************//
264 // MULTIPLICATION //
265 //****************//
266
267
268
269Tenseur 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
292Tenseur operator* (const Tenseur& t, double x) {
293 return x * t ;
294}
295
296Tenseur operator*(int m, const Tenseur& t) {
297 return double(m) * t ;
298}
299
300
301Tenseur operator* (const Tenseur& t, int m) {
302 return double(m) * t ;
303}
304
305
306 //**********//
307 // DIVISION //
308 //**********//
309
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
356Tenseur 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
384Tenseur 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
401Tenseur operator/ (const Tenseur& t, int m) {
402
403 return t / double(m) ;
404}
405
406
407Tenseur operator/ (int m, const Tenseur& t) {
408
409 return double(m) / t ;
410}
411
412
413
414
415
416}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:699
double get_poids() const
Returns the weight.
Definition tenseur.h:738
int get_valence() const
Returns the valence.
Definition tenseur.h:710
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition tenseur.h:745
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Cmp operator-(const Cmp &)
- Cmp
Definition cmp_arithm.C:108
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition cmp_arithm.C:457
Cmp operator+(const Cmp &)
Definition cmp_arithm.C:104
Lorene prototypes.
Definition app_hor.h:64