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