From e0b8a997bc1cbb959cb739c37ad14344f7f49873 Mon Sep 17 00:00:00 2001
From: diener <diener@4ec1db94-0e4f-0410-ada3-8bed251432c9>
Date: Tue, 9 Oct 2007 11:00:41 +0000
Subject: Add the option to use user defined regions for the elliptic smoothing
 procedure. Works by setting the admbase variables to zero in these regions
 before doing the smoothing. Needs some testing.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@12 4ec1db94-0e4f-0410-ada3-8bed251432c9
---
 param.ccl    |  4 ++++
 schedule.ccl |  7 ++++++
 src/cg.F90   | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 86 insertions(+)

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
-- 
cgit v1.2.3