Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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
if (nf.ne.3) write(*,*) 'error: octree_interpolate_three, nf = ', nf
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
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