From 9dc2b759e766024564a7d20c1d316569bcbc8d43 Mon Sep 17 00:00:00 2001 From: Douglas Guptill <douglas.guptill@dal.ca> Date: Mon, 20 Jul 2009 04:41:37 +0000 Subject: [PATCH] another specialization of octree_interpolate_many --- OCTREE/octree_interpolate_five.f90 | 132 +++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 OCTREE/octree_interpolate_five.f90 diff --git a/OCTREE/octree_interpolate_five.f90 b/OCTREE/octree_interpolate_five.f90 new file mode 100644 index 00000000..021cba63 --- /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 + -- GitLab