try 2 specialized octree_interpolate_ subroutines

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
 

OCTREE/octree_interpolate_three.f90

+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
+

OCTREE/octree_interpolate_three_derivative.f90

+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
+

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