LORENE
map_af_poisson2d.C
1 /*
2  * Method of the class Map_af for the resolution of the 2-D Poisson
3  * equation.
4  *
5  * (see file map.h for documentation).
6  */
7 
8 /*
9  * Copyright (c) 2000-2001 Eric Gourgoulhon
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 char map_af_poisson2d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $" ;
31 
32 /*
33  * $Id: map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $
34  * $Log: map_af_poisson2d.C,v $
35  * Revision 1.6 2014/10/13 08:53:02 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.5 2012/09/04 14:53:28 j_novak
39  * Replacement of the FORTRAN version of huntm by a C one.
40  *
41  * Revision 1.4 2012/08/12 17:48:36 p_cerda
42  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
43  *
44  * Revision 1.3 2002/09/09 13:54:20 e_gourgoulhon
45  *
46  * Change of name of the Fortran subroutines
47  * poisson2d -> poiss2d
48  * poisson2di -> poiss2di
49  * to avoid any clash with Map::poisson2d and Map::poisson2di
50  *
51  * Revision 1.2 2002/09/09 13:00:39 e_gourgoulhon
52  * Modification of declaration of Fortran 77 prototypes for
53  * a better portability (in particular on IBM AIX systems):
54  * All Fortran subroutine names are now written F77_* and are
55  * defined in the new file C++/Include/proto_f77.h.
56  *
57  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
58  * LORENE
59  *
60  * Revision 2.1 2000/10/11 15:15:26 eric
61  * 1ere version operationnelle.
62  *
63  * Revision 2.0 2000/10/09 13:47:10 eric
64  * *** empty log message ***
65  *
66  *
67  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $
68  *
69  */
70 
71 // Header Lorene:
72 #include "map.h"
73 #include "cmp.h"
74 #include "param.h"
75 #include "proto_f77.h"
76 
77 //*****************************************************************************
78 
79 namespace Lorene {
80 
81 void Map_af::poisson2d(const Cmp& source_mat, const Cmp& source_quad,
82  Param& par, Cmp& uu) const {
83 
84  assert(source_mat.get_etat() != ETATNONDEF) ;
85  assert(source_quad.get_etat() != ETATNONDEF) ;
86  assert(source_mat.get_mp()->get_mg() == mg) ;
87  assert(source_quad.get_mp()->get_mg() == mg) ;
88  assert(uu.get_mp()->get_mg() == mg) ;
89 
90  assert( source_quad.check_dzpuis(4) ) ;
91 
92  int mpsymm = uu.get_mp()->get_mg()->get_type_t();
93 
94  double& lambda = par.get_double_mod(0) ;
95 
96  // Special case of a vanishing source
97  // ----------------------------------
98 
99  if ( (source_mat.get_etat() == ETATZERO)
100  && (source_quad.get_etat() == ETATZERO) ) {
101 
102  uu = 0 ;
103  lambda = 1 ;
104  return ;
105  }
106 
107  // ---------------------------------------------------------------------
108  // Preparation of the parameters for the Fortran subroutines F77_poisson2d
109  // and F77_poisson2di
110  // ---------------------------------------------------------------------
111 
112  int nz = mg->get_nzone() ;
113  int np1 = 1 ; // Axisymmetry enforced
114  int nt = mg->get_nt(0) ;
115  int nt2 = 0 ;
116 
117  switch ( mpsymm ){
118  case SYM: {
119  nt2 = 2*nt - 1 ; // Number of points for the theta sampling
120  break; // in [0,Pi], instead of [0,Pi/2]
121  }
122  case NONSYM: {
123  nt2 = nt;
124  break;
125  }
126  }
127 
128 
129  // Array NDL
130  // ---------
131  int* ndl = new int[nz+4] ;
132  ndl[0] = nz ;
133  for (int l=0; l<nz; l++) {
134  ndl[1+l] = mg->get_nr(l) ;
135  }
136  ndl[1+nz] = nt2 ;
137  ndl[2+nz] = np1 ;
138  ndl[3+nz] = nz ;
139 
140  // Array INDD
141  // ----------
142  int* indd = new int[nz] ;
143  for (int l=0; l<nz; l++) {
144  switch ( mg->get_type_r(l) ) {
145  case RARE : {
146  indd[l] = 0 ;
147  break ;
148  }
149  case FIN : {
150  indd[l] = 1 ;
151  break ;
152  }
153  case UNSURR : {
154  indd[l] = 2 ;
155  break ;
156  }
157  default : {
158  cout << "Map_af::poisson2d: unknown type_r !" << endl ;
159  abort() ;
160  break ;
161  }
162  }
163  }
164 
165  // Parameters NDR, NDT, NDP and NDZ
166  // --------------------------------
167  int nrmax = 0 ;
168  for (int l=0; l<nz ; l++) {
169  nrmax = ( ndl[1+l] > nrmax ) ? ndl[1+l] : nrmax ;
170  }
171  int ndr = nrmax + 5 ; // Le +5 est impose par les routines de resolution
172  // de l'equation de Poisson du type gr2p3s_
173  int ndt = nt2 + 2 ;
174  int ndp = np1 + 2 ;
175  int ndz = nz ;
176 
177  // Array ERRE
178  // ----------
179 
180  double* erre = new double [ndz*ndr] ;
181 
182  for (int l=0; l<nz; l++) {
183  for (int i=0; i<ndl[1+l]; i++) {
184  double xr = mg->get_grille3d(l)->x[i] ;
185  erre[ ndr*l + i ] = alpha[l] * xr + beta[l] ;
186  }
187  }
188 
189  // Arrays containing the data
190  // --------------------------
191 
192  int ndrt = ndr*ndt ;
193  int ndrtp = ndr*ndt*ndp ;
194  int taille = ndrtp*ndz ;
195 
196  double* tsou_m = new double[ taille ] ;
197  double* tsou_q = new double[ taille ] ;
198  double* tuu = new double[ taille ] ;
199 
200  // Initialisation to zero :
201  for (int i=0; i<taille; i++) {
202  tsou_m[i] = 0 ;
203  tsou_q[i] = 0 ;
204  tuu[i] = 0 ;
205  }
206 
207  // Copy of source_mat into tsou_m
208  // ------------------------------
209  const Valeur& va_m = source_mat.va ;
210  assert(va_m.get_etat() == ETATQCQ) ;
211  va_m.coef_i() ;
212  const Mtbl* s_m = va_m.c ;
213  assert(s_m->get_etat() == ETATQCQ) ;
214 
215  Base_val base_s = va_m.base ;
216 
217  for (int l=0; l<nz; l++) {
218  int nr = mg->get_nr(l) ;
219  int nrt = nr*nt ;
220  if (s_m->t[l]->get_etat() == ETATZERO) {
221  for (int k=0; k<np1; k++) {
222  for (int j=0; j<nt; j++) {
223  for (int i=0; i<nr; i++) {
224  tsou_m[ndrtp*l + ndrt*k + ndr*j + i] = 0 ;
225  // point symetrique par rapport au plan theta = pi/2 :
226  if ( mpsymm == SYM ) tsou_m[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = 0 ;
227  }
228  }
229  }
230 
231  }
232  else {
233  assert( s_m->t[l]->get_etat() == ETATQCQ ) ;
234  for (int k=0; k<np1; k++) {
235  for (int j=0; j<nt; j++) {
236  for (int i=0; i<nr; i++) {
237  double xx = s_m->t[l]->t[nrt*k + nr*j + i] ;
238  tsou_m[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
239  // point symetrique par rapport au plan theta = pi/2 :
240  if ( mpsymm == SYM ) tsou_m[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;
241  }
242  }
243  }
244  } // End of case etat != ETATZERO
245  }
246 
247  // Copy of source_quad into tsou_q
248  // -------------------------------
249 
250  if (source_quad.get_etat() != ETATZERO) {
251 
252  const Valeur& va_q = source_quad.va ;
253  assert(va_q.get_etat() == ETATQCQ) ;
254  va_q.coef_i() ;
255  const Mtbl* s_q = va_q.c ;
256 
257  assert( va_q.base == base_s ) ;
258 
259  assert(s_q->get_etat() == ETATQCQ) ;
260 
261  for (int l=0; l<nz; l++) {
262  int nr = mg->get_nr(l) ;
263  int nrt = nr*nt ;
264  if (s_q->t[l]->get_etat() == ETATZERO) {
265  for (int k=0; k<np1; k++) {
266  for (int j=0; j<nt; j++) {
267  for (int i=0; i<nr; i++) {
268  tsou_q[ndrtp*l + ndrt*k + ndr*j + i] = 0 ;
269  // point symetrique par rapport au plan theta = pi/2 :
270  if ( mpsymm == SYM ) tsou_q[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = 0 ;
271  }
272  }
273  }
274 
275  }
276  else {
277  assert( s_q->t[l]->get_etat() == ETATQCQ ) ;
278  for (int k=0; k<np1; k++) {
279  for (int j=0; j<nt; j++) {
280  for (int i=0; i<nr; i++) {
281  double xx = s_q->t[l]->t[nrt*k + nr*j + i] ;
282  tsou_q[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
283  // point symetrique par rapport au plan theta = pi/2 :
284  if ( mpsymm == SYM ) tsou_q[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;
285  }
286  }
287  }
288  } // End of case s_q->t[l]->etat != ETATZERO
289  }
290 
291  } // End of case source_quad.etat != ETATZERO
292 
293  //-----------------------------------------------------------
294  // Call of the Fortran subroutine poisson2d_ or poisson2di_
295  //-----------------------------------------------------------
296 
297  int base_t = (va_m.base).get_base_t(0) ;
298 
299  Base_val base_uu(nz) ; // Output spectral bases
300 
301  switch (base_t) {
302 
303  case T_COS :
304  case T_COS_P : {
305 
306  double lambda0 ;
307 
308  F77_poiss2d(ndl, &ndr, &ndt, &ndp, indd, erre, tsou_m, tsou_q,
309  &lambda0, tuu) ;
310 
311  base_uu = base_s ; // output bases = input bases
312 
313  lambda = lambda0 ;
314  break ;
315  }
316  case T_SIN :
317  case T_SIN_I : {
318 
319  double* tsou = new double[taille] ;
320  for (int i=0; i<taille; i++) {
321  tsou[i] = tsou_m[i] + tsou_q[i] ;
322  }
323 
324  F77_poiss2di(ndl, &ndr, &ndt, &ndp, indd, erre, tsou, tuu) ;
325 
326  base_uu = base_s ; // output bases = input bases
327 
328  lambda = double(1) ;
329 
330  delete [] tsou ;
331  break ;
332  }
333 
334  default : {
335  cout << "Map_af::poisson2d : unkown theta basis !" << endl ;
336  cout << " basis : " << hex << base_t << endl ;
337  abort() ;
338  break ;
339  }
340  }
341 
342  //-------------------------------
343  // Copy of tuu into uu
344  //-------------------------------
345 
346  uu.allocate_all() ;
347  (uu.va).set_etat_c_qcq() ; // To make sure that the coefficients are
348  // deleted
349 
350  for (int l=0; l<nz; l++) {
351  int nr = mg->get_nr(l) ;
352  for (int k=0; k<mg->get_np(l); k++) {
353  for (int j=0; j<nt; j++) {
354  for (int i=0; i<nr; i++) {
355  uu.set(l, k, j, i) = tuu[ndrtp*l + ndr*j + i] ;
356  }
357  }
358  }
359  }
360 
361  (uu.va).set_base( base_uu ) ; // Bases for spectral expansions
362 
363  uu.set_dzpuis(0) ;
364 
365  // Cleaning
366  // --------
367 
368  delete [] ndl ;
369  delete [] indd ;
370  delete [] erre ;
371  delete [] tsou_m ;
372  delete [] tsou_q ;
373  delete [] tuu ;
374 
375 
376 }
377 
378 
379 }
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Cmp::check_dzpuis
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: cmp.C:715
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
T_COS
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
T_SIN
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid
Definition: valeur.h:299
Lorene::Map_af::beta
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2035
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
Lorene::Mg3d::get_type_t
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
Lorene::Map_af::poisson2d
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const
Computes the solution of a 2-D Poisson equation.
Definition: map_af_poisson2d.C:81
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Mtbl::get_etat
int get_etat() const
Gives the logical state.
Definition: mtbl.h:277
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
Lorene::Mg3d::get_grille3d
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:500
Lorene::Param::get_double_mod
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:498
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition: mtbl.h:132
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Cmp::allocate_all
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: cmp.C:323
Lorene::Grille3d::x
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:209
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition: valeur_coef_i.C:140
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Lorene::Map_af::alpha
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2033
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
Lorene::Map::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676