LORENE
interpol_bifluid.C
1 /*
2  * Hermite interpolation functions for 2-fluid EoS.
3  *
4  */
5 
6 /*
7  * Copyright (c) 2014 Aurelien Sourie
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 as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 
29 
30 /*
31  * $Id: interpol_bifluid.C,v 1.2 2016/12/05 16:18:11 j_novak Exp $
32  * $Log: interpol_bifluid.C,v $
33  * Revision 1.2 2016/12/05 16:18:11 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.1 2015/06/15 15:08:22 j_novak
37  * New file interpol_bifluid for interpolation of 2-fluid EoSs
38  *
39  *
40  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/interpol_bifluid.C,v 1.2 2016/12/05 16:18:11 j_novak Exp $
41  *
42  */
43 
44 #include<cmath>
45 
46 // Headers Lorene
47 #include "tbl.h"
48 namespace Lorene {
49 
50 void interpol_herm(const Tbl&, const Tbl&, const Tbl&, double, int&, double&,
51  double&) ;
52 
53 /*interpolation for functions of 3 variables : hermite interpolation on the 2 last variables and linear interpolation on the first one*/
54 
55 /*
56 xtab/x refers to delta_car (logdelta_car)
57 ytab/y refers to mu_n (logent1)
58 ztab/z refers to mu_p (logent2)
59 ftab/f refers to psi (logp) or alpha (dlpsddelta_car or logalpha)
60 dfdytab/dfdy refers to dpsi/dmu_n (dlpsdlent1) or dalpha/dmu_n (d2lpsdlent1ddelta_car or dlalphadlent1)
61 dfdztab/dfdz refers to dpsi/dmu_p (dlpsdlent2) or dalpha/dmu_p (dl2psdlent2ddelta_car or dlalphadlent2)
62 d2fdytabdztab refers to d2psi/dmu_ndmu_p (d2lpsdlent1dlent2) or d2alpha/dmu_ndmu_p (d3lpsdlent1dlent2ddelta_car or d2lalphadlent1dlent2)
63 */
64 
65 /*this routine provides the interpolated values of f, dfdy and dfdz at point (x,y,z) via the use of adapted tables*/
66 
67  void interpol_mixed_3d(const Tbl& xtab, const Tbl& ytab, const Tbl& ztab, const Tbl& ftab,
68  const Tbl& dfdytab, const Tbl& dfdztab, const Tbl& d2fdydztab,
69  double x, double y, double z, double& f, double& dfdy, double& dfdz)
70 {
71  assert(ytab.dim == xtab.dim) ;
72  assert(ztab.dim == xtab.dim) ;
73  assert(ftab.dim == xtab.dim) ;
74  assert(dfdytab.dim == xtab.dim) ;
75  assert(dfdztab.dim == xtab.dim) ;
76  assert(d2fdydztab.dim == xtab.dim) ;
77 
78  int nbp1, nbp2, nbp3;
79  nbp1 = xtab.get_dim(2) ; // \Delta^{2}
80  nbp2 = xtab.get_dim(1) ; // \mu_n
81  nbp3 = xtab.get_dim(0) ; // \mu_p
82 
83 
84  /* we can put these tests directly in the code (before calling interpol_mixed_3d)
85  assert(x >= xtab(0,0,0)) ;
86  assert(x <= xtab(nbp1,0,0)) ;
87  assert(y >= ytab(0,0,0)) ;
88  assert(y <= ytab(0,nbp2,0)) ;
89  assert(z >= ztab(0,0,0)) ;
90  assert(z <= ztab(0,0,nbp3)) ;
91  */
92 
93  int i_near = 0 ;
94  int j_near = 0 ;
95  int k_near = 0 ;
96 
97 
98  // look for the positions of (x,y,z) in the tables
99  while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
100  i_near++ ;
101  }
102  if (i_near != 0) {
103  i_near -- ;
104  }
105 
106  while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
107  j_near++ ;
108  }
109  if (j_near != 0) {
110  j_near -- ;
111  }
112 
113  while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
114  k_near++ ;
115  }
116  if (k_near != 0) {
117  k_near-- ;
118  }
119 
120  int i1 = i_near + 1 ;
121  int j1 = j_near + 1 ;
122  int k1 = k_near + 1 ;
123 
124 /*
125 bool hum1 = (( ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) ) == (ytab(i_near, j1, k1) - ytab(i_near, j_near, k1) ) ) ;
126 bool hum2 = ((ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) ) == ( ztab(i_near, j1, k1) - ztab(i_near, j1, k_near) )) ;
127 if (hum1 == false ) {
128 cout << setprecision(16);
129 cout << ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) << " " << ytab(i_near, j1, k1) - ytab(i_near, j_near, k1)<< endl ;
130 cout << j_near << " " << k_near << endl ;
131 }
132 
133 if (hum2 == false ) {
134 cout << setprecision(16);
135 cout << ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) << " " << ztab(i_near, j1, k1) - ztab(i_near, j1, k_near)<< endl ;
136 cout << j_near << " " << k_near << endl ;
137 }*/
138 
139  double dy_i_near = ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) ;
140  double dz_i_near = ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) ;
141  double u_i_near = (y - ytab(i_near, j_near, k_near)) / dy_i_near ;
142  double v_i_near = (z - ztab(i_near, j_near, k_near)) / dz_i_near ;
143 /*
144 //lineaire
145  double dy_anal = 1. ; // log(2.) ;
146  double dz_anal = 1. ; //log(2.) ;
147 //en log
148 cout << y << " " << z << endl ;
149  double u_i_anal = y /log(2.) ; //(y - 1.) / 1.
150  double v_i_anal = z /log(2.) ; //(z - 1.) / 1.
151 
152 cout << setprecision(16);
153 cout <<"interpol " << u_i_near <<" " << u_i_anal<< " "<< v_i_near << " " << v_i_anal << endl ;
154 */
155  double u2_i_near = u_i_near * u_i_near ;
156  double v2_i_near = v_i_near * v_i_near ;
157  double u3_i_near = u2_i_near * u_i_near ;
158  double v3_i_near = v2_i_near * v_i_near ;
159 
160  /*
161 //lineaire
162 double u2_i_anal = (y - 1. ) * (y - 1. ) ; //
163  double v2_i_anal = (z - 1.) * (z - 1.) ; //
164  double u3_i_anal = (y - 1. ) * (y - 1. )* (y - 1. ) ; //
165  double v3_i_anal = (z - 1.) * (z - 1.)* (z - 1.) ; //
166 */
167 
168  double psi0_u_i_near = double(2)*u3_i_near - double(3)*u2_i_near + double(1) ;
169  double psi0_1mu_i_near = -double(2)*u3_i_near + double(3)*u2_i_near ;
170  double psi1_u_i_near = u3_i_near - double(2)*u2_i_near + u_i_near ;
171  double psi1_1mu_i_near = -u3_i_near + u2_i_near ;
172  double psi0_v_i_near = double(2)*v3_i_near - double(3)*v2_i_near + double(1) ;
173  double psi0_1mv_i_near = -double(2)*v3_i_near + double(3)*v2_i_near ;
174  double psi1_v_i_near = v3_i_near - double(2)*v2_i_near + v_i_near ;
175  double psi1_1mv_i_near = -v3_i_near + v2_i_near ;
176 
177  double f_i_near ;
178 /*// lineaire
179  double psi0_u_i_anal = double(2)*(y - 1. ) * (y - 1. )* (y - 1. ) - double(3)*(y - 1. ) *(y - 1. ) + double(1) ;
180  double psi0_1mu_i_anal = -double(2)* (y - 1. ) * (y - 1. )* (y - 1. ) + double(3)*(y - 1. ) * (y - 1. );
181  double psi1_u_i_anal = (y - 1. ) * (y - 1. )* (y - 1. ) - double(2)*(y - 1. ) * (y - 1.) + (y - 1. ) ;
182  double psi1_1mu_i_anal = -(y - 1. ) * (y - 1. )* (y - 1. ) +(y - 1. ) * (y - 1. ) ;
183 //en log
184  double psi0_u_i_anal = double(2)*y /log(2.) * y /log(2.) * y /log(2.) - double(3)* y /log(2.) * y /log(2.) + double(1) ;
185  double psi0_1mu_i_anal = -double(2)*y /log(2.)*y /log(2.)*y /log(2.) + double(3)*y /log(2.)*y /log(2.);
186  double psi1_u_i_anal = y /log(2.)*y /log(2.)*y /log(2.) - double(2)* y /log(2.) * y /log(2.) + y /log(2.);
187  double psi1_1mu_i_anal = -y /log(2.)*y /log(2.)*y /log(2.) +y /log(2.)*y /log(2.) ;
188 // lineaire
189  double psi0_v_i_anal = double(2)*(z - 1. ) * (z - 1. )* (z - 1. ) - double(3)*(z - 1. ) *(z - 1. ) + double(1) ;
190  double psi0_1mv_i_anal = -double(2)* (z - 1. ) * (z - 1. )* (z - 1. ) + double(3)*(z - 1. ) * (z - 1. );
191  double psi1_v_i_anal = (z - 1. ) * (z - 1. )* (z - 1. ) - double(2)*(z - 1. ) * (z - 1.) + (z - 1. ) ;
192  double psi1_1mv_i_anal = -(z - 1. ) * (z - 1. )* (z - 1. ) +(z - 1. ) * (z - 1. ) ;
193 
194 */
195 
196  f_i_near = ftab(i_near, j_near, k_near) * psi0_u_i_near * psi0_v_i_near
197  + ftab(i_near, j1, k_near) * psi0_1mu_i_near * psi0_v_i_near
198  + ftab(i_near, j_near, k1) * psi0_u_i_near * psi0_1mv_i_near
199  + ftab(i_near, j1, k1) * psi0_1mu_i_near * psi0_1mv_i_near ;
200 
201  f_i_near += (dfdytab(i_near, j_near, k_near) * psi1_u_i_near * psi0_v_i_near
202  - dfdytab(i_near, j1, k_near) * psi1_1mu_i_near * psi0_v_i_near
203  + dfdytab(i_near, j_near, k1) * psi1_u_i_near * psi0_1mv_i_near
204  - dfdytab(i_near, j1, k1) * psi1_1mu_i_near * psi0_1mv_i_near) * dy_i_near ;
205 
206  f_i_near += (dfdztab(i_near, j_near, k_near) * psi0_u_i_near * psi1_v_i_near
207  + dfdztab(i_near, j1, k_near) * psi0_1mu_i_near * psi1_v_i_near
208  - dfdztab(i_near, j_near, k1) * psi0_u_i_near * psi1_1mv_i_near
209  - dfdztab(i_near, j1, k1) * psi0_1mu_i_near * psi1_1mv_i_near) * dz_i_near ;
210 
211  f_i_near += (d2fdydztab(i_near, j_near, k_near) * psi1_u_i_near * psi1_v_i_near
212  - d2fdydztab(i_near, j1, k_near) * psi1_1mu_i_near * psi1_v_i_near
213  - d2fdydztab(i_near, j_near, k1) * psi1_u_i_near * psi1_1mv_i_near
214  + d2fdydztab(i_near, j1, k1) * psi1_1mu_i_near * psi1_1mv_i_near) * dy_i_near * dz_i_near ;
215 
216  double dpsi0_u_i_near = 6.*(u2_i_near - u_i_near) ;
217  double dpsi0_1mu_i_near = 6.*(u2_i_near - u_i_near) ;
218  double dpsi1_u_i_near = 3.*u2_i_near - 4.*u_i_near + 1. ;
219  double dpsi1_1mu_i_near = 3.*u2_i_near - 2.*u_i_near ;
220 
221  double dfdy_i_near;
222 
223  dfdy_i_near = (ftab(i_near, j_near, k_near) * dpsi0_u_i_near * psi0_v_i_near
224  - ftab(i_near, j1, k_near) * dpsi0_1mu_i_near * psi0_v_i_near
225  + ftab(i_near, j_near, k1) * dpsi0_u_i_near * psi0_1mv_i_near
226  - ftab(i_near, j1, k1) * dpsi0_1mu_i_near * psi0_1mv_i_near ) / dy_i_near;
227 
228  dfdy_i_near += (dfdytab(i_near, j_near, k_near) * dpsi1_u_i_near * psi0_v_i_near
229  + dfdytab(i_near, j1, k_near) * dpsi1_1mu_i_near * psi0_v_i_near
230  + dfdytab(i_near, j_near, k1) * dpsi1_u_i_near * psi0_1mv_i_near
231  + dfdytab(i_near, j1, k1) * dpsi1_1mu_i_near * psi0_1mv_i_near) ;
232 
233  dfdy_i_near += (dfdztab(i_near, j_near, k_near) * dpsi0_u_i_near * psi1_v_i_near
234  - dfdztab(i_near, j1, k_near) * dpsi0_1mu_i_near * psi1_v_i_near
235  - dfdztab(i_near, j_near, k1) * dpsi0_u_i_near * psi1_1mv_i_near
236  + dfdztab(i_near, j1, k1) * dpsi0_1mu_i_near * psi1_1mv_i_near) * dz_i_near /dy_i_near ;
237 
238  dfdy_i_near += (d2fdydztab(i_near, j_near, k_near) * dpsi1_u_i_near * psi1_v_i_near
239  + d2fdydztab(i_near, j1, k_near) * dpsi1_1mu_i_near * psi1_v_i_near
240  - d2fdydztab(i_near, j_near, k1) * dpsi1_u_i_near * psi1_1mv_i_near
241  - d2fdydztab(i_near, j1, k1) * dpsi1_1mu_i_near * psi1_1mv_i_near) * dz_i_near ;
242 
243  double dpsi0_v_i_near = 6.*(v2_i_near - v_i_near) ;
244  double dpsi0_1mv_i_near = 6.*(v2_i_near - v_i_near) ;
245  double dpsi1_v_i_near= 3.*v2_i_near - 4.*v_i_near + 1. ;
246  double dpsi1_1mv_i_near= 3.*v2_i_near - 2.*v_i_near ;
247 
248  double dfdz_i_near;
249 
250  dfdz_i_near = (ftab(i_near, j_near, k_near) * psi0_u_i_near * dpsi0_v_i_near
251  + ftab(i_near, j1, k_near) * psi0_1mu_i_near * dpsi0_v_i_near
252  - ftab(i_near, j_near, k1) * psi0_u_i_near * dpsi0_1mv_i_near
253  - ftab(i_near, j1, k1) * psi0_1mu_i_near * dpsi0_1mv_i_near) / dz_i_near ;
254 
255  dfdz_i_near += (dfdytab(i_near, j_near, k_near) * psi1_u_i_near * dpsi0_v_i_near
256  - dfdytab(i_near, j1, k_near) * psi1_1mu_i_near * dpsi0_v_i_near
257  - dfdytab(i_near, j_near, k1) * psi1_u_i_near * dpsi0_1mv_i_near
258  + dfdytab(i_near, j1, k1) * psi1_1mu_i_near * dpsi0_1mv_i_near) * dy_i_near / dz_i_near ;
259 
260  dfdz_i_near += (dfdztab(i_near, j_near, k_near) * psi0_u_i_near * dpsi1_v_i_near
261  + dfdztab(i_near, j1, k_near) * psi0_1mu_i_near * dpsi1_v_i_near
262  + dfdztab(i_near, j_near, k1) * psi0_u_i_near * dpsi1_1mv_i_near
263  + dfdztab(i_near, j1, k1) * psi0_1mu_i_near * dpsi1_1mv_i_near) ;
264 
265  dfdz_i_near += (d2fdydztab(i_near, j_near, k_near) * psi1_u_i_near* dpsi1_v_i_near
266  - d2fdydztab(i_near, j1, k_near) * psi1_1mu_i_near * dpsi1_v_i_near
267  + d2fdydztab(i_near, j_near, k1) * psi1_u_i_near* dpsi1_1mv_i_near
268  - d2fdydztab(i_near, j1, k1) * psi1_1mu_i_near * dpsi1_1mv_i_near) * dy_i_near;
269 
270  // Hermite interpolation on the slice i1
271  // cette recherche est inutile car meme couverture du plan mu1, mu2 pour des delta differents
272  j_near = 0;
273  k_near = 0;
274  while ( ( ytab(i1,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
275  j_near++ ;
276  }
277  if (j_near != 0) {
278  j_near -- ;
279  }
280 
281  while ( ( ztab( i1, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
282  k_near++ ;
283  }
284  if (k_near != 0) {
285  k_near-- ;
286  }
287 
288  j1 = j_near + 1 ;
289  k1 = k_near + 1 ;
290 
291  double dy_i1 = ytab(i1, j1, k_near) - ytab(i1, j_near, k_near) ;
292  double dz_i1 = ztab(i1, j_near, k1) - ztab(i1, j_near, k_near) ;
293 
294  double u_i1 = (y - ytab(i1, j_near, k_near)) / dy_i1 ;
295  double v_i1 = (z - ztab(i1, j_near, k_near)) / dz_i1 ;
296 
297  double u2_i1 = u_i1 * u_i1 ;
298  double v2_i1 = v_i1 * v_i1 ;
299  double u3_i1 = u2_i1 * u_i1 ;
300  double v3_i1 = v2_i1 * v_i1 ;
301 
302  double psi0_u_i1 = 2.*u3_i1 - 3.*u2_i1 + 1. ;
303  double psi0_1mu_i1 = -2.*u3_i1 + 3.*u2_i1 ;
304  double psi1_u_i1 = u3_i1 - 2.*u2_i1 + u_i1 ;
305  double psi1_1mu_i1 = -u3_i1 + u2_i1 ;
306 
307  double psi0_v_i1 = 2.*v3_i1 - 3.*v2_i1 + 1. ;
308  double psi0_1mv_i1 = -2.*v3_i1 + 3.*v2_i1 ;
309  double psi1_v_i1 = v3_i1 - 2.*v2_i1 + v_i1 ;
310  double psi1_1mv_i1 = -v3_i1 + v2_i1 ;
311 
312  double f_i1;
313 
314  f_i1 = ftab(i1, j_near, k_near) * psi0_u_i1 * psi0_v_i1
315  + ftab(i1, j1, k_near) * psi0_1mu_i1 * psi0_v_i1
316  + ftab(i1, j_near, k1) * psi0_u_i1 * psi0_1mv_i1
317  + ftab(i1, j1, k1) * psi0_1mu_i1 * psi0_1mv_i1 ;
318 
319  f_i1 += (dfdytab(i1, j_near, k_near) * psi1_u_i1 * psi0_v_i1
320  - dfdytab(i1, j1, k_near) * psi1_1mu_i1 * psi0_v_i1
321  + dfdytab(i1, j_near, k1) * psi1_u_i1 * psi0_1mv_i1
322  - dfdytab(i1, j1, k1) * psi1_1mu_i1 * psi0_1mv_i1) * dy_i1 ;
323 
324  f_i1 += (dfdztab(i1, j_near, k_near) * psi0_u_i1 * psi1_v_i1
325  + dfdztab(i1, j1, k_near) * psi0_1mu_i1 * psi1_v_i1
326  - dfdztab(i1, j_near, k1) * psi0_u_i1 * psi1_1mv_i1
327  - dfdztab(i1, j1, k1) * psi0_1mu_i1 * psi1_1mv_i1) * dz_i1 ;
328 
329  f_i1 += (d2fdydztab(i1, j_near, k_near) * psi1_u_i1 * psi1_v_i1
330  - d2fdydztab(i1, j1, k_near) * psi1_1mu_i1 * psi1_v_i1
331  - d2fdydztab(i1, j_near, k1) * psi1_u_i1 * psi1_1mv_i1
332  + d2fdydztab(i1, j1, k1) * psi1_1mu_i1 * psi1_1mv_i1) * dy_i1 * dz_i1 ;
333 
334  double dpsi0_u_i1 = 6.*(u2_i1 - u_i1) ;
335  double dpsi0_1mu_i1 = 6.*(u2_i1 - u_i1) ;
336  double dpsi1_u_i1 = 3.*u2_i1 - 4.*u_i1 + 1. ;
337  double dpsi1_1mu_i1 = 3.*u2_i1 - 2.*u_i1 ;
338 
339  double dfdy_i1;
340 
341  dfdy_i1 = (ftab(i1, j_near, k_near) * dpsi0_u_i1 * psi0_v_i1
342  - ftab(i1, j1, k_near) * dpsi0_1mu_i1 * psi0_v_i1
343  + ftab(i1, j_near, k1) * dpsi0_u_i1 * psi0_1mv_i1
344  - ftab(i1, j1, k1) * dpsi0_1mu_i1 * psi0_1mv_i1 ) / dy_i1;
345 
346  dfdy_i1 += (dfdytab(i1, j_near, k_near) * dpsi1_u_i1 * psi0_v_i1
347  + dfdytab(i1, j1, k_near) * dpsi1_1mu_i1 * psi0_v_i1
348  + dfdytab(i1, j_near, k1) * dpsi1_u_i1 * psi0_1mv_i1
349  + dfdytab(i1, j1, k1) * dpsi1_1mu_i1 * psi0_1mv_i1) ;
350 
351  dfdy_i1 += (dfdztab(i1, j_near, k_near) * dpsi0_u_i1 * psi1_v_i1
352  - dfdztab(i1, j1, k_near) * dpsi0_1mu_i1 * psi1_v_i1
353  - dfdztab(i1, j_near, k1) * dpsi0_u_i1 * psi1_1mv_i1
354  + dfdztab(i1, j1, k1) * dpsi0_1mu_i1 * psi1_1mv_i1) * dz_i1 /dy_i1 ;
355 
356  dfdy_i1 += (d2fdydztab(i1, j_near, k_near) * dpsi1_u_i1 * psi1_v_i1
357  + d2fdydztab(i1, j1, k_near) * dpsi1_1mu_i1 * psi1_v_i1
358  - d2fdydztab(i1, j_near, k1) * dpsi1_u_i1 * psi1_1mv_i1
359  - d2fdydztab(i1, j1, k1) * dpsi1_1mu_i1 * psi1_1mv_i1) * dz_i1 ;
360 
361  double dpsi0_v_i1 = 6.*(v2_i1 - v_i1) ;
362  double dpsi0_1mv_i1 = 6.*(v2_i1 - v_i1) ;
363  double dpsi1_v_i1= 3.*v2_i1 - 4.*v_i1 + 1. ;
364  double dpsi1_1mv_i1= 3.*v2_i1 - 2.*v_i1 ;
365 
366  double dfdz_i1;
367 
368  dfdz_i1 = (ftab(i1, j_near, k_near) * psi0_u_i1 * dpsi0_v_i1
369  + ftab(i1, j1, k_near) * psi0_1mu_i1 * dpsi0_v_i1
370  - ftab(i1, j_near, k1) * psi0_u_i1 * dpsi0_1mv_i1
371  - ftab(i1, j1, k1) * psi0_1mu_i1 * dpsi0_1mv_i1) / dz_i1 ;
372 
373  dfdz_i1 += (dfdytab(i1, j_near, k_near) * psi1_u_i1 * dpsi0_v_i1
374  - dfdytab(i1, j1, k_near) * psi1_1mu_i1 * dpsi0_v_i1
375  - dfdytab(i1, j_near, k1) * psi1_u_i1 * dpsi0_1mv_i1
376  + dfdytab(i1, j1, k1) * psi1_1mu_i1 * dpsi0_1mv_i1) * dy_i1 / dz_i1 ;
377 
378  dfdz_i1 += (dfdztab(i1, j_near, k_near) * psi0_u_i1 * dpsi1_v_i1
379  + dfdztab(i1, j1, k_near) * psi0_1mu_i1 * dpsi1_v_i1
380  + dfdztab(i1, j_near, k1) * psi0_u_i1 * dpsi1_1mv_i1
381  + dfdztab(i1, j1, k1) * psi0_1mu_i1 * dpsi1_1mv_i1) ;
382 
383  dfdz_i1 += (d2fdydztab(i1, j_near, k_near) * psi1_u_i1* dpsi1_v_i1
384  - d2fdydztab(i1, j1, k_near) * psi1_1mu_i1 * dpsi1_v_i1
385  + d2fdydztab(i1, j_near, k1) * psi1_u_i1* dpsi1_1mv_i1
386  - d2fdydztab(i1, j1, k1) * psi1_1mu_i1 * dpsi1_1mv_i1) * dy_i1;
387 
388  /* linear interpolation on the first variable */
389 
390  double x1 = xtab(i_near, 0, 0) ;
391  double x2 = xtab(i1, 0, 0) ;
392 
393  double x12 = x1-x2 ;
394 
395  //for f
396  double y1 = f_i_near;
397  double y2 = f_i1;
398  double a = (y1-y2)/x12 ;
399  double b = (x1*y2-y1*x2)/x12 ;
400 
401  f = x*a+b ;
402  /*cout << " x1 " << x1 << endl ;
403  cout << " x2 " << x2 << endl ;
404  cout << " x " << x << endl ;
405  cout << " y1 " << y1 << endl ;
406  cout << " y2 " << y2 << endl ;
407  cout << " y " <<f << endl ;*/
408  //for df/dy
409  double y1_y = dfdy_i_near;
410  double y2_y = dfdy_i1;
411  double a_y = (y1_y-y2_y)/x12 ;
412  double b_y = (x1*y2_y-y1_y*x2)/x12 ;
413 
414  dfdy = x*a_y+b_y ;
415 
416  //for df/dz
417  double y1_z = dfdz_i_near;
418  double y2_z = dfdz_i1;
419  double a_z = (y1_z-y2_z)/x12 ;
420  double b_z = (x1*y2_z-y1_z*x2)/x12 ;
421 
422  dfdz = x*a_z+b_z ;
423  /*cout << "i_near " << i_near << endl;
424  cout << "j_near " << j_near << endl;
425  cout << "k_near " << k_near << endl;
426  abort () ;*/
427  //cout << dfdy_i_near << " " <<dfdy_i1 << " " << dfdz_i_near << " " << dfdz_i1<< endl;
428  return;
429 
430 }
431 
432 
433 
434 void interpol_mixed_3d_mod(const Tbl& xtab, const Tbl& ytab, const Tbl& ztab, const Tbl& ftab,
435  const Tbl& dfdytab, const Tbl& dfdztab,
436  double x, double y, double z, double& f, double& dfdy, double& dfdz)
437 {
438  assert(ytab.dim == xtab.dim) ;
439  assert(ztab.dim == xtab.dim) ;
440  assert(ftab.dim == xtab.dim) ;
441  assert(dfdytab.dim == xtab.dim) ;
442  assert(dfdztab.dim == xtab.dim) ;
443 
444  int nbp1, nbp2, nbp3;
445  nbp1 = xtab.get_dim(2) ; // \Delta^{2}
446  nbp2 = xtab.get_dim(1) ; // \mu_n
447  nbp3 = xtab.get_dim(0) ; // \mu_p
448 
449 
450  /* we can put these tests directly in the code (before calling interpol_mixed_3d)
451  assert(x >= xtab(0,0,0)) ;
452  assert(x <= xtab(nbp1,0,0)) ;
453  assert(y >= ytab(0,0,0)) ;
454  assert(y <= ytab(0,nbp2,0)) ;
455  assert(z >= ztab(0,0,0)) ;
456  assert(z <= ztab(0,0,nbp3)) ;
457  */
458 
459  int i_near = 0 ;
460  int j_near = 0 ;
461  int k_near = 0 ;
462 
463  // look for the positions of (x,y,z) in the tables
464  while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
465  i_near++ ;
466  }
467  if (i_near != 0) {
468  i_near -- ;
469  }
470 
471  while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
472  j_near++ ;
473  }
474  if (j_near != 0) {
475  j_near -- ;
476  }
477 
478  while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
479  k_near++ ;
480  }
481  if (k_near != 0) {
482  k_near-- ;
483  }
484 
485  int i1 = i_near + 1 ;
486  int j1 = j_near + 1 ;
487  int k1 = k_near + 1 ;
488 
489 /*
490  bool hum1 = (( ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) ) == (ytab(i_near, j1, k1) - ytab(i_near, j_near, k1) ) ) ;
491  bool hum2 = ((ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) ) == ( ztab(i_near, j1, k1) - ztab(i_near, j1, k_near) )) ;
492  if (hum1 == false ) {
493  cout << setprecision(16);
494  cout << ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) << " " << ytab(i_near, j1, k1) - ytab(i_near, j_near, k1)<< endl ;
495  cout << j_near << " " << k_near << endl ;
496 }
497 
498 if (hum2 == false ) {
499 cout << setprecision(16);
500 cout << ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) << " " << ztab(i_near, j1, k1) - ztab(i_near, j1, k_near)<< endl ;
501 cout << j_near << " " << k_near << endl ;
502 }*/
503 
504  double dy_i_near = ytab(i_near, j1, k_near) - ytab(i_near, j_near, k_near) ;
505  double dz_i_near = ztab(i_near, j_near, k1) - ztab(i_near, j_near, k_near) ;
506  double u_i_near = (y - ytab(i_near, j_near, k_near)) / dy_i_near ;
507  double v_i_near = (z - ztab(i_near, j_near, k_near)) / dz_i_near ;
508  /*
509  //lineaire
510  double dy_anal = 1. ; // log(2.) ;
511  double dz_anal = 1. ; //log(2.) ;
512  //en log
513 cout << y << " " << z << endl ;
514 double u_i_anal = y /log(2.) ; //(y - 1.) / 1.
515 double v_i_anal = z /log(2.) ; //(z - 1.) / 1.
516 
517 cout << setprecision(16);
518 cout <<"interpol " << u_i_near <<" " << u_i_anal<< " "<< v_i_near << " " << v_i_anal << endl ;
519  */
520 
521  double u2_i_near = u_i_near * u_i_near ;
522  double v2_i_near = v_i_near * v_i_near ;
523  double u3_i_near = u2_i_near * u_i_near ;
524  double v3_i_near = v2_i_near * v_i_near ;
525 
526  /*
527 //lineaire
528 double u2_i_anal = (y - 1. ) * (y - 1. ) ; //
529  double v2_i_anal = (z - 1.) * (z - 1.) ; //
530  double u3_i_anal = (y - 1. ) * (y - 1. )* (y - 1. ) ; //
531  double v3_i_anal = (z - 1.) * (z - 1.)* (z - 1.) ; //
532 
533 cout << " interpol2" << u2_i_near << " " << u2_i_anal << " " << u3_i_near << " " << u3_i_anal << endl ;
534 */
535 
536  double psi0_u_i_near = double(2)*u3_i_near - double(3)*u2_i_near + double(1) ;
537  double psi0_1mu_i_near = -double(2)*u3_i_near + double(3)*u2_i_near ;
538  double psi1_u_i_near = u3_i_near - double(2)*u2_i_near + u_i_near ;
539  double psi1_1mu_i_near = -u3_i_near + u2_i_near ;
540  double psi0_v_i_near = double(2)*v3_i_near - double(3)*v2_i_near + double(1) ;
541  double psi0_1mv_i_near = -double(2)*v3_i_near + double(3)*v2_i_near ;
542  double psi1_v_i_near = v3_i_near - double(2)*v2_i_near + v_i_near ;
543  double psi1_1mv_i_near = -v3_i_near + v2_i_near ;
544 
545  double f_i_near ;
546 /*// lineaire
547  double psi0_u_i_anal = double(2)*(y - 1. ) * (y - 1. )* (y - 1. ) - double(3)*(y - 1. ) *(y - 1. ) + double(1) ;
548  double psi0_1mu_i_anal = -double(2)* (y - 1. ) * (y - 1. )* (y - 1. ) + double(3)*(y - 1. ) * (y - 1. );
549  double psi1_u_i_anal = (y - 1. ) * (y - 1. )* (y - 1. ) - double(2)*(y - 1. ) * (y - 1.) + (y - 1. ) ;
550  double psi1_1mu_i_anal = -(y - 1. ) * (y - 1. )* (y - 1. ) +(y - 1. ) * (y - 1. ) ;
551 //en log
552  double psi0_u_i_anal = double(2)*y /log(2.) * y /log(2.) * y /log(2.) - double(3)* y /log(2.) * y /log(2.) + double(1) ;
553  double psi0_1mu_i_anal = -double(2)*y /log(2.)*y /log(2.)*y /log(2.) + double(3)*y /log(2.)*y /log(2.);
554  double psi1_u_i_anal = y /log(2.)*y /log(2.)*y /log(2.) - double(2)* y /log(2.) * y /log(2.) + y /log(2.);
555  double psi1_1mu_i_anal = -y /log(2.)*y /log(2.)*y /log(2.) +y /log(2.)*y /log(2.) ;
556 // lineaire
557  double psi0_v_i_anal = double(2)*(z - 1. ) * (z - 1. )* (z - 1. ) - double(3)*(z - 1. ) *(z - 1. ) + double(1) ;
558  double psi0_1mv_i_anal = -double(2)* (z - 1. ) * (z - 1. )* (z - 1. ) + double(3)*(z - 1. ) * (z - 1. );
559  double psi1_v_i_anal = (z - 1. ) * (z - 1. )* (z - 1. ) - double(2)*(z - 1. ) * (z - 1.) + (z - 1. ) ;
560  double psi1_1mv_i_anal = -(z - 1. ) * (z - 1. )* (z - 1. ) +(z - 1. ) * (z - 1. ) ;
561 
562 cout << " Interpol3 " << psi0_u_i_near << " " << psi0_u_i_anal << " " << psi0_1mu_i_near << " " << psi0_1mu_i_anal << " " <<
563 psi1_u_i_near << " " << psi1_u_i_anal << " " << psi1_1mu_i_near << " " << psi1_1mu_i_anal << endl;
564 
565 cout << " Interpol4 " << psi0_v_i_near << " " << psi0_v_i_anal << " " << psi0_1mv_i_near << " " << psi0_1mv_i_anal << " " <<
566 psi1_v_i_near << " " << psi1_v_i_anal << " " << psi1_1mv_i_near << " " << psi1_1mv_i_anal << endl;
567 */
568 
569  f_i_near = ftab(i_near, j_near, k_near) * psi0_u_i_near * psi0_v_i_near
570  + ftab(i_near, j1, k_near) * psi0_1mu_i_near * psi0_v_i_near
571  + ftab(i_near, j_near, k1) * psi0_u_i_near * psi0_1mv_i_near
572  + ftab(i_near, j1, k1) * psi0_1mu_i_near * psi0_1mv_i_near ;
573 
574  f_i_near += (dfdytab(i_near, j_near, k_near) * psi1_u_i_near * psi0_v_i_near
575  - dfdytab(i_near, j1, k_near) * psi1_1mu_i_near * psi0_v_i_near
576  + dfdytab(i_near, j_near, k1) * psi1_u_i_near * psi0_1mv_i_near
577  - dfdytab(i_near, j1, k1) * psi1_1mu_i_near * psi0_1mv_i_near) * dy_i_near ;
578 
579  f_i_near += (dfdztab(i_near, j_near, k_near) * psi0_u_i_near * psi1_v_i_near
580  + dfdztab(i_near, j1, k_near) * psi0_1mu_i_near * psi1_v_i_near
581  - dfdztab(i_near, j_near, k1) * psi0_u_i_near * psi1_1mv_i_near
582  - dfdztab(i_near, j1, k1) * psi0_1mu_i_near * psi1_1mv_i_near) * dz_i_near ;
583 
584  double dpsi0_u_i_near = 6.*(u2_i_near - u_i_near) ;
585  double dpsi0_1mu_i_near = 6.*(u2_i_near - u_i_near) ;
586  double dpsi1_u_i_near = 3.*u2_i_near - 4.*u_i_near + 1. ;
587  double dpsi1_1mu_i_near = 3.*u2_i_near - 2.*u_i_near ;
588 
589  double dfdy_i_near;
590 
591  dfdy_i_near = (ftab(i_near, j_near, k_near) * dpsi0_u_i_near * psi0_v_i_near
592  - ftab(i_near, j1, k_near) * dpsi0_1mu_i_near * psi0_v_i_near
593  + ftab(i_near, j_near, k1) * dpsi0_u_i_near * psi0_1mv_i_near
594  - ftab(i_near, j1, k1) * dpsi0_1mu_i_near * psi0_1mv_i_near ) / dy_i_near;
595 
596  dfdy_i_near += (dfdytab(i_near, j_near, k_near) * dpsi1_u_i_near * psi0_v_i_near
597  + dfdytab(i_near, j1, k_near) * dpsi1_1mu_i_near * psi0_v_i_near
598  + dfdytab(i_near, j_near, k1) * dpsi1_u_i_near * psi0_1mv_i_near
599  + dfdytab(i_near, j1, k1) * dpsi1_1mu_i_near * psi0_1mv_i_near) ;
600 
601  dfdy_i_near += (dfdztab(i_near, j_near, k_near) * dpsi0_u_i_near * psi1_v_i_near
602  - dfdztab(i_near, j1, k_near) * dpsi0_1mu_i_near * psi1_v_i_near
603  - dfdztab(i_near, j_near, k1) * dpsi0_u_i_near * psi1_1mv_i_near
604  + dfdztab(i_near, j1, k1) * dpsi0_1mu_i_near * psi1_1mv_i_near) * dz_i_near /dy_i_near ;
605 
606 
607 
608  double dpsi0_v_i_near = 6.*(v2_i_near - v_i_near) ;
609  double dpsi0_1mv_i_near = 6.*(v2_i_near - v_i_near) ;
610  double dpsi1_v_i_near= 3.*v2_i_near - 4.*v_i_near + 1. ;
611  double dpsi1_1mv_i_near= 3.*v2_i_near - 2.*v_i_near ;
612 
613  double dfdz_i_near;
614 
615  dfdz_i_near = (ftab(i_near, j_near, k_near) * psi0_u_i_near * dpsi0_v_i_near
616  + ftab(i_near, j1, k_near) * psi0_1mu_i_near * dpsi0_v_i_near
617  - ftab(i_near, j_near, k1) * psi0_u_i_near * dpsi0_1mv_i_near
618  - ftab(i_near, j1, k1) * psi0_1mu_i_near * dpsi0_1mv_i_near) / dz_i_near ;
619 
620  dfdz_i_near += (dfdytab(i_near, j_near, k_near) * psi1_u_i_near * dpsi0_v_i_near
621  - dfdytab(i_near, j1, k_near) * psi1_1mu_i_near * dpsi0_v_i_near
622  - dfdytab(i_near, j_near, k1) * psi1_u_i_near * dpsi0_1mv_i_near
623  + dfdytab(i_near, j1, k1) * psi1_1mu_i_near * dpsi0_1mv_i_near) * dy_i_near / dz_i_near ;
624 
625  dfdz_i_near += (dfdztab(i_near, j_near, k_near) * psi0_u_i_near * dpsi1_v_i_near
626  + dfdztab(i_near, j1, k_near) * psi0_1mu_i_near * dpsi1_v_i_near
627  + dfdztab(i_near, j_near, k1) * psi0_u_i_near * dpsi1_1mv_i_near
628  + dfdztab(i_near, j1, k1) * psi0_1mu_i_near * dpsi1_1mv_i_near) ;
629 
630  // Hermite interpolation on the slice i1
631  // cette recherche est inutile car meme couverture du plan mu1, mu2 pour des delta differents
632  j_near = 0;
633  k_near = 0;
634  while ( ( ytab(i1,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
635  j_near++ ;
636  }
637  if (j_near != 0) {
638  j_near -- ;
639  }
640 
641  while ( ( ztab( i1, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
642  k_near++ ;
643  }
644  if (k_near != 0) {
645  k_near-- ;
646  }
647 
648  j1 = j_near + 1 ;
649  k1 = k_near + 1 ;
650 
651  double dy_i1 = ytab(i1, j1, k_near) - ytab(i1, j_near, k_near) ;
652  double dz_i1 = ztab(i1, j_near, k1) - ztab(i1, j_near, k_near) ;
653 
654  double u_i1 = (y - ytab(i1, j_near, k_near)) / dy_i1 ;
655  double v_i1 = (z - ztab(i1, j_near, k_near)) / dz_i1 ;
656 
657  double u2_i1 = u_i1 * u_i1 ;
658  double v2_i1 = v_i1 * v_i1 ;
659  double u3_i1 = u2_i1 * u_i1 ;
660  double v3_i1 = v2_i1 * v_i1 ;
661 
662  double psi0_u_i1 = 2.*u3_i1 - 3.*u2_i1 + 1. ;
663  double psi0_1mu_i1 = -2.*u3_i1 + 3.*u2_i1 ;
664  double psi1_u_i1 = u3_i1 - 2.*u2_i1 + u_i1 ;
665  double psi1_1mu_i1 = -u3_i1 + u2_i1 ;
666 
667  double psi0_v_i1 = 2.*v3_i1 - 3.*v2_i1 + 1. ;
668  double psi0_1mv_i1 = -2.*v3_i1 + 3.*v2_i1 ;
669  double psi1_v_i1 = v3_i1 - 2.*v2_i1 + v_i1 ;
670  double psi1_1mv_i1 = -v3_i1 + v2_i1 ;
671 
672  double f_i1;
673 
674  f_i1 = ftab(i1, j_near, k_near) * psi0_u_i1 * psi0_v_i1
675  + ftab(i1, j1, k_near) * psi0_1mu_i1 * psi0_v_i1
676  + ftab(i1, j_near, k1) * psi0_u_i1 * psi0_1mv_i1
677  + ftab(i1, j1, k1) * psi0_1mu_i1 * psi0_1mv_i1 ;
678 
679  f_i1 += (dfdytab(i1, j_near, k_near) * psi1_u_i1 * psi0_v_i1
680  - dfdytab(i1, j1, k_near) * psi1_1mu_i1 * psi0_v_i1
681  + dfdytab(i1, j_near, k1) * psi1_u_i1 * psi0_1mv_i1
682  - dfdytab(i1, j1, k1) * psi1_1mu_i1 * psi0_1mv_i1) * dy_i1 ;
683 
684  f_i1 += (dfdztab(i1, j_near, k_near) * psi0_u_i1 * psi1_v_i1
685  + dfdztab(i1, j1, k_near) * psi0_1mu_i1 * psi1_v_i1
686  - dfdztab(i1, j_near, k1) * psi0_u_i1 * psi1_1mv_i1
687  - dfdztab(i1, j1, k1) * psi0_1mu_i1 * psi1_1mv_i1) * dz_i1 ;
688 
689  double dpsi0_u_i1 = 6.*(u2_i1 - u_i1) ;
690  double dpsi0_1mu_i1 = 6.*(u2_i1 - u_i1) ;
691  double dpsi1_u_i1 = 3.*u2_i1 - 4.*u_i1 + 1. ;
692  double dpsi1_1mu_i1 = 3.*u2_i1 - 2.*u_i1 ;
693 
694  double dfdy_i1;
695 
696  dfdy_i1 = (ftab(i1, j_near, k_near) * dpsi0_u_i1 * psi0_v_i1
697  - ftab(i1, j1, k_near) * dpsi0_1mu_i1 * psi0_v_i1
698  + ftab(i1, j_near, k1) * dpsi0_u_i1 * psi0_1mv_i1
699  - ftab(i1, j1, k1) * dpsi0_1mu_i1 * psi0_1mv_i1 ) / dy_i1;
700 
701  dfdy_i1 += (dfdytab(i1, j_near, k_near) * dpsi1_u_i1 * psi0_v_i1
702  + dfdytab(i1, j1, k_near) * dpsi1_1mu_i1 * psi0_v_i1
703  + dfdytab(i1, j_near, k1) * dpsi1_u_i1 * psi0_1mv_i1
704  + dfdytab(i1, j1, k1) * dpsi1_1mu_i1 * psi0_1mv_i1) ;
705 
706  dfdy_i1 += (dfdztab(i1, j_near, k_near) * dpsi0_u_i1 * psi1_v_i1
707  - dfdztab(i1, j1, k_near) * dpsi0_1mu_i1 * psi1_v_i1
708  - dfdztab(i1, j_near, k1) * dpsi0_u_i1 * psi1_1mv_i1
709  + dfdztab(i1, j1, k1) * dpsi0_1mu_i1 * psi1_1mv_i1) * dz_i1 /dy_i1 ;
710 
711  double dpsi0_v_i1 = 6.*(v2_i1 - v_i1) ;
712  double dpsi0_1mv_i1 = 6.*(v2_i1 - v_i1) ;
713  double dpsi1_v_i1= 3.*v2_i1 - 4.*v_i1 + 1. ;
714  double dpsi1_1mv_i1= 3.*v2_i1 - 2.*v_i1 ;
715 
716  double dfdz_i1;
717 
718  dfdz_i1 = (ftab(i1, j_near, k_near) * psi0_u_i1 * dpsi0_v_i1
719  + ftab(i1, j1, k_near) * psi0_1mu_i1 * dpsi0_v_i1
720  - ftab(i1, j_near, k1) * psi0_u_i1 * dpsi0_1mv_i1
721  - ftab(i1, j1, k1) * psi0_1mu_i1 * dpsi0_1mv_i1) / dz_i1 ;
722 
723  dfdz_i1 += (dfdytab(i1, j_near, k_near) * psi1_u_i1 * dpsi0_v_i1
724  - dfdytab(i1, j1, k_near) * psi1_1mu_i1 * dpsi0_v_i1
725  - dfdytab(i1, j_near, k1) * psi1_u_i1 * dpsi0_1mv_i1
726  + dfdytab(i1, j1, k1) * psi1_1mu_i1 * dpsi0_1mv_i1) * dy_i1 / dz_i1 ;
727 
728  dfdz_i1 += (dfdztab(i1, j_near, k_near) * psi0_u_i1 * dpsi1_v_i1
729  + dfdztab(i1, j1, k_near) * psi0_1mu_i1 * dpsi1_v_i1
730  + dfdztab(i1, j_near, k1) * psi0_u_i1 * dpsi1_1mv_i1
731  + dfdztab(i1, j1, k1) * psi0_1mu_i1 * dpsi1_1mv_i1) ;
732 
733  /* linear interpolation on the first variable */
734 
735  double x1 = xtab(i_near, 0, 0) ;
736  double x2 = xtab(i1, 0, 0) ;
737 
738  double x12 = x1-x2 ;
739 
740  //for f
741  double y1 = f_i_near;
742  double y2 = f_i1;
743  double a = (y1-y2)/x12 ;
744  double b = (x1*y2-y1*x2)/x12 ;
745 
746  f = x*a+b ;
747  // cout << " x1 " << x1 << " x2 " << x2 << " x " << x << " y1 " << y1 << " y2 " << y2 << " y " <<f << endl ;
748  //for df/dy
749  double y1_y = dfdy_i_near;
750  double y2_y = dfdy_i1;
751  double a_y = (y1_y-y2_y)/x12 ;
752  double b_y = (x1*y2_y-y1_y*x2)/x12 ;
753 
754  dfdy = x*a_y+b_y ;
755 
756  //for df/dz
757  double y1_z = dfdz_i_near;
758  double y2_z = dfdz_i1;
759  double a_z = (y1_z-y2_z)/x12 ;
760  double b_z = (x1*y2_z-y1_z*x2)/x12 ;
761 
762  dfdz = x*a_z+b_z ;
763  /*cout << "i_near " << i_near << endl;
764  cout << "j_near " << j_near << endl;
765  cout << "k_near " << k_near << endl;
766  abort () ;*/
767  return;
768 
769 }
770 
771 
772 void interpol_herm_2d_new_avec( double y, double z,
773  double mu1_11, double mu1_21, double mu2_11, double mu2_12,
774  double p_11, double p_21,double p_12, double p_22,
775  double n1_11, double n1_21, double n1_12,double n1_22,
776  double n2_11, double n2_21,double n2_12, double n2_22,
777  double cross_11, double cross_21, double cross_12, double cross_22,
778  double& f, double& dfdy, double& dfdz) {
779 
780 
781  double dy = mu1_21 - mu1_11 ;
782  double dz = mu2_12 - mu2_11;
783 
784  double u = (y - mu1_11) / dy ;
785  double v = (z - mu2_11) / dz ;
786 
787  double u2 = u*u ; double v2 = v*v ;
788  double u3 = u2*u ; double v3 = v2*v ;
789 
790  double psi0_u = 2.*u3 - 3.*u2 + 1. ;
791  double psi0_1mu = -2.*u3 + 3.*u2 ;
792  double psi1_u = u3 - 2.*u2 + u ;
793  double psi1_1mu = -u3 + u2 ;
794 
795  double psi0_v = 2.*v3 - 3.*v2 + 1. ;
796  double psi0_1mv = -2.*v3 + 3.*v2 ;
797  double psi1_v = v3 - 2.*v2 + v ;
798  double psi1_1mv = -v3 + v2 ;
799 
800  f = p_11 * psi0_u * psi0_v
801  + p_21 * psi0_1mu * psi0_v
802  + p_12 * psi0_u * psi0_1mv
803  + p_22 * psi0_1mu * psi0_1mv ;
804 
805  f += (n1_11 * psi1_u * psi0_v
806  - n1_21 * psi1_1mu * psi0_v
807  + n1_12* psi1_u * psi0_1mv
808  - n1_22 * psi1_1mu * psi0_1mv) * dy ;
809 
810  f += (n2_11 * psi0_u * psi1_v
811  + n2_21 * psi0_1mu * psi1_v
812  - n2_12 * psi0_u * psi1_1mv
813  - n2_22* psi0_1mu * psi1_1mv) * dz ;
814 
815  f += (cross_11 * psi1_u * psi1_v
816  - cross_21 * psi1_1mu * psi1_v
817  - cross_12 * psi1_u * psi1_1mv
818  + cross_22 * psi1_1mu * psi1_1mv) * dy * dz ;
819 
820  double dpsi0_u = 6.*(u2 - u) ;
821  double dpsi0_1mu = 6.*(u2 - u) ;
822  double dpsi1_u = 3.*u2 - 4.*u + 1. ;
823  double dpsi1_1mu = 3.*u2 - 2.*u ;
824 
825  dfdy = (p_11 * dpsi0_u * psi0_v
826  - p_21 * dpsi0_1mu * psi0_v
827  + p_12 * dpsi0_u * psi0_1mv
828  - p_22 * dpsi0_1mu * psi0_1mv ) / dy;
829 
830  dfdy += (n1_11* dpsi1_u * psi0_v
831  + n1_21 * dpsi1_1mu * psi0_v
832  + n1_12 * dpsi1_u * psi0_1mv
833  + n1_22 * dpsi1_1mu * psi0_1mv) ;
834 
835  dfdy += (n2_11 * dpsi0_u * psi1_v
836  - n2_21 * dpsi0_1mu * psi1_v
837  - n2_12 * dpsi0_u * psi1_1mv
838  + n2_22 * dpsi0_1mu * psi1_1mv) * dz /dy ;
839 
840  dfdy += (cross_11 * dpsi1_u * psi1_v
841  + cross_21* dpsi1_1mu * psi1_v
842  - cross_12 * dpsi1_u * psi1_1mv
843  - cross_22 * dpsi1_1mu * psi1_1mv) * dz ;
844 
845  double dpsi0_v = 6.*(v2 - v) ;
846  double dpsi0_1mv = 6.*(v2 - v) ;
847  double dpsi1_v = 3.*v2 - 4.*v + 1. ;
848  double dpsi1_1mv = 3.*v2 - 2.*v ;
849 
850  dfdz = (p_11* psi0_u * dpsi0_v
851  + p_21 * psi0_1mu * dpsi0_v
852  - p_12 * psi0_u * dpsi0_1mv
853  - p_22 * psi0_1mu * dpsi0_1mv) / dz ;
854 
855  dfdz += (n1_11 * psi1_u * dpsi0_v
856  - n1_21 * psi1_1mu * dpsi0_v
857  - n1_12 * psi1_u * dpsi0_1mv
858  + n1_22 * psi1_1mu * dpsi0_1mv) * dy / dz ;
859 
860  dfdz += (n2_11 * psi0_u * dpsi1_v
861  + n2_21 * psi0_1mu * dpsi1_v
862  + n2_12 * psi0_u * dpsi1_1mv
863  + n2_22 * psi0_1mu * dpsi1_1mv) ;
864 
865  dfdz += (cross_11 * psi1_u * dpsi1_v
866  - cross_21 * psi1_1mu * dpsi1_v
867  + cross_12 * psi1_u * dpsi1_1mv
868  - cross_22 * psi1_1mu * dpsi1_1mv) * dy ;
869 
870  return ;
871 }
872 
873 void interpol_herm_2d_new_sans( double y, double z,
874  double mu1_11, double mu1_21, double mu2_11, double mu2_12,
875  double p_11, double p_21,double p_12, double p_22,
876  double n1_11, double n1_21, double n1_12,double n1_22,
877  double n2_11, double n2_21,double n2_12, double n2_22,
878  double& f, double& dfdy, double& dfdz) {
879 
880  double dy = mu1_21 - mu1_11 ;
881  double dz = mu2_12 - mu2_11;
882 
883  double u = (y - mu1_11) / dy ;
884  double v = (z - mu2_11) / dz ;
885 
886  double u2 = u*u ; double v2 = v*v ;
887  double u3 = u2*u ; double v3 = v2*v ;
888 
889  double psi0_u = 2.*u3 - 3.*u2 + 1. ;
890  double psi0_1mu = -2.*u3 + 3.*u2 ;
891  double psi1_u = u3 - 2.*u2 + u ;
892  double psi1_1mu = -u3 + u2 ;
893 
894  double psi0_v = 2.*v3 - 3.*v2 + 1. ;
895  double psi0_1mv = -2.*v3 + 3.*v2 ;
896  double psi1_v = v3 - 2.*v2 + v ;
897  double psi1_1mv = -v3 + v2 ;
898 
899  f = p_11 * psi0_u * psi0_v
900  + p_21 * psi0_1mu * psi0_v
901  + p_12 * psi0_u * psi0_1mv
902  + p_22 * psi0_1mu * psi0_1mv ;
903 
904  f += (n1_11 * psi1_u * psi0_v
905  - n1_21 * psi1_1mu * psi0_v
906  + n1_12* psi1_u * psi0_1mv
907  - n1_22 * psi1_1mu * psi0_1mv) * dy ;
908 
909  f += (n2_11 * psi0_u * psi1_v
910  + n2_21 * psi0_1mu * psi1_v
911  - n2_12 * psi0_u * psi1_1mv
912  - n2_22* psi0_1mu * psi1_1mv) * dz ;
913 
914  double dpsi0_u = 6.*(u2 - u) ;
915  double dpsi0_1mu = 6.*(u2 - u) ;
916  double dpsi1_u = 3.*u2 - 4.*u + 1. ;
917  double dpsi1_1mu = 3.*u2 - 2.*u ;
918 
919  dfdy = (p_11 * dpsi0_u * psi0_v
920  - p_21 * dpsi0_1mu * psi0_v
921  + p_12 * dpsi0_u * psi0_1mv
922  - p_22 * dpsi0_1mu * psi0_1mv ) / dy;
923 
924  dfdy += (n1_11* dpsi1_u * psi0_v
925  + n1_21 * dpsi1_1mu * psi0_v
926  + n1_12 * dpsi1_u * psi0_1mv
927  + n1_22 * dpsi1_1mu * psi0_1mv) ;
928 
929  dfdy += (n2_11 * dpsi0_u * psi1_v
930  - n2_21 * dpsi0_1mu * psi1_v
931  - n2_12 * dpsi0_u * psi1_1mv
932  + n2_22 * dpsi0_1mu * psi1_1mv) * dz /dy ;
933 
934  double dpsi0_v = 6.*(v2 - v) ;
935  double dpsi0_1mv = 6.*(v2 - v) ;
936  double dpsi1_v = 3.*v2 - 4.*v + 1. ;
937  double dpsi1_1mv = 3.*v2 - 2.*v ;
938 
939 
940  dfdz = (p_11* psi0_u * dpsi0_v
941  + p_21 * psi0_1mu * dpsi0_v
942  - p_12 * psi0_u * dpsi0_1mv
943  - p_22 * psi0_1mu * dpsi0_1mv) / dz ;
944 
945  dfdz += (n1_11 * psi1_u * dpsi0_v
946  - n1_21 * psi1_1mu * dpsi0_v
947  - n1_12 * psi1_u * dpsi0_1mv
948  + n1_22 * psi1_1mu * dpsi0_1mv) * dy / dz ;
949 
950  dfdz += (n2_11 * psi0_u * dpsi1_v
951  + n2_21 * psi0_1mu * dpsi1_v
952  + n2_12 * psi0_u * dpsi1_1mv
953  + n2_22 * psi0_1mu * dpsi1_1mv) ;
954 
955  return ;
956 }
957 
958 /*interpolation for functions of 3 variables : hermite interpolation on the 2 last variables and linear interpolation on the first one*/
959 
960 /*
961 xtab/x refers to delta_car (logdelta_car)
962 ytab/y refers to mu_n (logent1)
963 ztab/z refers to mu_p (logent2)
964 ftab/f refers to psi (logp) or alpha (dlpsddelta_car or logalpha)
965 dfdytab/dfdy refers to dpsi/dmu_n (dlpsdlent1) or dalpha/dmu_n (d2lpsdlent1ddelta_car or dlalphadlent1)
966 dfdztab/dfdz refers to dpsi/dmu_p (dlpsdlent2) or dalpha/dmu_p (dl2psdlent2ddelta_car or dlalphadlent2)
967 d2fdytabdztab refers to d2psi/dmu_ndmu_p (d2lpsdlent1dlent2) or d2alpha/dmu_ndmu_p (d3lpsdlent1dlent2ddelta_car or d2lalphadlent1dlent2)
968 */
969 
970 /*this routine provides the interpolated values of f, dfdy and dfdz at point (x,y,z) via the use of adapted tables*/
971 
972  void interpol_mixed_3d_new(double m_1, double m_2, const Tbl& xtab, const Tbl& ytab,
973  const Tbl& ztab, const Tbl& ftab, const Tbl& dfdytab,
974  const Tbl& dfdztab, const Tbl& d2fdydztab, const Tbl&
975  dlpsddelta_car, const Tbl& d2lpsdlent1ddelta_car, const
976  Tbl& d2lpsdlent2ddelta_car, const Tbl& mu2_P, const Tbl&
977  n_p_P, const Tbl& press_P, const Tbl& mu1_N, const Tbl&
978  n_n_N, const Tbl& press_N, const Tbl& delta_car_n0, const
979  Tbl& mu1_n0, const Tbl& mu2_n0, const Tbl& delta_car_p0,
980  const Tbl& mu1_p0, const Tbl& mu2_p0, double x, double y,
981  double z, double& f, double& dfdy, double& dfdz, double& alpha)
982  {
983  assert(ytab.dim == xtab.dim) ;
984  assert(ztab.dim == xtab.dim) ;
985  assert(ftab.dim == xtab.dim) ;
986  assert(dfdytab.dim == xtab.dim) ;
987  assert(dfdztab.dim == xtab.dim) ;
988  assert(d2fdydztab.dim == xtab.dim) ;
989 
990  int nbp1, nbp2, nbp3;
991  nbp1 = xtab.get_dim(2) ; // \Delta^{2}
992  nbp2 = xtab.get_dim(1) ; // \mu_n
993  nbp3 = xtab.get_dim(0) ; // \mu_p
994 
995  int i_near = 0 ;
996  int j_near = 0 ;
997  int k_near = 0 ;
998 
999  // look for the positions of (x,y,z) in the tables
1000  while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
1001  i_near++ ;
1002  }
1003  if (i_near != 0) {
1004  i_near -- ;
1005  }
1006 
1007  while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
1008  j_near++ ;
1009  }
1010  if (j_near != 0) {
1011  j_near -- ;
1012  }
1013 
1014  while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
1015  k_near++ ;
1016  }
1017  if (k_near != 0) {
1018  k_near-- ;
1019  }
1020 
1021  int i1 = i_near + 1 ;
1022  int j1 = j_near + 1 ;
1023  int k1 = k_near + 1 ;
1024 
1025  // Remarque : cette recherche implique qu'on est au moins supérieure a la premiere valeur (il faudrait mettre un abort sinon!)
1026  // Pour vérifier que la recherche se déroule comme prévue :
1027  // cout << " DEBUG mode : " << endl ;
1028  if ( ( xtab( i_near, j_near, k_near) > x ) || (x > xtab( i1, j_near, k_near) ) ) {
1029  cout << "mauvais positionnement de x dans xtab " << endl ;
1030  cout << xtab( i_near, j_near, k_near) << " " << x << " "
1031  << xtab( i1, j_near, k_near) << endl;
1032  abort();
1033  }
1034 
1035  if ( ( ytab( i_near, j_near, k_near) > y ) || (y > ytab( i1, j1, k_near) ) ) {
1036  cout << "mauvais positionnement de y dans ytab " << endl ;
1037  cout << ytab( i_near, j_near, k_near) * 1.009000285 * 1.66e-27 * 3e8 * 3e8
1038  /(1.6e-19 * 1e6) << " " << y << " " << ytab( i1, j1, k_near)
1039  * 1.009000285 * 1.66e-27 * 3e8 * 3e8 /(1.6e-19 * 1e6) << endl;
1040  abort();
1041  }
1042 
1043  if ( ( ztab( i_near, j_near, k_near) > z ) || ( z > ztab( i1, j_near, k1) ) ){
1044  cout << "mauvais positionnement de z dans ztab " << endl ;
1045  cout << ztab( i_near, j_near, k_near) << " " << z << " "
1046  << ztab( i1, j_near, k1) << endl;
1047  abort();
1048  }
1049 
1050  double f_i_near =0. ;
1051  double dfdy_i_near =0. ;
1052  double dfdz_i_near =0. ;
1053  double alpha_i_near = 0.;
1054  double f_i1 =0. ;
1055  double dfdy_i1 =0. ;
1056  double dfdz_i1 =0. ;
1057  double alpha_i1 = 0.;
1058 
1059  int n_deltaN = delta_car_n0.get_dim(1) ;
1060  int n_mu1N = delta_car_n0.get_dim(0) ;
1061  int n_deltaP = delta_car_p0.get_dim(1) ;
1062  int n_mu2P = delta_car_p0.get_dim(0) ;
1063 
1064  //-----------------------------------------
1065  //----------------TRANCHE i --------------
1066  //-----------------------------------------
1067 
1068  //REPERAGE DES ZONES:
1069  int Placei = 0 ; // 0 = 2fluides, 1 = fluideN seul, 2 = fluideP seul
1070 
1071  int i_nearN_i = 0;
1072  int j_nearN_i = 0;
1073  int i_nearP_i = 0;
1074  int j_nearP_i = 0;
1075 
1076 
1077  if ( y > m_1 ) // Dans cette zone : soit deux fluides, soit fluideN seul (ie np=0)
1078  {
1079  // on enregistre l'indice correpondant au delta_car le plus proche de xtab(i_near, j_near, k_near)
1080  while ( ( (delta_car_n0)(i_nearN_i,0) <= xtab(i_near, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i ) ) {
1081  i_nearN_i++ ;
1082  }
1083  if (i_nearN_i != 0) {
1084  i_nearN_i -- ;
1085  }
1086  // on localise maintenant ou se situe mu_n dans la table de mun
1087  while ( ( (mu1_n0)(i_nearN_i,j_nearN_i) <= y ) && ( ( n_mu1N-1 ) > j_nearN_i ) ) {
1088  j_nearN_i++ ;
1089  }
1090  if (j_nearN_i != 0) {
1091  j_nearN_i -- ;
1092  }
1093 
1094 
1095  // Pour vérifier le positionnement :
1096  //cout << " DEBUG mode = " << (delta_car_n0)(i_nearN_i,0) << " " << xtab( i_near, j_near, k_near) << endl ;
1097 // if ( (delta_car_n0)(i_nearN_i,0) != xtab( i_near, j_near, k_near) ) {
1098 // cout << "c'est un test : " << endl;
1099 // abort();
1100 // }
1101 // if ( ( (delta_car_n0)(i_nearN_i,0) > xtab(i_near, j_near, k_near) ) || (xtab(i_near, j_near, k_near) > (delta_car_n0)(i_nearN_i+1,0) ) ) {
1102 // cout << "mauvais positionnement de delta_car_i dans delta_car_n0 (courbe limite np = 0) " << endl ;
1103 // cout << (delta_car_n0)(i_nearN_i,0) << " " << xtab(i_near, j_near, k_near) << " " << (delta_car_n0)(i_nearN_i+1,0) << endl;
1104 // abort();
1105 // }
1106 // if ( ( (mu1_n0)(i_nearN_i,j_nearN_i) > y ) || (y > (mu1_n0)(i_nearN_i,j_nearN_i+1) ) ) {
1107 // cout << "mauvais positionnement demu_n dans mu1_n0 (courbe limite np = 0) " << endl ;
1108 // cout << (mu1_n0)(i_nearN_i,j_nearN_i) << " " << y << " " << (mu1_n0)(i_nearN_i,j_nearN_i+1) << endl;
1109 // abort();
1110 // }
1111 
1112 
1113  // on regarde alors ou on est par rapport a la courbe np = 0.
1114  double aN_i, bN_i;
1115  aN_i = ((mu2_n0)(i_nearN_i,j_nearN_i+1) - (mu2_n0)(i_nearN_i,j_nearN_i) ) / ((mu1_n0)(i_nearN_i,j_nearN_i+1) - (mu1_n0)(i_nearN_i,j_nearN_i) ) ;
1116  bN_i = (mu2_n0)(i_nearN_i,j_nearN_i) - aN_i * (mu1_n0)(i_nearN_i,j_nearN_i) ;
1117  double zN_i = aN_i * y + bN_i ;
1118 
1119  if (zN_i < z) { // Zone ou deux fluides
1120  Placei = 0;
1121  }
1122  else { // Zone ou fluide N seul
1123  Placei = 1 ;
1124  }
1125 // cout << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1126 // z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << m_1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<
1127 // m_2 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<endl;
1128 // cout << " Placei = " << Placei<< endl ;
1129 
1130  }
1131 
1132  else //if ( y <= m_1) // Dans cette zone : soit deux fluides, soit fluideP seul (ie nn=0), soit zéro fluide !!!!
1133  {
1134 
1135  //cout << y << " " << m_1 << endl ;
1136 
1137  if ( z <= m_2) {
1138  Placei = 3 ; // NO fluid at all !
1139  }
1140  else {
1141 
1142  // on enregistre l'indice correpondant au delta_car le plus proche de xx
1143  while ( ( (delta_car_p0)(i_nearP_i,0) <= xtab(i_near, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i ) ) {
1144  i_nearP_i++ ;
1145  }
1146  if (i_nearP_i != 0) {
1147  i_nearP_i -- ;
1148  }
1149 
1150  // on localise maintenant ou se situe z dans la table de mup
1151  while ( ( (mu2_p0)(i_nearP_i,j_nearP_i) <= z ) && ( ( n_mu2P-1 ) > j_nearP_i ) ) {
1152  j_nearP_i++ ;
1153  }
1154  if (j_nearP_i != 0) {
1155  j_nearP_i -- ;
1156  }
1157 
1158  // Pour vérifier le positionnement :
1159  // cout << " DEBUG mode = " << (delta_car_p0)(i_nearP_i,0) << " " << xtab( i_near, j_near, k_near) << endl ;
1160 // if ( (delta_car_p0)(i_nearP_i,0) != xtab( i_near, j_near, k_near) ) {
1161 // cout << "c'est un test : " << endl;
1162 // abort();
1163 // }
1164 // if ( ( (delta_car_p0)(i_nearP_i,0) > xtab(i_near, j_near, k_near) ) || (xtab(i_near, j_near, k_near) > (delta_car_p0)(i_nearP_i+1,0) ) ) {
1165 // cout << "mauvais positionnement de delta_car_i dans delta_car_p0 (courbe limite nn = 0) " << endl ;
1166 // cout << (delta_car_p0)(i_nearP_i,0) << " " << xtab(i_near, j_near, k_near) << " " << (delta_car_p0)(i_nearP_i+1,0) << endl;
1167 // abort();
1168 // }
1169 // if ( ( (mu2_p0)(i_nearP_i,j_nearP_i) > y ) || (y > (mu2_p0)(i_nearP_i,j_nearP_i+1) ) ) {
1170 // cout << "mauvais positionnement demu_p dans mu2_p0 (courbe limite nn = 0) " << endl ;
1171 // cout << (mu2_p0)(i_nearP_i,j_nearP_i) << " " << y << " " << (mu2_p0)(i_nearP_i,j_nearP_i+1) << endl;
1172 // abort();
1173 // }
1174 
1175 
1176  // on regarde alors ou on est par rapport a la droite nn = 0.
1177  double aP_i, bP_i;
1178  aP_i = ( (mu2_p0)(i_nearP_i,j_nearP_i+1) - (mu2_p0)(i_nearP_i,j_nearP_i) ) / ( (mu1_p0)(i_nearP_i,j_nearP_i+1) - (mu1_p0)(i_nearP_i,j_nearP_i) ) ;
1179  bP_i = (mu2_p0)(i_nearP_i,j_nearP_i) - aP_i * (mu1_p0)(i_nearP_i,j_nearP_i) ;
1180  double yP_i = (z- bP_i) /aP_i ;
1181 
1182 
1183  if (yP_i < y) { // Zone ou deux fluides
1184  Placei = 0;
1185  }
1186  else { // Zone ou fluide 2 seul
1187  Placei = 2 ;
1188  }
1189 
1190  // cout << yP_i * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1191  // z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << endl;
1192 
1193  }
1194 // cout << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1195 // z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << m_1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<
1196 // m_2 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<endl;
1197 // cout << " Placei = " << Placei<< endl ;
1198 
1199 }
1200 //FIN DU REPERAGE -----------------------------------------------
1201 
1202  // traitement au cas par cas
1203 
1204  if (Placei == 3 ) { // NO FLUID
1205  f_i_near = 0. ;
1206  dfdy_i_near = 0. ;
1207  dfdz_i_near = 0. ;
1208  alpha_i_near = 0.;
1209  }
1210  else if (Placei == 1 ) {
1211  //cout << "fluide N seul" << endl;
1212  alpha_i_near = 0.;
1213  dfdz_i_near = 0. ;
1214  int i = 0;
1215  interpol_herm(mu1_N, press_N, n_n_N,
1216  y, i, f_i_near , dfdy_i_near ) ;
1217  if (f_i_near < 0.) {
1218  cout << " INTERPOLATION FLUID N --> negative pressure " << endl ;
1219  abort();
1220  // f_i_near = 0.;
1221  }
1222  if (dfdy_i_near < 0.) {
1223  cout << " INTERPOLATION FLUID N --> negative density " << endl ;
1224  abort();
1225  // dfdy_i_near = 0.;
1226  }
1227  }
1228  else if (Placei == 2 ) {
1229  //cout << "fluide P seul" << endl;
1230  alpha_i_near = 0.;
1231  dfdy_i_near = 0. ;
1232  int i =0;
1233  interpol_herm( mu2_P, press_P, n_p_P,
1234  z, i, f_i_near, dfdz_i_near) ;
1235  if (f_i_near < 0.) {
1236  cout << " INTERPOLATION FLUID P --> negative pressure " << endl ;
1237  abort();
1238  // f_i_near = 0.;
1239  }
1240  if (dfdz_i_near < 0.) {
1241  cout << " INTERPOLATION FLUID P --> negative density " << endl ;
1242  abort();
1243  // dfdz_i_near = 0.;
1244  }
1245  }
1246  else if (Placei == 0 ) {
1247  /*
1248  // les deux fluides sont presents
1249  // on regarde alors si on a des termes nuls montrant qu'on est a la surface
1250 
1251  if ( (dfdytab( i_near, j_near, k_near) > 0. ) && (dfdztab( i_near, j_near, k_near) > 0. ) &&
1252  (dfdytab( i_near, j1, k_near) > 0. ) && (dfdztab( i_near, j1, k_near) > 0. ) &&
1253  (dfdytab( i_near, j_near, k1) > 0. ) && (dfdztab( i_near, j_near, k1) > 0. ) &&
1254  (dfdytab( i_near, j1, k1) > 0. ) && (dfdztab( i_near, j1, k1) > 0. ) ) {
1255 
1256  //cout << "DEBUG : LOIN DES COURBES LIMITES" << endl;
1257  interpol_mixed_3d(xtab, ytab, ztab, ftab,
1258  dfdytab, dfdztab, d2fdydztab,
1259  xtab(i_near, j_near, k_near), y, z, f_i_near, dfdy_i_near , dfdz_i_near ) ;
1260 
1261  //cout << "f_i_near " << f_i_near << " dfdy_i_near " << dfdy_i_near << " dfdz_i_near " << dfdz_i_near << endl ;
1262 
1263  if (f_i_near < 0.) {
1264 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1265  cout << " INTERPOLATION 2 FLUIDS --> negative pressure " << endl ;
1266 // abort();
1267  }
1268  if (dfdy_i_near < 0.) {
1269 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1270  cout << " INTERPOLATION 2 FLUIDS --> negative density for fluid N " << endl ;
1271 // abort();
1272  }
1273  if (dfdz_i_near < 0.) {
1274 
1275 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1276  cout << " INTERPOLATION 2 FLUIDS --> negative density for fluid P " << endl ;
1277 // abort();
1278  }
1279 
1280  double der1 = 0., der2 = 0.;
1281  interpol_mixed_3d_mod(xtab, ytab, ztab, dlpsddelta_car,
1282  d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car,
1283  xtab(i_near, j_near, k_near), y, z, alpha_i_near, der1 , der2 ) ;
1284 
1285  alpha_i_near = - alpha_i_near ;
1286  //cout << " alpha_i_near" << alpha_i_near << endl ;
1287 
1288  }
1289 
1290  else if ( (dfdytab( i_near, j_near, k_near) == 0. ) && (dfdztab( i_near, j_near, k_near) == 0. ) &&
1291  (dfdytab( i_near, j1, k_near) > 0. ) && (dfdztab( i_near, j1, k_near) == 0. ) &&
1292  (dfdytab( i_near, j_near, k1) == 0. ) && (dfdztab( i_near, j_near, k1) > 0. ) &&
1293  (dfdytab( i_near, j1, k1) > 0. ) && (dfdztab( i_near, j1, k1) > 0. ) ) {
1294 
1295 
1296  //cout << "Croisement des courbes limites - Tranche i" << endl;
1297 // // ********** point pris en haut a gauche
1298  double aP_inf, bP_inf;
1299  aP_inf = ( mu2_p0(i_nearP_i,j_nearP_i+1) - mu2_p0(i_nearP_i,j_nearP_i) ) /
1300  ( mu1_p0(i_nearP_i,j_nearP_i+1) - mu1_p0(i_nearP_i,j_nearP_i) ) ;
1301  bP_inf = mu2_p0(i_nearP_i,j_nearP_i) - aP_inf * mu1_p0(i_nearP_i,j_nearP_i) ;
1302 
1303  double mu2_nul_inf_Left = ztab(i_near, j1, k1) ;
1304  double mu1_nul_inf_Left = (mu2_nul_inf_Left - bP_inf) / aP_inf ;
1305 //
1306 // // ********** point pris en bas a droite
1307 //
1308  double aN_inf, bN_inf;
1309  aN_inf = (mu2_n0(i_nearN_i,j_nearN_i +1) - mu2_n0(i_nearN_i,j_nearN_i ) ) /
1310  (mu1_n0(i_nearN_i,j_nearN_i +1) - mu1_n0(i_nearN_i,j_nearN_i ) ) ;
1311  bN_inf = mu2_n0(i_nearN_i,j_nearN_i ) - aN_inf * mu1_n0(i_nearN_i,j_nearN_i ) ;
1312 
1313 
1314  double mu1_nul_inf_Right= ytab(i_near, j1, k1) ;
1315  double mu2_nul_inf_Right = aN_inf * mu1_nul_inf_Right + bN_inf ;
1316 //
1317  double mu1_11, mu1_21, mu2_11, mu2_12 ;
1318  double p_11,p_21,p_12,p_22 ;
1319  double n1_11,n1_21,n1_12, n1_22, n2_11 ,n2_21, n2_12, n2_22 ;
1320  double cross_11 , cross_21, cross_12, cross_22 ;
1321  double dp2sddelta_car_11,dp2sddelta_car_21,dp2sddelta_car_12, dp2sddelta_car_22 ;
1322  double d2psdent1ddelta_car_11 ,d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22 ;
1323  double d2psdent2ddelta_car_11,d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22 ;
1324 
1325 
1326  mu1_11 = mu1_nul_inf_Left ;
1327  mu1_21 = mu1_nul_inf_Right ;
1328  mu2_11 = mu2_nul_inf_Right ;
1329  mu2_12 = mu2_nul_inf_Left ;
1330  p_21 = ftab(i_near,j1, k_near) ;
1331  p_12 = ftab(i_near,j_near, k1) ;
1332  p_22 = ftab(i_near,j1, k1) ;
1333  n1_21 = dfdytab(i_near,j1, k_near) ;
1334  n1_12 = dfdytab(i_near,j_near, k1) ;
1335  n1_22 = dfdytab(i_near,j1, k1) ;
1336  n2_21 = dfdztab(i_near,j1, k_near) ;
1337  n2_12 = dfdztab(i_near,j_near, k1) ;
1338  n2_22 = dfdztab(i_near,j1, k1) ;
1339  cross_21 = d2fdydztab(i_near,j1, k_near) ;
1340  cross_12 = d2fdydztab(i_near,j_near, k1) ;
1341  cross_22 = d2fdydztab(i_near,j1, k1) ;
1342  dp2sddelta_car_21 = dlpsddelta_car(i_near,j1, k_near) ;
1343  dp2sddelta_car_12 = dlpsddelta_car(i_near,j_near, k1) ;
1344  dp2sddelta_car_22 = dlpsddelta_car(i_near,j1, k1) ;
1345  d2psdent1ddelta_car_21 = d2lpsdlent1ddelta_car(i_near,j1, k_near) ;
1346  d2psdent1ddelta_car_12 = d2lpsdlent1ddelta_car(i_near,j_near, k1) ;
1347  d2psdent1ddelta_car_22 = d2lpsdlent1ddelta_car(i_near,j1, k1) ;
1348  d2psdent2ddelta_car_21 = d2lpsdlent2ddelta_car(i_near,j1, k_near) ;
1349  d2psdent2ddelta_car_12 = d2lpsdlent2ddelta_car(i_near,j_near, k1) ;
1350  d2psdent2ddelta_car_22 = d2lpsdlent2ddelta_car(i_near,j1, k1) ;
1351  p_11 = 0.;
1352  n1_11 = 0. ;
1353  n2_11 = 0.;
1354  cross_11 = 0.;
1355  dp2sddelta_car_11 = 0. ;
1356  d2psdent1ddelta_car_11 = 0. ;
1357  d2psdent2ddelta_car_11 = 0. ;
1358 
1359 // // changmeent
1360  interpol_herm_2d_new_avec( y, z,
1361  mu1_11, mu1_21, mu2_11,mu2_12,
1362  p_11,p_21, p_12, p_22,
1363  n1_11, n1_21, n1_12, n1_22,
1364  n2_11, n2_21, n2_12, n2_22,
1365  cross_11, cross_21, cross_12, cross_22,
1366  f_i_near, dfdy_i_near, dfdz_i_near) ;
1367 
1368 
1369  // cout << "f_i_near " << f_i_near << " dfdy_i_near " << dfdy_i_near << " dfdz_i_near " << dfdz_i_near << endl ;
1370 
1371  double der1= 0., der2=0.;
1372  interpol_herm_2d_new_sans( y, z,
1373  mu1_11, mu1_21, mu2_11, mu2_12,
1374  dp2sddelta_car_11, dp2sddelta_car_21, dp2sddelta_car_12, dp2sddelta_car_22,
1375  d2psdent1ddelta_car_11, d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22,
1376  d2psdent2ddelta_car_11, d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22,
1377  alpha_i_near, der1, der2) ;
1378  alpha_i_near = - alpha_i_near ;
1379 
1380  if (f_i_near < 0.) {
1381 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1382  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative pressure " << endl ;
1383 // abort();
1384  }
1385  if (dfdy_i_near < 0.) {
1386 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1387  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative density for fluid N " << endl ;
1388 // abort();
1389  }
1390  if (dfdz_i_near < 0.) {
1391 
1392 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1393  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative density for fluid P " << endl ;
1394 // abort();
1395  }
1396 
1397 
1398  //cout << " alpha_i_near" << alpha_i_near << endl ;
1399  }
1400  else if ( (dfdytab( i_near, j_near, k_near) == 0. ) && (dfdztab( i_near, j_near, k_near) > 0. ) &&
1401  (dfdytab( i_near, j_near, k1) == 0. ) && (dfdztab( i_near, j_near, k1) > 0. ) &&
1402  (dfdytab( i_near, j1, k_near) > 0. ) && (dfdztab( i_near, j1, k_near) > 0. ) &&
1403  (dfdytab( i_near, j1, k1) > 0. ) && (dfdztab( i_near, j1, k1) > 0. ) ) {
1404 //
1405  //cout << "cas Fluide de Neutrons seul / Courbe limite - Tranche i" << endl;
1407  double aP_inf, bP_inf;
1408  aP_inf = ( mu2_p0(i_nearP_i,j_nearP_i+1) - mu2_p0(i_nearP_i,j_nearP_i) ) /
1409  ( mu1_p0(i_nearP_i,j_nearP_i+1) - mu1_p0(i_nearP_i,j_nearP_i) ) ;
1410  bP_inf = mu2_p0(i_nearP_i,j_nearP_i) - aP_inf * mu1_p0(i_nearP_i,j_nearP_i) ;
1411 
1412  double mu2_nul_inf = ztab(i_near, j1, k1) ;
1413  double mu1_nul_inf = (mu2_nul_inf - bP_inf)/aP_inf;
1414  //cout << mu1_nul_inf * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " << mu2_nul_inf* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << endl;
1415  double mu1_11, mu1_21, mu2_11, mu2_12 ;
1416  double p_11,p_21,p_12,p_22 ;
1417  double n1_11,n1_21,n1_12, n1_22, n2_11 ,n2_21, n2_12, n2_22 ;
1418  double cross_11 , cross_21, cross_12, cross_22 ;
1419  double dp2sddelta_car_11,dp2sddelta_car_21,dp2sddelta_car_12, dp2sddelta_car_22 ;
1420  double d2psdent1ddelta_car_11 ,d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22 ;
1421  double d2psdent2ddelta_car_11,d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22 ;
1423 
1424  //cout << " aP_inf = " << aP_inf << endl ;
1425  //cout << " bP_inf = " << bP_inf << endl ;
1426  //cout << " mu2_nul_inf " << mu2_nul_inf << endl ;
1427  //cout << " mu1_nul_inf " << mu1_nul_inf << endl ;
1428 
1429  mu1_11 = mu1_nul_inf ;
1430  mu1_21 = ytab(i_near,j1, k_near) ;
1431  mu2_11 = ztab(i_near,j1, k_near) ;
1432  mu2_12 = mu2_nul_inf ;
1433  p_21 = ftab(i_near,j1, k_near) ;
1434  p_12 = ftab(i_near,j_near, k1) ;
1435  p_22 = ftab(i_near,j1, k1) ;
1436  n1_21 = dfdytab(i_near,j1, k_near) ;
1437  n1_12 = dfdytab(i_near,j_near, k1) ;
1438  n1_22 = dfdytab(i_near,j1, k1) ;
1439  n2_21 = dfdztab(i_near,j1, k_near) ;
1440  n2_12 = dfdztab(i_near,j_near, k1) ;
1441  n2_22 = dfdztab(i_near,j1, k1) ;
1442  cross_21 = d2fdydztab(i_near,j1, k_near) ;
1443  cross_12 = d2fdydztab(i_near,j_near, k1) ;
1444  cross_22 = d2fdydztab(i_near,j1, k1) ;
1445  p_11 = ftab(i_near,j_near, k_near) ;
1446  n1_11 = 0. ; //dfdytab(i_near,j_near, k_near) ;
1447  //cout << "n1_11 " << n1_11 << endl ;
1448  n2_11 = dfdztab(i_near,j_near, k_near) ;
1449  cross_11 = 0.; //d2fdydztab(i_near,j_near, k_near) ;
1450  //cout << "cross_11 " << cross_11 << endl ;
1451  dp2sddelta_car_21 = dlpsddelta_car(i_near,j1, k_near) ;
1452  dp2sddelta_car_12 = dlpsddelta_car(i_near,j_near, k1) ;
1453  dp2sddelta_car_22 = dlpsddelta_car(i_near,j1, k1) ;
1454  d2psdent1ddelta_car_21 = d2lpsdlent1ddelta_car(i_near,j1, k_near) ;
1455  d2psdent1ddelta_car_12 = d2lpsdlent1ddelta_car(i_near,j_near, k1) ;
1456  d2psdent1ddelta_car_22 = d2lpsdlent1ddelta_car(i_near,j1, k1) ;
1457  d2psdent2ddelta_car_21 = d2lpsdlent2ddelta_car(i_near,j1, k_near) ;
1458  d2psdent2ddelta_car_12 = d2lpsdlent2ddelta_car(i_near,j_near, k1) ;
1459  d2psdent2ddelta_car_22 = d2lpsdlent2ddelta_car(i_near,j1, k1) ;
1460  dp2sddelta_car_11 = 0. ;//dlpsddelta_car(i_near,j_near, k_near) ;
1461  //cout << dp2sddelta_car_11 << endl ;
1462  d2psdent1ddelta_car_11 =0. ;// d2lpsdlent1ddelta_car(i_near,j_near, k_near) ;
1463  //cout << d2psdent1ddelta_car_11 << endl ;
1464  d2psdent2ddelta_car_11 = 0. ; //d2lpsdlent2ddelta_car(i_near,j_near, k_near) ;
1465  //cout << d2psdent2ddelta_car_11 << endl ;
1466 
1467  interpol_herm_2d_new_avec( y, z,
1468  mu1_11, mu1_21, mu2_11,mu2_12,
1469  p_11,p_21, p_12, p_22,
1470  n1_11, n1_21, n1_12, n1_22,
1471  n2_11, n2_21, n2_12, n2_22,
1472  cross_11, cross_21, cross_12, cross_22,
1473  f_i_near, dfdy_i_near, dfdz_i_near) ;
1474 
1475  //cout << "f_i_near " << f_i_near << " dfdy_i_near " << dfdy_i_near << " dfdz_i_near " << dfdz_i_near << endl ;
1476 //
1477  double der1= 0., der2=0.;
1478 
1479  //cout << "y " << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " z " << z * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< endl ;
1480  interpol_herm_2d_new_sans( y, z,
1481  mu1_11, mu1_21, mu2_11, mu2_12,
1482  dp2sddelta_car_11, dp2sddelta_car_21, dp2sddelta_car_12, dp2sddelta_car_22,
1483  d2psdent1ddelta_car_11, d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22,
1484  d2psdent2ddelta_car_11, d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22,
1485  alpha_i_near, der1, der2) ;
1486  alpha_i_near = - alpha_i_near ;
1487  //cout << " alpha_i_near " << alpha_i_near << endl ;
1488 
1489  if (f_i_near < 0.) {
1490 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1491  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative pressure " << endl ;
1492 // abort();
1493  }
1494  if (dfdy_i_near < 0.) {
1495 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1496  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative density for fluid N " << endl ;
1497 // abort();
1498  }
1499  if (dfdz_i_near < 0.) {
1500 
1501 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1502  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative density for fluid P " << endl ;
1503 // abort();
1504  }
1505 
1506 //
1507  }
1508  else {
1509 
1510 */
1511 
1512 //cout << "TRAITEMENT DE LA SURFACE PAS ENCORE REALISE - Tranche i" << endl;
1513  interpol_mixed_3d(xtab, ytab, ztab, ftab,
1514  dfdytab, dfdztab, d2fdydztab,
1515  xtab(i_near, j_near, k_near), y, z, f_i_near, dfdy_i_near , dfdz_i_near ) ;
1516 
1517  //cout << "f_i_near " << f_i_near << " dfdy_i_near " << dfdy_i_near << " dfdz_i_near " << dfdz_i_near << endl ;
1518 
1519  double der1 = 0., der2 = 0.;
1520  interpol_mixed_3d_mod(xtab, ytab, ztab, dlpsddelta_car,
1521  d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car,
1522  xtab(i_near, j_near, k_near), y, z, alpha_i_near, der1 , der2 ) ;
1523  // cout << "alpha " << alpha << endl;
1524  alpha_i_near = - alpha_i_near ;
1525 
1526 /* if (f_i_near < 0.) {
1527 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1528  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative pressure " << endl ;
1529 // abort();
1530  }
1531  if (dfdy_i_near < 0.) {
1532 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1533  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative density for fluid N " << endl ;
1534 // abort();
1535  }
1536  if (dfdz_i_near < 0.) {
1537 
1538 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1539  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative density for fluid P " << endl ;
1540 // abort();
1541  }
1542 */
1543 
1544  //cout << " alpha_i_near" << alpha_i_near << endl ;
1545 // //cout << "surface non 6" << endl;
1546 // // Zone en plein milieu de la table a deux fluides
1547 //
1548 
1549 
1550 /*}*/
1551 
1552 
1553 }
1554 
1555 
1556 
1557  //-----------------------------------------
1558  //--------------TRANCHE i+1 --------------
1559  //-----------------------------------------
1560 
1561 
1562 
1563 
1564  //REPERAGE DES ZONES:
1565  int Placei1 = 0 ; // 0 = 2fluides, 1 = fluideN seul, 2 = fluideP seul, 3 = 0 fluide !
1566 
1567  int i_nearN_i1 = 0;
1568  int j_nearN_i1 = 0;
1569  int i_nearP_i1 = 0;
1570  int j_nearP_i1 = 0;
1571 
1572  if ( y > m_1 ) // Dans cette zone : soit deux fluides, soit fluideN seul (ie np=0)
1573  {
1574 
1575  // on enregistre l'indice correpondant au delta_car le plus proche de xtab(i_near, j_near, k_near)
1576  while ( ( (delta_car_n0)(i_nearN_i1,0) <= xtab(i_near+1, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i1 ) ) {
1577  i_nearN_i1++ ;
1578  }
1579  if (i_nearN_i1 != 0) {
1580  i_nearN_i1 -- ;
1581  }
1582  // on localise maintenant ou se situe mu_n dans la table de mun
1583  while ( ( (mu1_n0)(i_nearN_i1,j_nearN_i1) <= y ) && ( ( n_mu1N-1 ) > j_nearN_i1 ) ) {
1584  j_nearN_i1++ ;
1585  }
1586  if (j_nearN_i1 != 0) {
1587  j_nearN_i1 -- ;
1588  }
1589 
1590  // Pour vérifier le positionnement :
1591 // cout << " DEBUG mode = " << (delta_car_n0)(i_nearN_i1,0) << " " << xtab( i_near+1, j_near, k_near) << endl ;
1592 // if ( (delta_car_n0)(i_nearN_i1,0) != xtab( i_near+1, j_near, k_near) ) {
1593 // cout << "c'est un test : " << endl;
1594 // abort();
1595 // }
1596 // if ( ( (delta_car_n0)(i_nearN_i1,0) > xtab(i_near+1, j_near, k_near) ) || (xtab(i_near+1, j_near, k_near) > (delta_car_n0)(i_nearN_i1+1,0) ) ) {
1597 // cout << "mauvais positionnement de delta_car_i+1 dans delta_car_n0 (courbe limite np = 0) " << endl ;
1598 // cout << (delta_car_n0)(i_nearN_i1,0) << " " << xtab(i_near+1, j_near, k_near) << " " << (delta_car_n0)(i_nearN_i1+1,0) << endl;
1599 // abort();
1600 // }
1601 // if ( ( (mu1_n0)(i_nearN_i1,j_nearN_i1) > y ) || (y > (mu1_n0)(i_nearN_i1,j_nearN_i1+1) ) ) {
1602 // cout << "mauvais positionnement demu_n dans mu1_n0 (courbe limite np = 0) " << endl ;
1603 // cout << (mu1_n0)(i_nearN_i1,j_nearN_i1) << " " << y << " " << (mu1_n0)(i_nearN_i1,j_nearN_i1+1) << endl;
1604 // abort();
1605 // }
1606 
1607  // on regarde alors ou on est par rapport a la courbe np = 0.
1608  double aN_i1, bN_i1;
1609  aN_i1 = ((mu2_n0)(i_nearN_i1,j_nearN_i1+1) - (mu2_n0)(i_nearN_i1,j_nearN_i1) ) / ((mu1_n0)(i_nearN_i1,j_nearN_i1+1) - (mu1_n0)(i_nearN_i1,j_nearN_i1) ) ;
1610  bN_i1 = (mu2_n0)(i_nearN_i1,j_nearN_i1) - aN_i1 * (mu1_n0)(i_nearN_i1,j_nearN_i1) ;
1611  double zN_i1 = aN_i1 * y + bN_i1 ;
1612 
1613  if (zN_i1 < z) { // Zone ou deux fluides
1614  Placei1 = 0;
1615  }
1616  else { // Zone ou fluide N seul
1617  Placei1 = 1 ;
1618  }
1619 
1620 // cout << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1621 // z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << m_1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<
1622 // m_2 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<endl;
1623 // cout << " Placei1 = " << Placei1<< endl ;
1624  }
1625 
1626  else //if ( y <= m_1) // Dans cette zone : soit deux fluides, soit fluideP seul (ie nn=0), soit zéro fluide !!!!
1627  {
1628  if ( z <= m_2) {
1629  Placei1 = 3 ; // NO fluid at all !
1630  }
1631  else {
1632 
1633 
1634  // on enregistre l'indice correpondant au delta_car le plus proche de xx
1635  while ( ( (delta_car_p0)(i_nearP_i1,0) <= xtab(i_near+1, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i1) ) {
1636  i_nearP_i1++ ;
1637  }
1638  if (i_nearP_i1 != 0) {
1639  i_nearP_i1 -- ;
1640  }
1641  // on localise maintenant ou se situe z dans la table de mup
1642  while ( ( (mu2_p0)(i_nearP_i1,j_nearP_i1) <= z ) && ( ( n_mu2P-1 ) > j_nearP_i1 ) ) {
1643  j_nearP_i1++ ;
1644  }
1645  if (j_nearP_i1 != 0) {
1646  j_nearP_i1 -- ;
1647  }
1648 
1649  // Pour vérifier le positionnement :
1650 // cout << " DEBUG mode = " << (delta_car_p0)(i_nearP_i1,0) << " " << xtab( i_near+1, j_near, k_near) << endl ;
1651 // if ( (delta_car_p0)(i_nearP_i1,0) != xtab( i_near+1, j_near, k_near) ) {
1652 // cout << "c'est un test : " << endl;
1653 // abort();
1654 // }
1655 // if ( ( (delta_car_p0)(i_nearP_i1,0) > xtab(i_near+1, j_near, k_near) ) || (xtab(i_near+1, j_near, k_near) > (delta_car_p0)(i_nearP_i1+1,0) ) ) {
1656 // cout << "mauvais positionnement de delta_car_i+1 dans delta_car_p0 (courbe limite nn = 0) " << endl ;
1657 // cout << (delta_car_p0)(i_nearP_i1,0) << " " << xtab(i_near+1, j_near, k_near) << " " << (delta_car_p0)(i_nearP_i1+1,0) << endl;
1658 // abort();
1659 // }
1660 // if ( ( (mu2_p0)(i_nearP_i1,j_nearP_i1) > y ) || (y > (mu2_p0)(i_nearP_i1,j_nearP_i1+1) ) ) {
1661 // cout << "mauvais positionnement demu_p dans mu2_p0 (courbe limite nn = 0) " << endl ;
1662 // cout << (mu2_p0)(i_nearP_i1,j_nearP_i1) << " " << y << " " << (mu2_p0)(i_nearP_i1,j_nearP_i1+1) << endl;
1663 // abort();
1664 // }
1665 
1666  // on regarde alors ou on est par rapport a la droite nn = 0.
1667  double aP_i1, bP_i1;
1668  aP_i1 = ( (mu2_p0)(i_nearP_i1,j_nearP_i1+1) - (mu2_p0)(i_nearP_i1,j_nearP_i1) ) / ( (mu1_p0)(i_nearP_i1,j_nearP_i1+1) - (mu1_p0)(i_nearP_i1,j_nearP_i1) ) ;
1669  bP_i1 = (mu2_p0)(i_nearP_i1,j_nearP_i1) - aP_i1 * (mu1_p0)(i_nearP_i1,j_nearP_i1) ;
1670  double yP_i1 = (z- bP_i1) /aP_i1 ;
1671 
1672  if (yP_i1 < y) { // Zone ou deux fluides
1673  Placei1 = 0;
1674  }
1675  else { // Zone ou fluide 2 seul
1676  Placei1 = 2 ;
1677  }
1678 
1679  //cout << yP_i1 * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1680  //z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << endl;
1681 
1682  }
1683 
1684 // cout << y * 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13)<< " " <<
1685 // z* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << m_1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<
1686 // m_2 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13)<<endl;
1687 // cout << " Placei1 = " << Placei1<< endl ;
1688 
1689 }
1690 
1691 //FIN DU REPERAGE -----------------------------------------------
1692 
1693 
1694  // traitement au cas par cas
1695  if (Placei1 == 3 ) { // NO FLUID
1696  f_i1 = 0. ;
1697  dfdy_i1 = 0. ;
1698  dfdz_i1 = 0. ;
1699  alpha_i1 = 0.;
1700  }
1701  else if (Placei1 == 1 ) {
1702  //cout << "fluideN seul" << endl;
1703  alpha_i1 = 0. ;
1704  dfdz_i1 = 0. ;
1705  int i =0;
1706  interpol_herm(mu1_N, press_N, n_n_N,
1707  y, i, f_i1 , dfdy_i1 ) ;
1708  if (f_i1 < 0.) {
1709  cout << " INTERPOLATION FLUID N i+1 --> negative pressure " << endl ;
1710  abort();
1711  // f_i1 = 0.;
1712  }
1713  if (dfdy_i1 < 0.) {
1714  cout << " INTERPOLATION FLUID N i+1--> negative density " << endl ;
1715  abort();
1716  // dfdy_i1 = 0.;
1717  }
1718  }
1719  else if (Placei1 == 2 ) {
1720  //cout << "fluideP seul" << endl;
1721  alpha_i1 = 0.;
1722  dfdy_i1 = 0. ;
1723  int i =0;
1724  interpol_herm( mu2_P, press_P, n_p_P,
1725  z, i, f_i1, dfdz_i1) ;
1726  if (f_i1 < 0.) {
1727  cout << " INTERPOLATION FLUID P i+1--> negative pressure " << endl ;
1728  abort();
1729  // f_i1 = 0.;
1730  }
1731  if (dfdz_i1 < 0.) {
1732  cout << " INTERPOLATION FLUID P i+1 --> negative density " << endl ;
1733  abort();
1734  // dfdz_i1= 0.;
1735  }
1736 }
1737 else if (Placei1 == 0 ) {
1738  /*
1739  // les deux fluides sont presents
1740  // on regarde alors si on a des termes nuls montrant qu'on est a la surface
1741 
1742  // ----------pour savoir si on est proche ou non d'une surface
1743 
1744 
1745  if ( (dfdytab( i1, j_near, k_near) > 0. ) && (dfdztab( i1, j_near, k_near) > 0. ) &&
1746  (dfdytab( i1, j1, k_near) > 0. ) && (dfdztab( i1, j1, k_near) > 0. ) &&
1747  (dfdytab( i1, j_near, k1) > 0. ) && (dfdztab( i1, j_near, k1) > 0. ) &&
1748  (dfdytab( i1, j1, k1) > 0. ) && (dfdztab( i1, j1, k1) > 0. ) ) {
1749 
1750  //cout << "DEBUG : LOIN DES COURBES LIMITES" << endl;
1751  interpol_mixed_3d(xtab, ytab, ztab, ftab,
1752  dfdytab, dfdztab, d2fdydztab,
1753  xtab(i1, j_near, k_near), y, z, f_i1, dfdy_i1 , dfdz_i1) ;
1754 
1755  if (f_i1 < 0.) {
1756 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1757  cout << " INTERPOLATION 2 FLUIDS i+1--> negative pressure " << endl ;
1758 // abort();
1759  }
1760  if (dfdy_i1 < 0.) {
1761 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1762  cout << " INTERPOLATION 2 FLUIDS i+1--> negative density for fluid N " << endl ;
1763 // abort();
1764  }
1765  if (dfdz_i1 < 0.) {
1766 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1767  cout << " INTERPOLATION 2 FLUIDS i+1--> negative density for fluid P " << endl ;
1768 // abort();
1769  }
1770 
1771  double der1 = 0., der2 = 0.;
1772  interpol_mixed_3d_mod(xtab, ytab, ztab, dlpsddelta_car,
1773  d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car,
1774  xtab(i1, j_near, k_near), y, z, alpha_i1, der1 , der2 ) ;
1775  alpha_i1 = - alpha_i1 ;
1776  }
1777 //
1778  else if ( (dfdytab( i1, j_near, k_near) == 0. ) && (dfdztab( i1, j_near, k_near) == 0. ) &&
1779  (dfdytab( i1, j1, k_near) > 0. ) && (dfdztab( i1, j1, k_near) == 0. ) &&
1780  (dfdytab( i1, j_near, k1) == 0. ) && (dfdztab( i1, j_near, k1) > 0. ) &&
1781  (dfdytab( i1, j1, k1) > 0. ) && (dfdztab( i1, j1, k1) > 0. ) ) {
1782 //
1783 //
1784  //cout << "Croisement des courbes limites - Tranche i + 1 " << endl;
1785 // // ********** point pris en haut a gauche
1786 //
1787 //
1788  double aP_sup, bP_sup;
1789  aP_sup = ( mu2_p0(i_nearP_i1,j_nearP_i1+1) - mu2_p0(i_nearP_i1,j_nearP_i1) ) /
1790  ( mu1_p0(i_nearP_i1,j_nearP_i1+1) - mu1_p0(i_nearP_i1,j_nearP_i1) ) ;
1791  bP_sup = mu2_p0(i_nearP_i1,j_nearP_i1) - aP_sup * mu1_p0(i_nearP_i1,j_nearP_i1) ;
1792 
1793  double mu2_nul_sup_Left = ztab(i1, j1, k1) ;
1794  double mu1_nul_sup_Left = (mu2_nul_sup_Left - bP_sup) / aP_sup ;
1795 
1796 
1797  // ********** point pris en bas a droite
1798 
1799  double aN_sup, bN_sup;
1800  aN_sup = (mu2_n0(i_nearN_i1,j_nearN_i1 +1) - mu2_n0(i_nearN_i1,j_nearN_i1) ) /
1801  (mu1_n0(i_nearN_i1,j_nearN_i1 +1) - mu1_n0(i_nearN_i1,j_nearN_i1) ) ;
1802  bN_sup = mu2_n0(i_nearN_i1,j_nearN_i1 ) - aN_sup * mu1_n0(i_nearN_i1,j_nearN_i1 ) ;
1803 
1804  double mu1_nul_sup_Right = ytab(i1, j1, k1) ;
1805  double mu2_nul_sup_Right= aN_sup * mu1_nul_sup_Right + bN_sup;
1806 
1807  double mu1_11, mu1_21, mu2_11, mu2_12 ;
1808  double p_11,p_21,p_12,p_22 ;
1809  double n1_11,n1_21,n1_12, n1_22, n2_11 ,n2_21, n2_12, n2_22 ;
1810  double cross_11 , cross_21, cross_12, cross_22 ;
1811  double dp2sddelta_car_11,dp2sddelta_car_21,dp2sddelta_car_12, dp2sddelta_car_22 ;
1812  double d2psdent1ddelta_car_11 ,d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22 ;
1813  double d2psdent2ddelta_car_11,d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22 ;
1814 
1815  mu1_11 = mu1_nul_sup_Left ;
1816  mu1_21 = mu1_nul_sup_Right ;
1817  mu2_11 = mu2_nul_sup_Right ;
1818  mu2_12 = mu2_nul_sup_Left ;
1819  p_21 = ftab(i1,j1, k_near) ;
1820  p_12 = ftab(i1,j_near, k1) ;
1821  p_22 = ftab(i1,j1, k1) ;
1822  n1_21 = dfdytab(i1,j1, k_near) ;
1823  n1_12 = dfdytab(i1,j_near, k1) ;
1824  n1_22 = dfdytab(i1,j1, k1) ;
1825  n2_21 = dfdztab(i1,j1, k_near) ;
1826  n2_12 = dfdztab(i1,j_near, k1) ;
1827  n2_22 = dfdztab(i1,j1, k1) ;
1828  cross_21 = d2fdydztab(i1,j1, k_near) ;
1829  cross_12 = d2fdydztab(i1,j_near, k1) ;
1830  cross_22 = d2fdydztab(i1,j1, k1) ;
1831  dp2sddelta_car_21 = dlpsddelta_car(i1,j1, k_near) ;
1832  dp2sddelta_car_12 = dlpsddelta_car(i1,j_near, k1) ;
1833  dp2sddelta_car_22 = dlpsddelta_car(i1,j1, k1) ;
1834  d2psdent1ddelta_car_21 = d2lpsdlent1ddelta_car(i1,j1, k_near) ;
1835  d2psdent1ddelta_car_12 = d2lpsdlent1ddelta_car(i1,j_near, k1) ;
1836  d2psdent1ddelta_car_22 = d2lpsdlent1ddelta_car(i1,j1, k1) ;
1837  d2psdent2ddelta_car_21 = d2lpsdlent2ddelta_car(i1,j1, k_near) ;
1838  d2psdent2ddelta_car_12 = d2lpsdlent2ddelta_car(i1,j_near, k1) ;
1839  d2psdent2ddelta_car_22 = d2lpsdlent2ddelta_car(i1,j1, k1) ;
1840  p_11 = 0.;
1841  n1_11 = 0. ;
1842  n2_11 = 0.;
1843  cross_11 = 0.;
1844  dp2sddelta_car_11 = 0. ;
1845  d2psdent1ddelta_car_11 = 0. ;
1846  d2psdent2ddelta_car_11 = 0. ;
1847 
1848 
1849  // changement
1850  interpol_herm_2d_new_avec( y, z,
1851  mu1_11, mu1_21, mu2_11, mu2_12,
1852  p_11,p_21, p_12, p_22,
1853  n1_11, n1_21, n1_12, n1_22,
1854  n2_11, n2_21, n2_12, n2_22,
1855  cross_11, cross_21, cross_12, cross_22,
1856  f_i1, dfdy_i1, dfdz_i1) ;
1857 
1858  double der1= 0., der2=0.;
1859  interpol_herm_2d_new_sans( y, z,
1860  mu1_11, mu1_21, mu2_11, mu2_12,
1861  dp2sddelta_car_11, dp2sddelta_car_21, dp2sddelta_car_12, dp2sddelta_car_22,
1862  d2psdent1ddelta_car_11, d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22,
1863  d2psdent2ddelta_car_11, d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22,
1864  alpha_i1, der1, der2) ;
1865  alpha_i1 = - alpha_i1 ;
1866 
1867  if (f_i1 < 0.) {
1868 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1869  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative pressure " << endl ;
1870 // abort();
1871  }
1872  if (dfdy_i1 < 0.) {
1873 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1874  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative density for fluid N " << endl ;
1875 // abort();
1876  }
1877  if (dfdz_i1 < 0.) {
1878 
1879 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1880  cout << " INTERPOLATION 2 FLUIDS CAS 2 --> negative density for fluid P " << endl ;
1881 // abort();
1882  }
1883 
1884 
1885 
1886  }
1887 //
1888  else if ( (dfdytab( i_near+1, j_near, k_near) == 0. ) && (dfdztab( i_near+1, j_near, k_near) > 0. ) &&
1889  (dfdytab( i_near+1, j_near, k1) == 0. ) && (dfdztab( i_near+1, j_near, k1) > 0. ) &&
1890  (dfdytab( i_near+1, j1, k_near) > 0. ) && (dfdztab( i_near+1, j1, k_near) > 0. ) &&
1891  (dfdytab( i_near+1, j1, k1) > 0. ) && (dfdztab( i_near+1, j1, k1) > 0. ) ) {
1892 
1893  //cout << "cas Fluide de Neutrons seul / Courbe limite - Tranche i +1 " << endl;
1894 
1895  double aP_sup, bP_sup;
1896  aP_sup = (mu2_p0(i_nearP_i+1,j_nearP_i+1) - mu2_p0(i_nearP_i+1,j_nearP_i) ) /
1897  ( mu1_p0(i_nearP_i+1,j_nearP_i+1) - mu1_p0(i_nearP_i+1,j_nearP_i) ) ;
1898  bP_sup = mu2_p0(i_nearP_i+1,j_nearP_i) - aP_sup * mu1_p0(i_nearP_i+1,j_nearP_i) ;
1899 
1900 
1901  double mu2_nul_sup = ztab(i1, j1, k1) ;
1902  double mu1_nul_sup = (mu2_nul_sup - bP_sup)/aP_sup;
1903  double mu1_11, mu1_21, mu2_11, mu2_12 ;
1904  double p_11,p_21,p_12,p_22 ;
1905  double n1_11,n1_21,n1_12, n1_22, n2_11 ,n2_21, n2_12, n2_22 ;
1906  double cross_11 , cross_21, cross_12, cross_22 ;
1907  double dp2sddelta_car_11,dp2sddelta_car_21,dp2sddelta_car_12, dp2sddelta_car_22 ;
1908  double d2psdent1ddelta_car_11 ,d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22 ;
1909  double d2psdent2ddelta_car_11,d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22 ;
1910 
1911 
1912  mu1_11 = mu1_nul_sup ;
1913  mu1_21 = ytab(i1,j1, k_near) ;
1914  mu2_11 = ztab(i1,j1, k_near) ;
1915  mu2_12 = mu2_nul_sup ;
1916  p_21 = ftab(i1,j1, k_near) ;
1917  p_12 = ftab(i1,j_near, k1) ;
1918  p_22 = ftab(i1,j1, k1) ;
1919  n1_21 = dfdytab(i1,j1, k_near) ;
1920  n1_12 = dfdytab(i1,j_near, k1) ;
1921  n1_22 = dfdytab(i1,j1, k1) ;
1922  n2_21 = dfdztab(i1,j1, k_near) ;
1923  n2_12 = dfdztab(i1,j_near, k1) ;
1924  n2_22 = dfdztab(i1,j1, k1) ;
1925  cross_21 = d2fdydztab(i1,j1, k_near) ;
1926  cross_12 = d2fdydztab(i1,j_near, k1) ;
1927  cross_22 = d2fdydztab(i1,j1, k1) ;
1928  p_11 = ftab(i1,j_near, k_near) ;
1929  n1_11 = 0.;// dfdytab(i1,j_near, k_near) ;
1930  n2_11 = dfdztab(i1,j_near, k_near) ;
1931  cross_11 = 0.; //d2fdydztab(i1,j_near, k_near) ;
1932  dp2sddelta_car_21 = dlpsddelta_car(i1,j1, k_near) ;
1933  dp2sddelta_car_12 = dlpsddelta_car(i1,j_near, k1) ;
1934  dp2sddelta_car_22 = dlpsddelta_car(i1,j1, k1) ;
1935  d2psdent1ddelta_car_21 = d2lpsdlent1ddelta_car(i1,j1, k_near) ;
1936  d2psdent1ddelta_car_12 = d2lpsdlent1ddelta_car(i1,j_near, k1) ;
1937  d2psdent1ddelta_car_22 = d2lpsdlent1ddelta_car(i1,j1, k1) ;
1938  d2psdent2ddelta_car_21 = d2lpsdlent2ddelta_car(i1,j1, k_near) ;
1939  d2psdent2ddelta_car_12 = d2lpsdlent2ddelta_car(i1,j_near, k1) ;
1940  d2psdent2ddelta_car_22 = d2lpsdlent2ddelta_car(i1,j1, k1) ;
1941  dp2sddelta_car_11 = 0.; //dlpsddelta_car(i1,j_near, k_near) ;
1942  d2psdent1ddelta_car_11 = 0.;//d2lpsdlent1ddelta_car(i1,j_near, k_near) ;
1943  d2psdent2ddelta_car_11 = 0.;//d2lpsdlent2ddelta_car(i1,j_near, k_near) ;
1944 
1945  interpol_herm_2d_new_avec( y, z,
1946  mu1_11, mu1_21, mu2_11,mu2_12,
1947  p_11,p_21, p_12, p_22,
1948  n1_11, n1_21, n1_12, n1_22,
1949  n2_11, n2_21, n2_12, n2_22,
1950  cross_11, cross_21, cross_12, cross_22,
1951  f_i1, dfdy_i1, dfdz_i1) ;
1952 
1953  double der1= 0., der2=0.;
1954  interpol_herm_2d_new_sans( y, z,
1955  mu1_11, mu1_21, mu2_11, mu2_12,
1956  dp2sddelta_car_11, dp2sddelta_car_21, dp2sddelta_car_12, dp2sddelta_car_22,
1957  d2psdent1ddelta_car_11, d2psdent1ddelta_car_21, d2psdent1ddelta_car_12, d2psdent1ddelta_car_22,
1958  d2psdent2ddelta_car_11, d2psdent2ddelta_car_21, d2psdent2ddelta_car_12, d2psdent2ddelta_car_22,
1959  alpha_i1, der1, der2) ;
1960  alpha_i1 = - alpha_i1 ;
1961 
1962  if (f_i1 < 0.) {
1963 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1964  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative pressure " << endl ;
1965 // abort();
1966  }
1967  if (dfdy_i1 < 0.) {
1968 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1969  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative density for fluid N " << endl ;
1970 // abort();
1971  }
1972  if (dfdz_i1 < 0.) {
1973 
1974 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1975  cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative density for fluid P " << endl ;
1976 // abort();
1977  }
1978 
1979 }
1980 
1981  else { */
1982 
1983 
1984 
1985  //cout << "TRAITEMENT DE LA SURFACE PAS ENCORE REALISE - Tranche i + 1 " << endl;
1986  interpol_mixed_3d(xtab, ytab, ztab, ftab,
1987  dfdytab, dfdztab, d2fdydztab,
1988  xtab(i1, j_near, k_near), y, z, f_i1, dfdy_i1 , dfdz_i1 ) ;
1989 
1990  double der1b = 0., der2b = 0.;
1991  interpol_mixed_3d_mod(xtab, ytab, ztab, dlpsddelta_car,
1992  d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car,
1993  xtab(i1, j_near, k_near), y, z, alpha_i1, der1b , der2b ) ;
1994  alpha_i1 = - alpha_i1 ;
1995  //cout << "surface non 6" << endl;
1996  // Zone en plein milieu de la table a deux fluides
1997 /*
1998  if (f_i1 < 0.) {
1999 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
2000  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative pressure " << endl ;
2001 // abort();
2002  }
2003  if (dfdy_i1 < 0.) {
2004 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
2005  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative density for fluid N " << endl ;
2006 // abort();
2007  }
2008  if (dfdz_i1 < 0.) {
2009 
2010 // cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
2011  cout << " INTERPOLATION 2 FLUIDS CAS NON TRAITE --> negative density for fluid P " << endl ;
2012 // abort();
2013  }
2014 */
2015 
2016  /*
2017  }*/
2018 //
2019 }
2020 
2021 
2022  //--------------------------------------------
2023  // ----- Linear Interpolation in Delta2 ------
2024  // -------------------------------------------
2025 
2026  double x1 = xtab(i_near, 0, 0) ;
2027  double x2 = xtab(i1, 0, 0) ;
2028  double x12 = x1-x2 ;
2029 
2030  //for f
2031  double y1 = f_i_near;
2032  double y2 = f_i1;
2033  double a = (y1-y2)/x12 ;
2034  double b = (x1*y2-y1*x2)/x12 ;
2035 
2036  f = x*a+b ;
2037 
2038  //for df/dy
2039  double y1_y = dfdy_i_near;
2040  double y2_y = dfdy_i1;
2041  double a_y = (y1_y-y2_y)/x12 ;
2042  double b_y = (x1*y2_y-y1_y*x2)/x12 ;
2043 
2044  dfdy = x*a_y+b_y ;
2045 
2046  //for df/dz
2047  double y1_z = dfdz_i_near;
2048  double y2_z = dfdz_i1;
2049  double a_z = (y1_z-y2_z)/x12 ;
2050  double b_z = (x1*y2_z-y1_z*x2)/x12 ;
2051 
2052  dfdz = x*a_z+b_z ;
2053 
2054  // for alpha
2055  double y1_alpha = alpha_i_near;
2056  double y2_alpha = alpha_i1;
2057  double a_alpha = (y1_alpha-y2_alpha)/x12 ;
2058  double b_alpha = (x1*y2_alpha-y1_alpha*x2)/x12 ;
2059 
2060  alpha = x*a_alpha+b_alpha ;
2061 
2062  // A vérifier mais on devrait pouvoir enlever ca :
2063  /*if (( y <= m_1)&& ( z <= m_2)) {
2064  alpha = 0. ;
2065  f = 0.;
2066  dfdy = 0.;
2067  dfdz = 0.;
2068  } */
2069 
2070  return;
2071 }
2072 
2073 
2074 } // End of namespace Lorene
2075 
Lorene
Lorene prototypes.
Definition: app_hor.h:67