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