!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
!use mpi
use threads
implicit none
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
integer err
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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)))
! 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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|