LORENE
tensor_arithm.C
1 /*
2  * Arithmetics functions for the Tensor class.
3  *
4  * These functions are not member functions of the Tensor class.
5  *
6  * (see file tensor.h for documentation).
7  *
8  */
9 
10 /*
11  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
12  * Copyright (c) 1999-2001 Philippe Grandclement (Tenseur version)
13  * Copyright (c) 2000-2001 Eric Gourgoulhon (Tenseur version)
14  *
15  * This file is part of LORENE.
16  *
17  * LORENE is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation; either version 2 of the License, or
20  * (at your option) any later version.
21  *
22  * LORENE is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with LORENE; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  *
31  */
32 
33 
34 
35 /*
36  * $Id: tensor_arithm.C,v 1.6 2016/12/05 16:18:17 j_novak Exp $
37  * $Log: tensor_arithm.C,v $
38  * Revision 1.6 2016/12/05 16:18:17 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.5 2014/10/13 08:53:44 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.4 2014/10/06 15:13:20 j_novak
45  * Modified #include directives to use c++ syntax.
46  *
47  * Revision 1.3 2004/01/08 09:24:11 e_gourgoulhon
48  * Added arithmetics common to Scalar and Tensor.
49  * Corrected treatment ETATUN in Tensor / Scalar.
50  *
51  * Revision 1.2 2003/10/01 13:04:44 e_gourgoulhon
52  * The method Tensor::get_mp() returns now a reference (and not
53  * a pointer) onto a mapping.
54  *
55  * Revision 1.1 2003/09/26 14:33:53 j_novak
56  * Arithmetic functions for the class Tensor
57  *
58  *
59  *
60  * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_arithm.C,v 1.6 2016/12/05 16:18:17 j_novak Exp $
61  *
62  */
63 
64 // Headers C
65 #include <cstdlib>
66 #include <cassert>
67 #include <cmath>
68 
69 // Headers Lorene
70 #include "tensor.h"
71 
72  //********************//
73  // OPERATEURS UNAIRES //
74  //********************//
75 
76 namespace Lorene {
77 Tensor operator+(const Tensor & t) {
78 
79  return t ;
80 
81 }
82 
83 Tensor operator-(const Tensor & t) {
84 
85  Tensor res(t.get_mp(), t.get_valence(), t.get_index_type(),
86  t.get_triad()) ;
87 
88  for (int i=0 ; i<res.get_n_comp() ; i++) {
89  Itbl ind (res.indices(i)) ;
90  res.set(ind) = -t(ind) ;
91  }
92  return res ;
93 
94 }
95 
96  //**********//
97  // ADDITION //
98  //**********//
99 
100 Tensor operator+(const Tensor & t1, const Tensor & t2) {
101 
102  assert (t1.get_valence() == t2.get_valence()) ;
103  assert (t1.get_mp() == t2.get_mp()) ;
104  if (t1.get_valence() != 0) {
105  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
106  }
107 
108  for (int i=0 ; i<t1.get_valence() ; i++)
109  assert(t1.get_index_type(i) == t2.get_index_type(i)) ;
110 
111  Tensor res(t1.get_mp(), t1.get_valence(), t1.get_index_type(),
112  t1.get_triad()) ;
113 
114  for (int i=0 ; i<res.get_n_comp() ; i++) {
115  Itbl ind (res.indices(i)) ;
116  res.set(ind) = t1(ind) + t2(ind) ;
117  }
118  return res ;
119 
120 }
121 
122 
123 Scalar operator+(const Tensor& t1, const Scalar& t2) {
124 
125  assert (t1.get_valence() == 0) ;
126  assert (t1.get_mp() == t2.get_mp()) ;
127 
128  return *(t1.cmp[0]) + t2 ;
129 
130 }
131 
132 Scalar operator+(const Scalar& t1, const Tensor& t2) {
133 
134  assert (t2.get_valence() == 0) ;
135  assert (t1.get_mp() == t2.get_mp()) ;
136 
137  return t1 + *(t2.cmp[0]) ;
138 
139 }
140 
141  //**************//
142  // SOUSTRACTION //
143  //**************//
144 
145 Tensor operator-(const Tensor & t1, const Tensor & t2) {
146 
147  assert (t1.get_valence() == t2.get_valence()) ;
148  assert (t1.get_mp() == t2.get_mp()) ;
149  if (t1.get_valence() != 0) {
150  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
151  }
152 
153  for (int i=0 ; i<t1.get_valence() ; i++)
154  assert(t1.get_index_type(i) == t2.get_index_type(i)) ;
155 
156  Tensor res(t1.get_mp(), t1.get_valence(), t1.get_index_type(),
157  t1.get_triad()) ;
158 
159  for (int i=0 ; i<res.get_n_comp() ; i++) {
160  Itbl ind (res.indices(i)) ;
161  res.set(ind) = t1(ind) - t2(ind) ;
162  }
163  return res ;
164 
165 }
166 
167 Scalar operator-(const Tensor& t1, const Scalar& t2) {
168 
169  assert (t1.get_valence() == 0) ;
170  assert (t1.get_mp() == t2.get_mp()) ;
171 
172  return *(t1.cmp[0]) - t2 ;
173 
174 }
175 
176 Scalar operator-(const Scalar& t1, const Tensor& t2) {
177 
178  assert (t2.get_valence() == 0) ;
179  assert (t1.get_mp() == t2.get_mp()) ;
180 
181  return t1 - *(t2.cmp[0]) ;
182 
183 }
184 
185 
186 
187  //****************//
188  // MULTIPLICATION //
189  //****************//
190 
191 Tensor operator*(const Scalar& t1, const Tensor& t2) {
192 
193  assert (&(t1.get_mp()) == &(t2.get_mp())) ;
194 
195  if (t1.get_etat() == ETATUN) return t2 ;
196 
197  Tensor res(t2.get_mp(), t2.get_valence(), t2.get_index_type(),
198  t2.get_triad()) ;
199 
200  for (int ic=0 ; ic<res.get_n_comp() ; ic++) {
201  Itbl ind = res.indices(ic) ;
202  res.set(ind) = t1 * t2(ind) ;
203  }
204 
205  return res ;
206 }
207 
208 
209 Tensor operator*(const Tensor& t2, const Scalar& t1) {
210 
211  return t1*t2 ;
212 }
213 
214 
215 
216 Tensor operator*(double x, const Tensor& t) {
217 
218  Tensor res(t.get_mp(), t.get_valence(), t.get_index_type(),
219  t.get_triad()) ;
220 
221  for (int i=0 ; i<res.get_n_comp() ; i++) {
222  Itbl ind (res.indices(i)) ;
223  res.set(ind) = x*t(ind) ;
224  }
225 
226  return res ;
227 
228 }
229 
230 
231 Tensor operator* (const Tensor& t, double x) {
232  return x * t ;
233 }
234 
235 Tensor operator*(int m, const Tensor& t) {
236  return double(m) * t ;
237 }
238 
239 
240 Tensor operator* (const Tensor& t, int m) {
241  return double(m) * t ;
242 }
243 
244 
245  //**********//
246  // DIVISION //
247  //**********//
248 
249 Tensor operator/(const Tensor& t1, const Scalar& s2) {
250 
251  // Protections
252  assert(s2.get_etat() != ETATNONDEF) ;
253  assert(t1.get_mp() == s2.get_mp()) ;
254 
255  // Cas particuliers
256  if (s2.get_etat() == ETATZERO) {
257  cout << "Division by 0 in Tensor / Scalar !" << endl ;
258  abort() ;
259  }
260 
261  if (s2.get_etat() == ETATUN) return t1 ;
262 
263  Tensor res(t1.get_mp(), t1.get_valence(), t1.get_index_type(),
264  t1.get_triad()) ;
265 
266  for (int i=0 ; i<res.get_n_comp() ; i++) {
267  Itbl ind (res.indices(i)) ;
268  res.set(ind) = t1(ind) / s2 ; // Scalar / Scalar
269  }
270  return res ;
271 
272 }
273 
274 
275 Tensor operator/ (const Tensor& t, double x) {
276 
277  if ( x == double(0) ) {
278  cout << "Division by 0 in Tensor / double !" << endl ;
279  abort() ;
280  }
281 
282  if (x == double(1))
283  return t ;
284  else {
285  Tensor res(t.get_mp(), t.get_valence(), t.get_index_type(),
286  t.get_triad()) ;
287 
288  for (int i=0 ; i<res.get_n_comp() ; i++) {
289  Itbl ind (res.indices(i)) ;
290  res.set(ind) = t(ind) / x ; // Scalar / double
291  }
292  return res ;
293  }
294 
295 }
296 
297 Tensor operator/ (const Tensor& t, int m) {
298 
299  return t / double(m) ;
300 }
301 
302 
303 
304 
305 
306 
307 
308 }
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Lorene prototypes.
Definition: app_hor.h:67
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition: cmp_arithm.C:460
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Cmp operator+(const Cmp &)
Definition: cmp_arithm.C:107
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:899
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
Tensor handling.
Definition: tensor.h:294
int get_valence() const
Returns the valence.
Definition: tensor.h:882
Cmp operator-(const Cmp &)
- Cmp
Definition: cmp_arithm.C:111
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874