LORENE
map_et_adapt.C
1/*
2 * Method of the class Map_et to compute the mapping parameters to let the
3 * domain boundaries coincide with isosurfaces of a given scalar field.
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_et_adapt_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $" ;
31
32/*
33 * $Id: map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $
34 * $Log: map_et_adapt.C,v $
35 * Revision 1.11 2014/10/13 08:53:03 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.10 2014/10/06 15:13:12 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.9 2010/01/31 16:58:39 e_gourgoulhon
42 * Back to rev. 1.7 (1.8 has been committed by mistake).
43 *
44 * Revision 1.8 2010/01/31 16:51:47 e_gourgoulhon
45 * File local_settings_linux is now called local_settings_linux_old (for
46 * it is quite old and not compatible with the gcc 4.x series).
47 *
48 * Revision 1.7 2006/09/05 13:39:45 p_grandclement
49 * update of the bin_ns_bh project
50 *
51 * Revision 1.6 2006/05/17 13:27:12 j_novak
52 * Removed strange character on line 534.
53 *
54 * Revision 1.5 2006/05/03 11:15:25 p_grandclement
55 * modif filtering
56 *
57 * Revision 1.4 2006/05/03 07:07:52 p_grandclement
58 * Petite correction
59 *
60 * Revision 1.3 2006/04/25 07:21:59 p_grandclement
61 * Various changes for the NS_BH project
62 *
63 * Revision 1.2 2003/10/03 15:58:48 j_novak
64 * Cleaning of some headers
65 *
66 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
67 * LORENE
68 *
69 * Revision 2.6 2000/11/10 15:19:300 eric
70 * Appel de Valeur::equipot_outward plutot que Valeur::equipot
71 * dans la determination des isopotentielles.
72 *
73 * Revision 2.5 2000/10/23 14:01:17 eric
74 * Changement des arguments (en autre, ajout de nz_search, qui est
75 * desormais a priori independant de nzadapt).
76 *
77 * Revision 2.4 2000/10/20 13:13:01 eric
78 * nzet --> nzadapt.
79 * nz_search = nzadapt --> nz_search = nzadapt + 1
80 *
81 * Revision 2.3 2000/02/17 19:00:53 eric
82 * nz_search = nzet et non plus nz_search = nz-1.
83 *
84 * Revision 2.2 2000/02/15 15:09:12 eric
85 * Changement du Param : fact_echelle est desormais passe en double_mod.
86 *
87 * Revision 2.1 2000/01/03 15:46:59 eric
88 * *** empty log message ***
89 *
90 * Revision 2.0 2000/01/03 13:30:59 eric
91 * *** empty log message ***
92 *
93 *
94 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_adapt.C,v 1.11 2014/10/13 08:53:03 j_novak Exp $
95 *
96 */
97
98// Headers C
99#include <cassert>
100
101// Headers Lorene
102#include "cmp.h"
103#include "itbl.h"
104#include "param.h"
105#include "proto.h"
106
107namespace Lorene {
108void Map_et::adapt(const Cmp& ent, const Param& par, int nbr_filtre) {
109
110 // Parameters of the computation
111 // -----------------------------
112
113 int nitermax = par.get_int(0) ;
114 int& niter = par.get_int_mod(0) ;
115 int nzadapt = par.get_int(1) ;
116 int nz_search = par.get_int(2) ;
117 int adapt_flag = par.get_int(3) ;
118 int j_bord = par.get_int(4) ; // index of theta_*
119 int k_bord = par.get_int(5) ; // index of phi_*
120 double precis = par.get_double(0) ;
121 double fact_lamu = par.get_double(1) ;
122 double fact_echelle = par.get_double(2) ;
123
124 const Tbl& ent_limit = par.get_tbl(0) ;
125
126 int nz = mg->get_nzone() ;
127 int nt = mg->get_nt(0) ;
128 int np = mg->get_np(0) ;
129
130
131 // Protections
132 // -----------
133 assert(ent.get_mp() == this) ;
134 assert(nzadapt < nz) ;
135 assert(nzadapt > 0) ;
136 assert(nz_search >= nzadapt) ;
137 for (int l=1; l<nz; l++) {
138 assert( mg->get_nt(l) == nt ) ;
139 assert( mg->get_np(l) == np ) ;
140 }
141 assert(ent_limit.get_ndim() == 1) ;
142 assert(ent_limit.get_taille() >= nzadapt) ;
143 assert(ent_limit.get_etat() == ETATQCQ) ;
144
145 // Initializations
146 // ---------------
147
148 niter = 0 ;
149 const double xi_max = 1 ;
150
151 //=======================================================================
152 // Special case of a fixed mapping : only a rescaling is performed
153 //=======================================================================
154
155 if ( (adapt_flag == 0) || (nzadapt == 0) ) {
156
157 homothetie(fact_echelle) ;
158
159 return ;
160 }
161
162 //=======================================================================
163 // General case
164 //=======================================================================
165
166 // New mapping parameters
167 // ----------------------
168
169 double* nalpha = new double[nz] ;
170 double* nbeta = new double[nz] ;
171
172 Itbl l_ext(np, nt) ;
173 Tbl x_ext(np, nt) ;
174
175 Tbl xtgg(np, nt, 1) ; // working array constructed from the isosurface
176 // equation
177 Tbl xtff(np, nt, 1) ; // idem
178
179 //----------------------------------------------------------------
180 // Adaptation of the nucleus
181 //----------------------------------------------------------------
182
183 double ent0 = ent_limit(0) ;
184
185 // Search for the equipotential surface ent = ent0
186 // ---> l = l_ext(theta, phi)
187 // xi = x_ext(theta, phi)
188 // -----------------------------------------------
189
190 int niter0 ;
191 (ent.va).equipot_outward(ent0, nz_search, precis, nitermax, niter0,
192 l_ext, x_ext) ;
193
194 niter = ( niter0 > niter ) ? niter0 : niter ;
195
196 // The new mapping must be such that the found isosurface corresponds
197 // to xi = 1
198
199 // Computation of r_ext(theta, phi) - r_eq ---> xtgg
200 // -------------------------------------------------
201
202 xtgg.set_etat_qcq() ;
203 assert(l_ext.get_etat() == ETATQCQ) ;
204
205 double r_eq = val_r_jk(0, xi_max, j_bord, k_bord) ;
206
207 double* pxtgg = xtgg.t ;
208 int* pl_ext = l_ext.t ;
209 double* px_ext = x_ext.t ;
210
211 for (int k=0; k<np; k++) {
212 for (int j=0; j<nt; j++) {
213
214 *pxtgg = val_r_jk(*pl_ext, *px_ext, j, k) - r_eq ;
215
216 // ... next one
217 pl_ext++ ;
218 px_ext++ ;
219 pxtgg++ ;
220 }
221 }
222
223 // Decomposition into even and odd harmonics in phi of xtgg = r_ext - r_eq :
224 // even phi harmonics --> xtgg
225 // odd phi harmonics --> xtff
226 // ------------------------------------------------------------------------
227
228 int base_p = ( ((ent.va).base).b[0] ) & MSQ_P ;
229
230 switch( base_p ) {
231
232 case P_COSSIN_P : { // case of only even harmonics in phi
233
234 xtff.set_etat_zero() ;
235 break ;
236 }
237
238 case P_COSSIN : { // general case: a Fourier transform must be
239 // done
240
241 xtff.set_etat_qcq() ;
242 double* pxtff = xtff.t ;
243 pxtgg = xtgg.t ;
244
245 assert(np>=4) ;
246 int deg[3] ;
247 int dimc[3] ;
248 deg[0] = np ;
249 deg[1] = nt ;
250 deg[2] = 1;
251 dimc[0] = np + 2 ;
252 dimc[1] = nt ;
253 dimc[2] = 1 ;
254 double* cf = new double[(np+2)*nt] ;
255 double* cf0 = new double[(np+2)*nt] ;
256 double* ff0 = new double[np*nt] ;
257
258 for (int i=0; i < np*nt; i++) {
259 cf[i] = *pxtgg ;
260 pxtgg++ ;
261 }
262
263 cfpcossin(deg, dimc, cf) ; // Fourier transform
264
265 // Even harmonics
266 // --------------
267 double* pcf0 = cf0 ;
268 double* pcf = cf ;
269 for (int k=0; k<np-1; k += 4) {
270 for(int j=0; j<2*nt; j++) {
271 *pcf0 = *pcf ;
272 pcf0++ ;
273 pcf++ ;
274 }
275 for(int j=0; j<2*nt; j++) {
276 *pcf0 = 0 ;
277 pcf0++ ;
278 pcf++ ;
279 }
280 }
281 if (np%4 == 0) { // particular case of np multiple of 4
282 for(int j=0; j<2*nt; j++) {
283 *pcf0 = *pcf ;
284 pcf0++ ;
285 pcf++ ;
286 }
287 }
288
289 cipcossin(deg, dimc, deg, cf0, ff0) ; // Inverse Fourier transform
290
291 pxtgg = xtgg.t ;
292 for (int i=0; i < np*nt; i++) {
293 *pxtgg = ff0[i] ;
294 pxtgg++ ;
295 }
296
297 // Odd harmonics
298 // -------------
299 pcf0 = cf0 ;
300 pcf = cf ;
301 for (int k=0; k<np-1; k += 4) {
302 for(int j=0; j<2*nt; j++) {
303 *pcf0 = 0 ;
304 pcf0++ ;
305 pcf++ ;
306 }
307 for(int j=0; j<2*nt; j++) {
308 *pcf0 = *pcf ;
309 pcf0++ ;
310 pcf++ ;
311 }
312 }
313 if (np%4 == 0) { // particular case of np multiple of 4
314 for(int j=0; j<2*nt; j++) {
315 *pcf0 = 0 ;
316 pcf0++ ;
317 }
318 }
319
320 cipcossin(deg, dimc, deg, cf0, ff0) ; // TF inverse
321
322 pxtff = xtff.t ;
323 for (int i=0; i < np*nt; i++) {
324 *pxtff = ff0[i] ;
325 pxtff++ ;
326 }
327
328 delete [] cf ;
329 delete [] cf0 ;
330 delete [] ff0 ;
331 break ;
332 }
333
334 default : {
335 cout << "Map_et::adapt: unknown phi basis !" << endl ;
336 cout << " base_p = " << base_p << endl ;
337 abort() ;
338 }
339 }
340
341 // Computation of lambda and mu in the nucleus :
342 // --------------------------------------------
343
344 double lambda = 0 ; // lambda is set to zero because F(theta, phi)
345 // must not have constant terms in the nucleus
346
347 double mu = - fact_lamu * min(xtgg) ;
348
349 // Computation of alpha and beta in the nucleus :
350 // --------------------------------------------
351
352 nalpha[0] = r_eq - lambda - mu ;
353 nbeta[0] = 0 ;
354
355 // Computation of F_0, G_0 and {\tilde F_1} :
356 // ------------------------------------------
357
358 Mtbl nff(ff.get_mg()) ;
359 Mtbl ngg(gg.get_mg()) ;
360 nff.set_etat_qcq() ;
361 ngg.set_etat_qcq() ;
362
363 *(nff.t[0]) = ( xtff + lambda ) / nalpha[0] ;
364 *(ngg.t[0]) = ( xtgg + mu ) / nalpha[0] ;
365 xtff += xtgg ;
366
367
368 //----------------------------------------------------------------
369 // Adaptation of the shells
370 //----------------------------------------------------------------
371
372 double r_eqlm1 = r_eq ;
373
374
375 // Loop on the shells where the adaptation must be performed
376
377 for (int l=1; l<nzadapt; l++) {
378
379 ent0 = ent_limit(l) ;
380
381 // Search for the equipotential surface ent = ent0
382 // ---> l = l_ext(theta, phi)
383 // xi = x_ext(theta, phi)
384 // -----------------------------------------------
385
386 (ent.va).equipot_outward(ent0, nz_search, precis, nitermax, niter0,
387 l_ext, x_ext) ;
388
389 niter = ( niter0 > niter ) ? niter0 : niter ;
390
391 // The new mapping must be such that the found isosurface corresponds
392 // to xi = 1
393
394 // Computation of r_ext(theta, phi) - r_eq ---> xtgg
395 // -------------------------------------------------
396
397 xtgg.set_etat_qcq() ;
398 assert(l_ext.get_etat() == ETATQCQ) ;
399
400 r_eq = val_r_jk(l, xi_max, j_bord, k_bord) ;
401
402 pxtgg = xtgg.t ;
403 pl_ext = l_ext.t ;
404 px_ext = x_ext.t ;
405
406 for (int k=0; k<np; k++) {
407 for (int j=0; j<nt; j++) {
408
409 *pxtgg = val_r_jk(*pl_ext, *px_ext, j, k) - r_eq ;
410
411 // ... next one
412 pl_ext++ ;
413 px_ext++ ;
414 pxtgg++ ;
415 }
416 }
417
418 // Computation of lambda and mu in domain no. l :
419 // --------------------------------------------
420
421 lambda = - fact_lamu * max(xtff) ;
422 mu = - fact_lamu * min(xtgg) ;
423
424 // Computation of alpha and beta in domain no. l :
425 // --------------------------------------------
426
427 nalpha[l] = .5 * ( r_eq - r_eqlm1 + lambda - mu ) ;
428 nbeta[l] = .5 * ( r_eq + r_eqlm1 - lambda - mu ) ;
429
430 // Computation of F_l, G_l and {\tilde F_(l+1)} :
431 // ------------------------------------------
432
433 *(nff.t[l]) = ( xtff + lambda ) / nalpha[l] ;
434 *(ngg.t[l]) = ( xtgg + mu ) / nalpha[l] ;
435 xtff = xtgg ;
436
437 r_eqlm1 = r_eq ;
438
439 } // end of the loop on the shells where the adaptation must be performed
440
441 //----------------------------------------------------------------
442 // Adaptation of the domain of index nzadapt
443 //----------------------------------------------------------------
444
445 if (mg->get_type_r(nzadapt) == UNSURR) {
446
447 // Case of a compactified domain
448 // -----------------------------
449
450 xtff = 1 / (xtff + r_eqlm1) - double(1) / r_eqlm1 ;
451
452 lambda = - fact_lamu * min(xtff) ;
453
454 nalpha[nzadapt] = .5 * ( lambda - double(1) / r_eqlm1 ) ;
455 nbeta[nzadapt] = - nalpha[nzadapt] ;
456
457 // Computation of F_nzadapt :
458 *(nff.t[nzadapt]) = ( xtff + lambda ) / nalpha[nzadapt] ;
459
460 }
461 else {
462 // Case of a non-compactified domain
463 // ----------------------------------
464
465 r_eq = val_r_jk(nzadapt, xi_max, j_bord, k_bord) ;
466
467 lambda = - fact_lamu * max(xtff) ;
468
469 nalpha[nzadapt] = .5 * ( r_eq - r_eqlm1 + lambda ) ;
470 nbeta[nzadapt] = .5 * ( r_eq + r_eqlm1 - lambda ) ;
471
472 // Computation of F_l :
473
474 *(nff.t[nzadapt]) = ( xtff + lambda ) / nalpha[nzadapt] ;
475
476 }
477
478 // In all cases, G_nzadapt is set to zero :
479 ngg.t[nzadapt]->set_etat_zero() ;
480
481 //----------------------------------------------------------------
482 // Values of alpha, beta, F and G in the most external domains
483 // where the mapping is unchanged
484 //----------------------------------------------------------------
485
486 for (int l=nzadapt+1; l<nz; l++) {
487
488 nalpha[l] = alpha[l] ;
489 nbeta[l] = beta[l] ;
490
491 }
492
493 if (ff.get_etat() == ETATZERO) {
494 for (int l=nzadapt+1; l<nz; l++) {
495 nff.t[l]->set_etat_zero() ;
496 }
497 }
498 else {
499 assert(ff.get_etat() == ETATQCQ) ;
500 assert(ff.c != 0x0) ;
501 assert(ff.c->get_etat() == ETATQCQ) ;
502 for (int l=nzadapt+1; l<nz; l++) {
503 *(nff.t[l]) = *(ff.c->t[l]) ;
504 }
505 }
506
507 if (gg.get_etat() == ETATZERO) {
508 for (int l=nzadapt+1; l<nz; l++) {
509 ngg.t[l]->set_etat_zero() ;
510 }
511 }
512 else {
513 assert(gg.get_etat() == ETATQCQ) ;
514 assert(gg.c != 0x0) ;
515 assert(gg.c->get_etat() == ETATQCQ) ;
516 for (int l=nzadapt+1; l<nz; l++) {
517 *(ngg.t[l]) = *(gg.c->t[l]) ;
518 }
519 }
520
521
522 //=============================================================================
523 // Assignment of the mapping parameters alpha, beta, ff and gg to
524 // the newly computed values
525 //=============================================================================
526
527 for (int l=0; l<nz; l++) {
528
529 if (mg->get_type_r(l) == UNSURR) {
530 alpha[l] = nalpha[l] / fact_echelle ;
531 beta[l] = nbeta[l] / fact_echelle ;
532 }
533 else {
534 alpha[l] = fact_echelle * nalpha[l] ;
535 beta[l] = fact_echelle * nbeta[l] ;
536 }
537
538 }
539
540 ff = nff ;
541 gg = ngg ;
542
543 ff.std_base_scal() ; // Standard spectral bases for F
544 gg.std_base_scal() ; // Standard spectral bases for G
545
546
547 if (nbr_filtre !=0) {
548 ff.coef() ;
549 gg.coef() ;
552 for (int l=0 ; l<nzadapt+1 ; l++)
553 for (int k=np-nbr_filtre ; k<np ; k++)
554 for (int j=0 ; j<nt ; j++) {
555 if (ff.c_cf->t[l]->get_etat() != ETATZERO)
556 ff.c_cf->set(l, k,j,0) = 0 ;
557
558 if (gg.c_cf->t[l]->get_etat() != ETATZERO)
559 gg.c_cf->set(l,k,j,0) = 0 ;
560 }
561 }
562
563 // The derived quantities must be reset
564 // ------------------------------------
565
566 reset_coord() ;
567
568
569 // Clean exit
570 // ----------
571
572 delete [] nalpha ;
573 delete [] nbeta ;
574
575}
576}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Basic integer array class.
Definition itbl.h:122
int get_etat() const
Gives the logical state.
Definition itbl.h:317
int * t
The array of int 's.
Definition itbl.h:132
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2760
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
virtual void adapt(const Cmp &ent, const Param &par, int nbr_filtre=0)
Adaptation of the mapping to a given scalar field.
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2819
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2758
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2826
virtual double val_r_jk(int l, double xi, int j, int k) const
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
virtual void reset_coord()
Resets all the member Coords.
Definition map_et.C:628
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
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_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Parameter storage.
Definition param.h:125
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
const Tbl & get_tbl(int position=0) const
Returns the reference of a Tbl stored in the list.
Definition param.C:567
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition param.C:430
Basic array class.
Definition tbl.h:161
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition tbl.h:400
int get_etat() const
Gives the logical state.
Definition tbl.h:394
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
double * t
The array of double.
Definition tbl.h:173
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
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
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define P_COSSIN
dev. standart
#define MSQ_P
Extraction de l'info sur Phi.
Lorene prototypes.
Definition app_hor.h:64