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
28char 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"
45namespace Lorene {
46
47void 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/*
53xtab/x refers to delta_car (logdelta_car)
54ytab/y refers to mu_n (logent1)
55ztab/z refers to mu_p (logent2)
56ftab/f refers to psi (logp) or alpha (dlpsddelta_car or logalpha)
57dfdytab/dfdy refers to dpsi/dmu_n (dlpsdlent1) or dalpha/dmu_n (d2lpsdlent1ddelta_car or dlalphadlent1)
58dfdztab/dfdz refers to dpsi/dmu_p (dlpsdlent2) or dalpha/dmu_p (dl2psdlent2ddelta_car or dlalphadlent2)
59d2fdytabdztab 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/*
122bool 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) ) ) ;
123bool 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) )) ;
124if (hum1 == false ) {
125cout << setprecision(16);
126cout << 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 ;
127cout << j_near << " " << k_near << endl ;
128}
129
130if (hum2 == false ) {
131cout << setprecision(16);
132cout << 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 ;
133cout << 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
145cout << 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
149cout << setprecision(16);
150cout <<"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
159double 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
431void 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
495if (hum2 == false ) {
496cout << setprecision(16);
497cout << 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 ;
498cout << 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
510cout << y << " " << z << endl ;
511double u_i_anal = y /log(2.) ; //(y - 1.) / 1.
512double v_i_anal = z /log(2.) ; //(z - 1.) / 1.
513
514cout << setprecision(16);
515cout <<"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
525double 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
530cout << " 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
559cout << " Interpol3 " << psi0_u_i_near << " " << psi0_u_i_anal << " " << psi0_1mu_i_near << " " << psi0_1mu_i_anal << " " <<
560psi1_u_i_near << " " << psi1_u_i_anal << " " << psi1_1mu_i_near << " " << psi1_1mu_i_anal << endl;
561
562cout << " Interpol4 " << psi0_v_i_near << " " << psi0_v_i_anal << " " << psi0_1mv_i_near << " " << psi0_1mv_i_anal << " " <<
563psi1_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
769void 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
870void 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/*
958xtab/x refers to delta_car (logdelta_car)
959ytab/y refers to mu_n (logent1)
960ztab/z refers to mu_p (logent2)
961ftab/f refers to psi (logp) or alpha (dlpsddelta_car or logalpha)
962dfdytab/dfdy refers to dpsi/dmu_n (dlpsdlent1) or dalpha/dmu_n (d2lpsdlent1ddelta_car or dlalphadlent1)
963dfdztab/dfdz refers to dpsi/dmu_p (dlpsdlent2) or dalpha/dmu_p (dl2psdlent2ddelta_car or dlalphadlent2)
964d2fdytabdztab 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}
1734else 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 prototypes.
Definition app_hor.h:64