Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CALCULATE_LSF Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine calculate_lsf (lsf,octree_lsf,noctree_lsf,icon_lsf,nleaves_lsf, &
na_lsf,xl,yl,zl,xln,yln,zln,icon,nelem,na, &
levelmax_oct,normaladvect)
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
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this subroutine calculates the lsf (level set function) on an octree from
! the position of a series of particles connected by a 2D triangulation
! the algorithm goes first through the triangles
! to calculate a 3D normal at the "centre" of each triangle (from the position of the
! corresponding particles). From that normal, the levelset function is
! estimated at the nodes of the leaf of the levelset octree containing the
! centre of the particles (this algorithm assumes that the density of the
! particles is larger than the density of the octree)
! in a second pass, the algorithm propagates the sign of the level set function
! by successive pass through the leaves of the octree. The algorithm converges
! very rapidly because each successive path is in reverse order from the previous
! lsf is the values of the lsf at the nodes (main output of the routine)
! octree_lsf is the octree on which the lsf is to be built known
! noctree_lsf is its size
! icon_lsf is the node-leaf connectivity of the octree
! nleaves_lsf is the number of leaves in the octree
! na_lsf is the number of nodes in the octree
! xl,yl,zl are the coordinates of the particles on the surface
! icon is the connectivity matrix defining the triangulation
! nelem is the number of triangles
! na is the number of particles defining the surface
! levelmax_oct is the maximum level of refinement of the octree
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
!use mpi
Dave Whipp
committed
include 'mpif.h'
double precision lsf(na_lsf)
integer octree_lsf(noctree_lsf)
integer icon_lsf(8,nleaves_lsf)
integer nleaves_lsf
integer na_lsf
double precision xl(na),yl(na),zl(na)
double precision xln(na),yln(na),zln(na)
integer icon(3,nelem)
integer nelem
integer na
integer levelmax_oct
logical normaladvect
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ierr,iproc,nproc,ii,jj,kk,iii,istart,ijk,iend,itest,i,j,k,leaf_lsf
integer i1,i2,j1,j2,k1,k2,it,inode,niter,iinc,num_lsf,level,loc,nleaf_test,err
integer,dimension(:),allocatable::leaf_test
integer,dimension(:),allocatable::done_el,done_nd,norm
double precision x0,y0,z0,dxyz,sum_lsf,test
double precision x0_loc,y0_loc,z0_loc
double precision x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3
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
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
double precision,dimension(:),allocatable::distmin
double precision,dimension(:),allocatable::xxln,yyln,zzln
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
! first pass using the centres of the triangles connecting
! the particles of the surface
! (note that we also take the opportunity to update the normal of
! the particles)
lsf=0.d0
allocate (done_el(nleaves_lsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc done_el in calculate_lsf$')
allocate (done_nd(na_lsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc done_nd in calculate_lsf$')
allocate (distmin(na_lsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc distmin in calculate_lsf$')
allocate (norm(na),stat=err) ; if (err.ne.0) call stop_run ('Error alloc norm in calculate_lsf$')
allocate (xxln(na),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xxln in calculate_lsf$')
allocate (yyln(na),stat=err) ; if (err.ne.0) call stop_run ('Error alloc yyln in calculate_lsf$')
allocate (zzln(na),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zzln in calculate_lsf$')
done_el=1
done_nd=1
distmin=1.d6
! NOTE by JEAN
! here we should not be recalculating the normals but using those that have been advected
! But I have checked for the ball that the normals are OK during the first step for the ball
if (normaladvect) then
xxln=xln
yyln=yln
zzln=zln
else
call compute_normals (na,xl,yl,zl,nelem,icon,xxln,yyln,zzln)
endif
! estimates signed distance to surface
do it=1,nelem
x1=xl(icon(1,it))
x2=xl(icon(2,it))
x3=xl(icon(3,it))
y1=yl(icon(1,it))
y2=yl(icon(2,it))
y3=yl(icon(3,it))
z1=zl(icon(1,it))
z2=zl(icon(2,it))
z3=zl(icon(3,it))
x1n=xxln(icon(1,it))
x2n=xxln(icon(2,it))
x3n=xxln(icon(3,it))
y1n=yyln(icon(1,it))
y2n=yyln(icon(2,it))
y3n=yyln(icon(3,it))
z1n=zzln(icon(1,it))
z2n=zzln(icon(2,it))
z3n=zzln(icon(3,it))
i1=int(minval(xl(icon(:,it)))/.5d0**levelmax_oct)
i2=int(maxval(xl(icon(:,it)))/.5d0**levelmax_oct)+1
j1=int(minval(yl(icon(:,it)))/.5d0**levelmax_oct)
j2=int(maxval(yl(icon(:,it)))/.5d0**levelmax_oct)+1
k1=int(minval(zl(icon(:,it)))/.5d0**levelmax_oct)
k2=int(maxval(zl(icon(:,it)))/.5d0**levelmax_oct)+1
allocate (leaf_test((k2-k1+1)*(j2-j1+1)*(i2-i1+1)),stat=err)
if (err.ne.0) call stop_run ('Error alloc leaf_test in calculate_lsf$')
nleaf_test=0
! do k=k1,k2-1
! z=(real(k)+.5d0)*(.5d0**levelmax_oct)
! do j=j1,j2-1
! y=(real(j)+.5d0)*(.5d0**levelmax_oct)
! do i=i1,i2-1
! x=(real(i)+.5d0)*(.5d0**levelmax_oct)
do k=k1,k2
z=dfloat(k)*(.5d0**levelmax_oct)-.25d0**levelmax_oct
do j=j1,j2
y=dfloat(j)*(.5d0**levelmax_oct)-.25d0**levelmax_oct
do i=i1,i2
x=dfloat(i)*(.5d0**levelmax_oct)-.25d0**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) then
! find octree leaf
call octree_find_leaf (octree_lsf,noctree_lsf,x,y,z,leaf_lsf,level,loc,x0,y0,z0,dxyz)
do itest=1,nleaf_test
if (leaf_lsf.eq.leaf_test(itest)) goto 1111
enddo
nleaf_test=nleaf_test+1
leaf_test(nleaf_test)=leaf_lsf
! compute potential lsf's
ijk=0
do kk=0,1
do jj=0,1
do ii=0,1
ijk=ijk+1
x0_loc=x0+ii*dxyz
y0_loc=y0+jj*dxyz
z0_loc=z0+kk*dxyz
inode=icon_lsf(ijk,leaf_lsf)
call distance_to_triangle &
(x1,y1,z1,x2,y2,z2,x3,y3,z3, &
x1n,y1n,z1n,x2n,y2n,z2n,x3n,y3n,z3n, &
x0_loc,y0_loc,z0_loc,test)
lsf8(ijk)=test
enddo
enddo
enddo
do ijk=1,8
inode=icon_lsf(ijk,leaf_lsf)
if (abs(lsf8(ijk)).lt.distmin(inode)) then
lsf(inode)=lsf8(ijk)
distmin(inode)=abs(lsf8(ijk))
done_nd(inode)=0
endif
enddo
1111 continue
endif
enddo
enddo
enddo
deallocate (leaf_test)
enddo
! now propagates the lsf into elements that are not cut by the surface
istart=1
iend=nleaves_lsf
iinc=1
niter=0
do while (sum(done_nd).ne.0)
do i=istart,iend,iinc
if (done_el(i).eq.1) then
do k=1,8
if (done_nd(icon_lsf(k,i)).eq.0) then
done_el(i)=0
! note that we average the value of the lsf around the nodes of the leaf
! to find the sign of the nodes that have not yet been assigned
! the reason for this is that there might be leaves where the surface is
! passing through but where no particle resides; these leaves would have
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
! lsf values of different sign on different nodes from the previous two passes
! this should never happen...
sum_lsf=0.d0
num_lsf=0
do kk=1,8
if (done_nd(icon_lsf(kk,i)).eq.0) then
num_lsf=num_lsf+1
sum_lsf=sum_lsf+lsf(icon_lsf(kk,i))
endif
enddo
sum_lsf=sum_lsf/num_lsf
do kk=1,8
if (done_nd(icon_lsf(kk,i)).eq.1) then
lsf(icon_lsf(kk,i))=sign(1.d0,sum_lsf)
done_nd(icon_lsf(kk,i))=0
endif
enddo
goto 1
endif
enddo
endif
1 continue
enddo
iii=istart
istart=iend
iend=iii
iinc=-iinc
niter=niter+1
if (niter.eq.100) then
print*,'Too many iterations'
print*,'in algoritm propagating lsf function'
print*,'number of nodes remaining in lsf build ',sum(done_nd)
call stop_run ('program terminated in calculate_lsf$')
endif
enddo
deallocate (done_el)
deallocate (done_nd)
deallocate (distmin)
deallocate (norm)
deallocate (xxln)
deallocate (yyln)
deallocate (zzln)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine distance_to_triangle (x1,y1,z1,x2,y2,z2,x3,y3,z3, &
x1n,y1n,z1n,x2n,y2n,z2n,x3n,y3n,z3n, &
x0,y0,z0,dist)
implicit none
double precision x0,x1,x2,x3
double precision y0,y1,y2,y3
double precision z0,z1,z2,z3
double precision x1n,y1n,z1n
double precision x2n,y2n,z2n
double precision x3n,y3n,z3n
double precision dist
double precision x(4),y(4),z(4)
double precision xxn(4),yyn(4),zzn(4)
double precision xn,yn,zn,xyzn
double precision xc,yc,zc
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
double precision xnp(3),ynp(3),znp(3),alphap(3),dd(3)
double precision xp,yp,zp,xpp,ypp,zpp
double precision xxxn,yyyn,zzzn
integer k
x(1)=x1;x(2)=x2;x(3)=x3;x(4)=x1
y(1)=y1;y(2)=y2;y(3)=y3;y(4)=y1
z(1)=z1;z(2)=z2;z(3)=z3;z(4)=z1
xxn(1)=x1n;xxn(2)=x2n;xxn(3)=x3n;xxn(4)=x1n
yyn(1)=y1n;yyn(2)=y2n;yyn(3)=y3n;yyn(4)=y1n
zzn(1)=z1n;zzn(2)=z2n;zzn(3)=z3n;zzn(4)=z1n
xc=(x1+x2+x3)/3.d0
yc=(y1+y2+y3)/3.d0
zc=(z1+z2+z3)/3.d0
xn=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
yn=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
zn=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
xyzn=sqrt(xn**2+yn**2+zn**2)
xn=xn/xyzn
yn=yn/xyzn
zn=zn/xyzn
alpha=xn*(x0-xc)+yn*(y0-yc)+zn*(z0-zc)
do k=1,3
xnp(k)=(y(k+1)-y(k))*zn-(z(k+1)-z(k))*yn
ynp(k)=(z(k+1)-z(k))*xn-(x(k+1)-x(k))*zn
znp(k)=(x(k+1)-x(k))*yn-(y(k+1)-y(k))*xn
xyzn=sqrt(xnp(k)**2+ynp(k)**2+znp(k)**2)
xnp(k)=xnp(k)/xyzn
ynp(k)=ynp(k)/xyzn
znp(k)=znp(k)/xyzn
alphap(k)=(x0-x(k))*xnp(k)+(y0-y(k))*ynp(k)+(z0-z(k))*znp(k)
enddo
if (alphap(1).le.0.d0 .and. alphap(2).le.0.d0 .and. alphap(3).le.0.d0) then
dist=alpha
return
endif
do k=1,3
dd(k)=1.d6
if (alphap(k).ge.0.d0) then
xp=x0-alpha*xn
yp=y0-alpha*yn
zp=z0-alpha*zn
xpp=xp-alphap(k)*xnp(k)
ypp=yp-alphap(k)*ynp(k)
zpp=zp-alphap(k)*znp(k)
beta=((xpp-x(k))*(x(k+1)-x(k))+(ypp-y(k))*(y(k+1)-y(k))+(zpp-z(k))*(z(k+1)-z(k)))/ &
((x(k+1)-x(k))*(x(k+1)-x(k))+(y(k+1)-y(k))*(y(k+1)-y(k))+(z(k+1)-z(k))*(z(k+1)-z(k)))
if (beta.lt.0.d0) then
dd(k)=sqrt((x0-x(k))**2+(y0-y(k))**2+(z0-z(k))**2)
alpha=xxn(k)*(x0-x(k))+yyn(k)*(y0-y(k))+zzn(k)*(z0-z(k))
dd(k)=sign(dd(k),alpha)
elseif (beta.gt.1.d0) then
dd(k)=sqrt((x0-x(k+1))**2+(y0-y(k+1))**2+(z0-z(k+1))**2)
alpha=xxn(k+1)*(x0-x(k+1))+yyn(k+1)*(y0-y(k+1))+zzn(k+1)*(z0-z(k+1))
dd(k)=sign(dd(k),alpha)
else
xxxn=xxn(k)+beta*(xxn(k+1)-xxn(k))
yyyn=yyn(k)+beta*(yyn(k+1)-yyn(k))
zzzn=zzn(k)+beta*(zzn(k+1)-zzn(k))
dd(k)=sqrt((x0-xpp)**2+(y0-ypp)**2+(z0-zpp)**2)
alpha=xxxn*(x0-xpp)+yyyn*(y0-ypp)+zzzn*(z0-zpp)
dd(k)=sign(dd(k),alpha)
endif
endif
enddo
dist=dd(1)
if (abs(dd(2)).lt.abs(dist)) dist=dd(2)
if (abs(dd(3)).lt.abs(dist)) dist=dd(3)
return
end