EHFinder_IsoSurface_optimized.F90 compiles right now and it should compute the area of the event horizon- given the appropriate level set function- still in testing...

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

1 files changed, 180 insertions, 43 deletions

diff --git a/src/EHFinder_IsoSurface_optimized.F90 b/src/EHFinder_IsoSurface_optimized.F90
index 0b8eaf6..efb32fa 100755
--- a/src/EHFinder_IsoSurface_optimized.F90
+++ b/src/EHFinder_IsoSurface_optimized.F90
@@ -21,7 +21,7 @@ subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
   CCTK_REAL, dimension(:), allocatable :: x_coord, y_coord, z_coord
   CCTK_INT, dimension(:), allocatable :: t_index
   CCTK_REAL, parameter :: iso_value = zero
-  CCTK_INT :: l
+  CCTK_INT :: l,contor
 
 ! Find index ranges for points excluding ghost zones and symmetry points for
 ! the 3D grid functions.
@@ -34,18 +34,155 @@ subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
   print*,Xf0,Xf1,Xf2,f_length
   print*,size(x)
 !  allocate(x_coord(5),y_coord(5),z_coord(5),t_index(20))
+
   do l = 1, eh_number_level_sets
 !    call triangularization ( n_coords, n_triangles, nx, ny, nz, x_coord, &
 !                             y_coord, z_coord, t_index, f(:,:,:,l), x, y, z )
-    call triangularization ( n_coords, n_triangles, nx, ny, nz, f(:,:,:,l), x, y, z )
-
+!    call triangularization ( n_coords, n_triangles, nx, ny, nz, t_index, f(:,:,:,l), x, y, z )
+     call triangularization(l,contor,nx,ny,nz,t_index,f(:,:,:,l),x,y,z)
   print*,'The surface consists of ',n_triangles/3,' triangles and ',n_coords,' coords'
   end do
-
 end subroutine EHFinder_IsoSurface
 
