First step in making isosurface area calculation work.

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

1 files changed, 45 insertions, 19 deletions

diff --git a/src/EHFinder_IsoSurface_optimized.F90 b/src/EHFinder_IsoSurface_optimized.F90
index 6b49d16..90ba5d6 100755
--- a/src/EHFinder_IsoSurface_optimized.F90
+++ b/src/EHFinder_IsoSurface_optimized.F90
@@ -33,11 +33,19 @@ subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
   allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize))
 
   do p = 1, eh_number_level_sets
-     call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, &
-                            t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:))
+    call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, &
+                           t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:))
+
+    print*, 'contor = ', t_index(1:contor)
+    print*,'The surface consists of', contor/3,'triangles and',l,'coordinates'
+
+    call interpolation_area(cctkGH, l, x_coord,y_coord,z_coord, contor, &
+                            t_index, eh_area(1,p) )
+
+    print*,'The area is ', eh_area(1,p)   
   end do
+
   
-  print*,'The surface consists of', contor/3,'triangles and',l,'coordinates'
 
 end subroutine EHFinder_IsoSurface
 
@@ -52,8 +60,8 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
 
   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, dimension(:), pointer :: x_coord,y_coord,z_coord
+  CCTK_INT, dimension(:), pointer :: t_index
   CCTK_REAL, intent(OUT) :: area
   character(len=200) :: area_interp
   character(len=128) :: warn_message
@@ -120,31 +128,48 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
   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_VarIndex ( in_array(1), 'admbase::gxx' )
+  call CCTK_VarIndex ( in_array(2), 'admbase::gxy' )
+  call CCTK_VarIndex ( in_array(3), 'admbase::gxz' )
+  call CCTK_VarIndex ( in_array(4), 'admbase::gyy' )
+  call CCTK_VarIndex ( in_array(5), 'admbase::gyz' )
+  call CCTK_VarIndex ( in_array(6), 'admbase::gzz' )
 
   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)
+  
+  call area_computation(x_coord, y_coord, z_coord, t_index, &
+                        metric_xx, metric_xy, metric_xz, &
+                        metric_yy, metric_yz, metric_zz, &
+                        npoints, l, area )
+
   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)
+  subroutine area_computation(x_coord, y_coord, z_coord, t_index, &
+                              g11, g12, g13, &
+                              g22, g23, g33, &
+                              npoints, l, area )
+
   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_INT,intent(in) :: npoints, l
+  CCTK_REAL,dimension(:),pointer :: x_coord, y_coord, z_coord
+  CCTK_INT, dimension(:),pointer :: t_index
+  CCTK_REAL, dimension(npoints), intent(in) :: 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
+
+  print*, npoints, l
+  print*, t_index(1:l)
+  print*, 'x = ',x_coord
+  print*, 'g11 = ',g11
+  do i=1,l,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)
@@ -180,7 +205,7 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
   end subroutine area_computation
 
   
-subroutine triangularization(n_coords,contor,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,grid_value,grid_x,grid_y,grid_z)
+subroutine triangularization(n_coords,n_triangles,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,x_coord,y_coord,z_coord,t_index,f(:,:,:,1),x(:,1,1),y(1,:,1),z(1,1,:))
 
   use IsoSurface_mod
@@ -193,8 +218,8 @@ subroutine triangularization(n_coords,contor,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_
   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,j,i,k,l,n_triangles
-  CCTK_INT,intent(out)::n_coords,contor
+  CCTK_INT::ii,jj,kk,index,sum,nsize,boolean_finder1,boolean_finder2,boolean_finder3,t,equality,index2,pp,j,i,k,l,contor
+  CCTK_INT,intent(out)::n_coords,n_triangles
   CCTK_REAL, dimension(:),pointer:: x_coord,y_coord,z_coord
   CCTK_INT, dimension (:),pointer::t_index
   CCTK_INT,dimension(:),allocatable::resize_t_index
@@ -795,4 +820,5 @@ integer, dimension (16,16):: b9(16,16)=reshape&
   close(10)
   close(11)
   close(12)
+  print*,'t_index = ',t_index(1:n_triangles)
   end subroutine triangularization