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

4 files changed, 289 insertions, 16 deletions

diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c
index 87680cd..b1417d6 100644
--- a/src/ApplyCartoon.c
+++ b/src/ApplyCartoon.c
@@ -178,7 +178,7 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS)
     {
       if (CCTK_Equals(prolongtype, "None")) 
       {
-        prolongmethod = PROLONG_LAGRANGE; /* But why? */
+        prolongmethod = PROLONG_NONE; /* But why? */
       } 
       else if (CCTK_Equals(prolongtype, "Lagrange")) 
       {
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index f74c30f..7653424 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -24,6 +24,7 @@
 #include "cctk_FortranString.h"
 
 #include "Cartoon2D_tensors.h"
+#include "SpaceMask.h"
 
 static const char *rcsid = "$Id$";
 
@@ -38,9 +39,11 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
                                     CCTK_REAL *f, int i, int di, 
                                     int ijk, CCTK_REAL x);
 CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
-			 CCTK_REAL y[], CCTK_REAL x);
+                            CCTK_REAL y[], 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 y[], CCTK_REAL x,
+                          CCTK_INT excised_points[]);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, 
                 int *tensortype, int *vi);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH,
@@ -49,6 +52,9 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
                           (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG);
 
 
+void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS);
+void SetupExcisionVars(CCTK_ARGUMENTS);
+
 /* set boundaries of a grid tensor assuming axisymmetry 
    - handles lower boundary in x
    - does not touch other boundaries
@@ -90,7 +96,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
   CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz;
   CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz;
 
-
   gi = CCTK_GroupIndexFromVarI (vi);
   assert (gi >= 0);
   n_vars = CCTK_NumVarsInGroupI (gi);
@@ -342,6 +347,59 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
   int lnx;
   int n, offset;
 
+  CCTK_REAL *old_mask;
+  CCTK_INT *new_mask;
+
+  CCTK_INT excised_points[6];
+
+  if (*excision_active < 0)
+  {
+    SetupExcisionVars((cGH*)cctkGH);
+  }
+
+  if (old_excision)
+  {
+    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
+  }
+  else
+  {
+    old_mask = NULL;
+  }
+  
+  if (new_excision)
+  {
+    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
+  }
+  else
+  {
+    new_mask = NULL;
+  }
+
+  /* 
+     If this boundary point is excised then we return 
+     without doing anything 
+  */
+
+  if (*excision_active)
+  {
+    if ((*excision_active == 1)||(*excision_active == 3))
+    {
+      if (old_mask[ijk] < 0.75)
+      {
+        return f[ijk];
+      }
+    }
+    if ((*excision_active == 2)||(*excision_active == 3))
+    {
+      if (SpaceMask_CheckStateBits(new_mask, (ijk),
+                                   *new_excision_field,
+                                   *new_excision_descriptor))
+      {
+        return f[ijk];
+      }
+    }
+  }  
+
   lnx = cctkGH->cctk_lsh[0];
     
   dx0 = CCTK_DELTA_SPACE(0);
@@ -369,10 +427,25 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
   x0 = lx0 + dx0 * (i-offset);
   for (n = 0; n <= order; n++) {
     y[n] = f[ijk-offset+n];
+    excised_points[n] = 0;
+    if (*excision_active)
+    {
+      if ((*excision_active == 1)||(*excision_active == 3))
+      {
+        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1;
+      }
+      if ((*excision_active == 2)||(*excision_active == 3))
+      {
+        excised_points[n] |= 
+          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
+                                    *new_excision_field,
+                                    *new_excision_descriptor));
+      }
+    }
   }
   
   /* call interpolator */
-  r = interpolate_local(order, x0, dx0, y, x);
+  r = interpolate_local(order, x0, dx0, y, x, excised_points);
   
   return r;
 }
@@ -391,6 +464,59 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   int lnx;
   int n, offset;
   int tmp_order = 2;
+
+  CCTK_REAL *old_mask;
+  CCTK_INT *new_mask;
+
+  CCTK_INT excised_points[11];
+
+  if (*excision_active < 0)
+  {
+    SetupExcisionVars((cGH*)cctkGH);
+  }
+
+  if (old_excision)
+  {
+    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
+  }
+  else
+  {
+    old_mask = NULL;
+  }
+  
+  if (new_excision)
+  {
+    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
+  }
+  else
+  {
+    new_mask = NULL;
+  }
+
+  /* 
+     If this boundary point is excised then we return 
+     without doing anything 
+  */
+
+  if (*excision_active)
+  {
+    if ((*excision_active == 1)||(*excision_active == 3))
+    {
+      if (old_mask[ijk] < 0.75)
+      {
+        return f[ijk];
+      }
+    }
+    if ((*excision_active == 2)||(*excision_active == 3))
+    {
+      if (SpaceMask_CheckStateBits(new_mask, (ijk),
+                                   *new_excision_field,
+                                   *new_excision_descriptor))
+      {
+        return f[ijk];
+      }
+    }
+  }
   
   lnx = cctkGH->cctk_lsh[0];
   
@@ -424,11 +550,92 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   for (n = 0; n < 2 * eno_order + 1; ++n) 
   {
     y[n] = f[ijk-offset+n];
+    excised_points[n] = 0;
+    if (*excision_active)
+    {
+      if ((*excision_active == 1)||(*excision_active == 3))
+      {
+        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0;
+      }
+      if ((*excision_active == 2)||(*excision_active == 3))
+      {
+        excised_points[n] |= 
+          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
+                                    *new_excision_field,
+                                    *new_excision_descriptor));
+      }
+    }
   }
   
   /* call interpolator */
 /*   r = interpolate_eno(tmp_order, x0, dx0, y, x); */
-  r = interpolate_eno(eno_order, x0, dx0, y, x);
+  r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points);
   
   return r;
 }
+
+void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  *excision_active = -1;
+  *old_mask_vi = -1;
+  *new_mask_vi = -1;
+  *new_excision_field = -1;
+  *new_excision_descriptor = -1;
+
+}
+
+void SetupExcisionVars(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (old_excision)
+  {
+    *old_mask_vi = CCTK_VarIndex(old_style_excision_var);
+  }
+  else
+  {
+    *old_mask_vi = -1;
+  }
+  if (new_excision)
+  {
+    *new_mask_vi = CCTK_VarIndex(new_style_excision_var);
+  }
+  else
+  {
+    *new_mask_vi = -1;
+  }
+
+  if (*old_mask_vi > -1)
+  {
+    if (*new_mask_vi > -1)
+    {
+      *excision_active = 3;
+      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
+      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
+                                                        new_mask_excised_name);
+    }
+    else
+    {
+      *excision_active = 1;
+    }
+  }
+  else
+  {
+    if (*new_mask_vi > -1)
+    {
+      *excision_active = 2;
+      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
+      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
+                                                        new_mask_excised_name);
+    }
+    else
+    {
+      *excision_active = 0;
+    }
+  }
+  
+}
diff --git a/src/SymInterp.c b/src/SymInterp.c
index 0efc291..24668ab 100644
--- a/src/SymInterp.c
+++ b/src/SymInterp.c
@@ -360,6 +360,13 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   for (m=0; m<N_output_arrays; ++m) {
     new_output_arrays[m] = malloc (N_interp_points * sizeof(CCTK_REAL));
     assert (new_output_arrays[m]);
+    {
+      CCTK_INT mm;
+      for (mm = 0; mm < N_interp_points; ++mm)
+      {
+        ((CCTK_REAL*)new_output_arrays[m])[mm] = 0;
+      }
+    }
   }
   
   
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]; */