Added support for figuring out the symmetry of the surface based on the

symmetries of the underlying grid and the location of the surface on the
grid. Has been tested for all octant, quadrant and bitant cases I could
think of in a simple case, but hasn't been tested in a production run.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@91 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 204 insertions, 9 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 41bbac4..d9be219 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -27,6 +27,9 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   CCTK_REAL, dimension(3) :: center_loc
   CCTK_INT :: nc_loc, nc
   CCTK_REAL :: initial_radius, ncinv, maxdr_loc, maxdr, maxf_loc, maxf
+  CCTK_REAL, dimension(3) :: crange_min_loc, crange_max_loc, &
+                             crange_min_glob, crange_max_glob
+  CCTK_REAL :: ltheta, utheta, lphi, uphi
 
   CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv, delta, rloc
   CCTK_POINTER, dimension(3) :: interp_coords
@@ -127,29 +130,219 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     ncinv = one / nc
     center = ncinv * center
 !    print*,'EHFinder_FindSurfaces : ', sn, center
-
-    ! For now symmetries are not handled properly
-    if ( symx ) center(1) = zero
-    if ( symy ) center(2) = zero
-    if ( symz ) center(3) = zero
+!    print*
   end if
 !  print*,'Center : ',center
 !  print*,'use_user_center : ',use_user_center
 
