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
48
49
50
51
52
53
54
55
56
57
c find_neighbours
subroutine find_neighbours (x,y,nn,nb,nnode,nbmax,
& points,vertices,neighbour,nodes,
& vis_tlist,vis_elist,add_tlist,nt,
& surface,newsurface,eps,
& xx,pp,aa,bb,surfscale,sides)
c subroutine to find the list of natural neighbours attached to each
c nodes. It computes the nn and nb lists but also the surfaces
c attached to the nodes.
c INPUT: x,y = x,y nodal coordinates
c nnode = number of nodes
c nbmax = maximum number of neighbour per node
c points = working array
c nodes = working array
c vis_tlist = working array
c vis_elist = working array
c add_tlist = working array
c newsurface= here it is just used as a working array
c eps = precision
c xx = working array
c pp = working array
c aa = working array
c bb = working array
c surfscale = average nodal surface
c OUTPUT: nn = neighbour array
c nb = number of neighbour for each node
c vertices = triangle list
c neighbour = neighbour list
c nt = number of triangle
c surface = nodal surface (voronoi cell surface area)
c sides = length of the sides of the voronoi cells
c subroutines called:
c NONE
common /vocal/ ivocal
real x(*),y(*),surface(*)
integer nn(nbmax,*)
integer nb(*)
real sides(nbmax,*)
real*8 points(2,*),eps
integer vertices(3,*)
integer neighbour(3,*)
integer nodes(*)
integer vis_tlist(*)
integer vis_elist(*)
integer add_tlist(*)
integer ccw
real xx(2),pp(2,nbmax),aa(nbmax,2),bb(nbmax)
real newsurface(*)
logical dummy1(2)
integer dummy2(2),bfirst
do i=1,nnode
points(1,i)=dble(x(i))
points(2,i)=dble(y(i))
nodes(i)=i
enddo
np=nnode
nv_max=nnode
c sorting by x
if (ivocal.eq.1) call debug ('indexx$',0)
if (ivocal.eq.1) call debug ('find_neighbours$',1)
c ensure initial triangle is in ccw order
if(ccw(points(1,nodes(1)),
& points(1,nodes(2)),
& points(1,nodes(3)),k).eq.-1)then
itemp = nodes(1)
nodes(1) = nodes(2)
nodes(2) = itemp
end if
c Call the routine that
c does the work
mode=0
eps=0.d0
if (ivocal.eq.1) call debug ('delaun$',0)
call delaun (points,np,neighbour,vertices,nt,2*np,
& vis_tlist,vis_elist,add_tlist,eps,nv_max,
& mode,dummy1,bfirst,itstart,dummy2)
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
if (ivocal.eq.1) call debug ('find_neighbours$',1)
nt_max=nnode*3
if(nt.gt.nt_max)then
write(6,*)' Error number of triangles calculated is larger'
write(6,*)' than array sizes.'
write(6,*)' Number of triangles = ',nt
write(6,*)' Maximum value (nt_max) = ',nt_max
stop
endif
c find neighbour list
call find_neighbour_list (nb,nn,nnode,nbmax,x,y,sides,
& nt,vertices,neighbour)
c find Voronoi cell surface area
do i=1,nnode
newsurface(i)=1.
enddo
if (ivocal.eq.1) call debug ('find_surface$',0)
call find_surface (nn,nb,surface,nbmax,nnode,
& x,y,xx,pp,aa,bb,newsurface,surfscale)
if (ivocal.eq.1) call debug ('find_neighbours$',1)
do i=1,nnode
nb(i)=nb(i)+1
if (nb(i).gt.nbmax) then
stop 'nbmax too small...'
endif
nn(nb(i),i)=i
enddo
return
end