Cleaned version.

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

1 files changed, 60 insertions, 39 deletions

diff --git a/src/EHFinder_IsoSurface_optimized.F90 b/src/EHFinder_IsoSurface_optimized.F90
index 1733e20..6b49d16 100755
--- a/src/EHFinder_IsoSurface_optimized.F90
+++ b/src/EHFinder_IsoSurface_optimized.F90
@@ -17,12 +17,10 @@ subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: n_coords, n_triangles
   CCTK_REAL, dimension(:), pointer :: x_coord, y_coord, z_coord
   CCTK_INT, dimension(:), pointer :: t_index
-  CCTK_REAL, parameter :: iso_value = zero
-  CCTK_INT :: l,contor
-
+  CCTK_INT :: l,contor,p,nsize
+  
 ! Find index ranges for points excluding ghost zones and symmetry points for
 ! the 3D grid functions.
 #include "include/physical_part.h"
@@ -30,18 +28,20 @@ subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)
   call CCTK_INFO ('Entered IsoSurface')
   
   nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3)
-  print*,nx,ny,nz
-  print*,size(f)
-  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(l,contor,nx,ny,nz,x_coord,y_coord,z_coord,t_index,f(:,:,:,l),x(:,1,1),y(1,:,1),z(1,1,:))
-  print*,'The surface consists of ',n_triangles/3,' triangles and ',n_coords,' coords'
+
+  nsize=1000
+  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,:))
   end do
+  
+  print*,'The surface consists of', contor/3,'triangles and',l,'coordinates'
+
 end subroutine EHFinder_IsoSurface
 
+
 subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area)
   
   use EHFinder_mod
@@ -135,7 +135,7 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
   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
@@ -179,11 +179,11 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a
   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)
+  
+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(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
-
   implicit none 
 !declaring variables that will be used in the program
   CCTK_INT,intent(in)::Nx,Ny,Nz
@@ -193,16 +193,17 @@ subroutine triangularization(l,contor,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,g
   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
-  CCTK_INT,intent(out)::l,contor
+  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_REAL, dimension(:),pointer:: x_coord,y_coord,z_coord
   CCTK_INT, dimension (:),pointer::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
 
-  CCTK_INT, dimension (16,30):: b1(16,30)=reshape&
+  integer, dimension (16,30):: b1(16,30)=reshape&
 ((/-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
@@ -238,7 +239,7 @@ subroutine triangularization(l,contor,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,g
 
 
 
-CCTK_INT, dimension (16,30):: b2(16,30)=reshape&
+integer, dimension (16,30):: b2(16,30)=reshape&
 ((/4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1,&
 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1,&
 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
@@ -272,7 +273,7 @@ CCTK_INT, dimension (16,30):: b2(16,30)=reshape&
 
 
 
-CCTK_INT, dimension (16,30):: b3(16,30)=reshape&
+integer, dimension (16,30):: b3(16,30)=reshape&
 ((/9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1,&
 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1,&
 11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1,&
@@ -306,7 +307,7 @@ CCTK_INT, dimension (16,30):: b3(16,30)=reshape&
 
 
 
-CCTK_INT, dimension (16,30):: b4(16,30)=reshape&
+integer, dimension (16,30):: b4(16,30)=reshape&
 ((/0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1,&
 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1,&
 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1,&
@@ -341,7 +342,7 @@ CCTK_INT, dimension (16,30):: b4(16,30)=reshape&
 
 
 
-CCTK_INT, dimension (16,30):: b5(16,30)=reshape&
+integer, dimension (16,30):: b5(16,30)=reshape&
 ((/2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1,&
 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1,&
 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1,&
@@ -375,7 +376,7 @@ CCTK_INT, dimension (16,30):: b5(16,30)=reshape&
 
 
 
-CCTK_INT, dimension (16,30):: b6(16,30)=reshape&
+integer, dimension (16,30):: b6(16,30)=reshape&
 ((/4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1,&
 10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1,&
 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1,&
@@ -408,7 +409,7 @@ CCTK_INT, dimension (16,30):: b6(16,30)=reshape&
 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1/),(/16,30/))
 
 
-CCTK_INT, dimension (16,30):: b7(16,30)=reshape&
+integer, dimension (16,30):: b7(16,30)=reshape&
 ((/1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1,&
 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1,&
 11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1,&
@@ -441,7 +442,7 @@ CCTK_INT, dimension (16,30):: b7(16,30)=reshape&
 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1/),(/16,30/))
 
 
-CCTK_INT, dimension (16,30):: b8(16,30)=reshape&
+integer, dimension (16,30):: b8(16,30)=reshape&
 ((/0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1,&
 10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1,&
 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1,&
@@ -475,7 +476,7 @@ CCTK_INT, dimension (16,30):: b8(16,30)=reshape&
 
 
 
-CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
+integer, dimension (16,16):: b9(16,16)=reshape&
 ((/9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1,&
 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1,&
@@ -565,10 +566,10 @@ 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=50000
+  nsize=1000
   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.
+  !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 t_index that tells which points build which triangles.
   l=0
   contor=0
   do k=1,Nz-1
@@ -752,20 +753,21 @@ CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
                 end if
                
 
-              if ((nsize-l)<10) then
-                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)
+              if ((nsize-l)<20) then
+                allocate(resize_coord_x(nsize),resize_coord_y(nsize), &
+                         resize_coord_z(nsize),resize_t_index(nsize*5))
+                resize_coord_x(1:l)=x_coord(1:l)
+                resize_coord_y(1:l)=y_coord(1:l)
+                resize_coord_z(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(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)
+                x_coord(1:l)=resize_coord_x(1:l)
+                y_coord(1:l)=resize_coord_y(1:l)
+                z_coord(1:l)=resize_coord_z(1:l)
                 t_index(1:contor)=resize_t_index(1:contor)
-                deallocate(resize_x_coord,resize_y_coord,resize_z_coord,resize_t_index)
+                deallocate(resize_coord_x,resize_coord_y,resize_coord_z,resize_t_index)
               end if
          
               end if
@@ -774,4 +776,23 @@ CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
         end do
       end do
     end do 
+!writing data to the files-temporary output method
+  open (unit=9,file='output_x.txt',status='replace',action='write')
+  open (unit=10,file='output_y.txt',status='replace',action='write')
+  open (unit=11,file='output_z.txt',status='replace',action='write')
+  open (unit=12,file='triangle_test.txt',status='replace',action='write')  
+  n_triangles=contor
+  n_coords=l
+  do contor=1,n_triangles
+    write(12,*)t_index(contor) 
+  end do
+  do l=1,n_coords
+    write(9,*)x_coord(l)
+    write(10,*)y_coord(l)
+    write(11,*)z_coord(l)
+  end do
+  close(9)
+  close(10)
+  close(11)
+  close(12)
   end subroutine triangularization