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