respect emask if requested by parameter

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

2 files changed, 60 insertions, 20 deletions

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,