Improving boundary treatment.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@81 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 44 insertions, 75 deletions

diff --git a/src/AHFinder_exp.F b/src/AHFinder_exp.F
index 12ec8f8..880fa94 100644
--- a/src/AHFinder_exp.F
+++ b/src/AHFinder_exp.F
@@ -323,87 +323,56 @@
          end do
       end do
 
-!     Boundaries.
-
-      do k=1,nz
-         do j=1,ny
-
-            ahf_exp(1 ,j,k) = ahf_exp(2   ,j,k)
-            ahf_exp(nx,j,k) = ahf_exp(nx-1,j,k)
-
-            ahfgradx(nx,j,k) = ahfgradx(nx-1,j,k)
-            ahfgrady(nx,j,k) = ahfgrady(nx-1,j,k)
-            ahfgradz(nx,j,k) = ahfgradz(nx-1,j,k)
-            ahfgradn(nx,j,k) = ahfgradn(nx-1,j,k)
-
-            ahfgrady(1 ,j,k) = ahfgrady(2,j,k)
-            ahfgradz(1 ,j,k) = ahfgradz(2,j,k)
-            ahfgradn(1 ,j,k) = ahfgradn(2,j,k)
-
-            if ((CCTK_Equals(domain,"octant").eq.1).or.
-     .           (CCTK_Equals(domain,"quadrant").eq.1)) then
-               ahfgradx(1,j,k) = -ahfgradx(2,j,k)
-            else
-               ahfgradx(1,j,k) = +ahfgradx(2,j,k)
-            end if
-
-         end do
-      end do
-
-      do k=1,nz
-         do i=1,nx
-
-            ahf_exp(i,1 ,k) = ahf_exp(i,2   ,k)
-            ahf_exp(i,ny,k) = ahf_exp(i,ny-1,k)
-
-            ahfgradx(i,ny,k) = ahfgradx(i,ny-1,k)
-            ahfgrady(i,ny,k) = ahfgrady(i,ny-1,k)
-            ahfgradz(i,ny,k) = ahfgradz(i,ny-1,k)
-            ahfgradn(i,ny,k) = ahfgradn(i,ny-1,k)
-
-            ahfgradx(i,1 ,k) = ahfgradx(i,2,k)
-            ahfgradz(i,1 ,k) = ahfgradz(i,2,k)
-            ahfgradn(i,1 ,k) = ahfgradn(i,2,k)
-
-            if ((CCTK_Equals(domain,"octant").eq.1).or.
-     .           (CCTK_Equals(domain,"quadrant").eq.1)) then
-               ahfgrady(i,1,k) = -ahfgrady(i,2,k)
-            else
-               ahfgrady(i,1,k) = +ahfgrady(i,2,k)
-            end if
-
-         end do
-      end do
-
-      do j=1,ny
-         do i=1,nx
-
-            ahf_exp(i,j,1 ) = ahf_exp(i,j,2   )
-            ahf_exp(i,j,nz) = ahf_exp(i,j,nz-1)
-
-            ahfgradx(i,j,nz) = ahfgradx(i,j,nz-1)
-            ahfgrady(i,j,nz) = ahfgrady(i,j,nz-1)
-            ahfgradz(i,j,nz) = ahfgradz(i,j,nz-1)
-            ahfgradn(i,j,nz) = ahfgradn(i,j,nz-1)
-
-            ahfgradx(i,j,1 ) = ahfgradx(i,j,2)
-            ahfgrady(i,j,1 ) = ahfgrady(i,j,2)
-            ahfgradn(i,j,1 ) = ahfgradn(i,j,2)
-
-            if (CCTK_Equals(domain,"octant").eq.1) then
-               ahfgradz(i,j,1) = -ahfgradz(i,j,2)
-            else
-               ahfgradz(i,j,1) = +ahfgradz(i,j,2)
-            end if
-
-         end do
-      end do
+!     Boundaries on x direction.
+
+      ahf_exp(1,:,:)  = 2.0D0*ahf_exp(2,:,:) - ahf_exp(3,:,:)
+      ahfgradx(1,:,:) = 2.0D0*ahfgradx(2,:,:) - ahfgradx(3,:,:)
+      ahfgrady(1,:,:) = 2.0D0*ahfgrady(2,:,:) - ahfgrady(3,:,:)
+      ahfgradz(1,:,:) = 2.0D0*ahfgradz(2,:,:) - ahfgradz(3,:,:)
+      ahfgradn(1,:,:) = 2.0D0*ahfgradn(2,:,:) - ahfgradn(3,:,:)
+
+      ahf_exp(nx,:,:)  = 2.0D0*ahf_exp(nx-1,:,:) - ahf_exp(nx-2,:,:)
+      ahfgradx(nx,:,:) = 2.0D0*ahfgradx(nx-1,:,:) - ahfgradx(nx-2,:,:)
+      ahfgrady(nx,:,:) = 2.0D0*ahfgrady(nx-1,:,:) - ahfgrady(nx-2,:,:)
+      ahfgradz(nx,:,:) = 2.0D0*ahfgradz(nx-1,:,:) - ahfgradz(nx-2,:,:)
+      ahfgradn(nx,:,:) = 2.0D0*ahfgradn(nx-1,:,:) - ahfgradn(nx-2,:,:)
+
+!     Boundaries on y direction.
+
+      ahf_exp(:,1,:)  = 2.0D0*ahf_exp(:,2,:) - ahf_exp(:,3,:)
+      ahfgradx(:,1,:) = 2.0D0*ahfgradx(:,2,:) - ahfgradx(:,3,:)
+      ahfgrady(:,1,:) = 2.0D0*ahfgrady(:,2,:) - ahfgrady(:,3,:)
+      ahfgradz(:,1,:) = 2.0D0*ahfgradz(:,2,:) - ahfgradz(:,3,:)
+      ahfgradn(:,1,:) = 2.0D0*ahfgradn(:,2,:) - ahfgradn(:,3,:)
+
+      ahf_exp(:,ny,:)  = 2.0D0*ahf_exp(:,ny-1,:) - ahf_exp(:,ny-2,:)
+      ahfgradx(:,ny,:) = 2.0D0*ahfgradx(:,ny-1,:) - ahfgradx(:,ny-2,:)
+      ahfgrady(:,ny,:) = 2.0D0*ahfgrady(:,ny-1,:) - ahfgrady(:,ny-2,:)
+      ahfgradz(:,ny,:) = 2.0D0*ahfgradz(:,ny-1,:) - ahfgradz(:,ny-2,:)
+      ahfgradn(:,ny,:) = 2.0D0*ahfgradn(:,ny-1,:) - ahfgradn(:,ny-2,:)
+
+!     Boundaries on z direction.
+
+      ahf_exp(:,:,1)  = 2.0D0*ahf_exp(:,:,2) - ahf_exp(:,:,3)
+      ahfgradx(:,:,1) = 2.0D0*ahfgradx(:,:,2) - ahfgradx(:,:,3)
+      ahfgrady(:,:,1) = 2.0D0*ahfgrady(:,:,2) - ahfgrady(:,:,3)
+      ahfgradz(:,:,1) = 2.0D0*ahfgradz(:,:,2) - ahfgradz(:,:,3)
+      ahfgradn(:,:,1) = 2.0D0*ahfgradn(:,:,2) - ahfgradn(:,:,3)
+
+      ahf_exp(:,:,nz)  = 2.0D0*ahf_exp(:,:,nz-1) - ahf_exp(:,:,nz-2)
+      ahfgradx(:,:,nz) = 2.0D0*ahfgradx(:,:,nz-1) - ahfgradx(:,:,nz-2)
+      ahfgrady(:,:,nz) = 2.0D0*ahfgrady(:,:,nz-1) - ahfgrady(:,:,nz-2)
+      ahfgradz(:,:,nz) = 2.0D0*ahfgradz(:,:,nz-1) - ahfgradz(:,:,nz-2)
+      ahfgradn(:,:,nz) = 2.0D0*ahfgradn(:,:,nz-1) - ahfgradn(:,:,nz-2)
 
 !     Synchronize.
 
       call CCTK_SyncGroup(cctkGH,"ahfinder::ahfgradient")
       call CCTK_SyncGroup(cctkGH,"ahfinder::ahfinderexp")
 
+      call CartSymBCGroup(ierr,cctkGH,"ahfinder::ahfgradient")
+      call CartSymBCGroup(ierr,cctkGH,"ahfinder::ahfinderexp")
+
 
 !     ***************
 !     ***   END   ***