Zeroth attempt at doing excision with Cartoon. Now requires the

Spacemask thorn from CactusEinstein.

Initial tests suggest perfect second order convergence but not quite
perfect fourth order convergence, at least at reasonable
resolutions. Convergence is still better than 3rd order, though.

Excision with the new spacemask has been added but not tested.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@97 eec4d7dc-71c2-46d6-addf-10296150bf52

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]; */