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 $" ;
47void interpol_herm(
const Tbl&,
const Tbl&,
const Tbl&,
double,
int&,
double&,
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)
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) ;
76 nbp1 = xtab.get_dim(2) ;
77 nbp2 = xtab.get_dim(1) ;
78 nbp3 = xtab.get_dim(0) ;
96 while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
103 while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
110 while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
117 int i1 = i_near + 1 ;
118 int j1 = j_near + 1 ;
119 int k1 = k_near + 1 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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) ;
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;
271 while ( ( ytab(i1,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
278 while ( ( ztab( i1, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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) ;
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;
387 double x1 = xtab(i_near, 0, 0) ;
388 double x2 = xtab(i1, 0, 0) ;
393 double y1 = f_i_near;
395 double a = (y1-y2)/x12 ;
396 double b = (x1*y2-y1*x2)/x12 ;
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 ;
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 ;
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)
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) ;
441 int nbp1, nbp2, nbp3;
442 nbp1 = xtab.get_dim(2) ;
443 nbp2 = xtab.get_dim(1) ;
444 nbp3 = xtab.get_dim(0) ;
461 while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
468 while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
475 while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
482 int i1 = i_near + 1 ;
483 int j1 = j_near + 1 ;
484 int k1 = k_near + 1 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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) ;
631 while ( ( ytab(i1,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
638 while ( ( ztab( i1, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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) ;
732 double x1 = xtab(i_near, 0, 0) ;
733 double x2 = xtab(i1, 0, 0) ;
738 double y1 = f_i_near;
740 double a = (y1-y2)/x12 ;
741 double b = (x1*y2-y1*x2)/x12 ;
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 ;
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 ;
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) {
778 double dy = mu1_21 - mu1_11 ;
779 double dz = mu2_12 - mu2_11;
781 double u = (y - mu1_11) / dy ;
782 double v = (z - mu2_11) / dz ;
784 double u2 = u*u ;
double v2 = v*v ;
785 double u3 = u2*u ;
double v3 = v2*v ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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) ;
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 ;
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) {
877 double dy = mu1_21 - mu1_11 ;
878 double dz = mu2_12 - mu2_11;
880 double u = (y - mu1_11) / dy ;
881 double v = (z - mu2_11) / dz ;
883 double u2 = u*u ;
double v2 = v*v ;
884 double u3 = u2*u ;
double v3 = v2*v ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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) ;
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)
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) ;
987 int nbp1, nbp2, nbp3;
988 nbp1 = xtab.get_dim(2) ;
989 nbp2 = xtab.get_dim(1) ;
990 nbp3 = xtab.get_dim(0) ;
997 while ( ( xtab(i_near,0,0) <= x ) && ( ( nbp1-1 ) > i_near ) ) {
1004 while ( ( ytab(i_near,j_near, 0) <= y ) && ( ( nbp2-1 ) > j_near ) ) {
1011 while ( ( ztab( i_near, j_near, k_near) <= z) && ( ( nbp3-1 ) > k_near ) ) {
1018 int i1 = i_near + 1 ;
1019 int j1 = j_near + 1 ;
1020 int k1 = k_near + 1 ;
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;
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;
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;
1047 double f_i_near =0. ;
1048 double dfdy_i_near =0. ;
1049 double dfdz_i_near =0. ;
1050 double alpha_i_near = 0.;
1052 double dfdy_i1 =0. ;
1053 double dfdz_i1 =0. ;
1054 double alpha_i1 = 0.;
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) ;
1077 while ( ( (delta_car_n0)(i_nearN_i,0) <= xtab(i_near, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i ) ) {
1080 if (i_nearN_i != 0) {
1084 while ( ( (mu1_n0)(i_nearN_i,j_nearN_i) <= y ) && ( ( n_mu1N-1 ) > j_nearN_i ) ) {
1087 if (j_nearN_i != 0) {
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 ;
1140 while ( ( (delta_car_p0)(i_nearP_i,0) <= xtab(i_near, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i ) ) {
1143 if (i_nearP_i != 0) {
1148 while ( ( (mu2_p0)(i_nearP_i,j_nearP_i) <= z ) && ( ( n_mu2P-1 ) > j_nearP_i ) ) {
1151 if (j_nearP_i != 0) {
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 ;
1207 else if (Placei == 1 ) {
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 ;
1219 if (dfdy_i_near < 0.) {
1220 cout <<
" INTERPOLATION FLUID N --> negative density " << endl ;
1225 else if (Placei == 2 ) {
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 ;
1237 if (dfdz_i_near < 0.) {
1238 cout <<
" INTERPOLATION FLUID P --> negative density " << endl ;
1243 else if (Placei == 0 ) {
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) ;
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 ;
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 ;
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 ;
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) ;
1472 //cout << "f_i_near " << f_i_near << " dfdy_i_near " << dfdy_i_near << " dfdz_i_near " << dfdz_i_near << endl ;
1474 double der1= 0., der2=0.;
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 ;
1486 if (f_i_near < 0.) {
1487// cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1488 cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative pressure " << endl ;
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 ;
1496 if (dfdz_i_near < 0.) {
1498// cout << " DEBUG : FAR FROM LIMIT CURVES " << endl;
1499 cout << " INTERPOLATION 2 FLUIDS CAS 3 --> negative density for fluid P " << 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 ) ;
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 ) ;
1521 alpha_i_near = - alpha_i_near ;
1573 while ( ( (delta_car_n0)(i_nearN_i1,0) <= xtab(i_near+1, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i1 ) ) {
1576 if (i_nearN_i1 != 0) {
1580 while ( ( (mu1_n0)(i_nearN_i1,j_nearN_i1) <= y ) && ( ( n_mu1N-1 ) > j_nearN_i1 ) ) {
1583 if (j_nearN_i1 != 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 ;
1632 while ( ( (delta_car_p0)(i_nearP_i1,0) <= xtab(i_near+1, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i1) ) {
1635 if (i_nearP_i1 != 0) {
1639 while ( ( (mu2_p0)(i_nearP_i1,j_nearP_i1) <= z ) && ( ( n_mu2P-1 ) > j_nearP_i1 ) ) {
1642 if (j_nearP_i1 != 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 ;
1692 if (Placei1 == 3 ) {
1698 else if (Placei1 == 1 ) {
1703 interpol_herm(mu1_N, press_N, n_n_N,
1704 y, i, f_i1 , dfdy_i1 ) ;
1706 cout <<
" INTERPOLATION FLUID N i+1 --> negative pressure " << endl ;
1711 cout <<
" INTERPOLATION FLUID N i+1--> negative density " << endl ;
1716 else if (Placei1 == 2 ) {
1721 interpol_herm( mu2_P, press_P, n_p_P,
1722 z, i, f_i1, dfdz_i1) ;
1724 cout <<
" INTERPOLATION FLUID P i+1--> negative pressure " << endl ;
1729 cout <<
" INTERPOLATION FLUID P i+1 --> negative density " << endl ;
1734else if (Placei1 == 0 ) {
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 ) ;
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 ;
2023 double x1 = xtab(i_near, 0, 0) ;
2024 double x2 = xtab(i1, 0, 0) ;
2025 double x12 = x1-x2 ;
2028 double y1 = f_i_near;
2030 double a = (y1-y2)/x12 ;
2031 double b = (x1*y2-y1*x2)/x12 ;
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 ;
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 ;
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 ;
2057 alpha = x*a_alpha+b_alpha ;