Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CALCULATE_OCTREE Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine embed_surface_in_octree (osolve,params,surface,is,istep,iter,threadinfo)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine embed a surface in a new octree, computes a lsf, and then
! refines osolve accordingly and interpolates the lsf on osolve
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
Dave Whipp
committed
use mpi
use threads
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
implicit none
type(octreesolve) osolve
type(parameters) params
type (octreelsf) olsf
type (sheet) surface
integer is
integer istep
integer iter
type (thread) threadinfo
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer,dimension(:),allocatable::levs
type (octreesolve) oint
integer err,ierr,iproc,nproc
integer ie,ilsf,k,i
double precision xxx,yyy,zzz,dist
integer ioctree_number_of_elements
external ioctree_number_of_elements
character*72 :: shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
if(iproc.eq.0) write(8,*) 'embedding surface',is,'/',params%ns
! we want to create the octree olsf in which the surface under consideration will be imbedded.
! we first create a uniform octree at level leveluniform_oct
olsf%noctree=params%noctreemax
allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
call heap (threadinfo,'olsf%octree', 'embed_surf...',size(olsf%octree),'int',+1)
call octree_init (olsf%octree,olsf%noctree)
call octree_create_uniform (olsf%octree,olsf%noctree,params%leveluniform_oct)
! we then improve the octree to best fit/represent the surface
call build_surface_octree (surface,olsf,params,threadinfo)
if (iproc.eq.0) write (8,*) olsf%nleaves,' leaves in lsf octree'
! we then calculate the connectivity and node coordinates of the surface octree
olsf%nnode=olsf%nleaves*3
allocate (olsf%icon(8,olsf%nleaves),stat=threadinfo%err)
call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',+1)
call octree_find_node_connectivity (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%x,olsf%y,olsf%z,olsf%nnode)
Douglas Guptill
committed
! olsf%nnode has been changed by octree_find_node_connectivity, so re-size x,y,z
call octreelsf_shrink_xyz(olsf, threadinfo)
write(threadinfo%Logunit,*) '[same] embed_surface_in_octree: olsf%nnode ',olsf%nnode
call flush(threadinfo%Logunit)
allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%lsf', 'embed_surf...',size(olsf%lsf),'dp',+1)
! we then calculate the lsf on the olsf octree
call calculate_lsf (olsf%lsf,olsf%octree,olsf%noctree, &
olsf%icon,olsf%nleaves,olsf%nnode, &
surface%x,surface%y,surface%z, &
surface%xn,surface%yn,surface%zn, &
surface%icon,surface%nt,surface%nsurface, &
params%levelmax_oct,params%normaladvect,is)
if (params%debug>=2) call output_octree_lsf (olsf,surface,is,istep,iter)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! we create and allocate the oint octree that will receive the current osolve octree.
! this will allow us to resize osolve so that the new surface can be embedded in it.
oint%noctree=osolve%noctree
oint%nleaves=osolve%nleaves
oint%nnode=osolve%nnode
oint%nlsf=osolve%nlsf
! write(threadinfo%mem_heap_unit,*) 'oint noctree,nleaves,nnode,nlsf are', oint%noctree,oint%nleaves,oint%nnode,oint%nlsf
allocate (oint%octree(params%noctreemax),stat=threadinfo%err)
call heap (threadinfo,'oint%octree', 'embed_surf...',size(oint%octree),'int',+1)
allocate (oint%icon(8,oint%nleaves),stat=threadinfo%err)
call heap (threadinfo,'oint%icon', 'embed_surf...',size(oint%icon),'int',+1)
allocate (oint%lsf(oint%nnode,oint%nlsf),stat=threadinfo%err)
call heap (threadinfo,'oint%lsf', 'embed_surf...',size(oint%lsf),'dp',+1)
allocate (oint%x(oint%nnode),stat=threadinfo%err)
call heap (threadinfo,'oint%x', 'embed_surf...',size(oint%x),'dp',+1)
allocate (oint%y(oint%nnode),stat=threadinfo%err)
call heap (threadinfo,'oint%y', 'embed_surf...',size(oint%y),'dp',+1)
allocate (oint%z(oint%nnode),stat=threadinfo%err)
call heap (threadinfo,'oint%z', 'embed_surf...',size(oint%z),'dp',+1)
oint%octree=osolve%octree
oint%icon=osolve%icon
if (size(oint%lsf).ne.size(osolve%lsf)) then
write(*,*) 'oint%lsf = osolve%lsf; sizes differ', size(oint%lsf), size(osolve%lsf)
endif
oint%lsf=osolve%lsf
oint%x=osolve%x
oint%y=osolve%y
oint%z=osolve%z
! having replicated osolve, we can deallocate it
call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','embed_surf...',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','embed_surf...',size(osolve%y),'dp',-1) ; deallocate (osolve%y)
call heap (threadinfo,'osolve%z','embed_surf...',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
call heap (threadinfo,'osolve%lsf','embed_surf...',size(osolve%lsf),'dp',-1) ; deallocate (osolve%lsf)
! we reinitialise it
call octree_init (osolve%octree,osolve%noctree)
osolve%octree=oint%octree
! we modify the osolve%octree by embedding the olsf in it
allocate (levs(olsf%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',+1)
call octree_find_element_level (olsf%octree,olsf%noctree,levs,olsf%nleaves)
do ie=1,olsf%nleaves
xxx=sum(olsf%x(olsf%icon(:,ie)))/8.d0
yyy=sum(olsf%y(olsf%icon(:,ie)))/8.d0
zzz=sum(olsf%z(olsf%icon(:,ie)))/8.d0
call octree_create_from_particles (osolve%octree,osolve%noctree,xxx,yyy,zzz,1,levs(ie))
enddo
call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',-1) ; deallocate (levs)
! we smoothen osolve
call octree_smoothen (osolve%octree,osolve%noctree)
! the number of leaves in osolve is computed by means of an octree function, and the
! number of nodes is over estimated
osolve%nleaves = ioctree_number_of_elements(osolve%octree,osolve%noctree)
osolve%nnode = osolve%nleaves*3
osolve%nlsf = params%ns
allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',+1)
allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%x', 'embed_surf...',size(osolve%x),'dp',+1)
allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%y', 'embed_surf...',size(osolve%y),'dp',+1)
allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%z', 'embed_surf...',size(osolve%z),'dp',+1)
! we build the icon array of osolve. the routine also outputs the real number of nodes osolve%nnode
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
osolve%icon,osolve%nleaves, &
osolve%x,osolve%y,osolve%z,osolve%nnode)
Douglas Guptill
committed
! osolve%nnode has been changed by octree_find_node_connectivity, so re-size x,y,z
call octreesolve_shrink_xyz(osolve, threadinfo)
! now that osolve%nnode is known we can allocate osolve%lsf
allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
call heap (threadinfo,'osolve%lsf', 'embed_surf...',size(osolve%lsf),'dp',+1)
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
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
! if we are dealing with the last surface, nodes are renumbered
if (params%renumber_nodes .and. is.eq.params%ns) then
call octree_renumber_nodes (osolve%icon,osolve%nleaves, &
osolve%x,osolve%y,osolve%z,osolve%nnode)
end if
! the lsf values of the previous is-1 surfaces are known on the oint octree which is no more
! the same as osolve, so we have to interpolates these lsfs on final octree osolve
do ilsf=1,is-1
do k=1,osolve%nnode
call octree_interpolate (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf),oint%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))
! call octree_interpolate3 (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf), &
! oint%x(1:oint%nnode),oint%y(1:oint%nnode),oint%z(1:oint%nnode),oint%nnode, &
! osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))
enddo
enddo
! finally the lsf corresponding to the current surface also has to be interpolated onto osolve
! since it was computed on a dedicated olsf
do k=1,osolve%nnode
call octree_interpolate (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,olsf%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
! 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, &
! osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
enddo
call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',-1) ; deallocate (olsf%x)
call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',-1) ; deallocate (olsf%y)
call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',-1) ; deallocate (olsf%z)
call heap (threadinfo,'olsf%lsf','embed_surf...',size(olsf%lsf),'dp',-1) ; deallocate (olsf%lsf)
call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',-1) ; deallocate (olsf%icon)
call heap (threadinfo,'olsf%octree','embed_surf...',size(olsf%octree),'int',-1) ; deallocate (olsf%octree)
call heap (threadinfo,'oint%x','embed_surf...',size(oint%x),'dp',-1) ; deallocate (oint%x)
call heap (threadinfo,'oint%y','embed_surf...',size(oint%y),'dp',-1) ; deallocate (oint%y)
call heap (threadinfo,'oint%z','embed_surf...',size(oint%z),'dp',-1) ; deallocate (oint%z)
call heap (threadinfo,'oint%lsf','embed_surf...',size(oint%lsf),'dp',-1) ; deallocate (oint%lsf)
call heap (threadinfo,'oint%icon','embed_surf...',size(oint%icon),'int',-1) ; deallocate (oint%icon)
call heap (threadinfo,'oint%octree','embed_surf...',size(oint%octree),'int',-1) ; deallocate (oint%octree)
if (params%debug>=1 .and. iproc==0) then
write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%nleaves
end if
write(threadinfo%Logunit,*) '[same] build_surface_octree: osolve%nleaves ', osolve%nleaves
call flush(threadinfo%Logunit)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|