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