!-------------------------------------------
! 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
phi=phi+h(k)*field8(icon(k,leaf))
enddo
f8=phi
if (nf.eq.8) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field9(icon(k,leaf))
enddo
f9=phi
if (nf.eq.9) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field10(icon(k,leaf))
enddo
f10=phi
if (nf.eq.10) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field11(icon(k,leaf))
enddo
f11=phi
if (nf.eq.11) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field12(icon(k,leaf))
enddo
f12=phi
if (nf.eq.12) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field13(icon(k,leaf))
enddo
f13=phi
if (nf.eq.13) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field14(icon(k,leaf))
enddo
f14=phi
if (nf.eq.14) goto 222
phi=0.d0
do k=1,8
phi=phi+h(k)*field15(icon(k,leaf))
enddo
f15=phi
if (nf.gt.15) stop 'too many field'
222 continue
!if (abs(r)-1.d0+abs(s)-1.d0+abs(t)-1.d0.lt.1.d-10) return
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_MANY_DERIVATIVE]=====!
!==============================================!
subroutine octree_interpolate_many_derivative &
(nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
field1,f1,df1x,df1y,df1z, &
field2,f2,df2x,df2y,df2z, &
field3,f3,df3x,df3y,df3z, &
field4,f4,df4x,df4y,df4z, &
field5,f5,df5x,df5y,df5z, &
field6,f6,df6x,df6y,df6z, &
field7,f7,df7x,df7y,df7z, &
field8,f8,df8x,df8y,df8z, &
field9,f9,df9x,df9y,df9z, &
field10,f10,df10x,df10y,df10z, &
field11,f11,df11x,df11y,df11z, &
field12,f12,df12x,df12y,df12z, &
field13,f13,df13x,df13y,df13z, &
field14,f14,df14x,df14y,df14z, &
field15,f15,df15x,df15y,df15z)
! This function returns the value of several fields (fieldi) known at the nodes
! of an octree by trilinear interpolation as well as their 3 spatial derivatives
! 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 oarguments 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
optional :: df2x,df3x,df4x,df5x,df6x,df7x,df8x,df9x
optional :: df10x,df11x,df12x,df13x,df14x,df15x
optional :: df2y,df3y,df4y,df5y,df6y,df7y,df8y,df9y
optional :: df10y,df11y,df12y,df13y,df14y,df15y
optional :: df2z,df3z,df4z,df5z,df6z,df7z,df8z,df9z
optional :: df10z,df11z,df12z,df13z,df14z,df15z
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 df1x,df2x,df3x,df4x,df5x,df6x,df7x,df8x,df9x
double precision df10x,df11x,df12x,df13x,df14x,df15x
double precision df1y,df2y,df3y,df4y,df5y,df6y,df7y,df8y,df9y
double precision df10y,df11y,df12y,df13y,df14y,df15y
double precision df1z,df2z,df3z,df4z,df5z,df6z,df7z,df8z,df9z
double precision df10z,df11z,df12z,df13z,df14z,df15z
double precision x,y,z,x0,y0,z0,dxyz,r,s,t,h(8),phi,xt,yt,zt
double precision phix,phiy,phiz
integer leaf,level,loc,k,iii,jjj,kkk,ii,jj,kk,ic,ijk
double precision dhdr(8),dhds(8),dhdt(8),xg(8),yg(8),zg(8),jcb(3,3),jcbi(3,3),volume
double precision dhdx(8),dhdy(8),dhdz(8)
! 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
dhdr(1)=-(1.d0-s)*(1.d0-t)/8.d0
dhdr(2)=(1.d0-s)*(1.d0-t)/8.d0
dhdr(3)=-(1.d0+s)*(1.d0-t)/8.d0
dhdr(4)=(1.d0+s)*(1.d0-t)/8.d0
dhdr(5)=-(1.d0-s)*(1.d0+t)/8.d0
dhdr(6)=(1.d0-s)*(1.d0+t)/8.d0
dhdr(7)=-(1.d0+s)*(1.d0+t)/8.d0
dhdr(8)=(1.d0+s)*(1.d0+t)/8.d0
dhds(1)=-(1.d0-r)*(1.d0-t)/8.d0
dhds(2)=-(1.d0+r)*(1.d0-t)/8.d0
dhds(3)=(1.d0-r)*(1.d0-t)/8.d0
dhds(4)=(1.d0+r)*(1.d0-t)/8.d0
dhds(5)=-(1.d0-r)*(1.d0+t)/8.d0
dhds(6)=-(1.d0+r)*(1.d0+t)/8.d0
dhds(7)=(1.d0-r)*(1.d0+t)/8.d0
dhds(8)=(1.d0+r)*(1.d0+t)/8.d0
dhdt(1)=-(1.d0-r)*(1.d0-s)/8.d0
dhdt(2)=-(1.d0+r)*(1.d0-s)/8.d0
dhdt(3)=-(1.d0-r)*(1.d0+s)/8.d0
dhdt(4)=-(1.d0+r)*(1.d0+s)/8.d0
dhdt(5)=(1.d0-r)*(1.d0-s)/8.d0
dhdt(6)=(1.d0+r)*(1.d0-s)/8.d0
dhdt(7)=(1.d0-r)*(1.d0+s)/8.d0
dhdt(8)=(1.d0+r)*(1.d0+s)/8.d0
ijk=0
do kk=0,1
do jj=0,1
do ii=0,1
ijk=ijk+1
xg(ijk)=x0+ii*dxyz
yg(ijk)=y0+jj*dxyz
zg(ijk)=z0+kk*dxyz
enddo
enddo
enddo
jcb=0.
do k=1,8
jcb(1,1)=jcb(1,1)+dhdr(k)*xg(k)
jcb(1,2)=jcb(1,2)+dhdr(k)*yg(k)
jcb(1,3)=jcb(1,3)+dhdr(k)*zg(k)
jcb(2,1)=jcb(2,1)+dhds(k)*xg(k)
jcb(2,2)=jcb(2,2)+dhds(k)*yg(k)
jcb(2,3)=jcb(2,3)+dhds(k)*zg(k)
jcb(3,1)=jcb(3,1)+dhdt(k)*xg(k)
jcb(3,2)=jcb(3,2)+dhdt(k)*yg(k)
jcb(3,3)=jcb(3,3)+dhdt(k)*zg(k)
enddo
volume=jcb(1,1)*jcb(2,2)*jcb(3,3)+jcb(1,2)*jcb(2,3)*jcb(3,1) &
+jcb(2,1)*jcb(3,2)*jcb(1,3) &
-jcb(1,3)*jcb(2,2)*jcb(3,1)-jcb(1,2)*jcb(2,1)*jcb(3,3) &
-jcb(2,3)*jcb(3,2)*jcb(1,1)
jcbi(1,1)=(jcb(2,2)*jcb(3,3)-jcb(2,3)*jcb(3,2))/volume
jcbi(2,1)=(jcb(2,3)*jcb(3,1)-jcb(2,1)*jcb(3,3))/volume
jcbi(3,1)=(jcb(2,1)*jcb(3,2)-jcb(2,2)*jcb(3,1))/volume
jcbi(1,2)=(jcb(1,3)*jcb(3,2)-jcb(1,2)*jcb(3,3))/volume
jcbi(2,2)=(jcb(1,1)*jcb(3,3)-jcb(1,3)*jcb(3,1))/volume
jcbi(3,2)=(jcb(1,2)*jcb(3,1)-jcb(1,1)*jcb(3,2))/volume
jcbi(1,3)=(jcb(1,2)*jcb(2,3)-jcb(1,3)*jcb(2,2))/volume
jcbi(2,3)=(jcb(1,3)*jcb(2,1)-jcb(1,1)*jcb(2,3))/volume
jcbi(3,3)=(jcb(1,1)*jcb(2,2)-jcb(1,2)*jcb(2,1))/volume
do k=1,8
dhdx(k)=jcbi(1,1)*dhdr(k)+jcbi(1,2)*dhds(k)+jcbi(1,3)*dhdt(k)
dhdy(k)=jcbi(2,1)*dhdr(k)+jcbi(2,2)*dhds(k)+jcbi(2,3)*dhdt(k)
dhdz(k)=jcbi(3,1)*dhdr(k)+jcbi(3,2)*dhds(k)+jcbi(3,3)*dhdt(k)
enddo
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field1(ic)
phix=phix+dhdx(k)*field1(ic)
phiy=phiy+dhdy(k)*field1(ic)
phiz=phiz+dhdz(k)*field1(ic)
enddo
f1=phi
df1x=phix
df1y=phiy
df1z=phiz
if (nf.eq.1) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field2(ic)
phix=phix+dhdx(k)*field2(ic)
phiy=phiy+dhdy(k)*field2(ic)
phiz=phiz+dhdz(k)*field2(ic)
enddo
f2=phi
df2x=phix
df2y=phiy
df2z=phiz
if (nf.eq.2) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field3(ic)
phix=phix+dhdx(k)*field3(ic)
phiy=phiy+dhdy(k)*field3(ic)
phiz=phiz+dhdz(k)*field3(ic)
enddo
f3=phi
df3x=phix
df3y=phiy
df3z=phiz
if (nf.eq.3) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field4(ic)
phix=phix+dhdx(k)*field4(ic)
phiy=phiy+dhdy(k)*field4(ic)
phiz=phiz+dhdz(k)*field4(ic)
enddo
f4=phi
df4x=phix
df4y=phiy
df4z=phiz
if (nf.eq.4) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field5(ic)
phix=phix+dhdx(k)*field5(ic)
phiy=phiy+dhdy(k)*field5(ic)
phiz=phiz+dhdz(k)*field5(ic)
enddo
f5=phi
df5x=phix
df5y=phiy
df5z=phiz
if (nf.eq.5) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field6(ic)
phix=phix+dhdx(k)*field6(ic)
phiy=phiy+dhdy(k)*field6(ic)
phiz=phiz+dhdz(k)*field6(ic)
enddo
f6=phi
df6x=phix
df6y=phiy
df6z=phiz
if (nf.eq.6) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field7(ic)
phix=phix+dhdx(k)*field7(ic)
phiy=phiy+dhdy(k)*field7(ic)
phiz=phiz+dhdz(k)*field7(ic)
enddo
f7=phi
df7x=phix
df7y=phiy
df7z=phiz
if (nf.eq.7) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field8(ic)
phix=phix+dhdx(k)*field8(ic)
phiy=phiy+dhdy(k)*field8(ic)
phiz=phiz+dhdz(k)*field8(ic)
enddo
f8=phi
df8x=phix
df8y=phiy
df8z=phiz
if (nf.eq.8) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field9(ic)
phix=phix+dhdx(k)*field9(ic)
phiy=phiy+dhdy(k)*field9(ic)
phiz=phiz+dhdz(k)*field9(ic)
enddo
f9=phi
df9x=phix
df9y=phiy
df9z=phiz
if (nf.eq.9) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field10(ic)
phix=phix+dhdx(k)*field10(ic)
phiy=phiy+dhdy(k)*field10(ic)
phiz=phiz+dhdz(k)*field10(ic)
enddo
f10=phi
df10x=phix
df10y=phiy
df10z=phiz
if (nf.eq.10) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field11(ic)
phix=phix+dhdx(k)*field11(ic)
phiy=phiy+dhdy(k)*field11(ic)
phiz=phiz+dhdz(k)*field11(ic)
enddo
f11=phi
df11x=phix
df11y=phiy
df11z=phiz
if (nf.eq.11) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field12(ic)
phix=phix+dhdx(k)*field12(ic)
phiy=phiy+dhdy(k)*field12(ic)
phiz=phiz+dhdz(k)*field12(ic)
enddo
f12=phi
df12x=phix
df12y=phiy
df12z=phiz
if (nf.eq.12) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field13(ic)
phix=phix+dhdx(k)*field13(ic)
phiy=phiy+dhdy(k)*field13(ic)
phiz=phiz+dhdz(k)*field13(ic)
enddo
f13=phi
df13x=phix
df13y=phiy
df13z=phiz
if (nf.eq.13) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field14(ic)
phix=phix+dhdx(k)*field14(ic)
phiy=phiy+dhdy(k)*field14(ic)
phiz=phiz+dhdz(k)*field14(ic)
enddo
f14=phi
df14x=phix
df14y=phiy
df14z=phiz
if (nf.eq.14) goto 222
phi=0.d0
phix=0.d0
phiy=0.d0
phiz=0.d0
do k=1,8
ic=icon(k,leaf)
phi=phi+h(k)*field15(ic)
phix=phix+dhdx(k)*field15(ic)
phiy=phiy+dhdy(k)*field15(ic)
phiz=phiz+dhdz(k)*field15(ic)
enddo
f15=phi
df15x=phix
df15y=phiy
df15z=phiz
if (nf.gt.15) stop 'too many field'
222 continue
!if (abs(r)-1.d0+abs(s)-1.d0+abs(t)-1.d0.lt.1.d-10) return
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
!-------------------------------------
Subroutine mrgrnk (XDONT, IRNGT, np)
!Subroutine mrgrnk (XDONT, IRNGT)
! __________________________________________________________
! MRGRNK = Merge-sort ranking of an array
! For performance reasons, the first 2 passes are taken
! out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
Integer XDONT(np),IRNGT(np)
! Integer, Dimension (:), Intent (In) :: XDONT
! Integer, Dimension (:), Intent (Out) :: IRNGT
! __________________________________________________________
Integer :: XVALA, XVALB
!
Integer, Dimension(:),allocatable :: JWRKT
! Integer, Dimension (SIZE(IRNGT)) :: JWRKT
Integer :: LMTNA, LMTNC, IRNG1, IRNG2
Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
Select Case (NVAL)
Case (:0)
Return
Case (1)
IRNGT (1) = 1
Return
Case Default
Continue
End Select
allocate (JWRKT(np))
!
! Fill-in the index array, creating ordered couples
!
Do IIND = 2, NVAL, 2
If (XDONT(IIND-1) <= XDONT(IIND)) Then
IRNGT (IIND-1) = IIND - 1
IRNGT (IIND) = IIND
Else
IRNGT (IIND-1) = IIND
IRNGT (IIND) = IIND - 1
End If
End Do
If (Modulo(NVAL, 2) /= 0) Then
IRNGT (NVAL) = NVAL
End If
!
! We will now have ordered subsets A - B - A - B - ...
! and merge A and B couples into C - C - ...
!
LMTNA = 2
LMTNC = 4
!
! First iteration. The length of the ordered subsets goes from 2 to 4
!
Do
If (NVAL <= 2) Exit
!
! Loop on merges of A and B into C
!
Do IWRKD = 0, NVAL - 1, 4
If ((IWRKD+4) > NVAL) Then
If ((IWRKD+2) >= NVAL) Exit
!
! 1 2 3
!
If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
!
! 1 3 2
!
If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
IRNG2 = IRNGT (IWRKD+2)
IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
IRNGT (IWRKD+3) = IRNG2
!
! 3 1 2
!
Else
IRNG1 = IRNGT (IWRKD+1)
IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
IRNGT (IWRKD+2) = IRNG1
End If
Exit
End If
!
! 1 2 3 4
!
If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
!
! 1 3 x x
!
If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
IRNG2 = IRNGT (IWRKD+2)
IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
! 1 3 2 4
IRNGT (IWRKD+3) = IRNG2
Else
! 1 3 4 2
IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
IRNGT (IWRKD+4) = IRNG2
End If
!
! 3 x x x
!
Else
IRNG1 = IRNGT (IWRKD+1)
IRNG2 = IRNGT (IWRKD+2)
IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
IRNGT (IWRKD+2) = IRNG1
If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
! 3 1 2 4
IRNGT (IWRKD+3) = IRNG2
Else
! 3 1 4 2
IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
IRNGT (IWRKD+4) = IRNG2
End If
Else
! 3 4 1 2
IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
IRNGT (IWRKD+3) = IRNG1
IRNGT (IWRKD+4) = IRNG2
End If
End If
End Do
!
! The Cs become As and Bs
!
LMTNA = 4
Exit
End Do
!
! Iteration loop. Each time, the length of the ordered subsets
! is doubled.
!
Do
If (LMTNA >= NVAL) Exit
IWRKF = 0
LMTNC = 2 * LMTNC
!
! Loop on merges of A and B into C
!
Do
IWRK = IWRKF
IWRKD = IWRKF + 1
JINDA = IWRKF + LMTNA
IWRKF = IWRKF + LMTNC
If (IWRKF >= NVAL) Then
If (JINDA >= NVAL) Exit
IWRKF = NVAL
End If
IINDA = 1
IINDB = JINDA + 1
!
! Shortcut for the case when the max of A is smaller
! than the min of B. This line may be activated when the
! initial set is already close to sorted.
!
! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
!
! One steps in the C subset, that we build in the final rank array
!
! Make a copy of the rank array for the merge iteration
!
JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
!
XVALA = XDONT (JWRKT(IINDA))
XVALB = XDONT (IRNGT(IINDB))
!
Do
IWRK = IWRK + 1
!
! We still have unprocessed values in both A and B
!
If (XVALA > XVALB) Then
IRNGT (IWRK) = IRNGT (IINDB)
IINDB = IINDB + 1
If (IINDB > IWRKF) Then
! Only A still with unprocessed values
IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
Exit
End If
XVALB = XDONT (IRNGT(IINDB))
Else
IRNGT (IWRK) = JWRKT (IINDA)
IINDA = IINDA + 1
If (IINDA > LMTNA) Exit! Only B still with unprocessed values
XVALA = XDONT (JWRKT(IINDA))
End If
!
End Do
End Do
!
! The Cs become As and Bs
!
LMTNA = 2 * LMTNA
End Do
!
deallocate (JWRKT)
Return
!
End Subroutine mrgrnk