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
30char 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
79namespace Lorene {
80
81void 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}
Bases of the spectral expansions.
Definition base_val.h:322
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition cmp.C:323
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
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
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
double * x
Array of values of at the nr collocation points.
Definition grilles.h:209
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2035
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
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
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Multi-domain array.
Definition mtbl.h:118
int get_etat() const
Gives the logical state.
Definition mtbl.h:277
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Parameter storage.
Definition param.h:125
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition param.C:498
int get_etat() const
Gives the logical state.
Definition tbl.h:394
double * t
The array of double.
Definition tbl.h:173
Values and coefficients of a (real-value) function.
Definition valeur.h:287
int get_etat() const
Returns the logical state.
Definition valeur.h:726
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
void coef_i() const
Computes the physical value of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
#define T_COS_P
dev. cos seulement, harmoniques paires
#define T_SIN_I
dev. sin seulement, harmoniques impaires
#define T_COS
dev. cos seulement
#define T_SIN
dev. sin seulement
Lorene prototypes.
Definition app_hor.h:64