OctreeBitPlus.f90

!-------------------------------------------

! NEW OCTREEBIT LIBRARY (WRITTEN 26-27/02/2004) BY JEAN BRAUN

! This version addresses problems that arose in the computation
! of the icon array for very large octrees. This new version is much
! more efficient and uses the same amount of memory (sometimes less).
! Many of the querry routines have been greatly improved by storing
! additional information in the octree structure. Data storage and
! navigation within the octree are based on the methods described by
! Frisken and Perry from Mitsubishi Electric Research Lab.
! The computation of the icon array is based on a three stage method:
! (1) Compute icon and a nodal position array disregarding the nodal connectivity
! (2) Sort the nodal arrays by ordering along x, y and z coordinates
! (3) Remove the spurious nodes and modify the icon array accordingly
! Because a fast sorting method is used (O(n log log n)), this approach
! is rather efficient.


!=======================!
!=====[OCTREE_INIT]=====!
!=======================!

subroutine octree_init (octree,noctree)

! This routine must be called before any other routine when an octree is created
! the basic structure of the octree is 
! octree(1)=maximum level (unit cube is level 0)
! octree(2)=number of leaves
! octree(3)=total length of octree
! for each cube in the octree (at location loc)
! octree(loc)=level
! octree(loc+1)=address of parent
! octree(loc+2 to loc+9)=address of children (if negative the child is a leaf and the
!         value is the leaf number in the sequence of leaves)

implicit none

integer noctree,octree(noctree)

integer levelmax,nleaves,length,loc,k

levelmax=1
nleaves=8
length=13

octree(1)=levelmax
octree(2)=nleaves
octree(3)=length

loc=4

octree(loc)=1
octree(loc+1)=0
  do k=1,8
  octree(loc+1+k)=-k
  enddo

return
end

!====================================!
!=====[FIND_INTEGER_COORDINATES]=====!
!====================================!

subroutine find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)

! returns the integer coordinates of point (x,y,z) in (ix,iy,iz)
! the integer coordinates are determined by the maximum level in the octree
! levelmax

! the integer coordinate is comprise between 1 and 2**levelmax
! and corresponds to the location of the leaf containing the
! point of given coordinates

implicit none

double precision x,y,z
integer ix,iy,iz,levelmax

double precision powermax

powermax=2.d0**levelmax
ix=int(x*powermax)
iy=int(y*powermax)
iz=int(z*powermax)

return
end

!=================================!
!=====[FIND_REAL_COORDINATES]=====!
!=================================!

subroutine find_real_coordinates (ix,iy,iz,x,y,z,levelmax)

! returns the real coordinates of a point of integer coordinates (ix,iy,iz)
! the real coordinates are determined by the maximum level in the octree

implicit none

double precision x,y,z
integer ix,iy,iz,levelmax

double precision powermax

powermax=2**levelmax
x=dfloat(ix)/powermax
y=dfloat(iy)/powermax
z=dfloat(iz)/powermax

return
end

!========================================!
!=====[OCTREE_CREATE_FROM_PARTICLE]=====!
!========================================!

subroutine octree_create_from_particle (octree,noctree,x,y,z,np,level)

