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