From 9ef994cea210cf0789d7fae512f2b3f8cd33d890 Mon Sep 17 00:00:00 2001 From: florin Date: Mon, 11 Oct 2004 22:12:09 +0000 Subject: 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 --- src/EHFinder_IsoSurface_optimized.F90 | 223 +++++++++++++++++++++++++++------- 1 file changed, 180 insertions(+), 43 deletions(-) (limited to 'src') 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 -- cgit v1.2.3