-  !subroutine triangularization(l,contor,Nx,Ny,Nz,coord_x,coord_z,coord_y,triangle_index,grid_value,grid_x,grid_y,grid_z)
-subroutine triangularization(l,contor,Nx,Ny,Nz,grid_value,grid_x,grid_y,grid_z)
+subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area)
+  
+  use EHFinder_mod
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+!  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS 
+
+  CCTK_POINTER, intent(IN) :: cctkGH
+  CCTK_INT, intent(IN) :: npoints, l
+  CCTK_REAL, intent(IN), dimension(npoints) :: x_coord,y_coord,z_coord
+  CCTK_INT, intent(IN), dimension(l) :: t_index
+  CCTK_REAL, intent(OUT) :: area
+  character(len=200) :: area_interp
+  character(len=128) :: warn_message
+  CCTK_INT :: area_interp_len, local_contor
+  character(len=7) :: area_order
+  CCTK_INT :: interp_handle, table_handle, coord_system_handle
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(6) :: out_array
+  CCTK_INT, dimension(6) :: in_array
+  CCTK_REAL, dimension(:), allocatable :: metric_xx, metric_yy, metric_zz, metric_xy, metric_yz, metric_xz
+  
+  CCTK_INT, dimension(6), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL /)
+  ! Convert the interpolator parameter into a fortran string.
+  call CCTK_FortranString ( area_interp_len, area_interpolator, &
+                                             area_interp )
+
+!allocation of metric arrays
+  allocate(metric_xx(npoints), metric_yy(npoints), metric_zz(npoints), metric_xy(npoints), metric_xz(npoints), metric_yz(npoints))
+
+! Get an interpolation handle.
+  call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) )
+
+  if ( interp_handle .lt. 0 ) then
+    warn_message = 'Cannot get handle for interpolation. '
+    warn_message = trim(warn_message)//'Forgot to activate an implementation '
+    warn_message = trim(warn_message)//'providing interpolation operators?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  end if
+
+! Write the interpolation order parameter into a string
+  write(area_order,'(a6,i1)') 'order=',area_interpolation_order
+
+! Create a table from the interpolation order string.
+  call Util_TableCreateFromString ( table_handle, area_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
+  end if
+
+! Get a handle for the coordinate system.
+  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
+  if ( coord_system_handle .lt. 0) then
+    warn_message = 'Cannot get handle for cart3d coordinate system. '
+    warn_message = trim(warn_message)//'Forgot to activate an implementation '
+    warn_message = trim(warn_message)//'providing coordinates?'
+    call CCTK_WARN( 0, trim(warn_message) )
+  endif
+
+! Get the pointers to the interpolation points.
+  interp_coords(1) = CCTK_PointerTo(x_coord)
+  interp_coords(2) = CCTK_PointerTo(y_coord)
+  interp_coords(3) = CCTK_PointerTo(z_coord)
+
+! Get the pointers to the interpolation return arrays.
+  out_array(1) = CCTK_PointerTo(metric_xx)
+  out_array(2) = CCTK_PointerTo(metric_xy)
+  out_array(3) = CCTK_PointerTo(metric_xz)
+  out_array(4) = CCTK_PointerTo(metric_yy)
+  out_array(5) = CCTK_PointerTo(metric_yz)
+  out_array(6) = CCTK_PointerTo(metric_zz)
+!  out_array(7) = CCTK_PointerTo(psii)
+! Get the indices for the grid functions to be interpolated.
+  call CCTK_VarIndex ( in_array(1), 'admbase::metric_xx' )
+  call CCTK_VarIndex ( in_array(2), 'admbase::metric_xy' )
+  call CCTK_VarIndex ( in_array(3), 'admbase::metric_xz' )
+  call CCTK_VarIndex ( in_array(4), 'admbase::metric_yy' )
+  call CCTK_VarIndex ( in_array(5), 'admbase::metric_yz' )
+  call CCTK_VarIndex ( in_array(6), 'admbase::metric_zz' )
+
+  call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 npoints, CCTK_VARIABLE_REAL, &
+                                 interp_coords, 6, in_array(1:6), &
+                                 6, out_types(1:6), out_array(1:6) )
+  call area_computation(x_coord,y_coord,z_coord,t_index,metric_xx,metric_xy,metric_xz,metric_yy,metric_yz,metric_zz,npoints,area,local_contor)
+  deallocate(metric_xx, metric_yy, metric_zz, metric_xz, metric_yz, metric_xy)
+  end subroutine interpolation_area
+
+  subroutine area_computation(x_coord,y_coord,z_coord,t_index,g11,g12,g13,g22,g23,g33,npoints,area,contor)
+  use EHFinder_mod
+  implicit none
+  CCTK_INT,intent(in)::npoints,contor
+  CCTK_REAL,dimension(npoints),intent(in)::x_coord,y_coord,z_coord,t_index,g11,g12,g13,g22,g23,g33
+  CCTK_REAL::ds1,ds2,ds3,ds,dx1,dx2,dx3,dx_1,dx_2,dx_3,average_g11,average_g12,average_g13,average_g22,average_g33,average_g23,inverse_g11,inverse_g12,inverse_g13,inverse_g22,inverse_g33,inverse_g23,g
+  CCTK_INT::i,point1,point2,point3
+  CCTK_REAL,intent(out)::area
+  do i=1,contor,3
+!computing the determinant of the locale averaged metric (the metric is equal with the average of the metric in the corners of the triangle
+!we also use point1, point2, point3 as references from t_index that gives the corners of the current triangle  
+    point1=t_index(i)
+    point2=t_index(i+1)
+    point3=t_index(i+2)
+    average_g11=(g11(point1)+g11(point2)+g11(point3))/3
+    average_g12=(g12(point1)+g12(point2)+g12(point3))/3
+    average_g13=(g13(point1)+g13(point2)+g13(point3))/3
+    average_g22=(g22(point1)+g22(point2)+g22(point3))/3
+    average_g23=(g23(point1)+g23(point2)+g23(point3))/3
+    average_g33=(g33(point1)+g33(point2)+g33(point3))/3
+    g=average_g11*(average_g22*average_g33-average_g23*average_g23)-average_g12*(average_g12*average_g33-average_g23*average_g13)+average_g13*(average_g12*average_g23-average_g22*average_g13)
+!computing inverse of the g-matrix
+    inverse_g11=g*(average_g22*average_g33-average_g23*average_g23)
+    inverse_g22=g*(average_g11*average_g33-average_g13*average_g13)
+    inverse_g33=g*(average_g11*average_g22-average_g12*average_g12)
+    inverse_g12=g*(average_g13*average_g23-average_g12*average_g33)
+    inverse_g13=g*(average_g12*average_g23-average_g13*average_g22)
+    inverse_g23=g*(average_g13*average_g12-average_g11-average_g23)
+!computing length elements and area of an elementar triangle
+    dx1=x_coord(point2)-x_coord(point1)
+    dx2=y_coord(point2)-y_coord(point1)
+    dx3=z_coord(point2)-z_coord(point1)
+    dx_1=x_coord(point3)-x_coord(point1)
+    dx_2=y_coord(point3)-y_coord(point1)
+    dx_3=z_coord(point3)-z_coord(point1)
+    ds1=(1/4)*(g**0.5)*(dx2*dx_3-dx3*dx_2)
+    ds2=(1/4)*(g**0.5)*(dx3*dx_1-dx1*dx_3)
+    ds3=(1/4)*(g**0.5)*(dx1*dx_2-dx2*dx_1)
+    ds=(1/2)*(inverse_g11*ds1*ds1+inverse_g22*ds2*ds2+inverse_g33*ds3*ds3+2*inverse_g12*ds1*ds2+2*inverse_g13*ds1*ds3+2*inverse_g23*ds2*ds3)
+    area=area+ds
+  end do
+  end subroutine area_computation
+
+subroutine triangularization(l,contor,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,grid_value,grid_x,grid_y,grid_z)
+!subroutine triangularization(l,contor,Nx,Ny,Nz,grid_value,grid_x,grid_y,grid_z)
 
   use IsoSurface_mod
 
@@ -58,12 +195,12 @@ subroutine triangularization(l,contor,Nx,Ny,Nz,grid_value,grid_x,grid_y,grid_z)
   CCTK_REAL,dimension(Nx),intent(in)::grid_x
   CCTK_REAL,dimension(Ny),intent(in)::grid_y
   CCTK_REAL,dimension(Nz),intent(in)::grid_z
-  CCTK_INT::ii,jj,kk,index,sum,nsize,boolean_finder1,boolean_finder2,boolean_finder3,t,equality,index2,pp
+  CCTK_INT::ii,jj,kk,index,sum,nsize,boolean_finder1,boolean_finder2,boolean_finder3,t,equality,index2,pp,j,i,k
   CCTK_INT,intent(out)::l,contor
-  CCTK_REAL, dimension(:),allocatable,intent(out):: coord_x,coord_y,coord_z
-  CCTK_INT, dimension (:),allocatable,intent(out)::triangle_index
-  CCTK_INT,dimension(:),allocatable::resize_triangle_index
-  CCTK_REAL,dimension(:),allocatable::resize_coord_x,resize_coord_y,resize_coord_z
+  CCTK_REAL, dimension(:),allocatable,intent(out):: x_coord,y_coord,z_coord
+  CCTK_INT, dimension (:),allocatable,intent(out)::t_index
+  CCTK_INT,dimension(:),allocatable::resize_t_index
+  CCTK_REAL,dimension(:),allocatable::resize_x_coord,resize_y_coord,resize_z_coord
 !we store table with triangles in a 2 dimensional constant array using reshape function; Fortran stores column-like and we need row-like
   CCTK_INT, dimension (256,16):: a
 
@@ -430,8 +567,8 @@ CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
 
   
   !allocating initial space for coordinate and triangle index arrays.These arrays will be extended if needed during the execution of the program using allocate and deallocate functions.
-  nsize=500000
-  allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize))
+  nsize=50000
+  allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize))
 
   !marching through the grid and computing "sum" for each individual cube. "l" is the contor for the vectors that keeps the coordinates of intersection points. we will also produce a vector called triangle_index that tells which points build which triangles.
   l=0
