From c20d275cf681b0d25da32f59d4638e9218320cf2 Mon Sep 17 00:00:00 2001
From: Douglas Guptill <douglas.guptill@dal.ca>
Date: Mon, 20 Jul 2009 01:16:12 +0000
Subject: [PATCH] try 2 specialized octree_interpolate_ subroutines
---
OCTREE/Makefile | 2 +
OCTREE/octree_interpolate_three.f90 | 116 +++++++++
.../octree_interpolate_three_derivative.f90 | 220 ++++++++++++++++++
move_surface.f90 | 7 +-
4 files changed, 343 insertions(+), 2 deletions(-)
create mode 100644 OCTREE/octree_interpolate_three.f90
create mode 100644 OCTREE/octree_interpolate_three_derivative.f90
diff --git a/OCTREE/Makefile b/OCTREE/Makefile
index 49aa4396..1fc13da2 100644
--- a/OCTREE/Makefile
+++ b/OCTREE/Makefile
@@ -9,6 +9,8 @@ graph_sloan.o \
isorti_sloan.o \
label_sloan.o \
number_sloan.o \
+octree_interpolate_three.o \
+octree_interpolate_three_derivative.o \
profil_sloan.o \
rootls_sloan.o
diff --git a/OCTREE/octree_interpolate_three.f90 b/OCTREE/octree_interpolate_three.f90
new file mode 100644
index 00000000..b037413e
--- /dev/null
+++ b/OCTREE/octree_interpolate_three.f90
@@ -0,0 +1,116 @@
+subroutine octree_interpolate_three(nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
+ field1,f1,field2,f2,field3,f3)
+
+! octree_interpolate_many specialized for 3
+! optional parameters removed.
+! Author: Douglas Guptill
+! Date: 2009-07-19
+
+! This function returns the value of three fields (fieldi) known at the nodes
+! of an octree by trilinear interpolation
+
+! nf is the number of fields being interpolate (must be 3)
+! 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
+
+implicit none
+
+
+integer noctree,octree(noctree),nleaves,icon(8,nleaves)
+integer nfield,nf
+double precision field1(nfield),field2(nfield),field3(nfield)
+double precision f1,f2,f3
+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
+
+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
+
diff --git a/OCTREE/octree_interpolate_three_derivative.f90 b/OCTREE/octree_interpolate_three_derivative.f90
new file mode 100644
index 00000000..5556ada7
--- /dev/null
+++ b/OCTREE/octree_interpolate_three_derivative.f90
@@ -0,0 +1,220 @@
+subroutine octree_interpolate_three_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)
+
+! octree_interpolate_many_derivative specialized for 3
+! optional parameters removed.
+! Author: Douglas Guptill
+! Date: 2009-07-19
+
+! 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 3)
+! 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
+
+
+implicit none
+
+integer noctree,octree(noctree),nleaves,icon(8,nleaves)
+integer nfield,nf
+double precision field1(nfield),field2(nfield),field3(nfield)
+double precision f1,f2,f3
+double precision df1x,df2x,df3x
+double precision df1y,df2y,df3y
+double precision df1z,df2z,df3z
+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
+
+
+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
+
diff --git a/move_surface.f90 b/move_surface.f90
index 7e31fa14..5762339d 100644
--- a/move_surface.f90
+++ b/move_surface.f90
@@ -125,7 +125,7 @@ do i=iproc+1,surface%nsurface,nproc
xi=(x0+x)/2.d0
yi=(y0+y)/2.d0
zi=(z0+z)/2.d0
- call octree_interpolate_many (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
+ call octree_interpolate_three (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi,ov%unode,u,ov%vnode,v,ov%wnode,w)
if (i==1 .or. i==nside .or. i==2*nside-1 .or. i==3*nside-2) then !corner points
@@ -152,7 +152,7 @@ do i=iproc+1,surface%nsurface,nproc
ybuf(i)=y
zbuf(i)=z
if (params%normaladvect) then
- call octree_interpolate_many_derivative (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
+ call octree_interpolate_three_derivative (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi, &
ov%unode,u,ux,uy,uz, &
ov%vnode,v,vx,vy,vz, &
@@ -180,6 +180,9 @@ do i=iproc+1,surface%nsurface,nproc
xnbuf(i)=xn0+(yn1*zln2-zn1*yln2+yln1*zn2-zln1*yn2)*dt+(yln1*zln2-zln1*yln2)*dt**2
ynbuf(i)=yn0+(zn1*xln2-xn1*zln2+zln1*xn2-xln1*zn2)*dt+(zln1*xln2-xln1*zln2)*dt**2
znbuf(i)=zn0+(xn1*yln2-yn1*xln2+xln1*yn2-yln1*xn2)*dt+(xln1*yln2-yln1*xln2)*dt**2
+ if (xnbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', xnbuf is large: ', xnbuf(i)
+ if (ynbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', ynbuf is large: ', ynbuf(i)
+ if (znbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', znbuf is large: ', znbuf(i)
xyzn=sqrt(xnbuf(i)**2+ynbuf(i)**2+znbuf(i)**2)
xnbuf(i)=xnbuf(i)/xyzn
ynbuf(i)=ynbuf(i)/xyzn
--
GitLab