LORENE
sym_tensor_trans.C
1 /*
2  * Methods of class Sym_tensor_trans
3  *
4  * (see file sym_tensor.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 as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: sym_tensor_trans.C,v 1.19 2016/12/05 16:18:17 j_novak Exp $
34  * $Log: sym_tensor_trans.C,v $
35  * Revision 1.19 2016/12/05 16:18:17 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.18 2014/10/13 08:53:43 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.17 2014/10/06 15:13:19 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.16 2008/12/03 10:20:00 j_novak
45  * Modified output.
46  *
47  * Revision 1.15 2006/09/15 08:48:01 j_novak
48  * Suppression of some messages in the optimized version.
49  *
50  * Revision 1.14 2005/01/03 12:37:08 j_novak
51  * Sym_tensor_trans::trace_from_det_one : modified the test on hijtt to
52  * be compatible with older compilers.
53  *
54  * Revision 1.13 2005/01/03 08:36:36 f_limousin
55  * Come back to the previous version.
56  *
57  * Revision 1.12 2005/01/03 08:13:26 f_limousin
58  * The first argument of the function trace_from_det_one(...) is now
59  * a Sym_tensor_trans instead of a Sym_tensor_tt (because of a
60  * compilation error with some compilators).
61  *
62  * Revision 1.11 2004/12/28 14:21:48 j_novak
63  * Added the method Sym_tensor_trans::trace_from_det_one
64  *
65  * Revision 1.10 2004/05/25 15:07:12 f_limousin
66  * Add parameters in argument of the function tt_part for the case
67  * of a Map_et.
68  *
69  * Revision 1.9 2004/03/30 14:01:19 j_novak
70  * Copy constructors and operator= now copy the "derived" members.
71  *
72  * Revision 1.8 2004/03/30 08:01:16 j_novak
73  * Better treatment of dzpuis in mutators.
74  *
75  * Revision 1.7 2004/03/29 16:13:07 j_novak
76  * New methods set_longit_trans and set_tt_trace .
77  *
78  * Revision 1.6 2004/03/03 13:22:14 j_novak
79  * The case where dzpuis = 0 is treated in tt_part().
80  *
81  * Revision 1.5 2004/02/18 18:48:39 e_gourgoulhon
82  * Method trace() renamed the_trace() to avoid any confusion with
83  * the new method Tensor::trace().
84  *
85  * Revision 1.4 2004/02/09 12:57:13 e_gourgoulhon
86  * First implementation of method tt_part().
87  *
88  * Revision 1.3 2004/01/04 20:52:45 e_gourgoulhon
89  * Added assignement (operator=) to a Tensor_sym.
90  *
91  * Revision 1.2 2003/10/28 21:24:52 e_gourgoulhon
92  * Added new methods trace() and tt_part().
93  *
94  * Revision 1.1 2003/10/27 10:50:54 e_gourgoulhon
95  * First version.
96  *
97  *
98  *
99  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans.C,v 1.19 2016/12/05 16:18:17 j_novak Exp $
100  *
101  */
102 
103 // Headers C
104 #include <cstdlib>
105 
106 // Headers Lorene
107 #include "metric.h"
108 #include "param.h"
109 
110  //--------------//
111  // Constructors //
112  //--------------//
113 
114 // Standard constructor
115 // --------------------
116 namespace Lorene {
118  const Metric& met)
119  : Sym_tensor(map, CON, triad_i),
120  met_div(&met) {
121 
122  set_der_0x0() ;
123 
124 }
125 
126 // Copy constructor
127 // ----------------
129  : Sym_tensor(source),
130  met_div(source.met_div) {
131 
132  set_der_0x0() ;
133 
134  if (source.p_trace != 0x0) p_trace = new Scalar( *(source.p_trace) ) ;
135  if (source.p_tt != 0x0) p_tt = new Sym_tensor_tt( *(source.p_tt) ) ;
136 
137 }
138 
139 
140 // Constructor from a file
141 // -----------------------
142 Sym_tensor_trans::Sym_tensor_trans(const Map& mapping, const Base_vect& triad_i,
143  const Metric& met, FILE* fd)
144  : Sym_tensor(mapping, triad_i, fd),
145  met_div(&met) {
146 
147  set_der_0x0() ;
148 }
149 
150  //--------------//
151  // Destructor //
152  //--------------//
153 
155 
156  Sym_tensor_trans::del_deriv() ; // in order not to follow the virtual aspect
157  // of del_deriv()
158 
159 }
160 
161 
162 
163  //-------------------//
164  // Memory managment //
165  //-------------------//
166 
168 
169  if (p_trace != 0x0) delete p_trace ;
170  if (p_tt != 0x0) delete p_tt ;
171 
172  set_der_0x0() ;
174 
175 }
176 
178 
179  p_trace = 0x0 ;
180  p_tt = 0x0 ;
181 
182 }
183 
184 
185  //--------------//
186  // Assignment //
187  //--------------//
188 
190 
191  // Assignment of quantities common to all the derived classes of Sym_tensor
192  Sym_tensor::operator=(source) ;
193 
194  // Assignment of proper quantities of class Sym_tensor_trans
195  assert(met_div == source.met_div) ;
196 
197  del_deriv() ;
198 
199  if (source.p_trace != 0x0) p_trace = new Scalar( *(source.p_trace) ) ;
200  if (source.p_tt != 0x0) p_tt = new Sym_tensor_tt( *(source.p_tt) ) ;
201 
202 }
203 
204 
206 
207  // Assignment of quantities common to all the derived classes of Sym_tensor
208  Sym_tensor::operator=(source) ;
209 
210  // The metric which was set by the constructor is kept
211 
212  del_deriv() ;
213 }
214 
215 
217 
218  // Assignment of quantities common to all the derived classes of Sym_tensor
219  Sym_tensor::operator=(source) ;
220 
221  // The metric which was set by the constructor is kept
222 
223  del_deriv() ;
224 }
225 
226 
227 
228 void Sym_tensor_trans::operator=(const Tensor& source) {
229 
230  // Assignment of quantities common to all the derived classes of Sym_tensor
231  Sym_tensor::operator=(source) ;
232 
233  // The metric which was set by the constructor is kept
234 
235  del_deriv() ;
236 }
237 
239  const Scalar& htrace, Param* par ) {
240 
241  assert (met_div == &htt.get_met_div() ) ;
242 
243  assert (htrace.check_dzpuis(4)) ;
244 
245  Scalar pot (*mp) ;
246  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
247  pot = htrace.poisson() ;
248  }
249  else {
250  pot = 0. ;
251  htrace.poisson(*par, pot) ;
252  }
253 
254  Sym_tensor tmp = (pot.derive_con(*met_div)).derive_con(*met_div) ;
255  tmp.dec_dzpuis() ;
256 
257  *this = htrace * met_div->con() ;
258  dec_dzpuis(2) ; // this has now dzpuis = 2
259  *this = htt + 0.5 * ( *this - tmp ) ;
260 
261  del_deriv() ;
262 
263  p_trace = new Scalar( htrace ) ;
264  p_tt = new Sym_tensor_tt( htt ) ;
265 
266 }
267 
268 
269  //-----------------------------//
270  // Computational methods //
271  //-----------------------------//
272 
274 
275  if (p_trace == 0x0) { // a new computation is necessary
276 
277  assert( (type_indice(0)==CON) && (type_indice(1)==CON) ) ;
278  p_trace = new Scalar( trace(*met_div) ) ;
279 
280  }
281 
282  return *p_trace ;
283 
284 }
285 
286 
288 
289  if (p_tt == 0x0) { // a new computation is necessary
290 
291  int dzp = the_trace().get_dzpuis() ;
292 
293  assert((dzp == 0) || (dzp == 4)) ;
294 
295  p_tt = new Sym_tensor_tt(*mp, *triad, *met_div) ;
296 
297  Scalar pot (*mp) ;
298  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
299  pot = the_trace().poisson() ;
300  }
301  else {
302  pot = 0. ;
303  the_trace().poisson(*par, pot) ;
304  }
305 
306  Sym_tensor tmp = (pot.derive_con(*met_div)).derive_con(*met_div) ;
307  (dzp == 4) ? tmp.inc_dzpuis() : tmp.dec_dzpuis(3) ; //## to be improved ?
308 
309  *p_tt = *this - 0.5 * ( the_trace() * met_div->con() - tmp ) ;
310 
311  }
312 
313  return *p_tt ;
314 
315 }
316 
317 
319  double precis, int it_max) {
320 
321 #ifndef NDEBUG
322  const Sym_tensor_trans* ptmp =
323  dynamic_cast<const Sym_tensor_trans*>(&hijtt) ;
324  assert (ptmp != 0x0) ;
325  assert (ptmp != this) ;
326  for (int i=0; i<hijtt.get_n_comp(); i++)
327  assert(hijtt.cmp[i]->check_dzpuis(2)) ;
328 #endif
329  assert( (precis > 0.) && (it_max > 0) ) ;
330  assert (met_div == &hijtt.get_met_div() ) ;
331 
332  Sym_tensor_trans& hij = *this ;
333  hij = hijtt ; //initialization
334 
335  // The trace h = f_{ij} h^{ij} :
336  Scalar htrace(*mp) ;
337 
338  // Value of h at previous step of the iterative procedure below :
339  Scalar htrace_prev(*mp) ;
340  htrace_prev.set_etat_zero() ; // initialisation to zero
341 
342  for (int it=0; it<=it_max; it++) {
343 
344  // Trace h from the condition det(f^{ij} + h^{ij}) = det f^{ij} :
345 
346  htrace = hij(1,1) * hij(2,3) * hij(2,3)
347  + hij(2,2) * hij(1,3) * hij(1,3) + hij(3,3) * hij(1,2) * hij(1,2)
348  - 2.* hij(1,2) * hij(1,3) * hij(2,3)
349  - hij(1,1) * hij(2,2) * hij(3,3) ;
350 
351  htrace.dec_dzpuis(2) ; // dzpuis: 6 --> 4
352 
353  htrace += hij(1,2) * hij(1,2) + hij(1,3) * hij(1,3)
354  + hij(2,3) * hij(2,3) - hij(1,1) * hij(2,2)
355  - hij(1,1) * hij(3,3) - hij(2,2) * hij(3,3) ;
356 
357  // New value of hij from htrace and hijtt (obtained by solving
358  // the Poisson equation for Phi) :
359 
360  set_tt_trace(hijtt, htrace) ;
361 
362  double diff = max(max(abs(htrace - htrace_prev))) ;
363 #ifndef NDEBUG
364  cout << "Sym_tensor_trans::trace_from_det_one : "
365  << "iteration : " << it << " convergence on trace(h): " << diff << endl ;
366 #endif
367  if (diff < precis) break ;
368  else htrace_prev = htrace ;
369 
370  if (it == it_max) {
371  cout << "Sym_tensor_trans::trace_from_det_one : convergence not reached \n" ;
372  cout << " for the required accuracy (" << precis << ") ! " << endl ;
373  abort() ;
374  }
375  }
376 
377 }
378 
379 
380 
381 }
Metric for tensor calculation.
Definition: metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition: sym_tensor.h:617
Lorene prototypes.
Definition: app_hor.h:67
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
Base class for coordinate mappings.
Definition: map.h:682
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.h:885
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
Scalar * p_trace
Trace with respect to the metric *met_div.
Definition: sym_tensor.h:620
void set_tt_trace(const Sym_tensor_tt &a, const Scalar &h, Param *par=0x0)
Assigns the derived members p_tt and p_trace and updates the components accordingly.
virtual void operator=(const Sym_tensor &a)
Assignment to another Sym_tensor.
Definition: sym_tensor.C:237
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
virtual void operator=(const Sym_tensor_trans &a)
Assignment to another Sym_tensor_trans.
Scalar trace() const
Trace on two different type indices for a valence 2 tensor.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:817
virtual void del_deriv() const
Deletes the derived quantities.
Sym_tensor_tt * p_tt
Traceless part with respect to the metric *met_div.
Definition: sym_tensor.h:623
void set_der_0x0() const
Sets the pointers on derived quantities to 0x0.
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a cova...
Definition: tensor.h:316
virtual void del_deriv() const
Deletes the derived quantities.
Definition: sym_tensor.C:289
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
const Sym_tensor_tt & tt_part(Param *par=0x0) const
Returns the transverse traceless part of the tensor, the trace being defined with respect to metric *...
virtual ~Sym_tensor_trans()
Destructor.
Parameter storage.
Definition: param.h:125
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
const Metric & get_met_div() const
Returns the metric with respect to which the divergence and the trace are defined.
Definition: sym_tensor.h:672
void trace_from_det_one(const Sym_tensor_tt &htt, double precis=1.e-14, int it_max=100)
Assigns the derived member p_tt and computes the trace so that *this + the flat metric has a determin...
Tensor handling.
Definition: tensor.h:294
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:611
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:825
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
Sym_tensor_trans(const Map &map, const Base_vect &triad_i, const Metric &met)
Standard constructor.
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Scalar & the_trace() const
Returns the trace of the tensor with respect to metric *met_div.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Transverse and traceless symmetric tensors of rank 2.
Definition: sym_tensor.h:933