From 6155a4231cf69c74670edecc6765beceb96b78f9 Mon Sep 17 00:00:00 2001 From: florin Date: Sat, 12 Mar 2005 01:42:44 +0000 Subject: 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 --- src/EHFinder_IsoSurface_optimized.F90 | 64 ++++++++++++++++++++++++----------- 1 file 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 -- cgit v1.2.3