@@ -572,65 +709,65 @@ CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
                   boolean_finder3=0
                   if (l>1) then
                     do t=1,l
-                      if ((abs(temp_x(1)-coord_x(t))<1e-10).and.(abs(temp_y(1)-coord_y(t))<1e-10).and.(abs(temp_z(1)-coord_z(t))<1e-10)) then
+                      if ((abs(temp_x(1)-x_coord(t))<1e-10).and.(abs(temp_y(1)-y_coord(t))<1e-10).and.(abs(temp_z(1)-z_coord(t))<1e-10)) then
                         boolean_finder1=1
                         contor=contor+1
-                        triangle_index(contor)=t-1
+                        t_index(contor)=t-1
                       end if
 
-                      if ((abs(temp_x(2)-coord_x(t))<1e-10).and.(abs(temp_y(2)-coord_y(t))<1e-10).and.(abs(temp_z(2)-coord_z(t))<1e-10)) then
+                      if ((abs(temp_x(2)-x_coord(t))<1e-10).and.(abs(temp_y(2)-y_coord(t))<1e-10).and.(abs(temp_z(2)-z_coord(t))<1e-10)) then
                         boolean_finder2=1
                         contor=contor+1
-                        triangle_index(contor)=t-1
+                        t_index(contor)=t-1
                       end if
