LORENE
sym_tensor_trans_pde.C
1 /*
2  * Functions to solve various PDEs for a divergence-free symmetric tensor.
3  *
4  * (see file sym_tensor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005-2006 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_trans_pde.C,v 1.17 2016/12/05 16:18:17 j_novak Exp $
32  * $Log: sym_tensor_trans_pde.C,v $
33  * Revision 1.17 2016/12/05 16:18:17 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.16 2014/10/13 08:53:43 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.15 2014/10/06 15:13:19 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.14 2010/10/11 10:38:34 j_novak
43  * *** empty log message ***
44  *
45  * Revision 1.13 2010/10/11 10:23:03 j_novak
46  * Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
47  *
48  * Revision 1.12 2006/09/05 15:38:45 j_novak
49  * The fuctions sol_Dirac... are in a seperate file, with new parameters to
50  * control the boundary conditions.
51  *
52  * Revision 1.11 2006/08/31 12:13:22 j_novak
53  * Added an argument of type Param to Sym_tensor_trans::sol_Dirac_A().
54  *
55  * Revision 1.10 2006/06/28 07:48:26 j_novak
56  * Better treatment of some null cases.
57  *
58  * Revision 1.9 2006/06/21 15:42:47 j_novak
59  * Minor changes.
60  *
61  * Revision 1.8 2006/06/20 12:07:15 j_novak
62  * Improved execution speed for sol_Dirac_tildeB...
63  *
64  * Revision 1.7 2006/06/14 10:04:21 j_novak
65  * New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
66  *
67  * Revision 1.6 2006/06/13 13:30:12 j_novak
68  * New members sol_Dirac_A and sol_Dirac_tildeB (see documentation).
69  *
70  * Revision 1.5 2006/06/12 13:37:23 j_novak
71  * Added bounds in l (multipolar momentum) for Sym_tensor_trans::solve_hrr.
72  *
73  * Revision 1.4 2005/11/28 14:45:17 j_novak
74  * Improved solution of the Poisson tensor equation in the case of a transverse
75  * tensor.
76  *
77  * Revision 1.3 2005/11/24 14:07:54 j_novak
78  * Use of Matrice::annule_hard()
79  *
80  * Revision 1.2 2005/11/24 09:24:25 j_novak
81  * Corrected some missing references.
82  *
83  * Revision 1.1 2005/09/16 13:58:11 j_novak
84  * New Poisson solver for a Sym_tensor_trans.
85  *
86  *
87  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_pde.C,v 1.17 2016/12/05 16:18:17 j_novak Exp $
88  *
89  */
90 
91 // C headers
92 #include <cassert>
93 #include <cmath>
94 
95 // Lorene headers
96 #include "tensor.h"
97 #include "diff.h"
98 #include "proto.h"
99 #include "param.h"
100 
101 namespace Lorene {
103 
104  // All this has a meaning only for spherical components...
105  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
106  //## ... and affine mapping, for the moment!
107  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
108  assert( mpaff!= 0x0) ;
109  Sym_tensor_trans resu(*mp, *triad, *met_div) ;
110 
111  const Mg3d& gri = *mp->get_mg() ;
112  int np = gri.get_np(0) ;
113  int nt = gri.get_nt(0) ;
114  assert (nt > 4) ;
115  if (np == 1) {
116  int nz = gri.get_nzone() ;
117  double* bornes = new double[nz+1] ;
118  const double* alp = mpaff->get_alpha() ;
119  const double* bet = mpaff->get_beta() ;
120  for (int lz=0; lz<nz; lz++) {
121  assert (gri.get_np(lz) == np) ;
122  assert (gri.get_nt(lz) == nt) ;
123  switch (gri.get_type_r(lz)) {
124  case RARE: {
125  bornes[lz] = bet[lz] ;
126  break ;
127  }
128  case FIN: {
129  bornes[lz] = bet[lz] - alp[lz] ;
130  break ;
131  }
132  case UNSURR: {
133  bornes[lz] = double(1) / ( bet[lz] - alp[lz] ) ;
134  break ;
135  }
136  default: {
137  cout << "Sym_tensor_trans::poisson() : problem with the grid!"
138  << endl ;
139  abort() ;
140  break ;
141  }
142  }
143  }
144  if (gri.get_type_r(nz-1) == UNSURR)
145  bornes[nz] = 1./(alp[nz-1] + bet[nz-1]) ;
146  else
147  bornes[nz] = alp[nz-1] + bet[nz-1] ;
148 
149  const Mg3d& gr2 = *gri.get_non_axi() ;
150  Map_af mp2(gr2, bornes) ;
151  int np2 = ( np > 3 ? np : 4 ) ;
152 
153  Sym_tensor sou_cart(mp2, CON, mp2.get_bvect_spher()) ;
154  for (int l=1; l<=3; l++)
155  for (int c=l; c<=3; c++) {
156  switch (this->operator()(l,c).get_etat() ) {
157  case ETATZERO: {
158  sou_cart.set(l,c).set_etat_zero() ;
159  break ;
160  }
161  case ETATUN: {
162  sou_cart.set(l,c).set_etat_one() ;
163  break ;
164  }
165  case ETATQCQ : {
166  sou_cart.set(l,c).allocate_all() ;
167  for (int lz=0; lz<nz; lz++)
168  for (int k=0; k<np2; k++)
169  for (int j=0; j<nt; j++)
170  for(int i=0; i<gr2.get_nr(lz); i++)
171  sou_cart.set(l,c).set_grid_point(lz, k, j, i)
172  = this->operator()(l,c).val_grid_point(lz, 0, j, i) ;
173  break ;
174  }
175  default: {
176  cout <<
177  "Sym_tensor_trans::poisson() : source in undefined state!"
178  << endl ;
179  abort() ;
180  break ;
181  }
182  }
183  sou_cart.set(l,c).set_dzpuis(this->operator()(l,c).get_dzpuis()) ;
184  }
185  sou_cart.std_spectral_base() ;
186  sou_cart.change_triad(mp2.get_bvect_cart()) ;
187  Sym_tensor res_cart(mp2, CON, mp2.get_bvect_cart()) ;
188  for (int i=1; i<=3; i++)
189  for(int j=i; j<=3; j++)
190  res_cart.set(i,j) = sou_cart(i,j).poisson() ;
191  res_cart.change_triad(mp2.get_bvect_spher()) ;
192  Scalar res_A(*mp) ; Scalar big_A = res_cart.compute_A() ;
193  Scalar res_B(*mp) ; Scalar big_B = res_cart.compute_tilde_B_tt() ;
194 
195  switch (big_A.get_etat() ) {
196  case ETATZERO: {
197  res_A.set_etat_zero() ;
198  break ;
199  }
200  case ETATUN : {
201  res_A.set_etat_one() ;
202  break ;
203  }
204  case ETATQCQ : {
205  res_A.allocate_all() ;
206  for (int lz=0; lz<nz; lz++)
207  for (int k=0; k<np; k++)
208  for (int j=0; j<nt; j++)
209  for(int i=0; i<gri.get_nr(lz); i++)
210  res_A.set_grid_point(lz, k, j, i)
211  = big_A.val_grid_point(lz, k, j, i) ;
212  break ;
213  }
214  default: {
215  cout <<
216  "Sym_tensor_trans::poisson() : res_A in undefined state!"
217  << endl ;
218  abort() ;
219  break ;
220  }
221  }
222  res_A.set_spectral_base(big_A.get_spectral_base()) ;
223  int dzA = big_A.get_dzpuis() ;
224  res_A.set_dzpuis(dzA) ;
225 
226  switch (big_B.get_etat() ) {
227  case ETATZERO: {
228  res_B.set_etat_zero() ;
229  break ;
230  }
231  case ETATUN : {
232  res_B.set_etat_one() ;
233  break ;
234  }
235  case ETATQCQ : {
236  res_B.allocate_all() ;
237  for (int lz=0; lz<nz; lz++)
238  for (int k=0; k<np; k++)
239  for (int j=0; j<nt; j++)
240  for(int i=0; i<gri.get_nr(lz); i++)
241  res_B.set_grid_point(lz, k, j, i)
242  = big_B.val_grid_point(lz, k, j, i) ;
243  break ;
244  }
245  default: {
246  cout <<
247  "Sym_tensor_trans::poisson() : res_B in undefined state!"
248  << endl ;
249  abort() ;
250  break ;
251  }
252  }
253  res_B.set_spectral_base(big_B.get_spectral_base()) ;
254  int dzB = big_B.get_dzpuis() ;
255  res_B.set_dzpuis(dzB) ;
256 
257  resu.set_AtBtt_det_one(res_A, res_B, h_guess) ;
258 
259  delete [] bornes ;
260  }
261  else {
262  assert (np >=4) ;
263  Sym_tensor_trans sou_cart = *this ;
264  sou_cart.change_triad(mp->get_bvect_cart()) ;
265 
266  Sym_tensor res_cart(*mp, CON, mp->get_bvect_cart()) ;
267  for (int i=1; i<=3; i++)
268  for(int j=i; j<=3; j++)
269  res_cart.set(i,j) = sou_cart(i,j).poisson() ;
270 
271  res_cart.change_triad(*triad) ;
272 
273  resu.set_AtBtt_det_one(res_cart.compute_A(), res_cart.compute_tilde_B_tt(), h_guess) ;
274 
275  }
276 #ifndef NDEBUG
277  Vector dive = resu.divergence(*met_div) ;
278  dive.dec_dzpuis(2) ;
279  maxabs(dive, "Sym_tensor_trans::poisson : divergence of the solution") ;
280 #endif
281  return resu ;
282 }
283 
284 
285 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
const Mg3d * get_non_axi() const
Returns the pointer on the grid which has at least 4 points in the direction and at least 5 in the ...
Definition: mg3d.C:784
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.
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
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
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: scalar.C:1003
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Tensor field of valence 1.
Definition: vector.h:188
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
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
void set_AtBtt_det_one(const Scalar &a_in, const Scalar &tbtt_in, const Scalar *h_prev=0x0, Param *par_bc=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assigns the derived member A and computes from its TT-part (see Sym_tensor::compute_tilde_B_tt() )...
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition: scalar.C:340
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:611
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:803
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Sym_tensor_trans poisson(const Scalar *h_guess=0x0) const
Computes the solution of a tensorial transverse Poisson equation with *this as a source: In particu...
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:807
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
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