+  ! Try to figure out the symmetries. First find the range in coordinates
+  ! for the points inside the current surface.
+  crange_min_loc(1) = minval ( x, mask = sc .eq. sn )
+  crange_min_loc(2) = minval ( y, mask = sc .eq. sn )
+  crange_min_loc(3) = minval ( z, mask = sc .eq. sn )
+  crange_max_loc(1) = maxval ( x, mask = sc .eq. sn )
+  crange_max_loc(2) = maxval ( y, mask = sc .eq. sn )
+  crange_max_loc(3) = maxval ( z, mask = sc .eq. sn )
+
+  call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
+                                      crange_min_loc, crange_min_glob, &
+                                      3, CCTK_VARIABLE_REAL )
+  call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
+                                      crange_max_loc, crange_max_glob, &
+                                      3, CCTK_VARIABLE_REAL )
+!  print*,crange_min_loc,crange_max_loc
+!  print*,crange_min_glob,crange_max_glob
+!  print*
+!  print*,'Symmetries : ', symx, symy, symz
+
+  ! Full mode. Modify if necessary.
+  ltheta = zero; utheta = pi
+  lphi = zero; uphi = two * pi
+  sym_factor = one
+  s_symx = .false.; s_symy = .false.; s_symz = .false.
+
+  ! Bitant with symmetry in the z-direction
+  if ( ( .not. symx ) .and. ( .not. symy ) .and. ( symz ) ) then
+    if ( crange_min_glob(3) .lt. zero ) then
+      ltheta = zero; utheta = half * pi
+      lphi = zero; uphi = two * pi
+      sym_factor = two
+      center(3) = zero
+      s_symz = .true.
+    end if 
+  end if
+
+  ! Bitant with symmetry in the y-direction
+  if ( ( .not. symx ) .and. ( symy ) .and. ( .not. symz ) ) then
+    if ( crange_min_glob(2) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = zero; uphi = pi
+      sym_factor = two
+      center(2) = zero
+      s_symy = .true.
+    end if 
+  end if
+
+  ! Bitant with symmetry in the x-direction
+  if ( ( symx ) .and. ( .not. symy ) .and. ( .not. symz ) ) then
+    if ( crange_min_glob(1) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = two
+      center(1) = zero
+      s_symx = .true.
+    end if 
+  end if
+
+  ! Quadrant in x-direction
+  if ( ( .not. symx ) .and. ( symy ) .and. ( symz ) ) then
+    if ( ( crange_min_glob(3) .lt. zero ) .and. &
+         ( crange_min_glob(2) .lt. zero ) ) then
+      ltheta = zero; utheta = half * pi
+      lphi = zero; uphi = pi
+      sym_factor = four
+      center(3) = zero; center(2) = zero
+      s_symz = .true.; s_symy = .true.
+    else if ( crange_min_glob(3) .lt. zero ) then
+      ltheta = zero; utheta = half * pi
+      lphi = 0; uphi = two * pi
+      sym_factor = two
+      center(3) = zero
+      s_symz = .true.
+    else if ( crange_min_glob(2) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = 0; uphi = pi
+      sym_factor = two
+      center(2) = zero
+      s_symy = .true.
+    end if
+  end if
+    
+  ! Quadrant in y-direction
+  if ( ( symx ) .and. ( .not. symy ) .and. ( symz ) ) then
+    if ( ( crange_min_glob(3) .lt. zero ) .and. &
+         ( crange_min_glob(1) .lt. zero ) ) then
+      ltheta = zero; utheta = half * pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = four
+      center(3) = zero; center(1) = zero
+      s_symz = .true.; s_symx = .true.
+    else if ( crange_min_glob(3) .lt. zero ) then
+      ltheta = zero; utheta = half * pi
+      lphi = 0; uphi = two * pi
+      sym_factor = two
+      center(3) = zero
+      s_symz = .true.
+    else if ( crange_min_glob(1) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = two
+      center(1) = zero
+      s_symx = .true.
+    end if
+  end if
+
+  ! Quadrant in z-direction
+  if ( ( symx ) .and. ( symy ) .and. ( .not. symz ) ) then
+    if ( ( crange_min_glob(2) .lt. zero ) .and. &
+         ( crange_min_glob(1) .lt. zero ) ) then
+      ltheta = zero; utheta = pi
+      lphi = zero; uphi = half * pi
+      sym_factor = four
+      center(2) = zero; center(1) = zero
+      s_symy = .true.; s_symx = .true.
+    else if ( crange_min_glob(2) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = 0; uphi = pi
+      sym_factor = two
+      center(2) = zero
+      s_symy = .true.
+    else if ( crange_min_glob(1) .lt. zero ) then
+      ltheta = zero; utheta = pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = two
+      center(1) = zero
+      s_symx = .true.
+    end if
+  end if
+    
+  ! Octant
+  if ( ( symx ) .and. ( symy ) .and. ( symz ) ) then
+    if ( ( crange_min_glob(3) .lt. zero ) .and. &
+         ( crange_min_glob(2) .lt. zero ) .and. &
+         ( crange_min_glob(1) .lt. zero ) ) then
+      ltheta = zero; utheta = half * pi
+      lphi = zero; uphi = half * pi
+      sym_factor = eight
+      center(3) = zero; center(2) = zero; center(1) = zero
+      s_symz = .true.; s_symy = .true.; s_symx = .true.
+    else if ( ( crange_min_glob(3) .lt. zero ) .and. &
+              ( crange_min_glob(2) .lt. zero ) ) then
+      ltheta = zero; utheta = half * pi
+      lphi = zero; uphi = pi
+      sym_factor = four
+      center(3) = zero; center(2) = zero
+      s_symz = .true.; s_symy = .true.
+    else if ( ( crange_min_glob(3) .lt. zero ) .and. &
+              ( crange_min_glob(1) .lt. zero ) ) then
+      ltheta = zero; utheta = half * pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = four
+      center(3) = zero; center(1) = zero
+      s_symz = .true.; s_symx = .true.
+    else if ( ( crange_min_glob(2) .lt. zero ) .and. &
+              ( crange_min_glob(1) .lt. zero ) ) then
+      ltheta = zero; utheta = pi
+      lphi = zero; uphi = half * pi
+      sym_factor = four
+      center(2) = zero; center(1) = zero
+      s_symy = .true.; s_symx = .true.
+    else if ( crange_min_glob(3) .lt. zero ) then 
+      ltheta = zero; utheta = half * pi
+      lphi = zero; uphi = two * pi
+      sym_factor = two
+      center(3) = zero
+      s_symz = .true.
+    else if ( crange_min_glob(2) .lt. zero ) then 
+      ltheta = zero; utheta = pi
+      lphi = zero; uphi = pi
+      sym_factor = two
+      center(2) = zero
+      s_symy = .true.
+    else if ( crange_min_glob(1) .lt. zero ) then 
+      ltheta = zero; utheta = pi
+      lphi = -half * pi; uphi = half * pi
+      sym_factor = two
+      center(1) = zero
+      s_symx = .true.
+    end if
+  end if
+
+!  print*
+!  print*,ltheta/pi,utheta/pi,lphi/pi,uphi/pi
+!  print*,sym_factor
+!  print*,center
+!  ! For now symmetries are not handled properly
+!  if ( symx ) center(1) = zero
+!  if ( symy ) center(2) = zero
+!  if ( symz ) center(3) = zero
+
   ! Find dtheta and dphi and initialise the theta and phi grid arrays.
   ! Here i + lbnd(1) - 1 is the global index for theta and
   !      j + lbnd(2) - 1 is the global index for phi.
-  dtheta = pi / ( ntheta - 1 )
-  dphi = two * pi / ( nphi - 1 )
+
+  dtheta = ( utheta - ltheta ) / ( ntheta - 1 )
+  dphi = ( uphi - lphi ) / ( nphi - 1 )
   dthetainv = one / dtheta
   dphiinv = one / dphi
   do i = 1, lsh(1)
-    ctheta(i,:) = dtheta * ( i + lbnd(1) - 1 )
+    ctheta(i,:) = ltheta + dtheta * ( i + lbnd(1) - 1 )
   end do
   do j = 1, lsh(2)
-    cphi(:,j) = dphi * ( j +lbnd(2) - 1 )
+    cphi(:,j) = lphi + dphi * ( j +lbnd(2) - 1 )
   end do
 
+!  print*,'Delta : ',dtheta,dphi
   ! Calculate sines and cosines and store them in grid arrays since they
   ! are expensive.
   sintheta = sin(ctheta)
@@ -621,6 +814,8 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
   ! but for some reason that does not work.
   call CCTK_SyncGroup ( ierr, cctkGH, "ehfinder::surface_index" )
 
+  call CartSymVN(ierr,cctkGH,'ehfinder::sc')
+
 end subroutine EHFinder_MarkSurfaces