Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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
!use mpi
include 'mpif.h'
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|