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 char 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 
107 namespace Lorene {
108 void 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() ;
550  ff.set_etat_cf_qcq() ;
551  gg.set_etat_cf_qcq() ;
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 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:898
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
#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:462
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
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Multi-domain array.
Definition: mtbl.h:118
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2758
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
#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:394
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:824
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_etat() const
Gives the logical state.
Definition: mtbl.h:277
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
const Tbl & get_tbl(int position=0) const
Returns the reference of a Tbl stored in the list.
Definition: param.C:567
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:400
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Definition: map.h:2819
double * t
The array of double.
Definition: tbl.h:173
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
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:628
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:676
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:108
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
Definition: map.h:2826
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2760
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:361
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
#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:474
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:461
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205
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...