1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 #include <math.h>
7 #include "ProData.h"
8 #include "F.h"
9 #include "MallocUtils.h"
10 #include "Message.h"
11
12 extern struct CurrentData Current ;
13
14 /* ------------------------------------------------------------------------ */
15 /* Interpolation */
16 /* ------------------------------------------------------------------------ */
17
F_InterpolationLinear(F_ARG)18 void F_InterpolationLinear(F_ARG)
19 {
20 int N, up, lo ;
21 double xp, yp = 0., *x, *y, a ;
22 struct FunctionActive * D ;
23
24 if (!Fct->Active) Fi_InitListXY (Fct, A, V) ;
25
26 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
27 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
28
29 xp = A->Val[0] ;
30
31 if (xp < x[0]) {
32 Message::Error("Bad argument for linear interpolation (%g < %g)", xp, x[0]) ;
33 }
34 else if (xp > x[N-1]) {
35 a = (y[N-1] - y[N-2]) / (x[N-1] - x[N-2]) ;
36 yp = y[N-1] + ( xp - x[N-1] ) * a ;
37 }
38 else {
39 up = 0 ; while (x[++up] < xp){} ; lo = up - 1 ;
40 a = (y[up] - y[lo]) / (x[up] - x[lo]) ;
41 yp = y[up] + ( xp - x[up] ) * a ;
42 }
43
44 // Using interpolation function also in multi-harmonic case
45 // Input is scalar, output is also scalar, rest of Val is set to zero...
46 // Not very elegant, but was already done like that for NbrHar=2...
47 // Should we add here a warning?
48 // Message::Warning("Function 'Interpolation' interpolates only real values, i.e. result is a real scalar");
49 V->Val[0] = yp;
50 if (Current.NbrHar != 1){
51 V->Val[MAX_DIM] = 0. ;
52 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
53 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
54 }
55
56 /*
57 if (Current.NbrHar == 1)
58 V->Val[0] = yp ;
59 else if (Current.NbrHar == 2) {
60 V->Val[0] = yp ;
61 V->Val[1] = 0. ;
62 }
63 else {
64 Message::Error("Function 'Interpolation' not valid for Complex");
65 }
66 */
67
68 V->Type = SCALAR ;
69 }
70
F_dInterpolationLinear(F_ARG)71 void F_dInterpolationLinear(F_ARG)
72 {
73 int N, up, lo ;
74 double xp, dyp = 0., *x, *y ;
75 struct FunctionActive * D ;
76
77 if (!Fct->Active) Fi_InitListXY (Fct, A, V) ;
78
79 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
80 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
81
82 xp = A->Val[0] ;
83
84 if (xp < x[0]) {
85 Message::Error("Bad argument for linear Interpolation (%g < %g)", xp, x[0]) ;
86 }
87 else if (xp > x[N-1]) {
88 dyp = (y[N-1] - y[N-2]) / (x[N-1] - x[N-2]) ;
89 }
90 else {
91 up = 0 ; while (x[++up] < xp){} ; lo = up - 1 ;
92 dyp = (y[up] - y[lo]) / (x[up] - x[lo]) ;
93 }
94
95 // Allowing interpolation of a real in multi-harmonic case (to change...)
96 V->Val[0] = dyp;
97 if (Current.NbrHar != 1){
98 V->Val[MAX_DIM] = 0. ;
99 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
100 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
101 }
102
103 /*
104 if (Current.NbrHar == 1)
105 V->Val[0] = dyp ;
106 else if (Current.NbrHar == 2) {
107 V->Val[0] = dyp ;
108 V->Val[1] = 0. ;
109 }
110 else {
111 Message::Error("Function 'dInterpolation' not valid for Complex");
112 }
113 */
114 V->Type = SCALAR ;
115 }
116
F_dInterpolationLinear2(F_ARG)117 void F_dInterpolationLinear2(F_ARG)
118 {
119 int N, up, lo ;
120 double xp, yp = 0., *x, *y, a ;
121 struct FunctionActive * D ;
122
123 if (!Fct->Active) {
124 Fi_InitListXY (Fct, A, V) ;
125 Fi_InitListXY2 (Fct, A, V) ;
126 }
127
128 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
129 x = D->Case.Interpolation.xc ; y = D->Case.Interpolation.yc ;
130
131 xp = A->Val[0] ;
132
133 if (xp < x[0]) {
134 Message::Error("Bad argument for linear interpolation (%g < %g)", xp, x[0]) ;
135 }
136 else if (xp > x[N-1]) {
137 a = (y[N-1] - y[N-2]) / (x[N-1] - x[N-2]) ;
138 yp = y[N-1] + ( xp - x[N-1] ) * a ;
139 }
140 else {
141 up = 0 ; while (x[++up] < xp){} ; lo = up - 1 ;
142 a = (y[up] - y[lo]) / (x[up] - x[lo]) ;
143 yp = y[up] + ( xp - x[up] ) * a ;
144 }
145
146 // Allowing interpolation of a real in multi-harmonic case (to change...)
147 V->Val[0] = yp;
148 if (Current.NbrHar != 1){
149 V->Val[MAX_DIM] = 0. ;
150 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
151 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
152 }
153 /*
154 if (Current.NbrHar == 1)
155 V->Val[0] = yp ;
156 else if (Current.NbrHar == 2) {
157 V->Val[0] = yp ;
158 V->Val[1] = 0. ;
159 }
160 else {
161 Message::Error("Function 'dInterpolation' not valid for Complex");
162 }
163 */
164
165 V->Type = SCALAR ;
166 }
167
F_InterpolationAkima(F_ARG)168 void F_InterpolationAkima(F_ARG)
169 {
170 // Third order interpolation with slope control
171 int N, up, lo ;
172 double xp, yp = 0., *x, *y, a, a2, a3 ;
173 struct FunctionActive * D ;
174
175 if (!Fct->Active) {
176 Fi_InitListXY (Fct, A, V) ;
177 Fi_InitAkima (Fct, A, V) ;
178 }
179
180 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
181 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
182
183 xp = A->Val[0] ;
184
185 if (xp < x[0]) {
186 Message::Error("Bad argument for linear interpolation (%g < %g)", xp, x[0]) ;
187 }
188 else if (xp > x[N-1]) {
189 a = (y[N-1] - y[N-2]) / (x[N-1] - x[N-2]) ;
190 yp = y[N-1] + ( xp - x[N-1] ) * a ;
191 }
192 else {
193 up = 0 ; while (x[++up] < xp){} ; lo = up - 1 ;
194 a = xp - x[lo] ; a2 = a*a ; a3 = a2*a ;
195 yp = y[lo]
196 + D->Case.Interpolation.bi[lo] * a
197 + D->Case.Interpolation.ci[lo] * a2
198 + D->Case.Interpolation.di[lo] * a3 ;
199 }
200
201 // Allowing interpolation of a real in multi-harmonic case (to change...)
202 V->Val[0] = yp;
203 if (Current.NbrHar != 1){
204 V->Val[MAX_DIM] = 0. ;
205 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
206 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
207 }
208
209 /*
210 if (Current.NbrHar == 1)
211 V->Val[0] = yp ;
212 else if (Current.NbrHar == 2) {
213 V->Val[0] = yp ;
214 V->Val[1] = 0. ;
215 }
216 else {
217 Message::Error("Function 'InterpolationAkima' not valid for Complex");
218 }
219 */
220
221 V->Type = SCALAR ;
222 }
223
F_dInterpolationAkima(F_ARG)224 void F_dInterpolationAkima(F_ARG)
225 {
226 int N, up, lo ;
227 double xp, dyp = 0., *x, *y, a, a2 ;
228 struct FunctionActive * D ;
229
230 if (!Fct->Active) {
231 Fi_InitListXY (Fct, A, V) ;
232 Fi_InitAkima (Fct, A, V) ;
233 }
234
235 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
236 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
237
238 xp = A->Val[0] ;
239
240 if (xp < x[0]) {
241 Message::Error("Bad argument for linear interpolation (%g < %g)", xp, x[0]) ;
242 }
243 else if (xp > x[N-1]) {
244 dyp = (y[N-1] - y[N-2]) / (x[N-1] - x[N-2]) ;
245 }
246 else {
247 up = 0 ; while (x[++up] < xp){} ; lo = up - 1 ;
248 a = xp - x[lo] ; a2 = a*a ;
249 dyp = D->Case.Interpolation.bi[lo]
250 + D->Case.Interpolation.ci[lo] * 2. * a
251 + D->Case.Interpolation.di[lo] * 3. * a2 ;
252 }
253
254
255 // Extension to multi-harmonic case (to change...)
256 V->Val[0] = dyp;
257 if (Current.NbrHar != 1){
258 V->Val[MAX_DIM] = 0. ;
259 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
260 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
261 }
262 /*
263 if (Current.NbrHar == 1)
264 V->Val[0] = dyp ;
265 else if (Current.NbrHar == 2) {
266 V->Val[0] = dyp ;
267 V->Val[1] = 0. ;
268 }
269 else {
270 Message::Error("Function 'dInterpolationAkima' not valid for Complex");
271 }
272 */
273 V->Type = SCALAR ;
274 }
275
276
Fi_InterpolationBilinear(double * x,double * y,double * M,int NL,int NC,double xp,double yp,double * zp)277 bool Fi_InterpolationBilinear(double *x, double *y, double *M, int NL, int NC,
278 double xp, double yp, double *zp)
279 {
280 double a11, a12, a21, a22;
281 int i, j;
282
283 // Interpolate point (xp,yp) in a regular grid
284 // x[i] <= xp < x[i+1]
285 // y[j] <= yp < y[j+1]
286
287 *zp = 0.0 ;
288
289 // When (xp,yp) lays outside the boundaries of the table:
290 // the nearest border is taken
291 if (xp < x[0]) xp = x[0];
292 else if (xp > x[NL-1]) xp = x[NL-1];
293 for (i=0 ; i<NL-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
294 i = (i >= NL) ? NL-1 : i;
295
296 if (yp < y[0]) yp = y[0];
297 else if (yp > y[NC-1]) yp = y[NC-1];
298 for (j=0 ; j<NC-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
299 j = (j >= NC) ? NC-1 : j;
300
301
302 a11 = M[ i + NL * j ];
303 a21 = M[(i+1) + NL * j ];
304 a12 = M[ i + NL * (j+1)];
305 a22 = M[(i+1) + NL * (j+1)];
306
307 *zp = 1/((x[i+1]-x[i])*(y[j+1]-y[j])) *
308 ( a11 * ( x[i+1]-xp) * ( y[j+1]-yp) +
309 a21 * (-x[i ]+xp) * ( y[j+1]-yp) +
310 a12 * ( x[i+1]-xp) * (-y[j ]+yp) +
311 a22 * (-x[i ]+xp) * (-y[j ]+yp) );
312
313 return true ;
314 }
315
316
Fi_dInterpolationBilinear(double * x,double * y,double * M,int NL,int NC,double xp,double yp,double * dzp_dx,double * dzp_dy)317 bool Fi_dInterpolationBilinear(double *x, double *y, double *M, int NL, int NC,
318 double xp, double yp, double *dzp_dx, double *dzp_dy)
319 {
320 double a11, a12, a21, a22;
321 int i, j;
322
323 // When (xp,yp) lays outside the boundaries of the table:
324 // the nearest border is taken
325 if (xp < x[0]) xp = x[0];
326 else if (xp > x[NL-1]) xp = x[NL-1];
327 for (i=0 ; i<NL-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
328 i = (i >= NL) ? NL-1 : i;
329
330 if (yp < y[0]) yp = y[0];
331 else if (yp > y[NC-1]) yp = y[NC-1];
332 for (j=0 ; j<NC-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
333 j = (j >= NC) ? NC-1 : j;
334
335
336 a11 = M[ i + NL * j ];
337 a21 = M[(i+1) + NL * j ];
338 a12 = M[ i + NL * (j+1)];
339 a22 = M[(i+1) + NL * (j+1)];
340
341 *dzp_dx = 1/((x[i+1]-x[i])*(y[j+1]-y[j])) *
342 ( (a21-a11) * ( y[j+1]-yp) +
343 (a22-a12) * (-y[j ]+yp) );
344
345 *dzp_dy = 1/((x[i+1]-x[i])*(y[j+1]-y[j])) *
346 ( (a12-a11) * ( x[i+1]-xp) +
347 (a22-a21) * (-x[i ]+xp) );
348
349 return true ;
350 }
351
Fi_InterpolationTrilinear(double * x,double * y,double * z,double * M,int NX,int NY,int NZ,double xp,double yp,double zp,double * vp)352 bool Fi_InterpolationTrilinear (double *x, double *y, double *z, double *M, int NX, int NY, int NZ,
353 double xp, double yp, double zp, double *vp)
354 {
355 /*
356 *
357 * a122 **************************** a222
358 * * | * *
359 * * | * *
360 * * | * *
361 * **************************** a212 *
362 * * a112 | * *
363 * * | * *
364 * * | * *
365 * * | * *
366 * * | * *
367 * * | * *
368 * * | * *
369 * * | * *
370 * * | * *
371 * * a121 -------------------*-------* a221
372 * * / * *
373 * * / * *
374 * * / * *
375 * ****************************
376 * a111 a211
377 */
378
379 double a111, a121, a211, a221, a112, a122, a212, a222;
380 int i, j, k;
381
382 // Interpolate point (xp,yp,zp) in a regular grid
383 // x[i] <= xp < x[i+1]
384 // y[j] <= yp < y[j+1]
385 // z[k] <= zp < z[k+1]
386
387 *vp = 0.0 ;
388
389 // When (xp,yp,zp) lays outside the boundaries of the table:
390 // the nearest border is taken
391 if (xp < x[0]) xp = x[0];
392 else if (xp > x[NX-1]) xp = x[NX-1];
393 for (i=0 ; i<NX-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
394 i = (i >= NX) ? NX-1 : i;
395
396 if (yp < y[0]) yp = y[0];
397 else if (yp > y[NY-1]) yp = y[NY-1];
398 for (j=0 ; j<NY-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
399 j = (j >= NY) ? NY-1 : j;
400
401 if (zp < z[0]) zp = z[0];
402 else if (zp > z[NZ-1]) zp = z[NZ-1];
403 for (k=0 ; k<NZ-1 ; ++k) if (z[k+1] >= zp && zp >= z[k]) break;
404 k = (k >= NZ) ? NZ-1 : k;
405
406
407 a111 = M[ i + NX * j + NX * NY * k ];
408 a211 = M[(1+i) + NX * j + NX * NY * k ];
409 a121 = M[ i + NX * (1+j) + NX * NY * k ];
410 a221 = M[(1+i) + NX * (1+j) + NX * NY * k ];
411
412 a112 = M[ i + NX * j + NX * NY * (k+1)];
413 a212 = M[(1+i) + NX * j + NX * NY * (k+1)];
414 a122 = M[ i + NX * (1+j) + NX * NY * (k+1)];
415 a222 = M[(1+i) + NX * (1+j) + NX * NY * (k+1)];
416
417 double xd, yd, zd;
418 xd = (xp-x[i])/(x[i+1]-x[i]);
419 yd = (yp-y[j])/(y[j+1]-y[j]);
420 zd = (zp-z[k])/(z[k+1]-z[k]);
421
422 double a11, a12, a21, a22;
423 a11 = a111*(1-xd) + a211*xd;
424 a12 = a112*(1-xd) + a212*xd;
425 a21 = a121*(1-xd) + a221*xd;
426 a22 = a122*(1-xd) + a222*xd;
427
428 double a1, a2;
429 a1 = a11*(1-yd) + a21*yd;
430 a2 = a12*(1-yd) + a22*yd;
431
432 *vp = a1*(1-zd) + a2*zd;
433
434 return true ;
435 }
436
Fi_dInterpolationTrilinear(double * x,double * y,double * z,double * M,int NX,int NY,int NZ,double xp,double yp,double zp,double * dvp_dx,double * dvp_dy,double * dvp_dz)437 bool Fi_dInterpolationTrilinear (double *x, double *y, double *z, double *M, int NX, int NY, int NZ,
438 double xp, double yp, double zp, double *dvp_dx, double *dvp_dy, double *dvp_dz)
439 {
440 /*
441 *
442 * a122 **************************** a222
443 * * | * *
444 * * | * *
445 * * | * *
446 * **************************** a212 *
447 * * a112 | * *
448 * * | * *
449 * * | * *
450 * * | * *
451 * * | * *
452 * * | * *
453 * * | * *
454 * * | * *
455 * * | * *
456 * * a121 -------------------*-------* a221
457 * * / * *
458 * * / * *
459 * * / * *
460 * ****************************
461 * a111 a211
462 */
463
464 double a111, a121, a211, a221, a112, a122, a212, a222;
465 int i, j, k;
466
467 // Interpolate point (xp,yp,zp) in a regular grid
468 // x[i] <= xp < x[i+1]
469 // y[j] <= yp < y[j+1]
470 // z[k] <= zp < z[k+1]
471
472 *dvp_dx = 0.0 ;
473 *dvp_dy = 0.0 ;
474 *dvp_dz = 0.0 ;
475
476 // When (xp,yp,zp) lays outside the boundaries of the table:
477 // the nearest border is taken
478 if (xp < x[0]) xp = x[0];
479 else if (xp > x[NX-1]) xp = x[NX-1];
480 for (i=0 ; i<NX-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
481 i = (i >= NX) ? NX-1 : i;
482
483 if (yp < y[0]) yp = y[0];
484 else if (yp > y[NY-1]) yp = y[NY-1];
485 for (j=0 ; j<NY-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
486 j = (j >= NY) ? NY-1 : j;
487
488 if (zp < z[0]) zp = z[0];
489 else if (zp > z[NZ-1]) zp = z[NZ-1];
490 for (k=0 ; k<NZ-1 ; ++k) if (z[k+1] >= zp && zp >= z[j]) break;
491 k = (k >= NZ) ? NZ-1 : k;
492
493
494 a111 = M[ i + NX * j + NX * NY * k ];
495 a211 = M[(1+i) + NX * j + NX * NY * k ];
496 a121 = M[ i + NX * (1+j) + NX * NY * k ];
497 a221 = M[(1+i) + NX * (1+j) + NX * NY * k ];
498
499 a112 = M[ i + NX * j + NX * NY * (k+1)];
500 a212 = M[(1+i) + NX * j + NX * NY * (k+1)];
501 a122 = M[ i + NX * (1+j) + NX * NY * (k+1)];
502 a222 = M[(1+i) + NX * (1+j) + NX * NY * (k+1)];
503
504 double xd, yd, zd;
505 xd = (xp-x[i])/(x[i+1]-x[i]);
506 yd = (yp-y[j])/(y[j+1]-y[j]);
507 zd = (zp-z[k])/(z[k+1]-z[k]);
508
509 double dxd, dyd, dzd;
510 dxd = 1./(x[i+1]-x[i]);
511 dyd = 1./(y[j+1]-y[j]);
512 dzd = 1./(z[k+1]-z[k]);
513
514 double a11, a12, a21, a22, dxa11, dxa12, dxa21, dxa22;
515 a11 = a111*(1-xd) + a211*xd;
516 a12 = a112*(1-xd) + a212*xd;
517 a21 = a121*(1-xd) + a221*xd;
518 a22 = a122*(1-xd) + a222*xd;
519 dxa11 = -a111*dxd + a211*dxd;
520 dxa12 = -a112*dxd + a212*dxd;
521 dxa21 = -a121*dxd + a221*dxd;
522 dxa22 = -a122*dxd + a222*dxd;
523
524 double a1, a2, dya1, dya2, dxa1, dxa2;
525 a1 = a11*(1-yd) + a21*yd;
526 a2 = a12*(1-yd) + a22*yd;
527 dya1 = -a11*dyd + a21*dyd;
528 dya2 = -a12*dyd + a22*dyd;
529 dxa1 = dxa11*(1-yd) + dxa21*yd;
530 dxa2 = dxa12*(1-yd) + dxa22*yd;
531
532 *dvp_dx = dxa1*(1-zd) + dxa2*zd;
533 *dvp_dy = dya1*(1-zd) + dya2*zd;
534 *dvp_dz = -a1*dzd + a2*dzd;
535
536 return true ;
537 }
538
F_InterpolationBilinear(F_ARG)539 void F_InterpolationBilinear(F_ARG)
540 {
541 /*
542 It performs a bilinear interpolation at point (xp,yp) based
543 on a two-dimensional table (sorted grid).
544
545 Input parameters:
546 NL Number of lines
547 NC Number of columns
548 x values (ascending order) linked to the NL lines of the table
549 y values (ascending order) linked to the NC columns of the table
550 M Matrix M(x,y) = M[x+NL*y]
551
552 xp x coordinate of interpolation point
553 yp y coordinate of interpolation point
554
555 R. Scorretti
556 */
557
558 int NL, NC;
559 double xp, yp, zp = 0., *x, *y, *M;
560 struct FunctionActive * D;
561
562 if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR)
563 Message::Error("Two Scalar arguments required!");
564
565 if (!Fct->Active) Fi_InitListMatrix (Fct, A, V) ;
566
567 D = Fct->Active ;
568 NL = D->Case.ListMatrix.NbrLines ;
569 NC = D->Case.ListMatrix.NbrColumns ;
570
571 x = D->Case.ListMatrix.x ;
572 y = D->Case.ListMatrix.y ;
573 M = D->Case.ListMatrix.data ;
574
575 xp = (A+0)->Val[0] ;
576 yp = (A+1)->Val[0] ;
577
578 bool IsInGrid = Fi_InterpolationBilinear (x, y, M, NL, NC, xp, yp, &zp);
579 if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g)", xp, yp) ;
580
581 V->Type = SCALAR ;
582 V->Val[0] = zp ;
583 }
584
F_dInterpolationBilinear(F_ARG)585 void F_dInterpolationBilinear(F_ARG)
586 {
587 /*
588 It delivers the derivative of the bilinear interpolation at point (xp, yp)
589 based on a two-dimensional table (sorted grid).
590
591 Input parameters:
592 NL Number of lines
593 NC Number of columns
594 x values (ascending order) linked to the NL lines of the table
595 y values (ascending order) linked to the NC columns of the table
596 M Matrix M(x,y) = M[x+NL*y]
597
598 xp x coordinate of interpolation point
599 yp y coordinate of interpolation point
600
601 */
602
603 int NL, NC;
604 double xp, yp, dzp_dx = 0., dzp_dy = 0., *x, *y, *M;
605 struct FunctionActive * D;
606
607 if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR)
608 Message::Error("Two Scalar arguments required!");
609
610 if (!Fct->Active) Fi_InitListMatrix (Fct, A, V) ;
611
612 D = Fct->Active ;
613 NL = D->Case.ListMatrix.NbrLines ;
614 NC = D->Case.ListMatrix.NbrColumns ;
615
616 x = D->Case.ListMatrix.x ;
617 y = D->Case.ListMatrix.y ;
618 M = D->Case.ListMatrix.data ;
619
620 xp = (A+0)->Val[0] ;
621 yp = (A+1)->Val[0] ;
622
623 bool IsInGrid = Fi_dInterpolationBilinear (x, y, M, NL, NC, xp, yp, &dzp_dx, &dzp_dy);
624 if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g)", xp, yp) ;
625
626 V->Type = VECTOR ;
627 V->Val[0] = dzp_dx ;
628 V->Val[1] = dzp_dy ;
629 V->Val[2] = 0. ;
630 }
631
F_InterpolationTrilinear(F_ARG)632 void F_InterpolationTrilinear(F_ARG)
633 {
634 /*
635 It performs a trilinear interpolation at point (xp,yp,zp) based
636 on a three-dimensional table (sorted grid).
637
638 Input parameters:
639 NX Number of lines X
640 NY Number of lines Y
641 NZ Number of lines Z
642 x values (ascending order) linked to the NX lines of the table
643 y values (ascending order) linked to the NY lines of the table
644 z values (ascending order) linked to the NZ lines of the table
645 M Matrix M(x,y,z) = M[x+NX*y+NX*NY*z]
646
647 xp x coordinate of interpolation point
648 yp y coordinate of interpolation point
649 zp z coordinate of interpolation point
650
651 A. Royer
652 */
653
654 int NX, NY, NZ;
655 double xp, yp, zp, vp=0., *x, *y, *z, *M;
656 struct FunctionActive * D;
657
658 if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
659 Message::Error("Three Scalar arguments required!");
660
661 if (!Fct->Active) Fi_InitListMatrix3D (Fct, A, V) ;
662
663 D = Fct->Active ;
664 NX = D->Case.ListMatrix3D.NbrLinesX ;
665 NY = D->Case.ListMatrix3D.NbrLinesY ;
666 NZ = D->Case.ListMatrix3D.NbrLinesZ ;
667
668 x = D->Case.ListMatrix3D.x ;
669 y = D->Case.ListMatrix3D.y ;
670 z = D->Case.ListMatrix3D.z ;
671 M = D->Case.ListMatrix3D.data ;
672
673 xp = (A+0)->Val[0] ;
674 yp = (A+1)->Val[0] ;
675 zp = (A+2)->Val[0] ;
676
677 bool IsInGrid = Fi_InterpolationTrilinear (x, y, z, M, NX, NY, NZ, xp, yp, zp, &vp);
678 if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g, zp=%g)", xp, yp, zp) ;
679
680 V->Type = SCALAR ;
681 V->Val[0] = vp ;
682 }
683
F_dInterpolationTrilinear(F_ARG)684 void F_dInterpolationTrilinear(F_ARG)
685 {
686 /*
687 It delivers the derivative of the bilinear interpolation at point (xp, yp, zp)
688 based on a three-dimensional table (sorted grid).
689
690 Input parameters:
691 NX Number of lines X
692 NY Number of lines Y
693 NZ Number of lines Z
694 x values (ascending order) linked to the NX lines of the table
695 y values (ascending order) linked to the NY lines of the table
696 z values (ascending order) linked to the NZ lines of the table
697 M Matrix M(x,y,z) = M[x+NX*y+NX*NY*z]
698
699 xp x coordinate of interpolation point
700 yp y coordinate of interpolation point
701 zp z coordinate of interpolation point
702
703 A. Royer
704 */
705
706 int NX, NY, NZ;
707 double xp, yp, zp, dvp_dx =0., dvp_dy =0., dvp_dz =0., *x, *y, *z, *M;
708 struct FunctionActive * D;
709
710 if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
711 Message::Error("Three Scalar arguments required!");
712
713 if (!Fct->Active) Fi_InitListMatrix3D (Fct, A, V) ;
714
715 D = Fct->Active ;
716 NX = D->Case.ListMatrix3D.NbrLinesX ;
717 NY = D->Case.ListMatrix3D.NbrLinesY ;
718 NZ = D->Case.ListMatrix3D.NbrLinesZ ;
719
720 x = D->Case.ListMatrix3D.x ;
721 y = D->Case.ListMatrix3D.y ;
722 z = D->Case.ListMatrix3D.z ;
723 M = D->Case.ListMatrix3D.data ;
724
725 xp = (A+0)->Val[0] ;
726 yp = (A+1)->Val[0] ;
727 zp = (A+2)->Val[0] ;
728
729 bool IsInGrid = Fi_dInterpolationTrilinear (x, y, z, M, NX, NY, NZ, xp, yp, zp, &dvp_dx, &dvp_dy, &dvp_dz);
730 if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g, zp=%g)", xp, yp, zp) ;
731
732 V->Type = VECTOR ;
733 V->Val[0] = dvp_dx ;
734 V->Val[1] = dvp_dy ;
735 V->Val[2] = dvp_dz ;
736 }
737
Fi_InitListMatrix(F_ARG)738 void Fi_InitListMatrix(F_ARG)
739 {
740 int i=0, k, NL, NC, sz ;
741 struct FunctionActive * D ;
742
743 /*
744 The original table structure:
745 | y(1) y(2) ... y(NC)
746 ------+--------------------------------------------
747 x(1) | data(1) data(NL+1) ... .
748 x(2) | data(2) data(NL+2) .
749 . . . .
750 . . .
751 x(NL) | data(NL) data(2*NL) ... data(NL*NC)
752
753 is furnished with the following format:
754 [ NL, NC, x(1..NL), y(1..NC), data(1..NL*NC) ]
755
756 R. Scorretti
757 */
758
759 D = Fct->Active =
760 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
761
762 NL = Fct->Para[i++];
763 NC = Fct->Para[i++];
764
765 sz = 2 + NL + NC + NL*NC ; // expected size of list matrix
766 if (Fct->NbrParameters != sz)
767 Message::Error("Bad size of input data (expected = %d ; found = %d). "
768 "List with format: x(NbrLines=%d), y(NbrColumns=%d), "
769 "matrix(NbrLines*NbrColumns=%d)",
770 sz, Fct->NbrParameters, NL, NC, NL*NC);
771
772 // Initialize structure and allocate memory
773 D->Case.ListMatrix.NbrLines = NL;
774 D->Case.ListMatrix.NbrColumns = NC;
775 D->Case.ListMatrix.x = (double *) malloc (sizeof(double)*NL);
776 D->Case.ListMatrix.y = (double *) malloc (sizeof(double)*NC);
777 D->Case.ListMatrix.data = (double *) malloc (sizeof(double)*NL*NC);
778
779 // Assign values
780 for (k=0 ; k<NL ; ++k) D->Case.ListMatrix.x[k] = Fct->Para[i++];
781 for (k=0 ; k<NC ; ++k) D->Case.ListMatrix.y[k] = Fct->Para[i++];
782 for (k=0 ; k<NL*NC ; ++k) D->Case.ListMatrix.data[k] = Fct->Para[i++];
783 }
784
Fi_InitListMatrix3D(F_ARG)785 void Fi_InitListMatrix3D(F_ARG)
786 {
787 int i=0, k, NX, NY, NZ, sz ;
788 struct FunctionActive * D ;
789
790 /*
791 The original table structure data(i,j,k) is furnished with the following format:
792 [ NX, NY, NZ, x(1..NX), y(1..NY), y(1..NZ), data(1..NX*NY*NZ) ]
793
794 A. Royer
795 */
796
797 D = Fct->Active =
798 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
799
800 NX = Fct->Para[i++];
801 NY = Fct->Para[i++];
802 NZ = Fct->Para[i++];
803
804 sz = 3 + NX + NY + NZ + NX*NY*NZ ; // expected size of list matrix
805 if (Fct->NbrParameters != sz)
806 Message::Error("Bad size of input data (expected = %d ; found = %d). "
807 "List with format: x(NbrLines=%d), y(NbrLines=%d), z(NbrLines=%d), "
808 "matrix(NbrLinesX*NbrLinesY*NbrLinesZ=%d)",
809 sz, Fct->NbrParameters, NX, NY, NZ, NX*NY*NZ);
810
811 // Initialize structure and allocate memory
812 D->Case.ListMatrix3D.NbrLinesX = NX;
813 D->Case.ListMatrix3D.NbrLinesY = NY;
814 D->Case.ListMatrix3D.NbrLinesZ = NZ;
815 D->Case.ListMatrix3D.x = (double *) malloc (sizeof(double)*NX);
816 D->Case.ListMatrix3D.y = (double *) malloc (sizeof(double)*NY);
817 D->Case.ListMatrix3D.z = (double *) malloc (sizeof(double)*NZ);
818 D->Case.ListMatrix3D.data = (double *) malloc (sizeof(double)*NX*NY*NZ);
819
820 // Assign values
821 for (k=0 ; k<NX ; ++k) D->Case.ListMatrix3D.x[k] = Fct->Para[i++];
822 for (k=0 ; k<NY ; ++k) D->Case.ListMatrix3D.y[k] = Fct->Para[i++];
823 for (k=0 ; k<NZ ; ++k) D->Case.ListMatrix3D.z[k] = Fct->Para[i++];
824 for (k=0 ; k<NX*NY*NZ ; ++k) D->Case.ListMatrix3D.data[k] = Fct->Para[i++];
825 }
826
Fi_InitListX(F_ARG)827 void Fi_InitListX(F_ARG)
828 {
829 int i, N ;
830 double *x ;
831 struct FunctionActive * D ;
832
833 D = Fct->Active =
834 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
835 N = D->Case.Interpolation.NbrPoint = Fct->NbrParameters ;
836 x = D->Case.Interpolation.x = (double *)Malloc(sizeof(double)*N) ;
837
838 for (i = 0 ; i < N ; i++)
839 x[i] = Fct->Para[i] ;
840 }
841
Fi_InitListXY(F_ARG)842 void Fi_InitListXY(F_ARG)
843 {
844 int i, N ;
845 double *x, *y ;
846 struct FunctionActive * D ;
847
848 D = Fct->Active =
849 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
850 N = D->Case.Interpolation.NbrPoint = Fct->NbrParameters / 2 ;
851 x = D->Case.Interpolation.x = (double *)Malloc(sizeof(double)*N) ;
852 y = D->Case.Interpolation.y = (double *)Malloc(sizeof(double)*N) ;
853
854 for (i = 0 ; i < N ; i++) {
855 x[i] = Fct->Para[i*2 ] ;
856 y[i] = Fct->Para[i*2+1] ;
857 }
858 }
859
Fi_InitListXY2(F_ARG)860 void Fi_InitListXY2(F_ARG)
861 {
862 int i, N ;
863 double *x, *y, *xc, *yc ;
864 struct FunctionActive * D ;
865
866 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
867 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
868 xc = D->Case.Interpolation.xc = (double *)Malloc(sizeof(double)*N) ;
869 yc = D->Case.Interpolation.yc = (double *)Malloc(sizeof(double)*N) ;
870
871 xc[0] = 0. ;
872 yc[0] = (x[1]*y[1]-x[0]*y[0]) / (x[1]*x[1]-x[0]*x[0]) ;
873
874 for (i = 1 ; i < N ; i++) {
875
876 xc[i] = 0.5 * (x[i]+x[i-1]) ;
877 yc[i] = (x[i]*y[i]-x[i-1]*y[i-1]) / (x[i]*x[i]-x[i-1]*x[i-1]) ;
878
879 /*
880 xc[i] = x[i] ;
881 yc[i] = (y[i]-y[i-1]) / (x[i]-x[i-1]) ;
882 */
883 }
884 }
885
Fi_InitAkima(F_ARG)886 void Fi_InitAkima(F_ARG)
887 {
888 int i, N ;
889 double *x, *y, *mi, *bi, *ci, *di, a ;
890 struct FunctionActive * D ;
891
892 D = Fct->Active ; N = D->Case.Interpolation.NbrPoint ;
893 x = D->Case.Interpolation.x ; y = D->Case.Interpolation.y ;
894 mi = D->Case.Interpolation.mi = (double *)Malloc(sizeof(double)*(N+4)) ;
895 mi += 2 ;
896 bi = D->Case.Interpolation.bi = (double *)Malloc(sizeof(double)*N) ;
897 ci = D->Case.Interpolation.ci = (double *)Malloc(sizeof(double)*N) ;
898 di = D->Case.Interpolation.di = (double *)Malloc(sizeof(double)*N) ;
899
900 for (i = 0 ; i < N-1 ; i++)
901 mi[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) ;
902 mi[N-1] = 2.*mi[N-2] - mi[N-3] ;
903 mi[N ] = 2.*mi[N-1] - mi[N-2] ;
904 mi[ -1] = 2.*mi[ 0] - mi[ 1] ;
905 mi[ -2] = 2.*mi[ -1] - mi[ 0] ;
906
907 for (i = 0 ; i < N ; i++)
908 if ( (a = fabs(mi[i+1]-mi[i]) + fabs(mi[i-1]-mi[i-2])) > 1.e-30 )
909 bi[i] =
910 ( fabs(mi[i+1]-mi[i]) * mi[i-1] + fabs(mi[i-1]-mi[i-2]) * mi[i] ) / a ;
911 else
912 bi[i] = (mi[i] + mi[i-1]) / 2. ;
913
914 for (i = 0 ; i < N-1 ; i++) {
915 a = (x[i+1]-x[i]) ;
916 ci[i] = ( 3.*mi[i] - 2.*bi[i] - bi[i+1] ) / a ;
917 di[i] = ( bi[i] + bi[i+1] - 2.*mi[i] ) / (a*a) ;
918 }
919 }
920
921
922 struct IntDouble { int Int; double Double; } ;
923 struct IntVector { int Int; double Double[3]; } ;
924
F_ValueFromIndex(F_ARG)925 void F_ValueFromIndex (F_ARG)
926 {
927 struct FunctionActive * D ;
928 struct IntDouble * IntDouble_P;
929
930 if (!Fct->Active) Fi_InitValueFromIndex (Fct, A, V) ;
931
932 D = Fct->Active ;
933
934 if (List_Nbr(D->Case.ValueFromIndex.Table)){
935 IntDouble_P = (struct IntDouble *)
936 List_PQuery(D->Case.ValueFromIndex.Table, &Current.NumEntity, fcmp_int);
937 if(IntDouble_P){
938 V->Val[0] = IntDouble_P->Double ;
939 V->Type = SCALAR ;
940 return;
941 }
942 }
943
944 Message::Warning("Unknown entity index %d in ValueFromIndex table",
945 Current.NumEntity);
946 V->Val[0] = 0. ;
947 V->Type = SCALAR ;
948 }
949
F_VectorFromIndex(F_ARG)950 void F_VectorFromIndex(F_ARG)
951 {
952 struct FunctionActive * D ;
953 struct IntVector * IntVector_P;
954
955 if (!Fct->Active) Fi_InitVectorFromIndex (Fct, A, V) ;
956
957 D = Fct->Active ;
958
959 if (List_Nbr(D->Case.ValueFromIndex.Table)){
960 IntVector_P = (struct IntVector *)
961 List_PQuery(D->Case.ValueFromIndex.Table, &Current.NumEntity, fcmp_int);
962 if (IntVector_P){
963 V->Val[0] = IntVector_P->Double[0] ;
964 V->Val[1] = IntVector_P->Double[1] ;
965 V->Val[2] = IntVector_P->Double[2] ;
966 V->Type = VECTOR;
967 return;
968 }
969 }
970 Message::Warning("Unknown entity index %d in VectorFromIndex table",
971 Current.NumEntity);
972 V->Val[0] = 0.;
973 V->Val[1] = 0.;
974 V->Val[2] = 0.;
975 V->Type = VECTOR;
976 }
977
Fi_InitValueFromIndex(F_ARG)978 void Fi_InitValueFromIndex(F_ARG)
979 {
980 int i, N ;
981 struct IntDouble IntDouble_s;
982 struct FunctionActive * D ;
983
984 if ((Fct->NbrParameters)){
985 N = (int)Fct->Para[0];
986 D = Fct->Active =
987 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
988 D->Case.ValueFromIndex.Table = List_Create(N + 1, 1, sizeof(struct IntDouble));
989 for (i = 0 ; i < N ; i++) {
990 IntDouble_s.Int = (int)(Fct->Para[i*2+1]+0.1);
991 IntDouble_s.Double = Fct->Para[i*2+2];
992 List_Add(D->Case.ValueFromIndex.Table, &IntDouble_s);
993 }
994 }
995 else{
996 D = Fct->Active =
997 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
998 D->Case.ValueFromIndex.Table = NULL;
999 }
1000 }
1001
Fi_InitVectorFromIndex(F_ARG)1002 void Fi_InitVectorFromIndex(F_ARG)
1003 {
1004 int i, N ;
1005 struct IntVector IntVector_s;
1006 struct FunctionActive * D ;
1007
1008 if ((Fct->NbrParameters)){
1009 N = (int)Fct->Para[0];
1010 D = Fct->Active =
1011 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
1012 D->Case.ValueFromIndex.Table = List_Create(N, 1, sizeof(struct IntVector[3]));
1013 for (i = 0 ; i < N ; i++) {
1014 IntVector_s.Int = (int)(Fct->Para[i*4+1]+0.1);
1015 IntVector_s.Double[0] = Fct->Para[i*4+2];
1016 IntVector_s.Double[1] = Fct->Para[i*4+3];
1017 IntVector_s.Double[2] = Fct->Para[i*4+4];
1018 List_Add(D->Case.ValueFromIndex.Table, &IntVector_s);
1019 }
1020 }
1021 else{
1022 D = Fct->Active =
1023 (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
1024 D->Case.ValueFromIndex.Table = NULL;
1025 }
1026 }
1027