3 files changed, 86 insertions, 0 deletions
diff --git a/param.ccl b/param.ccl
index 6d93bd3..0b8c6b8 100644
--- a/param.ccl
+++ b/param.ccl
@@ -142,3 +142,7 @@ REAL smoothing_eps "CG smoothing stop criteria"
 {
   (0.0:* :: ""
 } 1e-6
+
+BOOLEAN use_user_regions "Use user defined regions for the smoothing regions"
+{
+} "no"
diff --git a/schedule.ccl b/schedule.ccl
index 54d402b..824d112 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -20,6 +20,13 @@ if (smooth_regions) {
       LANG: Fortran
     } "Smooth regions"
   } else {
+    
+    if (use_user_regions) {
+      SCHEDULE NoExcision_Set_Zero IN ADMBase_PostInitial BEFORE NoExcision_CGSmoothing
+      {
+        LANG: Fortran
+      } "set variablse to zero in used defined regions"
+    }
 
     SCHEDULE GROUP NoExcision_CGSmoothing IN ADMBase_PostInitial
     {
diff --git a/src/cg.F90 b/src/cg.F90
index ec4866e..c3dd6c8 100644
--- a/src/cg.F90
+++ b/src/cg.F90
@@ -645,3 +645,78 @@ subroutine NoExcision_CGApplySym(CCTK_ARGUMENTS)
   end if
 
 end subroutine NoExcision_CGApplySym
+
+
+subroutine NoExcision_Set_Zero(CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz
+  CCTK_REAL, dimension(:,:,:), allocatable :: dist2
+  integer :: n
+  CCTK_INT :: my_level, n_levels
+
+  call NoExcision_levelinfo ( cctkGH, my_level, n_levels )
+
+  if ( my_level == n_levels-1 ) then
+
+    allocate ( dist2(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) )
+
+    do n = 1, num_regions
+
+      cx = centre_x(n)
+      cy = centre_y(n)
+      cz = centre_z(n)
+
+      if (CCTK_EQUALS(region_shape(n), "sphere")) then
+
+        radx = radius(n)
+        rady = radius(n)
+        radz = radius(n)
+
+      else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+
+        radx = radius_x(n)
+        rady = radius_y(n)
+        radz = radius_z(n)
+
+      else
+
+        call CCTK_WARN (0, "internal error")
+
+      end if 
+
+      dist2 =   ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+           &  + ((z - cz) / radz)**2
+
+      where ( dist2 <= one )
+
+        gxx = zero
+        gxy = zero
+        gxz = zero
+        gyy = zero
+        gyz = zero
+        gzz = zero
+        kxx = zero
+        kxy = zero
+        kxz = zero
+        kyy = zero
+        kyz = zero
+        kzz = zero
+        alp = zero
+        betax = zero
+        betay = zero
+        betaz = zero
+
+      end where
+
+    end do
+
+    deallocate ( dist2 )
+
+  end if
+end subroutine NoExcision_Set_Zero