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
26char 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
53namespace Lorene {
54void radial_smoothing(double* , const double* , int , double) ;
55//****************************************************************************
56
57void 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
95void 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}
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
Definition tbl_val.h:110
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
Lorene prototypes.
Definition app_hor.h:64