-                      if ((abs(temp_x(3)-coord_x(t))<1e-10).and.(abs(temp_y(3)-coord_y(t))<1e-10).and.(abs(temp_z(3)-coord_z(t))<1e-10)) then
+                      if ((abs(temp_x(3)-x_coord(t))<1e-10).and.(abs(temp_y(3)-y_coord(t))<1e-10).and.(abs(temp_z(3)-z_coord(t))<1e-10)) then
                         boolean_finder3=1
                         contor=contor+1
-                        triangle_index(contor)=t-1
+                        t_index(contor)=t-1
                       end if
                     end do
                   end if
                   if (boolean_finder1==0) then
                     l=l+1
-                    coord_x(l)=temp_x(1)
-                    coord_y(l)=temp_y(1)
-                    coord_z(l)=temp_z(1)
+                    x_coord(l)=temp_x(1)
+                    y_coord(l)=temp_y(1)
+                    z_coord(l)=temp_z(1)
                     contor=contor+1
-                    triangle_index(contor)=l-1
+                    t_index(contor)=l-1
                   end if
                   if (boolean_finder2==0) then
                     l=l+1
-                    coord_x(l)=temp_x(2)
-                    coord_y(l)=temp_y(2)
-                    coord_z(l)=temp_z(2)
+                    x_coord(l)=temp_x(2)
+                    y_coord(l)=temp_y(2)
+                    z_coord(l)=temp_z(2)
                     contor=contor+1
-                    triangle_index(contor)=l-1
+                    t_index(contor)=l-1
                   end if
                   if (boolean_finder3==0) then
                     l=l+1
-                    coord_x(l)=temp_x(3)
-                    coord_y(l)=temp_y(3)
-                    coord_z(l)=temp_z(3)
+                    x_coord(l)=temp_x(3)
+                    y_coord(l)=temp_y(3)
+                    z_coord(l)=temp_z(3)
                     contor=contor+1
-                    triangle_index(contor)=l-1
+                    t_index(contor)=l-1
                   end if
                 end if
                
 
               if ((nsize-l)<10) then
-                allocate(resize_coord_x(nsize),resize_coord_y(nsize),resize_coord_z(nsize),resize_triangle_index(nsize*5))
-                resize_coord_x(1:l)=coord_x(1:l)
-                resize_coord_y(1:l)=coord_y(1:l)
-                resize_coord_z(1:l)=coord_z(1:l)
-                resize_triangle_index(1:contor)=triangle_index(1:contor)
-                deallocate(coord_x,coord_y,coord_z,triangle_index)
+                allocate(resize_x_coord(nsize),resize_y_coord(nsize),resize_z_coord(nsize),resize_t_index(5*nsize))
+                resize_x_coord(1:l)=x_coord(1:l)
+                resize_y_coord(1:l)=y_coord(1:l)
+                resize_z_coord(1:l)=z_coord(1:l)
+                resize_t_index(1:contor)=t_index(1:contor)
+                deallocate(x_coord,y_coord,z_coord,t_index)
                 nsize=nsize*2
-                allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize*5))
-                coord_x(1:l)=resize_coord_x(1:l)
-                coord_y(1:l)=resize_coord_y(1:l)
-                coord_z(1:l)=resize_coord_z(1:l)
-                triangle_index(1:contor)=resize_triangle_index(1:contor)
-                deallocate(resize_coord_x,resize_coord_y,resize_coord_z,resize_triangle_index)
+                allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(nsize*5))
+                x_coord(1:l)=resize_x_coord(1:l)
+                y_coord(1:l)=resize_y_coord(1:l)
+                z_coord(1:l)=resize_z_coord(1:l)
+                t_index(1:contor)=resize_t_index(1:contor)
+                deallocate(resize_x_coord,resize_y_coord,resize_z_coord,resize_t_index)
               end if
          
               end if