OctreeBitPlus.f90

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. &