! updates the octree by creating a leaf at point (x,y,z)
! of level level
! if the leaf (or a cube of smaller level) exists, the routine has no effect on the octree
! note that x,y,z must belong to [0,1[

implicit none

integer noctree,octree(noctree),np
double precision x,y,z
integer level
double precision xp,yp,zp

integer levelin,ix,iy,iz,levelmax,loc,ip

levelmax=octree(1)

xp=x
yp=y
zp=z
levelin=0
loc=4
  if (xp.eq.1.d0) xp=1.d0-1.d-20
  if (yp.eq.1.d0) yp=1.d0-1.d-20
  if (zp.eq.1.d0) zp=1.d0-1.d-20
  if (xp.eq.0.d0) xp=1.d-20
  if (yp.eq.0.d0) yp=1.d-20
  if (zp.eq.0.d0) zp=1.d-20
  if (xp*(xp-1.d0).lt.0.d0 .and. yp*(yp-1.d0).lt.0.d0 .and. zp*(zp-1.d0).lt.0.d0) &
    call update (octree,noctree,xp,yp,zp,level,levelin,loc)

return
end

!========================================!
!=====[OCTREE_CREATE_FROM_PARTICLES]=====!
!========================================!

subroutine octree_create_from_particles (octree,noctree,x,y,z,np,level)

! updates the octree by creating a leaf at points (x(1:np),y(1:np),z(1:np))
! of level level(1:np)
! if the leaf (or a cube of smaller level) exists, the routine has no effect on the octree
! note that x,y,z must belong to [0,1[

implicit none

integer noctree,octree(noctree),np
double precision x(np),y(np),z(np)
integer level(np)
double precision xp,yp,zp

integer levelin,ix,iy,iz,levelmax,loc,ip

levelmax=octree(1)

do ip=1,np
xp=x(ip)
yp=y(ip)
zp=z(ip)
levelin=0
loc=4
  if (xp.eq.1.d0) xp=1.d0-1.d-20
  if (yp.eq.1.d0) yp=1.d0-1.d-20
  if (zp.eq.1.d0) zp=1.d0-1.d-20
  if (xp.eq.0.d0) xp=1.d-20
  if (yp.eq.0.d0) yp=1.d-20
  if (zp.eq.0.d0) zp=1.d-20
  if (xp*(xp-1.d0).lt.0.d0 .and. yp*(yp-1.d0).lt.0.d0 .and. zp*(zp-1.d0).lt.0.d0) &
    call update (octree,noctree,xp,yp,zp,level(ip),levelin,loc)
enddo

return
end

!==================!
!=====[UPDATE]=====!
!==================!

recursive subroutine update (octree,noctree,x,y,z,level,levelin,loc)

! DO NOT USE
! internal routine called by OctreeUpdate
! this routine divides a cube of the octree in 8
! if the requested level is greater than0
! and if the requested level is stricly greater than levelin
! loc is the location in the octree of the current cube (to be divided)

implicit none

integer noctree,octree(noctree),ix,iy,iz,levelin,level,loc
double precision x,y,z
integer ibits_jean
external ibits_jean

integer levelmax,ibitx,ibity,ibitz,ichild,iaddress,newaddress

if (level.eq.0) return

if (levelin+1.eq.level) return

levelmax=octree(1)
call find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)
ibitx=ibits_jean(ix,levelmax-levelin-1)
ibity=ibits_jean(iy,levelmax-levelin-1)
ibitz=ibits_jean(iz,levelmax-levelin-1)

ichild=loc+2+ibitz*4+ibity*2+ibitx

if (ichild.gt.octree(3)) then
print*,ichild,' requested ',octree(3),' available'
stop 'octree needs to grow...'
endif

iaddress=octree(ichild)

  if (iaddress.lt.0) then
!adds a level
  newaddress=octree(3)+1
! updates current level
  octree(ichild)=newaddress
! updates octree length
  octree(3)=octree(3)+10
! updates number of leaves
  octree(2)=octree(2)+7
! updates maximum level
  if (levelin+2.gt.octree(1)) octree(1)=levelin+2
! creates new division
  octree(newaddress)=levelin+2
  octree(newaddress+1)=ichild
  octree(newaddress+2)=iaddress
  octree(newaddress+3)=-(octree(2)-6)
  octree(newaddress+4)=-(octree(2)-5)
  octree(newaddress+5)=-(octree(2)-4)
  octree(newaddress+6)=-(octree(2)-3)
  octree(newaddress+7)=-(octree(2)-2)
  octree(newaddress+8)=-(octree(2)-1)
  octree(newaddress+9)=-(octree(2))
  call update (octree,noctree,x,y,z,level,levelin+1,newaddress)
  else
  call update (octree,noctree,x,y,z,level,levelin+1,iaddress)
  endif

return
end


!=================!
!=====[IBITS]=====!
!=================!

integer function ibits_jean(i,ipos)

!internal function, do not modify

integer j
j=i/2**ipos
ibits_jean=j-(j/2)*2
return
end

!============================!
!=====[OCTREE_FIND_LEAF]=====!
!============================!

subroutine octree_find_leaf (octree,noctree,x,y,z,leaf,level,loc,x0,y0,z0,dxyz)

! given an octree of size noctree
! this routine returns the leaf number in which a point (x,y,z) resides
! the level of the leaf (0 is unit cube)
! the location in the octree of the part describing the parent of the leaf (loc)
! the centroid of the leaf (x0,y0,z0) and its size (dxyz)

implicit none

integer noctree,octree(noctree)
double precision x,y,z,x0,y0,z0,dxyz
integer leaf,level,loc

integer ix,iy,iz,levelin,locin,levelmax

levelmax=octree(1)
call find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)

levelin=0
locin=4
leaf=0
call find_leaf (octree,noctree,ix,iy,iz,levelin,locin,leaf,level,loc)

call find_integer_coordinates (x,y,z,ix,iy,iz,level)
call find_real_coordinates (ix,iy,iz,x0,y0,z0,level)
dxyz=1.d0/2.d0**level

return
end

!=====================!
!=====[FIND_LEAF]=====!
!=====================!

recursive subroutine find_leaf (octree,noctree,ix,iy,iz,levelin,locin,leaf,level,loc)

! DO NOT USE
! internal routine used by octree_find_leaf

implicit none

integer noctree,octree(noctree),ix,iy,iz,levelin,leaf,level,loc,locin
integer levelmax,ibitx,ibity,ibitz,ichild,iaddress
integer ibits_jean
external ibits_jean

levelmax=octree(1)
!ibitx=ibits(ix,levelmax-levelin-1,1)
!ibity=ibits(iy,levelmax-levelin-1,1)
!ibitz=ibits(iz,levelmax-levelin-1,1)
ibitx=ibits_jean(ix,levelmax-levelin-1)
ibity=ibits_jean(iy,levelmax-levelin-1)
ibitz=ibits_jean(iz,levelmax-levelin-1)

ichild=locin+2+ibitz*4+ibity*2+ibitx

if (ichild.gt.octree(3)) stop 'octree needs to grow...'
iaddress=octree(ichild)

  if (iaddress.lt.0) then
  leaf=-iaddress
  level=levelin+1
  loc=locin
  return
  else
  call find_leaf (octree,noctree,ix,iy,iz,levelin+1,iaddress,leaf,level,loc)
  endif

return
end

!===============================!
!=====[OCTREETHROUGHLEAVES]=====!
!===============================!

subroutine OctreeThroughLeaves (octree,noctree)

! DO NOT USE
! this is a general routine that simply goes through the leaves
! this routine is used as a template to build other routines.
! it should not be used

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz

loc=4
ix=0
iy=0
iz=0
call throughleaves (octree,noctree,loc,ix,iy,iz)

return
end

!=========================!
!=====[THROUGHLEAVES]=====!
!=========================!

recursive subroutine throughleaves (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    else
    call throughleaves (octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!===========================!
!=====[OCTREE_SMOOTHEN]=====!
!===========================!

subroutine octree_smoothen (octree,noctree)

! as it names indicates this routine smoothens the octree
! ie it ensures that no two adjacent leaves are more than one level apart

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz,nleaves,nleaves0
integer ioctree_number_of_elements

nleaves=ioctree_number_of_elements (octree,noctree)

do while (nleaves.ne.nleaves0)
nleaves0=nleaves
loc=4
ix=0
iy=0
iz=0
call smooth (octree,noctree,loc,ix,iy,iz)
nleaves=ioctree_number_of_elements (octree,noctree)
enddo

return
end

!==================!
!=====[SMOOTH]=====!
!==================!

recursive subroutine smooth (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer k,idx,idy,idz,ip
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp
double precision xp,yp,zp

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)

! first check that `right' neighbours are not down by more than 1 level
    ipmax=2**levelmax
    ip=2**level
    izp=izn+ipower
    iyp=iyn+ipower
    ixp=ixn+ipower
      if (izp.lt.ipmax) then
      zp=dfloat(izp)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (iyp.lt.ipmax) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyp)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (ixp.lt.ipmax) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixp)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif

! second check if the 'left' neighbours are not down by more than one level
    izp=izn-1
    iyp=iyn-1
    ixp=ixn-1
      if (izp.ge.0) then
      zp=dfloat(izp)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (iyp.ge.0) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyp)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (ixp.ge.0) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixp)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif

    else
    call smooth(octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!=================================!
!=====[OCTREE_SUPER_SMOOTHEN]=====!
!=================================!

subroutine octree_super_smoothen (octree,noctree)

! as it names indicates this routine further smoothens the octree
! ie it ensures that for any given leaf and its six closest neighbouring
! leaves, ther is no difference in level larger than 1

! Note that contrary to octree_smoothen, octree_super_smoothen can
! be used iteratively to increase the smoothness of the octree

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz,nleaves,nleaves0
integer ioctree_number_of_elements

nleaves=ioctree_number_of_elements (octree,noctree)

do while (nleaves.ne.nleaves0)
nleaves0=nleaves
loc=4
ix=0
iy=0
iz=0
call super_smooth (octree,noctree,loc,ix,iy,iz,nleaves)
nleaves=ioctree_number_of_elements (octree,noctree)
enddo

call octree_smoothen (octree,noctree)

return
end

!========================!
!=====[SUPER_SMOOTH]=====!
!========================!

recursive subroutine super_smooth (octree,noctree,loc,ix,iy,iz,nleaves0)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer k,idx,idy,idz,ip
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp
double precision xp,yp,zp
integer levelout(6),levelset,leaf,locin,nleaves0

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)

! first find `right' neighbours levels
    ipmax=2**levelmax
    ip=2**level
    izp=izn+ipower
    iyp=iyn+ipower
    ixp=ixn+ipower
    levelout=level
      if (izp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyn,izp,levelin,locin,leaf,levelout(1),locp)
      if (leaf.gt.nleaves0) levelout(1)=level
      endif
      if (iyp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyp,izn,levelin,locin,leaf,levelout(2),locp)
      if (leaf.gt.nleaves0) levelout(2)=level
      endif
      if (ixp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixp,iyn,izn,levelin,locin,leaf,levelout(3),locp)
      if (leaf.gt.nleaves0) levelout(3)=level
      endif

! second find 'left' neighbours levels
    izp=izn-1
    iyp=iyn-1
    ixp=ixn-1
      if (izp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyn,izp,levelin,locin,leaf,levelout(4),locp)
      if (leaf.gt.nleaves0) levelout(4)=level
      endif
      if (iyp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyp,izn,levelin,locin,leaf,levelout(5),locp)
      if (leaf.gt.nleaves0) levelout(5)=level
      endif
      if (ixp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixp,iyn,izn,levelin,locin,leaf,levelout(6),locp)
      if (leaf.gt.nleaves0) levelout(6)=level
      endif

! if the difference in level between the leaf and any of its neighbours is greater
! than 1, all leaves are set to the maximum level minus one at least

      if (maxval(levelout)-minval(levelout).gt.1) then
      levelset=maxval(levelout)-1
! now check that `right' neighbours are at least at this level
      izp=izn+ipower
      iyp=iyn+ipower
      ixp=ixn+ipower
        if (izp.lt.ipmax) then
        zp=dfloat(izp)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (iyp.lt.ipmax) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyp)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (ixp.lt.ipmax) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixp)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
  ! now check that 'left' neighbours are at least at this level
      izp=izn-1
      iyp=iyn-1
      ixp=ixn-1
        if (izp.ge.0) then
        zp=dfloat(izp)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (iyp.ge.0) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyp)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (ixp.ge.0) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixp)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
  
      endif

    else
    call super_smooth(octree,noctree,locn,ixn,iyn,izn,nleaves0)
    endif
  enddo
  enddo
  enddo

return
end

!===============================!
!=====[OCTREE_SHOW_AS_MESH]=====!
!===============================!

subroutine octree_show_as_mesh (octree,noctree)

! routine to show the octree as a mesh in vrml
! produces a VRML file called mesh.wrl
! this is mostly used for debugging purposes

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz

open (7,file='mesh.wrl',status='unknown')
write (7,'(a)') '#VRML V2.0 utf8'
write (7,'(a)') 'Transform { children ['
write (7,'(a)') 'NavigationInfo { '
write (7,'(a)') 'type ["EXAMINE"]'
write (7,'(a)') 'headlight FALSE}'
write (7,'(a)') 'Background{groundColor 1 1 1 skyColor 1 1 1}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  .8 1 .5}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  -.8 -1 -.5}'
write (7,'(a)') 'Transform { children Viewpoint {'
write (7,'(a)') '      description "Starting"'
write (7,'(a)') '      fieldOfView 1'
write (7,'(a,3f12.8,a)') '      position    ',-.41885125637054443, -.8311104774475098, 1.5406757593154907
write (7,'(a,4f12.8,a)') '      orientation ', .7352051138877869, -.10698327422142029, -.669348955154419, 1.3260198831558228,'}}'
write (7,'(a,a)') 'DEF Node0 Shape{geometry Sphere{radius 0.0075', &
                ' }appearance Appearance{material Material{diffuseColor 1 0 0}}}'
write (7,'(a,3f10.6,a)') 'Transform{translation',0.,0.,0., &
                         ' children [USE Node0]}'
loc=4
ix=0
iy=0
iz=0
call show (octree,noctree,loc,ix,iy,iz)
write (7,'(a)') ']}'
close (7)

return
end

!================!
!=====[SHOW]=====!
!================!

recursive subroutine show (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! to show the octree as a mesh

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k
double precision x1,y1,z1,x2,y2,z2,dxyz

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    call find_real_coordinates (ixn,iyn,izn,x1,y1,z1,levelmax)
    dxyz=1.d0/(2.d0**level)
    x2=x1+dxyz
    y2=y1+dxyz
    z2=z1+dxyz
    write (7,'(a,24f10.3,a,a,a)') &
         'Shape { geometry IndexedLineSet { coord Coordinate { point [', &
          x1,y1,z1,x2,y1,z1,x2,y2,z1,x1,y2,z1,x1,y1,z2,x2,y1,z2,x2,y2,z2,x1,y2,z2, &
          ']} coordIndex [0 1 2 3 0 -1 4 5 6 7 4 -1 0 4 -1 1 5 -1 2 6 -1 3 7 -1', &
          ']}appearance Appearance { material Material { emissiveColor 0 0 0}}}'
    else
    call show (octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!=====================================!
!=====[OCTREE_FIND_ELEMENT_LEVEL]=====!
!=====================================!

subroutine octree_find_element_level (octree,noctree,levs,nleaves)

! routine to return the level of each leaf
! the result is returned in the array levs of dimension nleaves

! the function ioctree_number_of_elements should be called first to
! find nleaves and dimension levs accordingly

implicit none

integer noctree,octree(noctree),nleaves,levs(nleaves)

integer loc,ix,iy,iz

loc=4
ix=0
iy=0
iz=0
call levels (octree,noctree,loc,ix,iy,iz,levs,nleaves)

return
end

!==================!
!=====[LEVELS]=====!
!==================!

recursive subroutine levels (octree,noctree,loc,ix,iy,iz,levs,nleaves)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz
integer nleaves,levs(nleaves)

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    if (-locn.gt.nleaves) stop 'array level needs to grow'
    levs(-locn)=level
    else
    call levels (octree,noctree,locn,ixn,iyn,izn,levs,nleaves)
    endif
  enddo
  enddo
  enddo

return
end

!======================================!
!=====[IOCTREE_NUMBER_OF_ELEMENTS]=====!
!======================================!

integer function ioctree_number_of_elements (octree,noctree)

! function that returns (in ioctree_number_of_elements) the number of leaves in the octree

implicit none

integer noctree,octree(noctree)

ioctree_number_of_elements=octree(2)

return
end

!=================================!
!=====[IOCTREE_MAXIMUM_LEVEL]=====!
!=================================!

integer function ioctree_maximum_level (octree,noctree)

! function that returns (in ioctree_maximum_level) the maximum level in the octree

implicit none

integer noctree,octree(noctree)

ioctree_maximum_level=octree(1)

return
end

!========================!
!=====[IOCTREE_SIZE]=====!
!========================!

integer function ioctree_size (octree,noctree)

! function that returns (in ioctree_size) the size of the octree

implicit none

integer noctree,octree(noctree)

ioctree_size=octree(3)

return
end

!=================================!
!=====[OCTREE_CREATE_UNIFORM]=====!
!=================================!

subroutine octree_create_uniform (octree,noctree,levelt)

! routine to generate a uniform octree down to level levelt

implicit none

integer noctree,octree(noctree),levelt

integer nl,levelin,loc,iz,ix,iy
double precision x,y,z

nl=2**levelt
do iz=0,nl-1
z=dfloat(iz)/nl
do iy=0,nl-1
y=dfloat(iy)/nl
do ix=0,nl-1
x=dfloat(ix)/nl
levelin=0
loc=4
call update (octree,noctree,x,y,z,levelt,levelin,loc)
enddo
enddo
enddo

return
end


!=================================!
!=====[OCTREE_RENUMBER_NODES]=====!
!=================================!

subroutine octree_renumber_nodes (icon,nleaves,xa,ya,za,na)

! This routine renumbers the nodes to minimize the maximum
! difference (in the least square sense) between the numbers
! of the nodes belonging to any given element.

! This is useful is one wishes to use the octree mesh for FE
! calculations.

! This uses SLOAN's algorithm and routines

! in output, the icon array is modified as well as the node
! location arrays (xa,ya,za) that have been reordered

implicit none

integer nleaves,icon(8,nleaves),na
double precision xa(na),ya(na),za(na)
integer,dimension(:),allocatable::npn,xnpn,adj,xadj,sort,working
integer,dimension(:,:),allocatable::jcon
double precision,dimension(:),allocatable::xyz
integer ie,k,inpn,iadj,oldpro,newpro

inpn=nleaves*8
iadj=nleaves*8*7

allocate (xnpn(nleaves+1),npn(inpn),adj(iadj),xadj(na+1))
allocate (sort(na),working(3*na+1))

xnpn(1)=1
  do ie=1,nleaves
  xnpn(ie+1)=xnpn(ie)+8-1
    do k=1,8
    npn(xnpn(ie)+k-1)=icon(k,ie)
    enddo
  enddo

call graph_sloan (na,nleaves,inpn,npn,xnpn,iadj,adj,xadj)

call label_sloan (na,xadj(na+1)-1,adj,xadj,sort,working,oldpro,newpro)

deallocate (xnpn,npn,adj,xadj,working)

allocate (jcon(8,nleaves),xyz(na))

  do ie=1,nleaves
    do k=1,8
    jcon(k,ie)=sort(icon(k,ie))
    enddo
  enddo
icon=jcon

  do k=1,na
  xyz(sort(k))=xa(k)
  enddo
xa=xyz

  do k=1,na
  xyz(sort(k))=ya(k)
  enddo
ya=xyz

  do k=1,na
  xyz(sort(k))=za(k)
  enddo
za=xyz

deallocate (sort,jcon,xyz)

return
end

!=========================================!
!=====[OCTREE_FIND_NODE_CONNECTIVITY]=====!
!=========================================!

subroutine octree_find_node_connectivity (octree,noctree,icon,nleaves,xa,ya,za,na)

! This routine finds the number (na) and locations (xa,ya,za)
! of the nodes of the octree
! It also computes the connectivity array between nodes and leaves
! (icon). Icon is dimensioned icon(8,nleaves) and contains the number
! of the 8 nodes connected by each leaf
! When calling this routine, na shold have the dimension used to dimension
! the coordinate arrys in the calling routine
! on return na contains the true dimension of these array (ie how many nodes 
! there are in the octree

implicit none

integer noctree,octree(noctree),nleaves,icon(8,nleaves),na
double precision xa(*),ya(*),za(*)
integer,dimension(:),allocatable::kx,ky,kz,isort,jsort,ksort
integer loc,ix,iy,iz,nk
integer k1,nnk,kk1,nnnk
integer i0,in,levelmax,npower,il,i,k,l,ii

allocate (kx(8*nleaves),ky(8*nleaves),kz(8*nleaves))

! first build a general/redondant icon array

loc=4
ix=0
iy=0
iz=0
nk=0
call iconfind (octree,noctree,loc,ix,iy,iz,icon,kx,ky,kz,nk)

if (nk.ne.8*nleaves) stop 'error in iconfind'

allocate (isort(nk),jsort(nk),ksort(nk))

! here we rank the nodes according to their x, y and then z coordinates

call mrgrnk (kx,isort,nk)
call sort (kx,isort,nk)
call sort (ky,isort,nk)
call sort (kz,isort,nk)

jsort=isort

k1=1
do i=1,nk
  if (kx(i).gt.kx(k1).or.i.eq.nk) then
  nnk=i-k1
    if (nnk.gt.1) then
    call mrgrnk (ky(k1),isort(k1),nnk)
    call sort (kx(k1),isort(k1),nnk)
    call sort (ky(k1),isort(k1),nnk)
    call sort (kz(k1),isort(k1),nnk)
      do l=k1,k1+nnk-1
      ksort(l)=jsort(k1+isort(l)-1)
      enddo
      do l=k1,k1+nnk-1
      jsort(l)=ksort(l)
      enddo
    kk1=k1
    do ii=k1,k1+nnk-1
      if (ky(ii).gt.ky(kk1).or.ii.eq.k1+nnk-1) then
      nnnk=ii-kk1
        if (nnnk.gt.1) then
        call mrgrnk (kz(kk1),isort(kk1),nnnk)
        call sort (kx(kk1),isort(kk1),nnnk)
        call sort (ky(kk1),isort(kk1),nnnk)
        call sort (kz(kk1),isort(kk1),nnnk)
          do l=kk1,kk1+nnnk-1
          ksort(l)=jsort(kk1+isort(l)-1)
          enddo
          do l=kk1,kk1+nnnk-1
          jsort(l)=ksort(l)
          enddo
        endif
      kk1=ii
      endif
    enddo
    endif
  k1=i
  endif
enddo

isort=jsort

do i=1,nk
jsort(isort(i))=i
enddo

! at this point isort(i) is the location of the point of rank i
! whereas jsort(i) is the rank of the point of location i

! we modify isort to remove reference to identical points

! we introduce ksort to reorder the node into consecutive numbers

i0=1
na=1
do i=1,nk
in=0
if (kz(i).eq.kz(i0)) then
if (ky(i).eq.ky(i0)) then
if (kx(i).eq.kx(i0)) then
isort(i)=isort(i0)
ksort(isort(i))=na
in=1
endif
endif
endif
  if (in.eq.0) then
  i0=i
  na=na+1
  isort(i)=isort(i0)
  ksort(isort(i))=na
  endif
enddo

!print*,'There are ',na,' nodes'

! we modify icon to represent the new node representation

do il=1,nleaves
do k=1,8
icon(k,il)=ksort(isort(jsort(icon(k,il))))
enddo
enddo

! we build coordinates arrays

deallocate (jsort)

allocate (jsort(na))

! we define the correspondence between global and reduced node set

do i=1,nk
jsort(ksort(isort(i)))=i
enddo

! we compute the geometry of the nodes

levelmax=octree(1)
npower=2**levelmax
do i=1,na
xa(i)=dfloat(kx(jsort(i)))/npower
ya(i)=dfloat(ky(jsort(i)))/npower
za(i)=dfloat(kz(jsort(i)))/npower
enddo

deallocate (kx,ky,kz,isort,jsort,ksort)

return
end

!================!
!=====[SORT]=====!
!================!

subroutine sort (k,is,n)

! DO NOT USE
! this routine is used to sort an array according to an sorting array

implicit none

integer n,k(n),is(n)
integer,dimension(:),allocatable::kk
integer i

allocate (kk(n))

do i=1,n
kk(i)=k(is(i))
enddo

do i=1,n
k(i)=kk(i)
enddo

deallocate(kk)

return
end

!====================!
!=====[ICONFIND]=====!
!====================!

recursive subroutine iconfind (octree,noctree,loc,ix,iy,iz,icon,kx,ky,kz,nk)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz
integer icon(8,*),kx(*),ky(*),kz(*),nk

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k,kkx,kky,kkz

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
      do kkz=0,1
      do kky=0,1
      do kkx=0,1
      nk=nk+1
      k=kkx+kky*2+kkz*4
      icon(k+1,-locn)=nk
      kx(nk)=ixn+kkx*ipower
      ky(nk)=iyn+kky*ipower
      kz(nk)=izn+kkz*ipower
      enddo
      enddo
      enddo
    else
    call iconfind (octree,noctree,locn,ixn,iyn,izn,icon,kx,ky,kz,nk)
    endif
  enddo
  enddo
  enddo

return
end

!============================!
!=====[OCTREE_SHOW_ICON]=====!
!============================!

subroutine octree_show_icon (icon,nelem,x,y,z,nnode)

! this routine creates a VRML file called icon.wrl
! that shows the octree as a mesh and the nodes as small spheres
! note that this routine does use icon
! note that it uses the octree definition for icon

! it creates a VRML file called icon.wrl

implicit none

integer nelem,icon(8,nelem),nnode
double precision x(nnode),y(nnode),z(nnode)
integer ie,i,kkk

open (7,file='icon.wrl',status='unknown')
write (7,'(a)') '#VRML V2.0 utf8'
write (7,'(a)') 'Transform { children ['
write (7,'(a)') 'NavigationInfo { '
write (7,'(a)') 'type ["EXAMINE"]'
write (7,'(a)') 'headlight FALSE}'
write (7,'(a)') 'Background{groundColor 1 1 1 skyColor 1 1 1}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  .8 1 .5}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  -.8 -1 -.5}'
write (7,'(a)') 'Transform { children Viewpoint {'
write (7,'(a)') '      description "Starting"'
write (7,'(a)') '      fieldOfView 1'
write (7,'(a,3f12.8,a)') '      position    ',-.41885125637054443, -.8311104774475098, 1.5406757593154907
write (7,'(a,4f12.8,a)') '      orientation ', .7352051138877869, -.10698327422142029, -.669348955154419, 1.3260198831558228,'}}'
write (7,'(a,a)') 'DEF Node0 Shape{geometry Sphere{radius 0.0075', &
                      ' }appearance Appearance{material Material{diffuseColor 1 0 0}}}'

do ie=1,nelem
write (7,'(a,24f10.3,a,a,a)') &
     'Shape { geometry IndexedLineSet { coord Coordinate { point [', &
      (x(icon(kkk,ie)),y(icon(kkk,ie)),z(icon(kkk,ie)),kkk=1,8), &
     ']} coordIndex [0 1 3 2 0 -1 4 5 7 6 4 -1 0 4 -1 1 5 -1 2 6 -1 3 7 -1', &
     ']}appearance Appearance { material Material { emissiveColor 0 0 0}}}'
enddo

do i=1,nnode
write (7,'(a,3f10.6,a)') 'Transform{translation',x(i),y(i),z(i), &
                         ' children [USE Node0]}'
enddo

write (7,'(a)') ']}'
close (7)

return
end

!=================================!
!=====[OCTREE_FIND_BAD_FACES]=====!
!=================================!

subroutine octree_find_bad_faces (octree,noctree,iface,nface,icon,nelem)

! returns the bad faces as an array (iface) of 9 nodes per face
! numbering used is different from earlier version of octreebit
! here it is:

! 4--7--3
! |  |  |
! 8--9--6
! |  |  |
! 1--5--2

! iface is the resulting bad face information iface(9,nface)
! nface is the number of bad faces found
! icon is the connectivity array
! of dimension nelem

implicit none

integer noctree,octree(noctree)
integer nface,iface(9,nface),nelem,icon(8,nelem)

integer loc,ix,iy,iz,mface

loc=4
ix=0
iy=0
iz=0
mface=nface
nface=0
call badface (octree,noctree,loc,ix,iy,iz,iface,nface,mface,icon,nelem)

return
end

!===================!
!=====[BADFACE]=====!
!===================!

recursive subroutine badface (octree,noctree,loc,ix,iy,iz,iface,nface,mface,icon,nelem)

! DO NOT USE
! internal routine used by octree_find_bad_faces

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz
integer mface,iface(9,mface),nelem,icon(8,nelem),nface

integer level,levelmax,ipower,ixn,iyn,izn,locn,ipowerne
integer k,idx,idy,idz,ip,ixpp,iypp,izpp
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp,locne,levelne,leaf
double precision xp,yp,zp,x0,y0,z0,dxyz

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)

! first check that 'right' neighbours are of a higher level
    ipmax=2**levelmax
    ip=2**level
    ixp=ixn+ipower
    iyp=iyn+ipower
    izp=izn+ipower
      if (izp.lt.ipmax) then
      zp=dble(izp)/ipmax
      yp=dble(iyn)/ipmax
      xp=dble(ixn)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(5,-locn)
        iface(2,nface)=icon(6,-locn)
        iface(3,nface)=icon(8,-locn)
        iface(4,nface)=icon(7,-locn)
        iface(5,nface)=icon(2,leaf)
        iface(8,nface)=icon(3,leaf)
        iface(9,nface)=icon(4,leaf)
        ixpp=ixn+ipowerne
        xp=dble(ixpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(4,leaf)
        iypp=iyn+ipowerne
        yp=dble(iypp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(3,leaf)
        endif
      endif
      if (iyp.lt.ipmax) then
      zp=dble(izn)/ipmax
      yp=dble(iyp)/ipmax
      xp=dble(ixn)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(3,-locn)
        iface(2,nface)=icon(7,-locn)
        iface(3,nface)=icon(8,-locn)
        iface(4,nface)=icon(4,-locn)
        iface(5,nface)=icon(5,leaf)
        iface(8,nface)=icon(2,leaf)
        iface(9,nface)=icon(6,leaf)
        izpp=izn+ipowerne
        zp=dble(izpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(6,leaf)
        ixpp=ixn+ipowerne
        xp=dble(ixpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(2,leaf)
        endif
      endif
      if (ixp.lt.ipmax) then
      zp=dble(izn)/ipmax
      yp=dble(iyn)/ipmax
      xp=dble(ixp)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(2,-locn)
        iface(2,nface)=icon(4,-locn)
        iface(3,nface)=icon(8,-locn)
        iface(4,nface)=icon(6,-locn)
        iface(5,nface)=icon(3,leaf)
        iface(8,nface)=icon(5,leaf)
        iface(9,nface)=icon(7,leaf)
        iypp=iyn+ipowerne
        yp=dble(iypp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(7,leaf)
        izpp=izn+ipowerne
        zp=dble(izpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(5,leaf)
        endif
      endif

! second check if the 'left' neighbours are of a higher level
    ixp=ixn-1
    iyp=iyn-1
    izp=izn-1
      if (izp.ge.0) then
      zp=dble(izp)/ipmax
      yp=dble(iyn)/ipmax
      xp=dble(ixn)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(1,-locn)
        iface(2,nface)=icon(2,-locn)
        iface(3,nface)=icon(4,-locn)
        iface(4,nface)=icon(3,-locn)
        iface(5,nface)=icon(6,leaf)
        iface(8,nface)=icon(7,leaf)
        iface(9,nface)=icon(8,leaf)
        ixpp=ixn+ipowerne
        xp=dble(ixpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(8,leaf)
        iypp=iyn+ipowerne
        yp=dble(iypp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(7,leaf)
        endif
      endif
      if (iyp.ge.0) then
      zp=dble(izn)/ipmax
      yp=dble(iyp)/ipmax
      xp=dble(ixn)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(1,-locn)
        iface(2,nface)=icon(2,-locn)
        iface(3,nface)=icon(6,-locn)
        iface(4,nface)=icon(5,-locn)
        iface(5,nface)=icon(4,leaf)
        iface(8,nface)=icon(7,leaf)
        iface(9,nface)=icon(8,leaf)
        izpp=izn+ipowerne
        zp=dble(izpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(8,leaf)
        ixpp=ixn+ipowerne
        xp=dble(ixpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(4,leaf)
        endif
      endif
      if (ixp.ge.0) then
      zp=dble(izn)/ipmax
      yp=dble(iyn)/ipmax
      xp=dble(ixp)/ipmax
      levelin=0
      locp=4
      call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        if (levelne.gt.level) then
        ipowerne=2**(levelmax-levelne)
        nface=nface+1
        if (nface.gt.mface) stop 'nface needs to grow'
        iface(1,nface)=icon(1,-locn)
        iface(2,nface)=icon(5,-locn)
        iface(3,nface)=icon(7,-locn)
        iface(4,nface)=icon(3,-locn)
        iface(5,nface)=icon(6,leaf)
        iface(8,nface)=icon(4,leaf)
        iface(9,nface)=icon(8,leaf)
        iypp=iyn+ipowerne
        yp=dble(iypp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(7,nface)=icon(8,leaf)
        izpp=izn+ipowerne
        zp=dble(izpp)/ipmax
        call octree_find_leaf (octree,noctree,xp,yp,zp,leaf,levelne,locne,x0,y0,z0,dxyz)
        iface(6,nface)=icon(6,leaf)
        endif
      endif

    else
    call badface (octree,noctree,locn,ixn,iyn,izn,iface,nface,mface,icon,nelem)
    endif
  enddo
  enddo
  enddo

return
end

!=================================!
!=====[OCTREE_SHOW_BAD_FACES]=====!
!=================================!

subroutine octree_show_bad_faces (octree,noctree,iface,nface,x,y,z,nnode)

! this routine creates a VRML file called bad_faces.wrl
! that shows the octree as produced by show_octree but it
! adds color to the bad faces to locate them

implicit none

integer noctree,octree(noctree)
integer nface,nnode,iface(9,nface)
double precision x(nnode),y(nnode),z(nnode)
integer loc,ix,iy,iz,i,k

open (7,file='bad_faces.wrl',status='unknown')
write (7,'(a)') '#VRML V2.0 utf8'
write (7,'(a)') 'Transform { children ['
write (7,'(a)') 'NavigationInfo { '
write (7,'(a)') 'type ["EXAMINE"]'
write (7,'(a)') 'headlight FALSE}'
write (7,'(a)') 'Background{groundColor 1 1 1 skyColor 1 1 1}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  .8 1 .5}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  -.8 -1 -.5}'
write (7,'(a)') 'Transform { children Viewpoint {'
write (7,'(a)') '      description "Starting"'
write (7,'(a)') '      fieldOfView 1'
write (7,'(a,3f12.8,a)') '      position    ',-.41885125637054443, -.8311104774475098, 1.5406757593154907
write (7,'(a,4f12.8,a)') '      orientation ', .7352051138877869, -.10698327422142029, -.669348955154419, 1.3260198831558228,'}}'

loc=4
ix=0
iy=0
iz=0
call show (octree,noctree,loc,ix,iy,iz)

do i=1,nface
write (7,'(a,27f10.3,a,a,a)') &
     'Shape { geometry IndexedLineSet { coord Coordinate { point [', &
     (x(iface(k,i)),y(iface(k,i)),z(iface(k,i)),k=1,9), &
     ']} coordIndex [8 1 -1 8 2 -1 8 3 -1 8 4 -1 8 5 -1 8 6 -1 8 7 -1 8 0 -1', &
     ']}appearance Appearance { material Material { emissiveColor 0 1 0}}}'
enddo

write (7,'(a)') ']}'
close (7)

return
end

!========================!
!=====[OCTREE_UNION]=====!
!========================!

subroutine octree_union (octree,noctree,octree1,noctree1,octree2,noctree2)

! NOTE: JEAN FOUND A BUG IN THIS ROUTINE ON JULY 4th 2006; HAS NOT BEEN FIXED!!!!
! IT SHOULD NOT BE USED.
! DOUAR HAS BEEN MODIFIED TO PERFORM THIS OPERATION FROM OTHER LOWER LEVEL ROUTINES
! IN OCTREEBITPLUS

! subroutine to calculate the union of two octrees (octree1,octree2)
! and store the result in a third octree (octree)
! it is recommended that the largest of the two octrees be octree1
! (as defined by their size ioctree_size (octree,noctree))
! on exit the two original octrees are left unchanged

implicit none

integer noctree,octree(noctree)
integer noctree1,octree1(noctree1)
integer noctree2,octree2(noctree2)

integer loc,ix,iy,iz

octree=octree1
noctree=noctree1

loc=4
ix=0
iy=0
iz=0
call unite (octree2,noctree2,loc,ix,iy,iz,octree,noctree)

return
end

!=================!
!=====[UNITE]=====!
!=================!

recursive subroutine unite (octree2,noctree2,loc,ix,iy,iz,octree,noctree)

! DO NOT USE
! internal routine
! called by octree_union

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz
integer noctree2,octree2(noctree2)

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer k,idx,idy,idz,ip
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp
double precision xp,yp,zp

level=octree2(loc)
levelmax=octree2(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree2(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)
    ip=2**level
    zp=dfloat(izn)/ip
    yp=dfloat(iyn)/ip
    xp=dfloat(ixn)/ip
    levelin=0
    locp=4
    call update (octree,noctree,xp,yp,zp,level,levelin,locp)
    else
    call unite (octree2,noctree2,locn,ixn,iyn,izn,octree,noctree)
    endif
  enddo
  enddo
  enddo

return
end


!==============================!
!=====[OCTREE_INTERPOLATE]=====!
!==============================!

subroutine octree_interpolate (octree,noctree,icon,nleaves,field,nfield,x,y,z,f)

! This function returns the value of a field (field) known at the nodes
! of an octree by trilinear interpolation

! icon is the connectivity matrix
! nleaves is the number of leaves in the octree
! field is the array of dimension nfield containing the field 
!     known at the nodes of the octree and to be interpolated
! x,y,z are the location of the point where the field is to be interpolated
! f is the resulting interpolated field

implicit none

integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nfield
double precision field(nfield),x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt,f
integer leaf,level,loc,k,iii,jjj,kkk

! function modified by JEAN BRAUN on September 26 2005
! to correct for an error in the logics that led to an interpolation
! from an octree to another identical octree with differences in the
! interpolated function. The reason for this problem was related to
! bad faces or hanging nodes. Indeed, for a hanging node it was very likely
! that the leaf that was detected as the loeaf in which the node resides
! was in fact a leave where the node was a hanging node (ie not one of the
! 4 corner nodes). This meant that the interpolated value was not equal
! to the "constrained" value imposed by the linear constraint at the 
! hanging node. To correct for this we first check if the node can
! be interpolated with r,s,t values that are equal to 1 or -1. If this is
! true than this value is chosen as this would correspond to a nodal value

xt=x
yt=y
zt=z

if (xt.lt.-1.e-11 .or. xt.gt.1.d0+1.d-11) return
if (yt.lt.-1.e-11 .or. yt.gt.1.d0+1.d-11) return
if (zt.lt.-1.e-11 .or. zt.gt.1.d0+1.d-11) return

if (x.lt.1.e-11) xt=1.e-11
if (x.gt.1.d0-1.d-11) xt=1.d0-1.d-11
if (y.lt.1.e-11) yt=1.e-11
if (y.gt.1.d0-1.d-11) yt=1.d0-1.d-11
if (z.lt.1.e-11) zt=1.e-11
if (z.gt.1.d0-1.d-11) zt=1.d0-1.d-11

do kkk=-1,1,2
do jjj=-1,1,2
do iii=-1,1,2

xt=x+iii*1.d-10
yt=y+jjj*1.d-10
zt=z+kkk*1.d-10

if (xt*(xt-1.d0).ge.0d0 .or. yt*(yt-1.d0).ge.0d0 .or. zt*(zt-1.d0).ge.0d0) goto 111

call octree_find_leaf (octree,noctree,xt,yt,zt,leaf,level,loc,x0,y0,z0,dxyz)

r=(x-x0)/dxyz*2.d0-1.d0
s=(y-y0)/dxyz*2.d0-1.d0
t=(z-z0)/dxyz*2.d0-1.d0
h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
phi=0.d0
  do k=1,8
  phi=phi+h(k)*field(icon(k,leaf))
  enddo
f=phi
if (abs(abs(r)-1.d0).lt.1.d-10 .and. abs(abs(s)-1.d0).lt.1.d-10 .and. abs(abs(t)-1.d0).lt.1.d-10) return

111 continue

enddo
enddo
enddo

return
end

!==============================!
!=====[OCTREE_INTERPOLATE]=====!
!==============================!

subroutine octree_interpolate3 (octree,noctree,icon,nleaves,field,nodex,nodey,nodez,nnode,x,y,z,f)

! This function returns the value of a field (field) known at the nodes
! of an octree by trilinear interpolation

! icon is the connectivity matrix
! nleaves is the number of leaves in the octree
! field is the array of dimension nfield containing the field 
!     known at the nodes of the octree and to be interpolated
! x,y,z are the location of the point where the field is to be interpolated
! f is the resulting interpolated field

implicit none

integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nnode,inode
double precision nodex(nnode),nodey(nnode),nodez(nnode)
double precision field(nnode),x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt,f
integer leaf,level,loc,k,iii,jjj,kkk

! function modified by JEAN BRAUN on September 26 2005
! to correct for an error in the logics that led to an interpolation
! from an octree to another identical octree with differences in the
! interpolated function. The reason for this problem was related to
! bad faces or hanging nodes. Indeed, for a hanging node it was very likely
! that the leaf that was detected as the loeaf in which the node resides
! was in fact a leave where the node was a hanging node (ie not one of the
! 4 corner nodes). This meant that the interpolated value was not equal
! to the "constrained" value imposed by the linear constraint at the 
! hanging node. To correct for this we first check if the node can
! be interpolated with r,s,t values that are equal to 1 or -1. If this is
! true than this value is chosen as this would correspond to a nodal value

! function modified by Cedric Thieulot on November 21st 2007
! to correct for an apparent little bug in the case of the interpolation of 
! a given field from an octree onto the same one.
! this modification simply checks whether the point on which we want to interpolate
! already exists in the octree structure. The first three do-loops explore the 8
! points distant from point (x,y,z) by a tiny distance in all three dimensions.
! After having checked that the predicted point falls within the unit cube, 
! the leaf in which the predicted point falls in is found, and we check whether
! any of the nodes of the leaf has the same coordinates as point (x,y,z).
! This slows down a bit the function but also insures and exact interpolation on 
! common nodes of both octrees which are many since they usually at least 
! share the nodes of a level 5 uniform octree, i.e. 35937 nodes.


!=====[CT]=====
do kkk=-1,1,2
   zt=z+kkk*1.d-8
   do jjj=-1,1,2
      yt=y+jjj*1.d-8
      do iii=-1,1,2
         xt=x+iii*1.d-8
         if (xt>0.d0 .and. xt<1.d0 .and. &
             yt>0.d0 .and. yt<1.d0 .and. &
             zt>0.d0 .and. zt<1.d0 ) then
            call octree_find_leaf (octree,noctree,xt,yt,zt,leaf,level,loc,x0,y0,z0,dxyz)
            do k=1,8
               inode=icon(k,leaf)
               if (abs(nodex(inode)-x)<1.d-10 .and. &
                   abs(nodey(inode)-y)<1.d-10 .and. &
                   abs(nodez(inode)-z)<1.d-10 ) then
                  f=field(inode)
                  return
               end if
            end do
         end if
      end do
   end do
end do
!==============


xt=x
yt=y
zt=z

do kkk=-1,1,2
   zt=z+kkk*1.d-10
   do jjj=-1,1,2
      yt=y+jjj*1.d-10
      do iii=-1,1,2
         xt=x+iii*1.d-10
         if (xt*(xt-1.d0).ge.0d0 .or. &
             yt*(yt-1.d0).ge.0d0 .or. &
             zt*(zt-1.d0).ge.0d0) goto 111
         call octree_find_leaf (octree,noctree,xt,yt,zt,leaf,level,loc,x0,y0,z0,dxyz)
         r=(x-x0)/dxyz*2.d0-1.d0
         s=(y-y0)/dxyz*2.d0-1.d0
         t=(z-z0)/dxyz*2.d0-1.d0
         h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
         h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
         h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
         h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
         h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
         h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
         h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
         h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
         phi=0.d0
         do k=1,8
            phi=phi+h(k)*field(icon(k,leaf))
         enddo
         f=phi
         if (abs(r)-1.d0+abs(s)-1.d0+abs(t)-1.d0.lt.1.d-10) return
         111 continue
      enddo
   enddo
enddo

return
end


!===================================!
!=====[OCTREE_INTERPOLATE_MANY]=====!
!===================================!

subroutine octree_interpolate_many (nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
                                    field1,f1,field2,f2,field3,f3, &
                                    field4,f4,field5,f5,field6,f6, &
                                    field7,f7,field8,f8,field9,f9, &
                                    field10,f10,field11,f11,field12,f12, &
                                    field13,f13,field14,f14,field15,f15)
! This function returns the value of several fields (fieldi) known at the nodes
! of an octree by trilinear interpolation

! nf is the number of fields being interpolate (must be comprised between 1 and 15)
! icon is the connectivity matrix
! nleaves is the number of leaves in the octree
! fieldi are the arrays of dimension nfield containing the fields
!     known at the nodes of the octree and to be interpolated
! x,y,z are the location of the point where the fields are to be interpolated
! fi are the resulting interpolated fields

! Note that the number of arguments to this routine depends on the number of
! fields to be interpolated (nf). This is why some of the arguments are declared
! as optional

implicit none

optional :: field2,f2,field3,f3,field4,f4,field5,f5,field6,f6
optional :: field7,f7,field8,f8,field9,f9,field10,f10
optional :: field11,f11,field12,f12,field13,f13,field14,f14,field15,f15

integer noctree,octree(noctree),nleaves,icon(8,nleaves)
integer nfield,nf
double precision field1(nfield),field2(nfield),field3(nfield),field4(nfield), &
                 field5(nfield),field6(nfield),field7(nfield),field8(nfield)
double precision field9(nfield),field10(nfield),field11(nfield),field12(nfield), &
                 field13(nfield),field14(nfield),field15(nfield)
double precision f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15
double precision x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt
integer leaf,level,loc,k,iii,jjj,kkk,ii

! function modified by JEAN BRAUN on September 26 2005
! to correct for an error in the logics that led to an interpolation
! from an octree to another identical octree with differences in the
! interpolated function. The reason for this problem was related to
! bad faces or hanging nodes. Indeed, for a hanging node it was very likely
! that the leaf that was detected as the loeaf in which the node resides
! was in fact a leave where the node was a hanging node (ie not one of the
! 4 corner nodes). This meant that the interpolated value was not equal
! to the "constrained" value imposed by the linear constraint at the
! hanging node. To correct for this we first check if the node can
! be interpolated with r,s,t values that are equal to 1 or -1. If this is
! true than this value is chosen as this would correspond to a nodal value

xt=x
yt=y
zt=z

if (xt.lt.-1.e-11 .or. xt.gt.1.d0+1.d-11) return
if (yt.lt.-1.e-11 .or. yt.gt.1.d0+1.d-11) return
if (zt.lt.-1.e-11 .or. zt.gt.1.d0+1.d-11) return

if (x.lt.1.e-11) xt=1.e-11
if (x.gt.1.d0-1.d-11) xt=1.d0-1.d-11
if (y.lt.1.e-11) yt=1.e-11
if (y.gt.1.d0-1.d-11) yt=1.d0-1.d-11
if (z.lt.1.e-11) zt=1.e-11
if (z.gt.1.d0-1.d-11) zt=1.d0-1.d-11

do kkk=-1,1,2
do jjj=-1,1,2
do iii=-1,1,2

xt=x+iii*1.d-10
yt=y+jjj*1.d-10
zt=z+kkk*1.d-10

if (xt*(xt-1.d0).ge.0d0 .or. yt*(yt-1.d0).ge.0d0 .or. zt*(zt-1.d0).ge.0d0) goto 111

call octree_find_leaf (octree,noctree,xt,yt,zt,leaf,level,loc,x0,y0,z0,dxyz)

r=(x-x0)/dxyz*2.d0-1.d0
s=(y-y0)/dxyz*2.d0-1.d0
t=(z-z0)/dxyz*2.d0-1.d0
h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field1(icon(k,leaf))
    enddo
  f1=phi
if (nf.eq.1) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field2(icon(k,leaf))
    enddo
  f2=phi
if (nf.eq.2) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field3(icon(k,leaf))
    enddo
  f3=phi
if (nf.eq.3) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field4(icon(k,leaf))
    enddo
  f4=phi
if (nf.eq.4) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field5(icon(k,leaf))
    enddo
  f5=phi
if (nf.eq.5) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field6(icon(k,leaf))
    enddo
  f6=phi
if (nf.eq.6) goto 222
  phi=0.d0
    do k=1,8
    phi=phi+h(k)*field7(icon(k,leaf))
    enddo
  f7=phi
if (nf.eq.7) goto 222
  phi=0.d0
    do k=1,8