1 files changed, 47 insertions, 0 deletions
diff --git a/src/SetMask.c b/src/SetMask.c
new file mode 100644
index 0000000..0e9cba9
--- /dev/null
+++ b/src/SetMask.c
@@ -0,0 +1,47 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "loopcontrol.h"
+#include "HydroBase.h"
+
+void SetMask_SphericalSurface (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT *mask = (CCTK_INT*) CCTK_VarDataPtr(cctkGH, 0, SetMask_MaskName);
+  if (!mask)
+    CCTK_WARN(0, "No such variable, or no storage enabled");
+
+  for (int smi = 0; smi < 10; smi++)
+  {
+    CCTK_INT sfi = SetMask_SurfaceIndex[smi];
+    if (sfi >= 0 && sf_active[sfi] && sf_valid[sfi] >= 0)
+    {
+      #pragma omp parallel
+      {
+        LC_LOOP3(setsurface, i,j,k, 0,0,0,
+                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2],
+                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
+        {
+          CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
+          CCTK_REAL dist2 = (sf_centroid_x[sfi]-x[i3D]) * (sf_centroid_x[sfi]-x[i3D]) +
+                            (sf_centroid_y[sfi]-y[i3D]) * (sf_centroid_y[sfi]-y[i3D]) +
+                            (sf_centroid_z[sfi]-z[i3D]) * (sf_centroid_z[sfi]-z[i3D]);
+          if (dist2 < SetMask_RadiusFactor[smi] * sf_min_radius[sfi] *
+                      SetMask_RadiusFactor[smi] * sf_min_radius[sfi])
+          {
+            mask[i3D] = HYDROEXCISION_EXCISED;
+          }
+          else if (SetMask_ResetAll)
+          {
+            mask[i3D] = HYDROEXCISION_NORMAL;
+          }
+        }
+        LC_ENDLOOP3(setsurface);
+      }
+    }
+  }
+}
+