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
!use mpi
Dave Whipp
committed
use threads
include 'mpif.h'
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 ierr,iproc,nproc
integer ie,ilsf,k
double precision xxx,yyy,zzz
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)
if (params%debug.gt.1) 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)
if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) 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, params)
Douglas Guptill
committed
if (params%debug.gt.1) then
write(threadinfo%Logunit,*) '[same] embed_surface_in_octree: olsf%nnode ',olsf%nnode
call flush(threadinfo%Logunit)
endif
allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) 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)
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)
if (params%debug.gt.1) call heap (threadinfo,'oint%octree', 'embed_surf...',size(oint%octree),'int',+1)
allocate (oint%icon(8,oint%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%icon', 'embed_surf...',size(oint%icon),'int',+1)
allocate (oint%lsf(oint%nnode,oint%nlsf),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%lsf', 'embed_surf...',size(oint%lsf),'dp',+1)
allocate (oint%x(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%x', 'embed_surf...',size(oint%x),'dp',+1)
allocate (oint%y(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%y', 'embed_surf...',size(oint%y),'dp',+1)
allocate (oint%z(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) 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
if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',-1)
deallocate (osolve%icon)
if (params%debug.gt.1) call heap (threadinfo,'osolve%x','embed_surf...',size(osolve%x),'dp',-1)
deallocate (osolve%x)
if (params%debug.gt.1) call heap (threadinfo,'osolve%y','embed_surf...',size(osolve%y),'dp',-1)
deallocate (osolve%y)
if (params%debug.gt.1) call heap (threadinfo,'osolve%z','embed_surf...',size(osolve%z),'dp',-1)
deallocate (osolve%z)
if (params%debug.gt.1) 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)
if (params%debug.gt.1) 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
if (params%debug.gt.1) 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)
if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',+1)
allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%x', 'embed_surf...',size(osolve%x),'dp',+1)
allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%y', 'embed_surf...',size(osolve%y),'dp',+1)
allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) 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, params)
! now that osolve%nnode is known we can allocate osolve%lsf
allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf', 'embed_surf...',size(osolve%lsf),'dp',+1)
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
! 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
if (params%debug.gt.1) call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',-1)
deallocate (olsf%x)
if (params%debug.gt.1) call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',-1)
deallocate (olsf%y)
if (params%debug.gt.1) call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',-1)
deallocate (olsf%z)
if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','embed_surf...',size(olsf%lsf),'dp',-1)
deallocate (olsf%lsf)
if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',-1)
deallocate (olsf%icon)
if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','embed_surf...',size(olsf%octree),'int',-1)
deallocate (olsf%octree)
if (params%debug.gt.1) call heap (threadinfo,'oint%x','embed_surf...',size(oint%x),'dp',-1)
deallocate (oint%x)
if (params%debug.gt.1) call heap (threadinfo,'oint%y','embed_surf...',size(oint%y),'dp',-1)
deallocate (oint%y)
if (params%debug.gt.1) call heap (threadinfo,'oint%z','embed_surf...',size(oint%z),'dp',-1)
deallocate (oint%z)
if (params%debug.gt.1) call heap (threadinfo,'oint%lsf','embed_surf...',size(oint%lsf),'dp',-1)
deallocate (oint%lsf)
if (params%debug.gt.1) call heap (threadinfo,'oint%icon','embed_surf...',size(oint%icon),'int',-1)
deallocate (oint%icon)
if (params%debug.gt.1) 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
if (params%debug.gt.1) then
write(threadinfo%Logunit,*) '[same] build_surface_octree: osolve%nleaves ', osolve%nleaves
call flush(threadinfo%Logunit)
endif
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|