LORENE
connection_fcart.C
1 /*
2  * Methods of class Connection_fcart.
3  *
4  * (see file connection.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: connection_fcart.C,v 1.15 2016/12/05 16:17:50 j_novak Exp $
32  * $Log: connection_fcart.C,v $
33  * Revision 1.15 2016/12/05 16:17:50 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.14 2014/10/13 08:52:50 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.13 2014/10/06 15:13:04 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.12 2004/01/28 13:25:40 j_novak
43  * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
44  * In the div/mult _r_dzpuis, there is no more default value.
45  *
46  * Revision 1.11 2004/01/04 21:00:50 e_gourgoulhon
47  * Better handling of tensor symmetries in methods p_derive_cov() and
48  * p_divergence() (thanks to the new class Tensor_sym).
49  *
50  * Revision 1.10 2004/01/01 11:24:04 e_gourgoulhon
51  * Full reorganization of method p_derive_cov: the main loop is now
52  * on the indices of the *output* tensor (to take into account
53  * symmetries in the input and output tensors).
54  *
55  * Revision 1.9 2003/12/27 14:59:52 e_gourgoulhon
56  * -- Method derive_cov() suppressed.
57  * -- Change of the position of the derivation index from the first one
58  * to the last one in methods p_derive_cov() and p_divergence().
59  *
60  * Revision 1.8 2003/10/17 13:46:15 j_novak
61  * The argument is now between 1 and 3 (instead of 0->2)
62  *
63  * Revision 1.7 2003/10/16 21:37:08 e_gourgoulhon
64  * Corrected deriv index in divergence.
65  *
66  * Revision 1.6 2003/10/16 15:26:03 e_gourgoulhon
67  * Suppressed unsued variable
68  *
69  * Revision 1.5 2003/10/16 14:21:36 j_novak
70  * The calculation of the divergence of a Tensor is now possible.
71  *
72  * Revision 1.4 2003/10/11 16:45:43 e_gourgoulhon
73  * Suppressed the call to Itbl::set_etat_qcq() after
74  * the construction of the Itbl's.
75  *
76  * Revision 1.3 2003/10/11 14:39:50 e_gourgoulhon
77  * Suppressed declaration of unusued arguments in some methods.
78  *
79  * Revision 1.2 2003/10/06 13:58:46 j_novak
80  * The memory management has been improved.
81  * Implementation of the covariant derivative with respect to the exact Tensor
82  * type.
83  *
84  * Revision 1.1 2003/10/03 14:11:48 e_gourgoulhon
85  * Methods of class Connection_fcart.
86  *
87  *
88  *
89  * $Header: /cvsroot/Lorene/C++/Source/Connection/connection_fcart.C,v 1.15 2016/12/05 16:17:50 j_novak Exp $
90  *
91  */
92 
93 // C++ headers
94 #include "headcpp.h"
95 
96 // C headers
97 #include <cstdlib>
98 
99 // Lorene headers
100 #include "connection.h"
101 
102 
103  //------------------------------//
104  // Constructors //
105  //------------------------------//
106 
107 
108 
109 // Contructor from a Cartesian flat-metric-orthonormal basis
110 
111 namespace Lorene {
113  : Connection_flat(mpi, bi) {
114 
115 }
116 
117 // Copy constructor
119  : Connection_flat(ci) {
120 
121 }
122 
123 
124  //----------------------------//
125  // Destructor //
126  //----------------------------//
127 
128 
130 
131 }
132 
133 
134  //--------------------------------//
135  // Mutators / assignment //
136  //--------------------------------//
137 
138 
140 
141  cout << "Connection_fcart::operator= : not implemented yet !" << endl ;
142  abort() ;
143 
144 }
145 
146 
147 
148  //-----------------------------//
149  // Computational methods //
150  //-----------------------------//
151 
152 // Covariant derivative, returning a pointer.
153 //-------------------------------------------
154 
156 
157  // Notations: suffix 0 in name <=> input tensor
158  // suffix 1 in name <=> output tensor
159 
160  int valence0 = uu.get_valence() ;
161  int valence1 = valence0 + 1 ;
162  int valence1m1 = valence1 - 1 ; // same as valence0, but introduced for
163  // the sake of clarity
164 
165  // Protections
166  // -----------
167  if (valence0 >= 1) {
168  assert(uu.get_triad() == triad) ;
169  }
170 
171  // Creation of the result (pointer)
172  // --------------------------------
173  Tensor* resu ;
174 
175  // If uu is a Scalar, the result is a vector
176  if (valence0 == 0)
177  resu = new Vector(*mp, COV, triad) ;
178  else {
179 
180  // Type of indices of the result :
181  Itbl tipe(valence1) ;
182  const Itbl& tipeuu = uu.get_index_type() ;
183  for (int id = 0; id<valence0; id++) {
184  tipe.set(id) = tipeuu(id) ; // First indices = same as uu
185  }
186  tipe.set(valence1m1) = COV ; // last index is the derivation index
187 
188  // if uu is a Tensor_sym, the result is also a Tensor_sym:
189  const Tensor* puu = &uu ;
190  const Tensor_sym* puus = dynamic_cast<const Tensor_sym*>(puu) ;
191  if ( puus != 0x0 ) { // the input tensor is symmetric
192  resu = new Tensor_sym(*mp, valence1, tipe, *triad,
193  puus->sym_index1(), puus->sym_index2()) ;
194  }
195  else {
196  resu = new Tensor(*mp, valence1, tipe, *triad) ; // no symmetry
197  }
198 
199  }
200 
201  int ncomp1 = resu->get_n_comp() ;
202 
203  Itbl ind1(valence1) ; // working Itbl to store the indices of resu
204  Itbl ind0(valence0) ; // working Itbl to store the indices of uu
205 
206  // Loop on all the components of the output tensor
207  // -----------------------------------------------
208  for (int ic=0; ic<ncomp1; ic++) {
209 
210  // indices corresponding to the component no. ic in the output tensor
211  ind1 = resu->indices(ic) ;
212 
213  // Component no. ic:
214  Scalar& cresu = resu->set(ind1) ;
215 
216  // Indices of the input tensor
217  for (int id = 0; id < valence0; id++) {
218  ind0.set(id) = ind1(id) ;
219  }
220 
221  // Value of last index (derivation index)
222  int k = ind1(valence1m1) ;
223 
224  // Partial derivation with respect to x^k:
225 
226  cresu = (uu(ind0)).deriv(k) ;
227 
228  }
229 
230  // C'est fini !
231  // -----------
232  return resu ;
233 
234 }
235 
236 
237 
238 // Divergence, returning a pointer.
239 //---------------------------------
240 
242 
243  // Notations: suffix 0 in name <=> input tensor
244  // suffix 1 in name <=> output tensor
245 
246  int valence0 = uu.get_valence() ;
247  int valence1 = valence0 - 1 ;
248 
249  // Protections
250  // -----------
251  assert (valence0 >= 1) ;
252  assert (uu.get_triad() == triad) ;
253 
254  // Last index must be contravariant:
255  assert (uu.get_index_type(valence0-1) == CON) ;
256 
257 
258  // Creation of the pointer on the result tensor
259  // --------------------------------------------
260  Tensor* resu ;
261 
262  if (valence0 == 1) // if u is a Vector, the result is a Scalar
263  resu = new Scalar(*mp) ;
264  else {
265 
266  // Type of indices of the result :
267  Itbl tipe(valence1) ;
268  const Itbl& tipeuu = uu.get_index_type() ;
269  for (int id = 0; id<valence1; id++) {
270  tipe.set(id) = tipeuu(id) ; // type of remaining indices =
271  } // same as uu indices
272 
273  if (valence0 == 2) { // if u is a rank 2 tensor, the result is a Vector
274  resu = new Vector(*mp, tipe(0), *triad) ;
275  }
276  else {
277  // if uu is a Tensor_sym, the result might be also a Tensor_sym:
278  const Tensor* puu = &uu ;
279  const Tensor_sym* puus = dynamic_cast<const Tensor_sym*>(puu) ;
280  if ( puus != 0x0 ) { // the input tensor is symmetric
281 
282  if (puus->sym_index2() != valence0 - 1) {
283 
284  // the symmetry is preserved:
285 
286  if (valence1 == 2) {
287  resu = new Sym_tensor(*mp, tipe, *triad) ;
288  }
289  else {
290  resu = new Tensor_sym(*mp, valence1, tipe, *triad,
291  puus->sym_index1(), puus->sym_index2()) ;
292  }
293  }
294  else { // the symmetry is lost:
295 
296  resu = new Tensor(*mp, valence1, tipe, *triad) ;
297  }
298  }
299  else { // no symmetry in the input tensor:
300  resu = new Tensor(*mp, valence1, tipe, *triad) ;
301  }
302  }
303  }
304 
305 
306  int ncomp1 = resu->get_n_comp() ;
307 
308  Itbl ind0(valence0) ; // working Itbl to store the indices of uu
309 
310  Itbl ind1(valence1) ; // working Itbl to store the indices of resu
311 
312  // Loop on all the components of the output tensor
313  for (int ic=0; ic<ncomp1; ic++) {
314 
315  ind1 = resu->indices(ic) ;
316  Scalar& cresu = resu->set(ind1) ;
317  cresu.set_etat_zero() ;
318 
319  for (int k=1; k<=3; k++) {
320 
321  // indices (ind1,k) in the input tensor
322  for (int id = 0; id < valence1; id++) {
323  ind0.set(id) = ind1(id) ;
324  }
325  ind0.set(valence0-1) = k ;
326 
327  cresu += uu(ind0).deriv(k) ; //Addition of dT^i/dx^i
328  }
329 
330  }
331 
332  // C'est fini !
333  // -----------
334  return resu ;
335 
336 }
337 
338 }
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
const Map *const mp
Reference mapping.
Definition: connection.h:119
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
Class Connection_flat.
Definition: connection.h:354
int sym_index1() const
Number of the first symmetric index (0<= id_sym1 < valence )
Definition: tensor.h:1162
Lorene prototypes.
Definition: app_hor.h:67
virtual Tensor * p_divergence(const Tensor &tens) const
Computes the divergence of a tensor (with respect to the current connection).
virtual Tensor * p_derive_cov(const Tensor &tens) const
Computes the covariant derivative of a tensor (with respect to the current connection).
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.h:885
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Definition: tensor.h:1167
Basic integer array class.
Definition: itbl.h:122
Tensor field of valence 1.
Definition: vector.h:188
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:899
Class Connection_fcart.
Definition: connection.h:546
void operator=(const Connection_fcart &)
Assignment to another Connection_fcart.
Tensor handling.
Definition: tensor.h:294
int get_valence() const
Returns the valence.
Definition: tensor.h:882
Cartesian vectorial bases (triads).
Definition: base_vect.h:201
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Connection_fcart(const Map &, const Base_vect_cart &)
Contructor from a Cartesian flat-metric-orthonormal basis.
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
const Base_vect *const triad
Triad with respect to which the connection coefficients are defined.
Definition: connection.h:124
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition: tensor.C:548
virtual ~Connection_fcart()
destructor