!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! EROSION Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine erosion (surface,olsf,is,zerosion)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine will need major improvement in the near future
! it is used to erode the surface "surface" according to a user supplied algorithm
! at the moment it simply "shaves" the top of the model at a given height
! surface are the surface/sheet object
! (the top free surface)
! olsf is a octreelsf object containig the geometry of the current velocity octree/object
! is is the number of the surface (0 playing a special role)
! zerosion is the level of erosion (between 0 and 1)
! Note that all surfaces are required because they can be affected by the erosion process
! for the top surface, we apply the erosion and on the octree
! we define an lsf that is simply the height of the eroded surface
! This allows us, for the other surfaces to find the height of the
! free/eroded surface at any point in space and check whether the
! vertical position of the point is above or below the free surface
! We then simply "erode" the following surfaces using this informatoin
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
type (sheet) surface
type (octreelsf) olsf
integer is
double precision zerosion
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ierr,iproc,nproc,i,nt,nh,nohalt_hull,loc,nnn,nnlist,ntrilist
integer ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,int_method
double precision,dimension(:),allocatable::field
double precision,dimension(:,:),allocatable::points,centres
integer,dimension(:,:),allocatable::vertices,neighbour
integer,dimension(:),allocatable::hulltriangles,lt_work,ln_work
integer,dimension(:),allocatable::vis_tlist,vis_elist,add_tlist
logical clockwise,out,nnext
double precision eps,zref
integer dmode,err
double precision x,y,f,df(2),ddf(3),v
double precision,dimension(:,:),allocatable::work_d1,work_d2,work_d3
double precision,dimension(:,:),allocatable::work_d4,work_d5,work_d6,work_d7
integer,dimension(:),allocatable::work_i1,work_i3
integer,dimension(:,:),allocatable::work_i2
real,dimension(:),allocatable::work_r2
real,dimension(:,:),allocatable::work_r1
double precision,dimension(:),allocatable::lsf
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
if (is.ne.1) then
do i=1,surface%nsurface
call octree_interpolate (olsf%octree,olsf%noctree,olsf%icon, &
olsf%nleaves,olsf%lsf,olsf%nnode, &
surface%x(i),surface%y(i),surface%z(i),f)
! call octree_interpolate3 (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf, &
! olsf%x(1:olsf%nnode),olsf%y(1:olsf%nnode),olsf%z(1:olsf%nnode),olsf%nnode, &
! surface%x(i),surface%y(i),surface%z(i),f)
surface%z(i)=min(f,surface%z(i))
enddo
return
endif
! first apply erosion to top surface
zref=zerosion
do i=1,surface%nsurface
surface%z(i)=zref
enddo
! nn_setup (calculates delaunay triangles and other info)
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 erode$')
allocate (field(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc field in erode$')
allocate (vertices(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vertices in erode$')
allocate (centres(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc centress in erode$')
allocate (neighbour(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc neighbour in erode$')
allocate (hulltriangles(nhmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in erode$')
allocate (vis_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in erode$')
allocate (vis_elist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_elist in erode$')
allocate (add_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc add_tlist in erode$')
allocate (lt_work(ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lt_work in erode$')
allocate (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in erode$')
points(1,:)=surface%x
points(2,:)=surface%y
field(:)=surface%z
dmode=-2
nmode=0
clockwise=.true.
eps=1.d-20
eps=0.d0
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)
allocate (work_d1(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d1 in erode$')
allocate (work_d2(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d2 in erode$')
allocate (work_d3(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d3 in erode$')
allocate (work_d4(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d4 in erode$')
allocate (work_d5(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d5 in erode$')
allocate (work_d6(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d6 in erode$')
allocate (work_d7(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d7 in erode$')
allocate (work_i1(nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_i1 in erode$')
allocate (work_i2(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_i2 in erode$')
allocate (work_i3(nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_i3 in erode$')
allocate (work_r1(nnpnmax,2),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_r1 in erode$')
allocate (work_r2(nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_r2 in erode$')
allocate (lsf(olsf%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf in erode$')
nnext=.true.
int_method=3
loc=1
lsf=0.d0
do i=iproc+1,olsf%nnode,nproc
x=olsf%x(i)
y=olsf%y(i)
call nn2D (x,y,points,vertices,neighbour,nnext,int_method,centres, &
hulltriangles,nh,loc,field,nnpnmax,work_r1,work_r2, &
work_d1,work_i1,work_i2,work_i3,work_d2,work_d3,work_d4, &
work_d5,work_d6,work_d7,lt_work,ln_work,out,f,df,ddf,v)
lsf(i)=f
enddo
olsf%lsf(1:olsf%nnode)=0.d0
call mpi_allreduce (lsf,olsf%lsf,olsf%nnode, &
mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
deallocate (points)
deallocate (field)
deallocate (vertices)
deallocate (centres)
deallocate (hulltriangles)
deallocate (vis_tlist)
deallocate (vis_elist)
deallocate (add_tlist)
deallocate (lt_work)
deallocate (ln_work)
deallocate (work_d1)
deallocate (work_d2)
deallocate (work_d3)
deallocate (work_d4)
deallocate (work_d5)
deallocate (work_d6)
deallocate (work_d7)
deallocate (work_r1)
deallocate (work_r2)
deallocate (work_i1)
deallocate (work_i2)
deallocate (work_i3)
deallocate (lsf)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|