another specialization of octree_interpolate_many

OCTREE/octree_interpolate_five.f90

+subroutine octree_interpolate_five(nf,octree,noctree,icon,nleaves,nfield,x,y,z, &
+     field1,f1,field2,f2,field3,f3,field4,f4,field5,f5)
+
+! octree_interpolate_many specialized for 5
+! optional parameters removed.
+! Author: Douglas Guptill
+! Date: 2009-07-20
+
+! This function returns the value of five fields (fieldi) known at the nodes
+! of an octree by trilinear interpolation
+
+! nf is the number of fields being interpolate (must be 5)
+! 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 field4(nfield),field5(nfield)
+double precision f1,f2,f3,f4,f5
+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
+
+if (nf.ne.5) write(*,*) 'error: octree_interpolate_five, nf = ', nf
+
+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
+
+         phi=0.d0
+         do k=1,8
+            phi=phi+h(k)*field2(icon(k,leaf))
+         enddo
+         f2=phi
+
+         phi=0.d0
+         do k=1,8
+            phi=phi+h(k)*field3(icon(k,leaf))
+         enddo
+         f3=phi
+
+         phi=0.d0
+         do k=1,8
+            phi=phi+h(k)*field4(icon(k,leaf))
+         enddo
+         f4=phi
+
+         phi=0.d0
+         do k=1,8
+            phi=phi+h(k)*field5(icon(k,leaf))
+         enddo
+         f5=phi
+
+
+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
+