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

4 files changed, 220 insertions, 15 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
 
 
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index e463a43..f454441 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -65,8 +65,8 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
   
-  dtheta = pi / ( ntheta - 1 )
-  dphi = two * pi / ( nphi - 1 )
+  dtheta = ctheta(2,1) - ctheta(1,1)
+  dphi = cphi(1,2) - cphi(1,1)
   dthetainv = one / dtheta
   dphiinv = one / dphi
   
@@ -207,7 +207,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   call CCTK_INFO ( 'Integrating area' )
 
   sn = surface_counter - integrate_counter
-  int_tmp = weights * da
+  int_tmp = sym_factor * weights * da
 
   call CCTK_VarIndex ( int_var, "ehfinder::int_tmp" )
   if ( int_var .lt. 0 ) then
@@ -272,18 +272,22 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
     Call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
   end if
 
-  int_tmp = interp_x * weights * da
+  int_tmp = sym_factor * interp_x * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
                      eh_centroid2_x(sn), 1, int_var )
 
-  int_tmp = interp_y * weights * da
+  int_tmp = sym_factor * interp_y * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
                      eh_centroid2_y(sn), 1, int_var )
 
-  int_tmp = interp_z * weights * da
+  int_tmp = sym_factor * interp_z * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
                      eh_centroid2_z(sn), 1, int_var )
 
+  if ( s_symx ) eh_centroid2_x(sn) = zero
+  if ( s_symy ) eh_centroid2_y(sn) = zero
+  if ( s_symz ) eh_centroid2_z(sn) = zero
+
 !  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
 !                     centroid_x, 1, int_var )
 !
diff --git a/src/EHFinder_SetSym.F90 b/src/EHFinder_SetSym.F90
index cc32b8f..817ec2e 100644
--- a/src/EHFinder_SetSym.F90
+++ b/src/EHFinder_SetSym.F90
@@ -23,6 +23,11 @@ subroutine EHFinder_SetSym(CCTK_ARGUMENTS)
     call CCTK_WARN(1,'Failed to register symmetry for the level set function')
   end if
 
+  call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::sc')
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN(1,'Failed to register symmetry for sc')
+  end if
+
 !  call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::rep_mask')
 !  if ( ierr .gt. 0 ) then
 !    call CCTK_WARN(1,'Failed to register symmetry for the level reparametrization mask')
diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90
index 0497e88..0588cbe 100644
--- a/src/EHFinder_mod.F90
+++ b/src/EHFinder_mod.F90
@@ -31,5 +31,6 @@ module EHFinder_mod
   CCTK_INT :: current_iteration
   logical :: mask_first = .true.
   logical :: symx, symy, symz
+  logical :: s_symx, s_symy, s_symz
   logical :: reparam_undone
 end module EHFinder_mod