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
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! COMPUTE_NORMALS Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine compute_normals (ns,x,y,z,nt,icon,xn,yn,zn)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! given a set of points, and a connectivity array of their triangulation,
! this routine computes the normal to each point through cross-products:
! at first one computes the normal to each triangle and add the vector to the
! three points that make the triangle
! then one loops over the points, and normalises the vectors.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
integer ns
double precision x(ns),y(ns),z(ns)
integer nt
integer icon(3,nt)
double precision xn(ns),yn(ns),zn(ns)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,i1,i2,i3,ij,k
integer iproc,nproc,ierr
integer, dimension(:),allocatable :: nb
double precision x1,x2,x3
double precision y1,y2,y3
double precision z1,z2,z3
double precision xne,yne,zne
double precision xyzn
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
allocate(nb(ns))
nb=0
xn=0.d0
yn=0.d0
zn=0.d0
do i=1,nt
i1=icon(1,i)
i2=icon(2,i)
i3=icon(3,i)
x1=x(i1)
x2=x(i2)
x3=x(i3)
y1=y(i1)
y2=y(i2)
y3=y(i3)
z1=z(i1)
z2=z(i2)
z3=z(i3)
xne=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
yne=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
zne=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
xyzn=sqrt(xne**2+yne**2+zne**2)
xne=xne/xyzn
yne=yne/xyzn
zne=zne/xyzn
do k=1,3
ij=icon(k,i)
xn(ij)=xn(ij)+xne
yn(ij)=yn(ij)+yne
zn(ij)=zn(ij)+zne
nb(ij)=nb(ij)+1
enddo
enddo
do i=1,ns
xyzn=sqrt(xn(i)**2+yn(i)**2+zn(i)**2)
xn(i)=xn(i)/xyzn
yn(i)=yn(i)/xyzn
zn(i)=zn(i)/xyzn
enddo
if (minval(nb).eq.0) then
do i=1,ns
if (nb(i).eq.0.and.iproc.eq.0) then
write (8,*) 'Particle ',i,' at ',x(i),y(i),z(i),' not connected'
end if
enddo
call stop_run ('error in compute_normals_from_triangles$')
else
if (iproc.eq.0) then
write(8,*) 'in compute_normals_from_triangles:'
write(8,*) 'minval(nb)=',minval(nb),' neighbours'
write(8,*) 'maxval(nb)=',maxval(nb),' neighbours'
end if
end if
deallocate(nb)
return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------