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
|
! $Header$
! Take a mask where the boundary has been been marked.
! Return information about the normals in these locations.
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
#include "maskvalues.h"
subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
! Arguments
! out: zero for success, nonzero for error
integer :: ierr
! in: array sizes for grid functions
! (you can pass in cctk_lsh(:) for these)
integer :: ni,nj,nk
! in: mask
CCTK_REAL :: mask(ni,nj,nk)
! out: normal directions to use for interpolation
CCTK_REAL :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk)
! Internal variables.
integer i,j,k,ii,jj,kk
CCTK_REAL sx,sy,sz,smag
CCTK_REAL vx,vy,vz
CCTK_REAL dot,dotmax
! Initialise direction arrays to zero.
dirx = 0
diry = 0
dirz = 0
! Loop over all grid points.
do k=2,nk-1
do j=2,nj-1
do i=2,ni-1
! Check if the point is on the boundary, if it
! is not forget about it.
if (mask(i,j,k)==MASK_BOUNDARY) then
! Initialise the vector (sx,sy,sz) to zero.
sx = 0.0D0
sy = 0.0D0
sz = 0.0D0
! Loop around nearest neighbours.
do kk=-1,1
do jj=-1,1
do ii=-1,1
! If a given neighbour is active, then add its
! relative coordinates to (sx,sy,sz).
if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
sx = sx + dble(ii)
sy = sy + dble(jj)
sz = sz + dble(kk)
end if
end do
end do
end do
! Find magnitude of (sx,sy,sz).
smag = sqrt(sx**2 + sy**2 + sz**2)
! Normalize (sx,sy,sz).
if (smag /= 0.0D0) then
sx = sx/smag
sy = sy/smag
sz = sz/smag
end if
! Now we have a unit vector pointing in the average
! direction on the unmasked neighbours. This is
! as close to the normal direction as we can get.
! What we need to do now is find that unmasked
! neighbour that is closest to this direction.
! For this I loop again over unmasked neighbours,
! and find the normalized dot product between
! (sx,sy,sz) and the direction of a given neighbour.
! The largest dot product will give us our best
! normal direction.
dotmax = 0.0D0
dirx(i,j,k) = 0.0D0
diry(i,j,k) = 0.0D0
dirz(i,j,k) = 0.0D0
do kk=-1,1
do jj=-1,1
do ii=-1,1
! If the point is not active, forget about it.
if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
! Find normalized dot product.
vx = dble(ii)
vy = dble(jj)
vz = dble(kk)
dot = (sx*vx + sy*vy + sz*vz) &
/ sqrt(vx**2 + vy**2 + vz**2)
! If the new product is larger than the
! largest so far, redifine the normal.
if (dot > dotmax) then
dotmax = dot
dirx(i,j,k) = vx
diry(i,j,k) = vy
dirz(i,j,k) = vz
end if
end if
end do
end do
end do
! Sanity check. If the check fails then
! something is very wrong.
ii = int(dirx(i,j,k))
jj = int(diry(i,j,k))
kk = int(dirz(i,j,k))
if (mask(i+ii,j+jj,k+kk) /= MASK_ACTIVE) then
call CCTK_WARN (0, "Mask boundary layer is too thick")
end if
end if
end do
end do
end do
! Success
ierr = 0
end subroutine excision_findnormals
|
