1 files changed, 150 insertions, 0 deletions
diff --git a/src/smooth.F90 b/src/smooth.F90
new file mode 100644
index 0000000..d6df156
--- /dev/null
+++ b/src/smooth.F90
@@ -0,0 +1,150 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Smooth (CCTK_ARGUMENTS)
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero=0
+  CCTK_REAL :: mask (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz
+  integer   :: n
+  integer   :: iter
+  integer   :: ierr
+  
+  do iter = 1, smoothing_iterations
+     
+     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
+        
+        mask = 1 - sqrt (((x - cx) / radx)**2 + &
+             &           ((y - cy) / rady)**2 + &
+             &           ((z - cz) / radz)**2)
+        where (mask <= 0)
+           mask = 0             ! outside
+        else where (mask >= smoothing_zone_width(n))
+           mask = 1             ! far inside
+        elsewhere
+           mask = mask / smoothing_zone_width(n) ! a bit inside
+        end where
+        
+        if (overwrite_geometry(n) /= 0) then
+           
+           if (conformal_state >= 1) then
+              call smooth (psi)
+           end if
+           if (conformal_state >= 2) then
+              call smooth (psix)
+              call smooth (psiy)
+              call smooth (psiz)
+           end if
+           if (conformal_state >= 3) then
+              call smooth (psixx)
+              call smooth (psixy)
+              call smooth (psixz)
+              call smooth (psiyy)
+              call smooth (psiyz)
+              call smooth (psizz)
+           end if
+           
+           call smooth (gxx)
+           call smooth (gxy)
+           call smooth (gxz)
+           call smooth (gyy)
+           call smooth (gyz)
+           call smooth (gzz)
+           call smooth (kxx)
+           call smooth (kxy)
+           call smooth (kxz)
+           call smooth (kyy)
+           call smooth (kyz)
+           call smooth (kzz)
+           
+        end if
+        
+        if (overwrite_lapse(n) /= 0) then
+           
+           call smooth (alp)
+           
+        end if
+        
+        if (overwrite_shift(n) /= 0) then
+           
+           if (shift_state /= 0) then
+              call smooth (betax)
+              call smooth (betay)
+              call smooth (betaz)
+           end if
+           
+        end if
+        
+     end do
+     
+     if (overwrite_geometry(n) /= 0) then
+        
+        if (conformal_state >= 1) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac")
+        end if
+        if (conformal_state >= 2) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_1derivs")
+        end if
+        if (conformal_state >= 3) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_2derivs")
+        end if
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::metric")
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::curv")
+        
+     end if
+     
+     if (overwrite_lapse(n) /= 0) then
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::lapse")
+        
+     end if
+     
+     if (overwrite_shift(n) /= 0) then
+        
+        if (shift_state /= 0) then
+           call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::shift")
+        end if
+        
+     end if
+     
+  end do
+  
+contains
+  
+  subroutine smooth (val)
+    CCTK_REAL, intent(inout) :: val(:,:,:)
+    
+    where (mask > 0)
+       val =  (1 - mask * smoothing_factor) * val &
+            & + (mask * smoothing_factor / 6) &
+            &   * (+ eoshift(val, +1, dim=1) + eoshift(val, -1, dim=1) &
+            &      + eoshift(val, +1, dim=2) + eoshift(val, -1, dim=2) &
+            &      + eoshift(val, +1, dim=3) + eoshift(val, -1, dim=3))
+    end where
+  end subroutine smooth
+  
+end subroutine NoExcision_Smooth