respect emask if requested by parameter

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Dissipation/trunk@15 850bcc8b-0e4f-0410-8c26-8d28fbf1eda9

5 files changed, 72 insertions, 23 deletions

diff --git a/interface.ccl b/interface.ccl
index 1ff5f1f..69313ae 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -3,7 +3,7 @@
 
 IMPLEMENTS: Dissipation
 
-INHERITS: grid SphericalSurface
+INHERITS: grid SphericalSurface SpaceMask
 
 CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
 REQUIRES FUNCTION MoLQueryEvolvedRHS
diff --git a/param.ccl b/param.ccl
index 29b8bc0..4ab5b61 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,6 +1,12 @@
 #: Parameter definitions for thorn Dissipation
 # $Header$
 
+shares:SpaceMask
+
+USES boolean use_mask
+
+restricted:
+
 BOOLEAN verbose "produce log output" STEERABLE=always
 {
 } "no"
@@ -56,7 +62,9 @@ INT horizon_number[100] "horizon number for extra dissipation in horizons -- AHF
   1:* :: "horizon number (from AHFinderDirect); starts from 1"
 } -1
 
-
+BOOLEAN respect_emask "respect excision mask" STEERABLE=always
+{
+} "no"
 
 # Extra dissipation at outer boundary
 
diff --git a/schedule.ccl b/schedule.ccl
index 649bcf5..c39b65e 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -13,9 +13,10 @@ SCHEDULE dissipation_basegrid AT basegrid
   LANG: C
 } "Ensure that there are enough ghost zones"
 
-SCHEDULE setup_epsdis AT prestep
+SCHEDULE setup_epsdis AT cctk_poststep after SphericalSurface_HasBeenSet
 {
   LANG: C
+  SYNC: epsdisA_group
 } "Setup spatially varying dissipation"
 
 SCHEDULE dissipation_add IN MoL_PostRHS
diff --git a/src/paramcheck.c b/src/paramcheck.c
index d3629fe..86b5c97 100644
--- a/src/paramcheck.c
+++ b/src/paramcheck.c
@@ -26,4 +26,8 @@ dissipation_paramcheck (CCTK_ARGUMENTS)
       CCTK_PARAMWARN ("The aliased function \"HorizonRadiusInDirection\" must be defined when the parameter \"extra_dissipation_in_horizons\" is set and one of the sources is AHFinderDirect");
     }
   }
+  if (respect_emask && !use_mask)
+  {
+    CCTK_PARAMWARN ("You can only respect the emask if it is used (use_mask)");
+  }
 }
diff --git a/src/setup_epsdis.c b/src/setup_epsdis.c
index 6239e6c..b17b8cc 100644
--- a/src/setup_epsdis.c
+++ b/src/setup_epsdis.c
@@ -22,7 +22,7 @@ setup_epsdis (CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS;
 
   int ni,nj,nk;
-  int i,j,k,s;
+  int i,j,k,s,l,m;
   int index,indexP;
   int ierr;
   int npts;
@@ -178,14 +178,32 @@ setup_epsdis (CCTK_ARGUMENTS)
       zmax=odz+maxrad;
 
       npts=0;
-      for (i=0;i<ni*nj*nk;i++) {
-        if (  x[i]<=xmax&&x[i]>=xmin
-            &&y[i]<=ymax&&y[i]>=ymin
-            &&z[i]<=zmax&&z[i]>=zmin)
-        {
-          npts++;
+      for (i=2;i<ni-2;i++)
+       for (j=2;j<nj-2;j++)
+        for (k=2;k<nk-2;k++) {
+          m=CCTK_GFINDEX3D(cctkGH,i,j,k);
+          if ((  x[m]<=xmax&&x[m]>=xmin
+               &&y[m]<=ymax&&y[m]>=ymin
+               &&z[m]<=zmax&&z[m]>=zmin) &&
+              (!respect_emask || (
+               (emask[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i+2,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j+2,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k+2)]>0.4)
+              ))
+             )
+          {
+            npts++;
+          }
         }
-      }
 
       xa=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
       ya=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
@@ -198,19 +216,37 @@ setup_epsdis (CCTK_ARGUMENTS)
         inds[i]=0;
       }
 
-      j=0;
-      for (i=0;i<ni*nj*nk;i++) {
-        if (  x[i]<=xmax&&x[i]>=xmin
-            &&y[i]<=ymax&&y[i]>=ymin
-            &&z[i]<=zmax&&z[i]>=zmin)
-        {
-          xa[j]=x[i];
-          ya[j]=y[i];
-          za[j]=z[i];
-          inds[j]=i;
-          j++;
+      l=0;
+      for (i=2;i<ni-2;i++)
+       for (j=2;j<nj-2;j++)
+        for (k=2;k<nk-2;k++) {
+          m=CCTK_GFINDEX3D(cctkGH,i,j,k);
+          if ((  x[m]<=xmax&&x[m]>=xmin
+               &&y[m]<=ymax&&y[m]>=ymin
+               &&z[m]<=zmax&&z[m]>=zmin) &&
+              (!respect_emask || (
+               (emask[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i+2,j,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j+2,k)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]>0.4) &&
+               (emask[CCTK_GFINDEX3D(cctkGH,i,j,k+2)]>0.4)
+              ))
+             )
+          {
+            xa[l]=x[m];
+            ya[l]=y[m];
+            za[l]=z[m];
+            inds[l]=m;
+            l++;
+          }
         }
-      }
 
       ierr=HorizonRadiusInDirection(horizon_number[s],
                                     npts,