!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! REMOVE_POINT May 2008 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine remove_point (iii,surface,flag0,surface0)
use definitions
implicit none
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine removes the node iii from the surface passed as argument.
! It especially ensures that the convex hull is kept intact, and repairs
! the triangulation (it will be Delaunay checked later)
! Bear in mind that if flag0=0, then the surface0 arrays are deallocated,
! so all operations on these arrays have to be toggled off.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
integer iii,flag0
type(sheet) surface,surface0
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
logical already_in_list,iii_found,jjj_found
double precision, dimension(:), allocatable :: new_x,new_y,new_z,new_r,new_s
double precision, dimension(:), allocatable :: new_xn,new_yn,new_zn
double precision, dimension(:), allocatable :: new_x0,new_y0,new_z0,new_r0,new_s0
double precision, dimension(:), allocatable :: new_xn0,new_yn0,new_zn0
double precision :: xiii,yiii,ziii,xjjj,yjjj,zjjj,dist,distt
integer, dimension(:,:), allocatable :: new_icon
integer,dimension(:), allocatable :: list_test
integer new_nt, new_nsurface
integer list(20),nb_in_list,iptcl,nb_test,jjj
integer tr_counter,i,j,k,kk,counter,nside
integer ierr,iproc,nproc
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
!nside=2**surface%levelt+1
nside=2**surface%levelt+2 !opla
if (iproc.eq.0) write(8,*) '=================================================='
if (iproc.eq.0) write(8,*) 'removing point iii=',iii
!===============================================================================
!=====[find how many triangles have point iii as a vertex]======================
!===============================================================================
nb_in_list=1
list(nb_in_list)=iii
tr_counter=0
do i=1,surface%nt
do k=1,3
if (surface%icon(k,i)==iii) then
tr_counter=tr_counter+1
do kk=1,3
iptcl=surface%icon(kk,i)
already_in_list=.false.
do j=1,nb_in_list
if (list(j)==iptcl) already_in_list=.true.
end do
if (.not.already_in_list) then
nb_in_list=nb_in_list+1
list(nb_in_list)=iptcl
end if
end do
end if
end do
end do
if(iproc.eq.0) write(8,*) 'nb of neighbouring triangles',tr_counter
if(iproc.eq.0) write(8,*) 'nb_in_list=',nb_in_list
if(iproc.eq.0) write(8,*) 'list of iii neighbours:',list(1:nb_in_list)
!===============================================================================
!=====[build list of neighbouring vertices that are not on convex hull]=========
!===============================================================================
nb_test=count(list(1:nb_in_list)>(4*nside-4))
allocate(list_test(nb_test))
list_test=-10
counter=0
do i=1,nb_in_list
j=list(i)
if(j>4*nside-4) then
counter=counter+1
list_test(counter)=j
end if
end do
if(iproc.eq.0) write(8,*) 'nb_test=',nb_test
if(iproc.eq.0) write(8,*) 'list test=',list_test
if(counter/=nb_test) call stop_run('pb in test')
!===============================================================================
!=====[find vertex closest to iii]==============================================
!===============================================================================
dist=10.d0
do i=1,nb_test
j=list_test(i)
xiii=surface%x(iii)
yiii=surface%y(iii)
ziii=surface%z(iii)
xjjj=surface%x(j)
yjjj=surface%y(j)
zjjj=surface%z(j)
distt=sqrt((xiii-xjjj)**2+(yiii-yjjj)**2+(ziii-zjjj)**2)
if (iii/=j .and. distt<dist) then
jjj=j
dist=distt
end if
end do
if (iproc.eq.0) write(8,*) 'closest point to iii is jjj=',jjj
!===============================================================================
!=====[carry out point removal]=================================================
!===============================================================================
new_nsurface=surface%nsurface-1
new_nt=surface%nt-2
allocate(new_icon(3,new_nt))
allocate(new_x(new_nsurface))
allocate(new_y(new_nsurface))
allocate(new_z(new_nsurface))
allocate(new_xn(new_nsurface))
allocate(new_yn(new_nsurface))
allocate(new_zn(new_nsurface))
allocate(new_r(new_nsurface))
allocate(new_s(new_nsurface))
if (flag0==1) then
allocate(new_x0(new_nsurface))
allocate(new_y0(new_nsurface))
allocate(new_z0(new_nsurface))
allocate(new_xn0(new_nsurface))
allocate(new_yn0(new_nsurface))
allocate(new_zn0(new_nsurface))
allocate(new_r0(new_nsurface))
allocate(new_s0(new_nsurface))
end if
!=====[update icon]=====
tr_counter=0
do i=1,surface%nt
iii_found=.false.
jjj_found=.false.
do k=1,3
if (surface%icon(k,i)==iii) iii_found=.true.
if (surface%icon(k,i)==jjj) jjj_found=.true.
end do
if (.not.(iii_found .and. jjj_found)) then
tr_counter=tr_counter+1
new_icon(1:3,tr_counter)=surface%icon(1:3,i)
do k=1,3
if (new_icon(k,tr_counter)==iii) new_icon(k,tr_counter)=jjj
if (new_icon(k,tr_counter)>iii) new_icon(k,tr_counter)=new_icon(k,tr_counter)-1
end do
end if
end do
!=====[update coordinates]=====
counter=0
do i=1,surface%nsurface
if (i/=iii) then
counter=counter+1
new_x (counter)=surface%x(i)
new_y (counter)=surface%y(i)
new_z (counter)=surface%z(i)
new_xn (counter)=surface%xn(i)
new_yn (counter)=surface%yn(i)
new_zn (counter)=surface%zn(i)
new_r (counter)=surface%r(i)
new_s (counter)=surface%s(i)
if (flag0==1) then
new_x0 (counter)=surface0%x(i)
new_y0 (counter)=surface0%y(i)
new_z0 (counter)=surface0%z(i)
new_xn0(counter)=surface0%xn(i)
new_yn0(counter)=surface0%yn(i)
new_zn0(counter)=surface0%zn(i)
new_r0 (counter)=surface0%r(i)
new_s0 (counter)=surface0%s(i)
end if
end if
end do
deallocate(surface%x)
deallocate(surface%y)
deallocate(surface%z)
deallocate(surface%xn)
deallocate(surface%yn)
deallocate(surface%zn)
deallocate(surface%icon)
deallocate(surface%r)
deallocate(surface%s)
if (flag0==1) then
deallocate(surface0%x)
deallocate(surface0%y)
deallocate(surface0%z)
deallocate(surface0%xn)
deallocate(surface0%yn)
deallocate(surface0%zn)
deallocate(surface0%r)
end if
surface%nt=new_nt
surface%nsurface=new_nsurface
allocate(surface%x(surface%nsurface))
allocate(surface%y(surface%nsurface))
allocate(surface%z(surface%nsurface))
allocate(surface%xn(surface%nsurface))
allocate(surface%yn(surface%nsurface))
allocate(surface%zn(surface%nsurface))
allocate(surface%icon(3,surface%nt))
allocate(surface%r(surface%nsurface))
allocate(surface%s(surface%nsurface))
if (flag0==1) then
surface0%nsurface=new_nsurface
allocate(surface0%x(surface0%nsurface))
allocate(surface0%y(surface0%nsurface))
allocate(surface0%z(surface0%nsurface))
allocate(surface0%xn(surface0%nsurface))
allocate(surface0%yn(surface0%nsurface))
allocate(surface0%zn(surface0%nsurface))
allocate(surface0%r(surface0%nsurface))
allocate(surface0%s(surface0%nsurface))
end if
surface%x=new_x
surface%y=new_y
surface%z=new_z
surface%xn=new_xn
surface%yn=new_yn
surface%zn=new_zn
surface%r=new_r
surface%s=new_s
surface%icon=new_icon
if (flag0==1) then
surface0%x=new_x0
surface0%y=new_y0
surface0%z=new_z0
surface0%xn=new_xn0
surface0%yn=new_yn0
surface0%zn=new_zn0
surface0%r=new_r0
surface0%s=new_s0
end if
deallocate(new_x)
deallocate(new_y)
deallocate(new_z)
deallocate(new_icon)
deallocate(new_xn)
deallocate(new_yn)
deallocate(new_zn)
deallocate(new_r)
deallocate(new_s)
if (flag0==1) then
deallocate(new_x0)
deallocate(new_y0)
deallocate(new_z0)
deallocate(new_xn0)
deallocate(new_yn0)
deallocate(new_zn0)
deallocate(new_r0)
deallocate(new_s0)
end if
deallocate(list_test)
end subroutine
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------