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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
|
! $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
integer :: ii,jj,kk
! 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)
! distinguish between faces, edges, and corners.
! look at the excised points, and distinguish the following cases,
! in this order:
! [trivial case]
! 1: if (none)
! case 0: interior
! [box shape cases]
! 2: if (only one corner)
! case I: corner
! 3: if (only points on one edge)
! case II: edge
! 4: if (only points on one face)
! case III: face
! 5: if (only points on two adjacent faces)
! case IIa: anti-edge
! 6: if (one points on three adjacent faces)
! case Ia: anti-corner
! [illegal cases]
! error
logical mask_interior (-1:+1,-1:+1,-1:+1, 1)
logical mask_corner (-1:+1,-1:+1,-1:+1, 8)
logical mask_edge (-1:+1,-1:+1,-1:+1, 12)
logical mask_face (-1:+1,-1:+1,-1:+1, 6)
logical mask_antiface (-1:+1,-1:+1,-1:+1, 6)
logical mask_antiedge (-1:+1,-1:+1,-1:+1, 12)
logical mask_anticorner(-1:+1,-1:+1,-1:+1, 8)
integer dir_interior (3, 1), num_interior
integer dir_corner (3, 8), num_corner
integer dir_edge (3, 12), num_edge
integer dir_face (3, 6), num_face
integer dir_antiface (3, 6), num_antiface
integer dir_antiedge (3, 12), num_antiedge
integer dir_anticorner (3, 8), num_anticorner
logical msk(-1:+1,-1:+1,-1:+1)
integer di,dj,dk
integer i,j,k
integer n
integer dir(3)
character(1000) message
!
! Create templates
!
num_interior = 0
num_corner = 0
num_edge = 0
num_face = 0
num_antiedge = 0
num_anticorner = 0
do dk=-1,+1
do dj=-1,+1
do di=-1,+1
select case (count( (/di,dj,dk/) /= 0 ))
case (0)
! interior (or error, which is not explicitely checked)
num_interior = num_interior + 1
dir_interior(:,num_interior) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_interior(i,j,k,num_interior) = &
any( (/di,dj,dk/)/=0 &
.and. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
case (1)
! face (or anti-face, which is the same thing)
num_face = num_face + 1
dir_face(:,num_face) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_face(i,j,k,num_face) = &
any( (/di,dj,dk/)/=0 &
.and. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
case (2)
! edge or anti-edge
num_edge = num_edge + 1
dir_edge(:,num_edge) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_edge(i,j,k,num_edge) = &
all( (/di,dj,dk/)==0 &
.or. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
num_antiedge = num_antiedge + 1
dir_antiedge(:,num_antiedge) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_antiedge(i,j,k,num_antiedge) = &
any( (/di,dj,dk/)/=0 &
.and. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
case (3)
! corner or anti-corner
num_corner = num_corner + 1
dir_corner(:,num_corner) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_corner(i,j,k,num_corner) = &
all( (/di,dj,dk/)==0 &
.or. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
num_anticorner = num_anticorner + 1
dir_anticorner(:,num_anticorner) = (/di,dj,dk/)
do k=-1,+1
do j=-1,+1
do i=-1,+1
mask_anticorner(i,j,k,num_anticorner) = &
any( (/di,dj,dk/)/=0 &
.and. (/i,j,k/)==(/di,dj,dk/) )
end do
end do
end do
case default
call CCTK_WARN (0, "internal error")
end select
end do
end do
end do
! sanity check
if ( num_interior /= 1 &
.or. num_corner /= 8 &
.or. num_edge /= 12 &
.or. num_face /= 6 &
.or. num_antiedge /= 12 &
.or. num_anticorner /= 8 ) then
call CCTK_WARN (0, "internal error")
end if
!
! Calculate normal directions
!
dirx = 0
diry = 0
dirz = 0
do k=2,nk-1
do j=2,nj-1
do i=2,ni-1
if (mask(i,j,k)==MASK_BOUNDARY) then
! make a snapshot
msk = mask(i-1:i+1,j-1:j+1,k-1:k+1) == MASK_ACTIVE
! for a template to match, all grid points have to be
! either active, or set in the template.
search: do
! 1: case 0: interior
do n=1,num_interior
if (all(msk .eqv. mask_interior(:,:,:,n))) then
dir = dir_interior(:,n)
exit search
end if
end do
! 2: case I: corner
do n=1,num_corner
if (all(msk .eqv. mask_corner(:,:,:,n))) then
dir = dir_corner(:,n)
exit search
end if
end do
! 3: case II: edge
do n=1,num_edge
if (all(msk .eqv. mask_edge(:,:,:,n))) then
dir = dir_edge(:,n)
exit search
end if
end do
! 4: case III: face
do n=1,num_face
if (all(msk .eqv. mask_face(:,:,:,n))) then
dir = dir_face(:,n)
exit search
end if
end do
! 5: case IIa: antiedge
do n=1,num_antiedge
if (all(msk .eqv. mask_antiedge(:,:,:,n))) then
dir = dir_antiedge(:,n)
exit search
end if
end do
! 6: case Ia: anticorner
do n=1,num_anticorner
if (all(msk .eqv. mask_anticorner(:,:,:,n))) then
dir = dir_anticorner(:,n)
exit search
end if
end do
! error: unclassified
print *,i,j,k,dir(1),dir(2),dir(3)
do kk=k-1,k+1
do jj=j-1,j+1
do ii=i-1,i+1
print*,ii-i,jj-j,kk-k,msk(ii-i,jj-j,kk-k)
end do
end do
end do
call CCTK_WARN (0, "Mask is not convex")
end do search
! dir(:) = -dir(:)
if (mask(i+dir(1), j+dir(2), k+dir(3)) /= MASK_ACTIVE) then
print *,i,j,k,dir(1),dir(2),dir(3)
do kk=k-1,k+1
do jj=j-1,j+1
do ii=i-1,i+1
print*,ii-i,jj-j,kk-k,msk(ii-i,jj-j,kk-k)
end do
end do
end do
call CCTK_WARN (0, "Mask boundary layer is too thick")
end if
dirx(i,j,k) = dir(1)
diry(i,j,k) = dir(2)
dirz(i,j,k) = dir(3)
end if
end do
end do
end do
! Success
ierr = 0
end subroutine excision_findnormals
|
