Skip to content
Snippets Groups Projects
Commit 17a6410f authored by Matthias Schmiddunser's avatar Matthias Schmiddunser
Browse files

Merge branch 'surface2' into douarw_rotsub

parents a2aa187e ec9275c8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/check_delaunay.f90
+ 50
− 76
View file @ 17a6410f
......@@ -75,7 +75,6 @@ integer nnmax,nflip,counter
integer ie,inode,kp,kpp
integer inodep,inodepp
integer iproc,nproc,ierr
integer checks
integer i1,i2,j1,j2,k,k12,k21,n1,n2,t1,t2,nn,j,jj,err
integer,dimension(:) ,allocatable :: nedgepernode
integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
......@@ -97,82 +96,57 @@ shift=' '
!------------------------
! build edge information
!------------------------
checksurf: do checks = 1,2
nsurfacen=surface%nsurface
ntriangle=surface%nt
nedge=nsurfacen+ntriangle-1
if (surface%closed.eq.1) nedge=nedge-1
nnmax=12 ! nnmax is maximum number of neighbours in triangulation (should be around 6 on average)
allocate (ed(nedge),stat=err)
allocate (nedgepernode(nsurfacen),stat=err)
allocate (nodenodenumber(nnmax,nsurfacen),stat=err)
allocate (nodeedgenumber(nnmax,nsurfacen),stat=err)
ed%t1=0
ed%t2=0
nedgepernode=0
iedge=0
do ie=1,ntriangle
do k=1,3
inode=surface%icon(k,ie)
kp=mod(k,3)+1
kpp=mod(k+1,3)+1
inodep=surface%icon(kp,ie)
inodepp=surface%icon(kpp,ie)
do j=1,nedgepernode(inodep)
if (nodenodenumber(j,inodep).eq.inode) then
jedge=nodeedgenumber(j,inodep)
ed(jedge)%m2=inodepp
ed(jedge)%t2=ie
goto 111
endif
enddo
iedge=iedge+1
if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
nedgepernode(inode)=nedgepernode(inode)+1
if (nedgepernode(inode).gt.nnmax) then
!retriangulate only on once
if (checks==1 .and. surface%closed/=1) then
if (iproc.eq.0) then
write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
write (8,*) 'retriangulate surface'
write (*,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
write (*,*) 'retriangulate surface'
end if
call retriangulate_surface (surface,is,istep,params)
!deallocate for next check run
deallocate (ed)
deallocate (nedgepernode)
deallocate (nodenodenumber)
deallocate (nodeedgenumber)
!redo check
cycle checksurf
else
!give up on second check
if (iproc.eq.0) then
write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
write (*,*) 'Node ',inode,' has ',nedgepernode(inode),'neighbors'
end if
call stop_run ('too many neighbours')
end if
endif
nodenodenumber(nedgepernode(inode),inode)=inodep
nodeedgenumber(nedgepernode(inode),inode)=iedge
ed(iedge)%n1=inode
ed(iedge)%n2=inodep
ed(iedge)%m1=inodepp
ed(iedge)%t1=ie
111 continue
enddo
enddo
!exit after first run if successfully completed
exit checksurf
end do checksurf
nsurfacen=surface%nsurface
ntriangle=surface%nt
nedge=nsurfacen+ntriangle-1
if (surface%closed.eq.1) nedge=nedge-1
nnmax=12 ! nnmax is maximum number of neighbours in triangulation (should be around 6 on average)
allocate (ed(nedge),stat=err)
allocate (nedgepernode(nsurfacen),stat=err)
allocate (nodenodenumber(nnmax,nsurfacen),stat=err)
allocate (nodeedgenumber(nnmax,nsurfacen),stat=err)
ed%t1=0
ed%t2=0
nedgepernode=0
iedge=0
do ie=1,ntriangle
do k=1,3
inode=surface%icon(k,ie)
kp=mod(k,3)+1
kpp=mod(k+1,3)+1
inodep=surface%icon(kp,ie)
inodepp=surface%icon(kpp,ie)
do j=1,nedgepernode(inodep)
if (nodenodenumber(j,inodep).eq.inode) then
jedge=nodeedgenumber(j,inodep)
ed(jedge)%m2=inodepp
ed(jedge)%t2=ie
goto 111
endif
enddo
iedge=iedge+1
if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
nedgepernode(inode)=nedgepernode(inode)+1
if (nedgepernode(inode).gt.nnmax) then
if (iproc.eq.0) then
write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
write (*,*) 'Node ',inode,' has ',nedgepernode(inode),'neighbors'
end if
call stop_run ('too many neighbours')
endif
nodenodenumber(nedgepernode(inode),inode)=inodep
nodeedgenumber(nedgepernode(inode),inode)=iedge
ed(iedge)%n1=inode
ed(iedge)%n2=inodep
ed(iedge)%m1=inodepp
ed(iedge)%t1=ie
111 continue
enddo
enddo
if (iedge/=nedge) then
if (iproc.eq.0) write (8,*) 'iedge=',iedge,'nedge=',nedge
......
src/move_surface.f90
+ 63
− 34
View file @ 17a6410f
......@@ -69,9 +69,10 @@ integer istep,is
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,iter,nedge,nrem
integer i,j,i1,iter,nedge,nrem,counter
integer err,nside
integer iproc,nproc,ierr
integer,dimension(:),allocatable :: remove_list
double precision u,v,w,rat,eps
double precision x0,y0,z0,xi,yi,zi,x,y,z
double precision,dimension(:),allocatable::xbuf,ybuf,zbuf
......@@ -81,7 +82,7 @@ double precision ux,uy,uz,vx,vy,vz,wx,wy,wz,xyzn
double precision xn0,yn0,zn0,xyn,theta,phi
double precision xn1,yn1,zn1,xn2,yn2,zn2
double precision xln1,yln1,zln1,xln2,yln2,zln2
logical remove
logical,dimension(:),allocatable :: remove
character*72 :: shift
!-------------------------------------------------------------------------------
......@@ -92,12 +93,12 @@ call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
allocate (xbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xbuf in move_surface$')
allocate (ybuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ybuf in move_surface$')
allocate (zbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zbuf in move_surface$')
allocate (xnbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xnbuf in move_surface$')
allocate (ynbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ynbuf in move_surface$')
allocate (znbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc znbuf in move_surface$')
allocate (xbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xbuf in move_surface')
allocate (ybuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ybuf in move_surface')
allocate (zbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zbuf in move_surface')
allocate (xnbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xnbuf in move_surface')
allocate (ynbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ynbuf in move_surface')
allocate (znbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc znbuf in move_surface')
xbuf=0.d0
ybuf=0.d0
......@@ -145,7 +146,7 @@ do i=iproc+1,surface%nsurface,nproc
ov%wnodecompact,w)
else
! We have a problem if we get here...
call stop_run ('Error: Unsupported value for flag1 in move_surface$')
call stop_run ('Error: Unsupported value for flag1 in move_surface')
endif
! This section is by Cedric and is very restrictive in the sense that it
......@@ -233,7 +234,7 @@ do i=iproc+1,surface%nsurface,nproc
vz,ov%wnodecompact,w,wx,wy,wz)
else
! We have a problem if we get here...
call stop_run ('Error: Unsupported value for flag1 in move_surface$')
call stop_run ('Error: Unsupported value for flag1 in move_surface')
endif
xn0=surface%xn(i)
......@@ -376,37 +377,65 @@ endif
! Again, disabled for closed surfaces
if (surface%closed.eq.0) then
nrem=0
remove=.true.
do while (remove)
remove=.false.
allocate (remove(surface%nsurface))
remove = .false.
!find all nodes outside of the domain
do i=4*nside-3,surface%nsurface
if (surface%x(i).lt.-eps .or. surface%x(i).gt.1.d0+eps .or. &
surface%y(i).lt.-eps .or. surface%y(i).gt.1.d0+eps .or. &
surface%z(i).lt.-eps .or. surface%z(i).gt.1.d0+eps) then
call remove_point(i,surface,flag0,surface0)
nrem=nrem+1
remove=.true.
exit
endif
if (surface%x(i).lt.-eps .or. surface%x(i).gt.1.d0+eps .or. &
surface%y(i).lt.-eps .or. surface%y(i).gt.1.d0+eps .or. &
surface%z(i).lt.-eps .or. surface%z(i).gt.1.d0+eps) then
remove(i)=.true.
endif
end do
end do
if (iproc.eq.0 .and. nrem/=0) write(*,'(a,i4,a)') shift//'removing',nrem,' particle from surface'
!create remove_list as in refine_surface
nrem=count(remove)
allocate (remove_list(nrem))
remove_list=0
call check_delaunay (params,surface,is,istep,'delaunay_aft_remove',0)
counter=0
do i=4*nside-3,surface%nsurface
if (remove(i)) then
counter = counter+1
remove_list(counter) = i
end if
end do
if (.not.params%normaladvect) then
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if
deallocate (remove)
if (flag0==0) then
nedge=surface%nsurface+surface%nt-1
call DoRuRe_surf_stats (params%doDoRuRe,istep,50,is,surface%nt,surface%nsurface,nedge,-nrem)
end if
!remove points one by one
do i=1,nrem
i1=remove_list(i)
call remove_point(i1,surface,flag0,surface0,remove_list,nrem)
!Modify remove_list to accout for removal of node from surface
remove_list(i) = 0
do j=i+1,nrem
if (remove_list(j) > i1) remove_list(j) = remove_list(j)-1
end do
end do
if (iproc.eq.0) write (8,*) 'Removal complete'
deallocate (remove_list)
if (iproc.eq.0 .and. params%debug>=1) &
write(*,'(a,i2,a,i8,a)') shift//'S.',is,':', nrem,' removed ptcls in refine_surface'
if (iproc.eq.0 .and. nrem/=0) write(*,'(a,i4,a)') shift//'removing',nrem,' particle from surface'
call check_delaunay (params,surface,is,istep,'delaunay_aft_remove',0)
if (.not.params%normaladvect) then
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if
if (flag0==0) then
nedge=surface%nsurface+surface%nt-1
call DoRuRe_surf_stats (params%doDoRuRe,istep,50,is,surface%nt,surface%nsurface,nedge,-nrem)
end if
if (iproc.eq.0) write (8,*) 'Removal after move complete'
endif
end subroutine move_surface
src/refine_surface.f90
+ 90
− 57
View file @ 17a6410f
......@@ -70,12 +70,12 @@ integer ref_count
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer naddp,nremp
integer naddp
integer nnmax
type (edge), dimension(:), allocatable :: ed
type (edge), dimension(:), allocatable :: edswap
integer,dimension(:) ,allocatable :: nedgepernode
integer,dimension(:) ,allocatable :: refine_list,remove_list
integer,dimension(:) ,allocatable :: refine_list,remove_list,remlistswap
integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
logical,dimension(:),allocatable::refine,remove,removep,blocked
double precision,dimension(:),allocatable::edlen,edlenp
......@@ -83,18 +83,19 @@ character*72 :: shift
integer ie,j,jedge,kp,kpp,k,inode,inodep,inodepp,iadd,irem
integer tr1,tr2
integer ntriangle
integer nedge
integer nedge,nside
integer nedgen,nsurfacen
double precision dist0,distmax,distmin,eps
double precision maxdist1,maxdist2
double precision, external :: dist
integer err,ierr,iproc,nproc,i,iedge,i1,i2
integer err,ierr,iproc,nproc,i,iedge,i0,i1,i2
integer nnodeint,nelemint,nedgeint,nelemmax,nnodemax,nedgemax
double precision,dimension(:,:),allocatable::memswap
integer,dimension(:,:),allocatable::iconswap
double precision xm,ym,zm,xyzn
logical :: rem
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
......@@ -108,13 +109,16 @@ eps=1.d-8
!------------------------
! build edge information
!------------------------
nsurfacen=surface%nsurface
ntriangle=surface%nt
nedge=nsurfacen+ntriangle-1
if (surface%closed.eq.1) nedge=nedge-1
nnmax=12 ! nnmax is maximum number of neighbours in triangulation (should be around 6 on average)
!the first nside nodes lie on the convex hull and must not be removed
nside=4*(2**surface%levelt+1)
allocate (ed(nedge),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'ed','refine_surface',size(ed),'int',+1)
allocate (nedgepernode(nsurfacen),stat=threadinfo%err)
......@@ -152,12 +156,12 @@ do ie=1,ntriangle !loop over all triangles
!check edge number and connectivity
iedge=iedge+1
if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2')
nedgepernode(inode)=nedgepernode(inode)+1
if (nedgepernode(inode).gt.nnmax) then
if (iproc.eq.0) write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
if (iproc.eq.0) write (8,'(a,3F10.6)') 'node position: ',surface%x(inode),surface%y(inode),surface%z(inode)
call stop_run ('too many neighbours$')
if (iproc.eq.0) write (*,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
call stop_run ('too many neighbours')
endif
!get edge information of left hand triangle and start end end points
......@@ -171,13 +175,12 @@ do ie=1,ntriangle !loop over all triangles
enddo
enddo
if (iedge/=nedge) call stop_run ('pb_b in refine_surface$')
if (iedge/=nedge) call stop_run ('pb_b in refine_surface')
!------------------------------------
! Calculate length for all edges
!------------------------------------
allocate (edlen(nedge),stat=err)
if (err/=0) call stop_run ('error allocating edlen in refine_surface')
allocate (edlenp(nedge),stat=err)
......@@ -251,6 +254,7 @@ if (params%remove_surf_pts) then
removep(iedge)=(dist0.gt.eps .and. dist0.le.distmin)
end do
call mpi_allreduce (removep,remove,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
deallocate (removep)
endif
......@@ -268,11 +272,12 @@ naddp=0
if (params%remove_surf_pts) then
nrem=count(remove)
allocate(remove_list(nrem))
allocate (remove_list(nrem))
remove_list=0
irem=0
nremp=0
else
nrem=0
irem=0
endif
do iedge=1,nedge
......@@ -283,32 +288,55 @@ do iedge=1,nedge
endif
if (params%remove_surf_pts) then
if (remove(iedge)) then
irem=irem+1
if (ed(iedge)%t1.eq.0 .or. ed(iedge)%t2.eq.0) nremp=nremp+1
!remove_list(irem)=iedge
remove_list(irem)=ed(iedge)%n1
if (ed(iedge)%n1 > nside) then
irem=irem+1
remove_list(irem) = ed(iedge)%n1
elseif (ed(iedge)%n2 > nside) then
irem=irem+1
remove_list(irem) = ed(iedge)%n2
end if
endif
endif
enddo
! Modify remove_list to account for removal of nodes from surface
if (params%remove_surf_pts) then
if (nrem.gt.0) then
do i=1,nrem
do j=i,nrem
if (remove_list(i).lt.remove_list(j)) remove_list(j)=remove_list(j)-1
enddo
enddo
endif
if (irem/=nrem) call stop_run ('pb_b in refine_surface')
deallocate (remove)
deallocate (removep)
endif
!if nodes part of the convex hull would have been removed,
!give warning and adjust remove_list to the correct length
if (irem/=nrem) then
if (iproc==0) then
write (*,'(a,i,a,i,a)') 'WARNING: ',nrem-irem, 'convex hull nodes would have been removed'
write (8,'(a,i,a,i,a)') 'WARNING: ',nrem-irem, 'convex hull nodes would have been removed'
end if
if (irem > 0) then
!adjust length of remove_list
allocate (remlistswap(irem))
remlistswap(:) = remove_list(1:irem)
deallocate (remove_list)
allocate (remove_list(irem))
remove_list = remlistswap
deallocate (remlistswap)
else
!if irem=0, allocate a list of length 0
deallocate (remove_list)
allocate (remove_list(0))
end if
nrem = irem
end if
!sort remove_list for faster checking in remove_point
if (nrem > 0) call iqsort_s (remove_list, nrem)
end if
if (iadd/=nadd) call stop_run ('pb_b in refine_surface')
deallocate (refine)
deallocate (blocked)
if (params%remove_surf_pts) deallocate (remove)
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i2,a,i8,a)') shift//'S.',is,':', nadd,' added ptcls in refine_surface'
......@@ -322,10 +350,10 @@ if (surface%closed.eq.1) nedgen=nedgen-1
! Not sure this is correct...
!nedgen=(surface%nsurface+(nadd-nrem))+(surface%nt+(((nadd-naddp)*2+naddp)-((nrem-nremp)*2+nremp))-1
allocate (edswap(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main$')
allocate (edswap(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main')
edswap(1:nedge)=ed(1:nedge)
deallocate (ed)
allocate (ed(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main$')
allocate (ed(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main')
ed(1:nedge)=edswap(1:nedge)
deallocate (edswap)
......@@ -338,10 +366,10 @@ nelemmax=surface%nt+(nadd-naddp)*2+naddp
nedgemax=nnodemax+nelemmax-1
if (surface%closed.eq.1) nedgemax=nedgemax-1
if (nedgemax.ne.nedgen) call stop_run ('error in counting number of edges in refine_surface$')
if (nedgemax.ne.nedgen) call stop_run ('error in counting number of edges in refine_surface')
allocate (memswap(surface%nsurface,8),stat=err) ; if (err.ne.0) call stop_run('Error alloc memswap in refine_surface$')
allocate (iconswap(3,surface%nt),stat=err) ; if (err.ne.0) call stop_run('Error alloc iconswap in refine_surface$')
allocate (memswap(surface%nsurface,8),stat=err) ; if (err.ne.0) call stop_run('Error alloc memswap in refine_surface')
allocate (iconswap(3,surface%nt),stat=err) ; if (err.ne.0) call stop_run('Error alloc iconswap in refine_surface')
memswap(:,1)=surface%x
memswap(:,2)=surface%y
......@@ -363,16 +391,16 @@ deallocate (surface%r)
deallocate (surface%s)
deallocate (surface%icon)
allocate (surface%x(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in refine_surface$')
allocate (surface%y(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in refine_surface$')
allocate (surface%z(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in refine_surface$')
allocate (surface%xn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xn in refine_surface$')
allocate (surface%yn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc yn in refine_surface$')
allocate (surface%zn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zn in refine_surface$')
allocate (surface%r(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc r in refine_surface$')
allocate (surface%s(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc s in refine_surface$')
allocate (surface%icon(3,nelemmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc icon in refine_surface$')
allocate (surface%x(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in refine_surface')
allocate (surface%y(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in refine_surface')
allocate (surface%z(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in refine_surface')
allocate (surface%xn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xn in refine_surface')
allocate (surface%yn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc yn in refine_surface')
allocate (surface%zn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zn in refine_surface')
allocate (surface%r(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc r in refine_surface')
allocate (surface%s(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc s in refine_surface')
allocate (surface%icon(3,nelemmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc icon in refine_surface')
surface%x(1:surface%nsurface)=memswap(:,1)
surface%y(1:surface%nsurface)=memswap(:,2)
surface%z(1:surface%nsurface)=memswap(:,3)
......@@ -389,7 +417,6 @@ deallocate (iconswap)
!---------------------
! updates the surface
!---------------------
nnodeint=surface%nsurface
nelemint=surface%nt
nedgeint=nedge
......@@ -414,9 +441,6 @@ do i=1,nadd
call update_icon (surface%icon,nelemmax,ed,nedgemax,iedge,nnodeint,nelemint,nedgeint)
enddo
if (iproc.eq.0 .and. params%debug>=1 .and. params%remove_surf_pts) &
write(*,'(a,i2,a,i8,a)') shift//'S.',is,':', nrem,' removed ptcls in refine_surface'
do i=1,size(ed%n1)
i1=ed(i)%n1
i2=ed(i)%n2
......@@ -433,20 +457,29 @@ enddo
surface%nsurface=nnodemax
surface%nt=nelemmax
if (params%remove_surf_pts) then
do i=1,nrem
i1=remove_list(i)
call remove_point(i1,surface,0,surface0)
end do
if (params%remove_surf_pts .and. irem>0) then
do i=1,irem
i1=remove_list(i)
call remove_point(i1,surface,0,surface0,remove_list,nrem)
!Modify remove_list to accout for removal of node from surface
remove_list(i) = 0
do j=i+1,nrem
if (remove_list(j) > i1) remove_list(j) = remove_list(j)-1
end do
end do
if (iproc.eq.0 .and. params%debug>=1) &
write(*,'(a,i2,a,i8,a)') shift//'S.',is,':', nrem,' removed ptcls in refine_surface'
endif
if (.not.params%normaladvect) then
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if
if (nnodeint.ne.nnodemax) call stop_run ('Error counting nodes in refine_surface$')
if (nelemint.ne.nelemmax) call stop_run ('Error counting elements in refine_surface$')
if (nnodeint.ne.nnodemax) call stop_run ('Error counting nodes in refine_surface')
if (nelemint.ne.nelemmax) call stop_run ('Error counting elements in refine_surface')
if (nrem.gt.0) then
nnodemax=surface%nsurface
......@@ -459,7 +492,7 @@ surface%nt=nelemmax
nedge=nedgemax
deallocate(refine_list)
if (params%remove_surf_pts) deallocate(remove_list)
if (params%remove_surf_pts) deallocate (remove_list)
if (params%debug.gt.1) call heap (threadinfo,'ed','refine_surface',size(ed),'int',-1)
deallocate (ed)
......
src/remove_point.f90
+ 26
− 15
View file @ 17a6410f
......@@ -16,7 +16,7 @@
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine remove_point (iii,surface,flag0,surface0)
subroutine remove_point (iii,surface,flag0,surface0,remove_list,nrem)
use definitions
!use mpi
......@@ -38,14 +38,15 @@ include 'mpif.h'
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
integer iii,flag0
integer iii,flag0,nrem
type(sheet) surface,surface0
integer,dimension(nrem) :: remove_list
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
logical already_in_list,iii_found,jjj_found
logical already_in_list,iii_found,jjj_found,rem
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
......@@ -131,18 +132,28 @@ if(counter/=nb_test) call stop_run('pb in test')
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
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
rem = .false.
!check remove_list if j will be removed in the future
do k=1,nrem
if (remove_list(k) >= j) then
if (remove_list(k) == j) rem = .true.
exit
end if
end do
if (.not. rem) then
jjj=j
dist=distt
end if
end if
end do
if (iproc.eq.0) write(8,*) 'closest point to iii is jjj=',jjj
......
src/retriangulate_surface.f90 deleted 100644 → 0
+ 0
− 151
View file @ a2aa187e
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_SURFACE Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine retriangulate_surface (surface,is,istep,params)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! If move_cloud in combination with regular remove_points failed (one
! node has too many neighbors), then this routine removes all points
! lying outside of the domain and recalculates the triangulation.
! WARNING: This does not work for closed surfaces
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
!use constants
use definitions
!use mpi
implicit none
include 'mpif.h'
type (sheet) :: surface
type (parameters) :: params
integer :: is, istep
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: nproc,iproc,ierr
integer :: err,ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,nh,nohalt_hull,loc
integer :: dmode,nt
integer :: nnn(1),nnlist(1),ntrilist(1)
integer,dimension(:), allocatable :: hulltriangles,vis_tlist,vis_elist
integer,dimension(:), allocatable :: add_tlist
integer,dimension(:,:),allocatable :: vertices,neighbour
double precision :: eps
double precision,dimension(:), allocatable :: field
double precision,dimension(:,:),allocatable :: points,centres
logical clockwise
logical*1,dimension(:),allocatable :: lt_work,ln_work
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
!if (params%debug>=2) call output_surf (surface,is,'faulty',istep,-1)
if (params%debug>=1) call output_surf (surface,is,'faulty',istep,-1)
!-------------------------------------------|
!-----building icon array-------------------|
!-------------------------------------------|
! we use the NN library to compute the triangulation
! of points for the surface
!
! Code taken from create_surfaces.f90
if (iproc==0) write (8,'(a)') 'Starting retriangulation'
ntmax=surface%nsurface*3
nhmax=surface%nsurface
npmax=surface%nsurface
nnpnmax=100
nmax=3*ntmax+npmax
nvmax=ntmax
allocate (points(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc points in retriangulate_surface')
allocate (field(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc field in retriangulate_surface')
allocate (vertices(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vertices in retriangulate_surface')
allocate (centres(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc centress in retriangulate_surface')
allocate (neighbour(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc neighbour in retriangulate_surface')
allocate (hulltriangles(nhmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in retriangulate_surface')
allocate (vis_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in retriangulate_surface')
allocate (vis_elist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_elist in retriangulate_surface')
allocate (add_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc add_tlist in iretriangulate_surface')
allocate (lt_work(ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lt_work in retriangulate_surface')
allocate (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in retriangulate_surface')
points(1,:)=surface%x
points(2,:)=surface%y
dmode=-2
nmode=0
loc=1
eps=1.0d-8
nohalt_hull=0
clockwise=.true.
call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
points,dmode,nmode,clockwise,field,nt,vertices, &
centres,neighbour,nh,hulltriangles,nohalt_hull, &
loc,nnn,nnlist,ntrilist, &
eps,nvmax,vis_tlist,vis_elist,add_tlist, &
lt_work,ln_work)
if (iproc.eq.0) then
write (*,'(a,i2,a,i8,a)') 'retriangulated surface ',is,' has ',nt,' triangles'
write (8,'(a,i2,a,i8,a)') 'retriangulated surface ',is,' has ',nt,' triangles'
end if
if (surface%nt /= nt) then
deallocate (surface%icon)
allocate (surface%icon(3,nt),stat=err); if (err/=0) call stop_run ('Error alloc icon in retriangulate_surface')
surface%nt = nt
end if
surface%icon(1:3,1:nt)=vertices(1:3,1:nt)
deallocate (points)
deallocate (field)
deallocate (vertices)
deallocate (centres)
deallocate (neighbour)
deallocate (hulltriangles)
deallocate (vis_tlist)
deallocate (vis_elist)
deallocate (add_tlist)
deallocate (lt_work)
deallocate (ln_work)
!----------------------------------------------|
!----computing the normals---------------------|
!----------------------------------------------|
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
!if (params%debug>=2) call output_surf (surface,is,'init',-1,-1)
if (params%debug>=1) call output_surf (surface,is,'redone',istep,-1)
end subroutine retriangulate_surface
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment