LORENE
sym_tensor_decomp.C
1 /*
2  * Methods transverse( ) and longit_pot( ) of class Sym_tensor
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 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 char sym_tensor_decomp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_decomp.C,v 1.14 2014/10/13 08:53:43 j_novak Exp $" ;
29 
30 /*
31  * $Id: sym_tensor_decomp.C,v 1.14 2014/10/13 08:53:43 j_novak Exp $
32  * $Log: sym_tensor_decomp.C,v $
33  * Revision 1.14 2014/10/13 08:53:43 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.13 2008/12/05 08:45:12 j_novak
37  * Modified dzpuis treatment.
38  *
39  * Revision 1.12 2008/12/03 10:20:00 j_novak
40  * Modified output.
41  *
42  * Revision 1.11 2004/06/17 06:56:42 e_gourgoulhon
43  * Simplified code for method transverse (use of Vector::ope_killing).
44  * Slight modif. of output in method longit_pot.
45  *
46  * Revision 1.10 2004/06/14 20:45:41 e_gourgoulhon
47  * Added argument method_poisson to methods longit_pot and transverse.
48  *
49  * Revision 1.9 2004/05/25 15:05:57 f_limousin
50  * Add parameters in argument of the functions transverse, longit_pot
51  * for the case of Map_et.
52  *
53  * Revision 1.8 2004/03/30 14:01:19 j_novak
54  * Copy constructors and operator= now copy the "derived" members.
55  *
56  * Revision 1.7 2004/03/30 08:01:16 j_novak
57  * Better treatment of dzpuis in mutators.
58  *
59  * Revision 1.6 2004/03/29 16:13:07 j_novak
60  * New methods set_longit_trans and set_tt_trace .
61  *
62  * Revision 1.5 2004/02/09 12:56:27 e_gourgoulhon
63  * Method longit_pot: added test of the vector Poisson equation.
64  *
65  * Revision 1.4 2004/02/02 09:18:11 e_gourgoulhon
66  * Method longit_pot: treatment of case divergence dzpuis = 5.
67  *
68  * Revision 1.3 2003/12/10 10:17:54 e_gourgoulhon
69  * First operational version.
70  *
71  * Revision 1.2 2003/11/27 16:01:47 e_gourgoulhon
72  * First implmentation.
73  *
74  * Revision 1.1 2003/11/26 21:57:03 e_gourgoulhon
75  * First version; not ready yet.
76  *
77  *
78  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_decomp.C,v 1.14 2014/10/13 08:53:43 j_novak Exp $
79  *
80  */
81 
82 
83 // Lorene headers
84 #include "metric.h"
85 #include "param.h"
86 
87 namespace Lorene {
89  const Sym_tensor_trans& ht ) {
90 
91  assert ( v_pot.get_index_type(0) == CON ) ;
92 
93  const Metric& metre = ht.get_met_div() ;
94 
95  *this = ht + v_pot.ope_killing(metre) ; // this has dzpuis = 2, if v_pot not 0
96  if ((*this)(1,1).get_dzpuis() == 2)
97  dec_dzpuis(2) ; // which is decreased so to add *this to a flat metric
98 
99  del_deriv() ;
100 
101  set_dependance(metre) ;
102  int jp = get_place_met(metre) ;
103  assert ((jp>=0) && (jp<N_MET_MAX)) ;
104 
105  p_transverse[jp] = new Sym_tensor_trans(ht) ;
106  p_longit_pot[jp] = new Vector( v_pot ) ;
107 
108 }
109 
111  int method_poisson) const {
112 
113  set_dependance(metre) ;
114  int jp = get_place_met(metre) ;
115  assert ((jp>=0) && (jp<N_MET_MAX)) ;
116 
117  if (p_transverse[jp] == 0x0) { // a new computation is necessary
118 
119  assert( (type_indice(0) == CON) && (type_indice(1) == CON) ) ;
120 
121  for (int ic=0; ic<n_comp; ic++) {
122  assert(cmp[ic]->check_dzpuis(4)) ; // dzpuis=4 is assumed
123  }
124 
125  const Vector& ww = longit_pot(metre, par, method_poisson) ;
126 
127  Sym_tensor lww = ww.ope_killing(metre) ; // D^i W^j + D^j W^i
128 
129  lww.inc_dzpuis(2) ;
130 
131  p_transverse[jp] = new Sym_tensor_trans(*mp, *triad, metre) ;
132 
133  *(p_transverse[jp]) = *this - lww ;
134 
135  }
136 
137  return *p_transverse[jp] ;
138 
139 
140 }
141 
142 
143 const Vector& Sym_tensor::longit_pot(const Metric& metre, Param* par,
144  int method_poisson) const {
145 
146  set_dependance(metre) ;
147  int jp = get_place_met(metre) ;
148  assert ((jp>=0) && (jp<N_MET_MAX)) ;
149 
150  if (p_longit_pot[jp] == 0x0) { // a new computation is necessary
151 
152  const Metric_flat* metf = dynamic_cast<const Metric_flat*>(&metre) ;
153  if (metf == 0x0) {
154  cout << "Sym_tensor::longit_pot : the case of a non flat metric"
155  << endl <<" is not treated yet !" << endl ;
156  abort() ;
157  }
158 
159  Vector hhh = divergence(metre) ;
160 
161 
162  // If dpzuis is 5, it should be decreased to 4 for the Poisson equation:
163  bool dzp5 = false ;
164  for (int i=1; i<=3; i++) {
165  dzp5 = dzp5 || hhh(i).check_dzpuis(5) ;
166  }
167  if (dzp5) hhh.dec_dzpuis() ;
168 
169  if (dynamic_cast<const Map_af*>(mp) != 0x0)
170  p_longit_pot[jp] = new Vector( hhh.poisson(double(1),
171  method_poisson) ) ;
172  else
173  p_longit_pot[jp] = new Vector( hhh.poisson(double(1), *par,
174  method_poisson) ) ;
175 
176 
177  // Test of resolution of the vector Poisson equation:
178  const Vector& ww = *(p_longit_pot[jp]) ;
179 
180  hhh.dec_dzpuis() ;
181 
182  Vector lapw = ww.derive_con(metre).divergence(metre)
183  + (ww.divergence(metre)).derive_con(metre) ;
184 #ifndef NDEBUG
185  cout << "## Sym_tensor::longit_pot : test of Poisson : \n" ;
186  cout <<
187  " Max absolute error in each domain on the vector Poisson equation:\n" ;
188  maxabs(lapw - hhh) ;
189 
190  int nz = mp->get_mg()->get_nzone() ; // total number of domains
191 
192  cout << " Relative error in each domain on the vector Poisson equation:\n" ;
193  for (int i=1; i<=3; i++){
194  cout << " Comp. " << i << " : " ;
195  for (int l=0; l<nz; l++){
196  cout << diffrel(lapw(i),hhh(i) )(l) << " " ;
197  }
198  cout << endl ;
199  }
200  cout << endl ;
201 #endif
202  }
203 
204  return *p_longit_pot[jp] ;
205 
206 
207 }
208 
209 
210 }
Affine radial mapping.
Definition: map.h:2027
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Flat metric for tensor calculation.
Definition: metric.h:261
Metric for tensor calculation.
Definition: metric.h:90
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Parameter storage.
Definition: param.h:125
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:608
const Metric & get_met_div() const
Returns the metric with respect to which the divergence and the trace are defined.
Definition: sym_tensor.h:669
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
void set_longit_trans(const Vector &v, const Sym_tensor_trans &a)
Assigns the derived members p_longit_pot and p_transverse and updates the components accordingly.
const Sym_tensor_trans & transverse(const Metric &gam, Param *par=0x0, int method_poisson=6) const
Computes the transverse part of the tensor with respect to a given metric, transverse meaning diverg...
const Vector & longit_pot(const Metric &gam, Param *par=0x0, int method_poisson=6) const
Computes the vector potential of longitudinal part of the tensor (see documentation of method transv...
Sym_tensor_trans * p_transverse[N_MET_MAX]
Array of the transverse part of the tensor with respect to various metrics, transverse meaning diver...
Definition: sym_tensor.h:239
virtual void del_deriv() const
Deletes the derived quantities.
Definition: sym_tensor.C:286
Vector * p_longit_pot[N_MET_MAX]
Array of the vector potential of the longitudinal part of the tensor with respect to various metrics ...
Definition: sym_tensor.h:246
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Tensor field of valence 1.
Definition: vector.h:188
Sym_tensor ope_killing(const Metric &gam) const
Computes the Killing operator associated with a given metric.
Definition: vector.C:438
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:886
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
int get_place_met(const Metric &) const
Returns the position of the pointer on metre in the array met_depend .
Definition: tensor.C:443
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
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:816
void set_dependance(const Metric &) const
To be used to describe the fact that the derivatives members have been calculated with met .
Definition: tensor.C:453
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
int n_comp
Number of stored components, depending on the symmetry.
Definition: tensor.h:312
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 Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
Definition: tensor.C:1055
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
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Lorene prototypes.
Definition: app_hor.h:64