Changed to use the newest version of MoL. Used its support for evolving

grid arrays to introduce tracking of the generators of the horizon along
with the horizon itself. Still needs some work on setting up the initial
generators (now it only does one) but first it needs to be tested.


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

6 files changed, 205 insertions, 20 deletions

diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
new file mode 100644
index 0000000..924588f
--- /dev/null
+++ b/src/EHFinder_Generator_Sources.F90
@@ -0,0 +1,157 @@
+! Calculation of the sources for the level set function.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i
+  CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+  CCTK_INT, dimension(1) :: lsh
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(14) :: out_arrays
+  CCTK_INT, dimension(12) :: in_arrays
+  CCTK_INT, dimension(14), parameter :: op_indices = (/ 0, 1, 2, 3, 4, &
+                                                        5, 6, 7, 8, 9, &
+                                                        10, 10, 10, 11 /), &
+                                        op_codes = (/ 0, 0, 0, 0, 0, &
+                                                      0, 0, 0, 0, 0, &
+                                                      1, 2, 3, 0 /)
+  CCTK_INT, dimension(14) :: out_types
+  CCTK_REAL :: alp2, dfux, dfuy, dfuz, factor
+  CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
+
+  out_types = CCTK_VARIABLE_REAL
+
+  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    end if
+!    print*,lsh 
+
+    call CCTK_InterpHandle ( interp_handle, &
+                                     "Lagrange polynomial interpolation" )
+    if ( interp_handle .lt. 0 ) then
+      call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+    end if
+    ! For now order=2 is hard wired.
+    call Util_TableCreateFromString ( table_handle, "order=3" )
+    if ( table_handle .lt. 0 ) then
+      call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+    end if
+
+    ! Get the 3D coordinate system handle.
+    call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+    if ( coord_system_handle .lt. 0) then
+      call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
+    endif
+
+    interp_coords(1) = CCTK_PointerTo(xg)
+    interp_coords(2) = CCTK_PointerTo(yg)
+    interp_coords(3) = CCTK_PointerTo(zg)
+
+    out_arrays(1) = CCTK_PointerTo(alpg)
+    out_arrays(2) = CCTK_PointerTo(betaxg)
+    out_arrays(3) = CCTK_PointerTo(betayg)
+    out_arrays(4) = CCTK_PointerTo(betazg)
+    out_arrays(5) = CCTK_PointerTo(gxxg)
+    out_arrays(6) = CCTK_PointerTo(gxyg)
+    out_arrays(7) = CCTK_PointerTo(gxzg)
+    out_arrays(8) = CCTK_PointerTo(gyyg)
+    out_arrays(9) = CCTK_PointerTo(gyzg)
+    out_arrays(10) = CCTK_PointerTo(gzzg)
+    out_arrays(11) = CCTK_PointerTo(dfxg)
+    out_arrays(12) = CCTK_PointerTo(dfyg)
+    out_arrays(13) = CCTK_PointerTo(dfzg)
+
+    call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
+    call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
+    call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
+    call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
+    call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
+    call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
+    call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
+    call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
+    call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
+    call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
+    call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" )
+
+    call Util_TableSetIntArray ( status, table_handle, 13, &
+                               op_indices(1:13), "operand_indices" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+    endif
+
+    call Util_TableSetIntArray ( status, table_handle, 13, &
+                               op_codes(1:13), "operation_codes" )
+    if ( status .lt. 0 ) then
+      call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+    endif
+
+
+    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 11, in_arrays(1:11), &
+                                 13, out_types(1:13), out_arrays(1:13) )
+
+    if ( status .lt. 0 ) then
+      call CCTK_INFO ( 'Interpolation failed.' )
+    end if
+    do i = 1, lsh(1)
+      alp2 = alpg(i)**2
+
+      guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
+      guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
+      guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
+
+      idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz )
+
+      guxx = idetg * guxx
+      guxy = idetg * guxy
+      guxz = idetg * guxz
+
+      guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
+      guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
+      guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
+      
+      dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
+      dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
+      dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
+      factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
+                                      dfuy * dfyg(i) + &
+                                      dfuz * dfzg(i) ) )
+      dxg(i) = - betaxg(i) + factor * dfux
+      dyg(i) = - betayg(i) + factor * dfuy
+      dzg(i) = - betazg(i) + factor * dfuz
+!      print*, alpg
+!      print*
+!      print*, betaxg, betayg, betazg
+!      print*
+!      print*, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg
+!      print*
+!      print*, dfxg, dfyg, dfzg
+!      print*
+!      print*, alp2
+!      print*
+!      print*, dfux, dfuy, dfuz
+!      print*
+!      print*, factor
+!      print*
+!      print*, dxg(i), dyg(i), dzg(i)
+!      stop
+    end do
+  end if
+  return
+end subroutine EHFinder_Generator_Sources
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index cf700f3..ebb7cb7 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -41,6 +41,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
     f = sqrt( ( x - translate_x )**2 + &
               ( y - translate_y )**2 + &
               ( z - translate_z )**2 ) - initial_rad
