diff --git a/OCTREE/octree_interpolate_five.f90 b/OCTREE/octree_interpolate_five.f90 new file mode 100644 index 0000000000000000000000000000000000000000..021cba63e6ed3a89f35558afde83a3816447f538 --- /dev/null +++ b/OCTREE/octree_interpolate_five.f90 @@ -0,0 +1,132 @@ +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 +