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
58
59
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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
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