build_surface_octree.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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 threads
use definitions

implicit none
     
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,dist
real cosd
integer ioctree_number_of_elements
external ioctree_number_of_elements

double precision, dimension(:), allocatable ::xn,yn,zn
double precision :: x1,x2,x3,y1,y2,y3,z1,z2,z3,xne,yne,zne,xyzn
integer i3,ij,err
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

INCLUDE 'mpif.h'

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

write(threadinfo%Logunit,*) '[same] build_surface_octree: olsf%octree(2) ',olsf%octree(2)

call flush(threadinfo%Logunit) 


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

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|