diff options
Diffstat (limited to 'src/interpolate.c')
| -rw-r--r-- | src/interpolate.c | 79 |
1 files changed, 69 insertions, 10 deletions
diff --git a/src/interpolate.c b/src/interpolate.c index 1d774a2..e207cb2 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -23,11 +23,13 @@ CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x); + CCTK_REAL x, + CCTK_INT excised_points[]); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x); + CCTK_REAL x, + CCTK_INT excised_points[]); /* prototypes for private functions defined in this file */ static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, @@ -75,14 +77,70 @@ static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x) + CCTK_REAL x, + CCTK_INT excised_points[]) { - switch (order) { - case 1: return interpolate_local_order1(x0, dx, y, x); - case 2: return interpolate_local_order2(x0, dx, y, x); - case 3: return interpolate_local_order3(x0, dx, y, x); - case 4: return interpolate_local_order4(x0, dx, y, x); - case 5: return interpolate_local_order5(x0, dx, y, x); + + /* + We assume that the excised region has one of the following 4 forms: + + o o o o o (No excision) + x o o o o (left excision) + o o o o x (right excision) + x o o o x (left and right excision) + + Anything else (e.g., x o x o x) will go horribly wrong. + */ + + int start, end, available_order; + + start = 0; + end = order; + + while ((start < order + 1)&&(excised_points[start])) + { + ++start; + } + while ((end > -1)&&(excised_points[end])) + { + --end; + } + + if (end < start) + { + /* The whole block is excised */ + return y[0]; + } + + /* The possible interpolation order that can be used */ + + available_order = end - start; + + if (available_order == 0) + { + /* There is only one non-excised point */ + return y[start]; + } + + if (dx * start + x0 > x) + { + /* We would have to extrapolate to the left */ + return y[start]; + } + + if (dx * end + x0 < x) + { + /* We would have to extrapolate to the right */ + return y[end]; + } + + switch (available_order) + { + case 1: return interpolate_local_order1(x0 + start * dx, dx, &y[start], x); + case 2: return interpolate_local_order2(x0 + start * dx, dx, &y[start], x); + case 3: return interpolate_local_order3(x0 + start * dx, dx, &y[start], x); + case 4: return interpolate_local_order4(x0 + start * dx, dx, &y[start], x); + case 5: return interpolate_local_order5(x0 + start * dx, dx, &y[start], x); default: CCTK_WARN(0, "Interpolation order not supported"); } } @@ -218,7 +276,8 @@ return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], - CCTK_REAL x) + CCTK_REAL x, + CCTK_INT excised_points[]) { /* CCTK_REAL diff[4]; */ |