From 09f472c557fd21350f5866c8152d4d6d1b5bfca2 Mon Sep 17 00:00:00 2001 From: florin Date: Fri, 11 Mar 2005 23:45:40 +0000 Subject: Cleaned version. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@195 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_IsoSurface_optimized.F90 | 99 +++++++++++++++++++++-------------- 1 file 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 -- cgit v1.2.3