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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! BUILD_SURFACE_OCTREE Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine build_surface_octree (surface,olsf,params,threadinfo)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this subroutine builds an octree to carry a level set function
! used to represent a surface
! surface is the surface represented by particles
! olsf is the octree that we must build (it has already been dimensioned to
! noctreemax in the main program)
! leveluniform_oct is the base level for the octree
! levelmax_oct is the maximum level allowed for leaves of the octree
! the octree resolution is based on one of a series of criteria
! defined by criterion (read in the input.txt file)
! if criterion=1 all triangles of the surface must be entirely
! comprised in leaf of maximum order
! if criterion=2 the octree is discretized to la level computed
! from the angle made by the normals; when the angle varies between
! 0 and angle max, the level varies from levelmin to levelmax
! ismooth determines whether additional smoothing is to be performed
! (ismooth is set in the input.txt file). Additional smoothing imposes
! that no leaf touches other leaves that are different in size by more than
! one level
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
Dave Whipp
committed
!use mpi
Dave Whipp
committed
use threads
include 'mpif.h'
type (sheet) surface
type (octreelsf) olsf
type (parameters) params
type (thread) threadinfo
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer level,iproc,nproc,ierr,i,j,k,seed
integer i1,i2,j1,j2,k1,k2,it,kp,ic,icp
double precision x,y,z,angle,anglemax
!double precision :: dist
real cosd
integer ioctree_number_of_elements
external ioctree_number_of_elements
double precision, dimension(:), allocatable ::xn,yn,zn
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
seed=56987
! criterion 1: all cells cut by any surface must be at surface_leveloct
if (surface%criterion.eq.1) then
do it=1,surface%nt
i1=int(minval(surface%x(surface%icon(:,it)))/.5d0**params%levelmax_oct)
i2=int(maxval(surface%x(surface%icon(:,it)))/.5d0**params%levelmax_oct)+1
j1=int(minval(surface%y(surface%icon(:,it)))/.5d0**params%levelmax_oct)
j2=int(maxval(surface%y(surface%icon(:,it)))/.5d0**params%levelmax_oct)+1
k1=int(minval(surface%z(surface%icon(:,it)))/.5d0**params%levelmax_oct)
k2=int(maxval(surface%z(surface%icon(:,it)))/.5d0**params%levelmax_oct)+1
do k=k1,k2-1
z=(real(k)+.5d0)*(.5d0**params%levelmax_oct)
do j=j1,j2-1
y=(real(j)+.5d0)*(.5d0**params%levelmax_oct)
do i=i1,i2-1
x=(real(i)+.5d0)*(.5d0**params%levelmax_oct)
if (x*(x-1.d0).le.0.d0 .and. &
y*(y-1.d0).le.0.d0 .and. &
z*(z-1.d0).le.0.d0) &
! call octree_create_from_particles (olsf%octree,olsf%noctree,x,y,z,1,params%levelmax_oct)
call octree_create_from_particles (olsf%octree,olsf%noctree,x,y,z,1,surface%leveloct)
enddo
enddo
enddo
enddo
! criterion 2: curvature based
elseif (surface%criterion.eq.2) then
allocate (xn(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xn in compute_normals$')
allocate (yn(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc yn in compute_normals$')
allocate (zn(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zn in compute_normals$')
if (params%normaladvect) then
xn(1:surface%nsurface)=surface%xn(1:surface%nsurface)
yn(1:surface%nsurface)=surface%yn(1:surface%nsurface)
zn(1:surface%nsurface)=surface%zn(1:surface%nsurface)
else
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z,surface%nt,surface%icon,xn,yn,zn)
endif
do i=1,surface%nt
level=params%leveluniform_oct
do k=1,3
kp=mod(k,3)+1
ic=surface%icon(k,i)
icp=surface%icon(kp,i)
cosd=xn(ic)*xn(icp)+yn(ic)*yn(icp)+zn(ic)*zn(icp)
! dist=sqrt((surface%x(ic)-surface%x(icp))**2 &
! +(surface%y(ic)-surface%y(icp))**2 &
! +(surface%z(ic)-surface%z(icp))**2)
! angle=acos(cosd)/dist
! anglemax=surface%anglemaxoctree/(0.5d0**surface%leveloct)
angle=acos(cosd)
anglemax=surface%anglemaxoctree
level=max(level,params%leveluniform_oct+ &
nint(min(angle/surface%anglemaxoctree,1.d0)*(params%levelmax_oct-params%leveluniform_oct)))
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
! level=max(level,params%leveluniform_oct+ nint(min(angle/surface%anglemaxoctree,1.d0)*(surface%leveloct-params%leveluniform_oct)))
! level=max(level,params%leveluniform_oct+ nint(min(angle/surface%anglemax,1.d0)*(surface%leveloct-params%leveluniform_oct)))
enddo
do k=1,3
ic=surface%icon(k,i)
x=surface%x(ic)
y=surface%y(ic)
z=surface%z(ic)
call octree_create_from_particles (olsf%octree,olsf%noctree,x,y,z,1,level)
enddo
enddo
deallocate (xn)
deallocate (yn)
deallocate (zn)
! criterion 3: all cells containing a node must be at surface%leveloct
elseif (surface%criterion.eq.3) then
do i=1,surface%nsurface
x=surface%x(i)
y=surface%y(i)
z=surface%z(i)
! call octree_create_from_particles(olsf%octree,olsf%noctree,x,y,z,1,params%levelmax_oct)
call octree_create_from_particles(olsf%octree,olsf%noctree,x,y,z,1,surface%leveloct)
enddo
else
call stop_run ('error in create_surface_octree; wrong criterion$')
endif
if (params%debug.gt.1) then
write(threadinfo%Logunit,*) '[same] build_surface_octree: olsf%octree(2) ',olsf%octree(2)
call flush(threadinfo%Logunit)
endif
call octree_smoothen (olsf%octree,olsf%noctree)
if (params%ismooth) then
if (iproc.eq.0) write (8,*) 'Before super_smooth ',ioctree_number_of_elements (olsf%octree,olsf%noctree)
call octree_super_smoothen (olsf%octree,olsf%noctree)
if (iproc.eq.0) write (8,*) 'After super_smooth ',ioctree_number_of_elements (olsf%octree,olsf%noctree)
endif
olsf%nleaves=ioctree_number_of_elements (olsf%octree,olsf%noctree)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|