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