+    if ( evolve_generators .gt. 0 ) then
+      xg(1) = initial_rad + translate_x
+      yg(1) = translate_y
+      zg(1) = translate_z
+    end if
   end if
 
   ! If an ellipsoid is requested...
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index 81745ea..48bd609 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -35,7 +35,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
                                                       CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL /)
   
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -212,7 +212,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   character(len=30) :: info_message
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -274,7 +274,7 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
   CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
@@ -382,7 +382,7 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS)
   CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha
   CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ
 
-! If finding of surface failed don't try to integrate but exit.
+! If finding of surface failed do not try to integrate but exit.
   if ( find_surface_status .lt. 0 ) then
     return
   endif
diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90
index 679c8fb..bccc870 100644
--- a/src/MoL_Init.F90
+++ b/src/MoL_Init.F90
@@ -10,6 +10,7 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
   implicit none
 
   CCTK_INT :: ierr, ierr_cum, varindex, rhsindex
+  CCTK_INT :: MoLRegisterEvolved
 
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
@@ -18,8 +19,28 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
   ierr_cum = 0
   call CCTK_VarIndex(varindex, 'ehfinder::f')
   call CCTK_VarIndex(rhsindex, 'ehfinder::sf')
-  call MoL_RegisterVar(ierr, varindex, rhsindex)
-  ierr_cum = ierr_cum + ierr
+!  call MoL_RegisterVar(ierr, varindex, rhsindex)
+!  ierr_cum = ierr_cum + ierr
+  ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+  if ( evolve_generators .gt. 0 ) then
+
+!    call CCTK_GroupIndex (varindex, 'ehfinder::generators')
+!    call CCTK_GroupIndex(rhsindex, 'ehfinder::generators_rhs')
+
+    call CCTK_VarIndex(varindex, 'ehfinder::xg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dxg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+    call CCTK_VarIndex(varindex, 'ehfinder::yg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dyg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+
+    call CCTK_VarIndex(varindex, 'ehfinder::zg')
+    call CCTK_VarIndex(rhsindex, 'ehfinder::dzg')
+    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
+  end if
+ 
 !  call CCTK_VarIndex(varindex, 'admbase::gxx')
 !  call MoL_RegisterPrimitive(ierr, varindex)
 !  ierr_cum = ierr_cum + ierr
diff --git a/src/make.code.defn b/src/make.code.defn
index c54b696..c8e8ce4 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -7,6 +7,7 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_Init.F90 \
        MoL_Init.F90 \
        EHFinder_Sources.F90 \
+       EHFinder_Generator_Sources.F90 \
        EHFinder_ReParametrize.F90 \
        EHFinder_Check.F90 \
        EHFinder_Sources2.F90 \
diff --git a/src/make.code.deps b/src/make.code.deps
index 1c09018..d8e8918 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -1,17 +1,18 @@
 # Dependencies file for thorn EHFinder.
 # $Header$
 
-EHFinder_mod.F90.o:		EHFinder_Constants.F90.o
-EHFinder_Init.F90.o:		EHFinder_mod.F90.o
-EHFinder_Sources.F90.o:		EHFinder_mod.F90.o
-EHFinder_ReParametrize.F90.o:	EHFinder_mod.F90.o
-EHFinder_Check.F90.o:		EHFinder_mod.F90.o
-EHFinder_Sources2.F90.o:	EHFinder_mod.F90.o
-EHFinder_Sources3.F90.o:	EHFinder_mod.F90.o
-EHFinder_SetMask.F90.o:		EHFinder_mod.F90.o
-EHFinder_SetSym.F90.o:		EHFinder_mod.F90.o
-EHFinder_ReadData.F90.o:	EHFinder_mod.F90.o
-EHFinder_Check.F90.o:		EHFinder_mod.F90.o
-EHFinder_Integrate.F90.o:	EHFinder_mod.F90.o
-EHFinder_Integrate2.F90.o:	EHFinder_mod.F90.o
-EHFinder_FindSurface.F90.o:	EHFinder_mod.F90.o
+EHFinder_mod.F90.o:			EHFinder_Constants.F90.o
+EHFinder_Init.F90.o:			EHFinder_mod.F90.o
+EHFinder_Sources.F90.o:			EHFinder_mod.F90.o
+EHFinder_Generator_Sources.F90.o:	EHFinder_mod.F90.o
+EHFinder_ReParametrize.F90.o:		EHFinder_mod.F90.o
+EHFinder_Check.F90.o:			EHFinder_mod.F90.o
+EHFinder_Sources2.F90.o:		EHFinder_mod.F90.o
+EHFinder_Sources3.F90.o:		EHFinder_mod.F90.o
+EHFinder_SetMask.F90.o:			EHFinder_mod.F90.o
+EHFinder_SetSym.F90.o:			EHFinder_mod.F90.o
+EHFinder_ReadData.F90.o:		EHFinder_mod.F90.o
+EHFinder_Check.F90.o:			EHFinder_mod.F90.o
+EHFinder_Integrate.F90.o:		EHFinder_mod.F90.o
+EHFinder_Integrate2.F90.o:		EHFinder_mod.F90.o
+EHFinder_FindSurface.F90.o:		EHFinder_mod.F90.o