LORENE
tbl_val_smooth.C
1 /*
2  * Smoothes the junction with an eventual atmosphere.
3  *
4  */
5 
6 /*
7  * Copyright (c) 2004 Jerome Novak
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 char tbl_val_smooth_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $" ;
27 
28 
29 /*
30  * $Id: tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $
31  * $Log: tbl_val_smooth.C,v $
32  * Revision 1.4 2014/10/13 08:53:49 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.3 2004/12/30 16:14:01 j_novak
36  * Changed the name of a shadowed variable.
37  *
38  * Revision 1.2 2004/12/03 13:24:01 j_novak
39  * Minor modif.
40  *
41  * Revision 1.1 2004/11/26 17:02:19 j_novak
42  * Added a function giving a smooth transition to the atmosphere.
43  *
44  *
45  * $Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $
46  *
47  */
48 
49 // Lorene headers
50 #include "tbl_val.h"
51 
52 //Local prototypes
53 namespace Lorene {
54 void radial_smoothing(double* , const double* , int , double) ;
55 //****************************************************************************
56 
57 void Tbl_val::smooth_atmosphere(double atmosphere_thr) {
58 
59  const Gval_spher* gspher = dynamic_cast<const Gval_spher*>(gval) ;
60  assert(gspher != 0x0) ;
61  int ndim = gspher->get_ndim() ;
62  int fant = gspher->get_fantome() ;
63  int nr = get_dim(0) + 2*fant;
64 
65  switch (ndim) {
66  case 1: {
67  radial_smoothing(t, gspher->zr->t, nr, atmosphere_thr) ;
68  break ;
69  }
70  case 2: {
71  int nt = get_dim(1) + 2*fant ;
72  for (int j=0; j<nt; j++)
73  radial_smoothing(t+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
74  break ;
75  }
76  case 3: {
77  int nt = get_dim(1) + 2*fant ;
78  int np = get_dim(2) + 2*fant ;
79  for (int j=0; j<nt; j++)
80  for (int k=0; k<np; k++)
81  radial_smoothing(t+k*nt*nr+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
82  break ;
83  }
84 
85  default: {
86  cerr << "Tbl_val::smooth_atmosphere : strange number of dimensions!"
87  << endl ;
88  abort() ;
89  break ;
90  }
91  }
92  return ;
93 }
94 
95 void radial_smoothing(double* tab, const double* rr, int n, double rho) {
96 
97  assert((tab!= 0x0)&&(rr!=0x0)) ;
98  assert (rho >= 0.) ;
99 
100  if (fabs(tab[n-1]) > rho) // no atmosphere here
101  return ;
102 
103  double* t = tab + (n-1) ;
104  int indice = -1 ;
105  bool atmos = true ;
106  bool jump = false ;
107  for (int i=0; ((i<n)&&(atmos)); i++) {
108  if (atmos) atmos = ( fabs(*t) < rho) ;
109  t-- ;
110  if (atmos) {
111  jump = ( fabs(*t) > rho ) ;
112  if (jump) // discontinuity found
113  indice = n - i - 2 ;
114  }
115  }
116  if (indice == -1) return ;
117  int np = 2*(n-indice-2)/3 ;
118  int nm = indice / 100 + 3 ;
119  assert(n > nm+np) ;
120  if (indice < n - np+1) { // enough points to interpolate
121 
122  // The inteprolation is done using a cubic polynomial
123  //---------------------------------------------------
124 
125  int ileft = indice - nm + 2 ;
126  int iright = indice + np - 1 ;
127  double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
128  ( rr[ileft -1] - rr[ileft]) ;
129  double der_l = ( alpha*(alpha+2.)*tab[ileft]
130  - (1.+alpha)*(1.+alpha)*tab[ileft-1]
131  + tab[ileft-2] ) /
132  ( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
133  double f_l = tab[ileft] ;
134  double f_r = tab[iright] ;
135  double tau = rr[ileft] - rr[iright] ;
136  double alp = der_l / (tau*tau) + 2.*(f_r - f_l)/(tau*tau*tau) ;
137  double bet = 0.5*(der_l + alp*tau*(rr[iright] - 3*rr[ileft])) / tau ;
138  for (int i=ileft; i<iright; i++) {
139  tab[i] = f_r + (alp*rr[i]+bet)*(rr[i] - rr[iright])*(rr[i] - rr[iright]);
140  }
141  }
142  else { // too few points to interpolate -> linear extrapolation ...
143  int ileft = indice ;
144  double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
145  ( rr[ileft -1] - rr[ileft]) ;
146  double der_l = ( alpha*(alpha+2.)*tab[ileft]
147  - (1.+alpha)*(1.+alpha)*tab[ileft-1]
148  + tab[ileft-2] ) /
149  ( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
150  for (int i=ileft; i<n; i++) {
151  tab[i] = tab[ileft] + (rr[i] - rr[ileft])*der_l ;
152  }
153  }
154  return ;
155 }
156 }
Lorene prototypes.
Definition: app_hor.h:64
int get_dim(int i) const
Gives the i th dimension (ie dim->dim[i] , without hidden cells)
Definition: tbl_val.h:485
double * t
The array of double at the nodes.
Definition: tbl_val.h:114
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
Definition: tbl_val.h:110