-
Notifications
You must be signed in to change notification settings - Fork 326
/
Div.jl
executable file
·126 lines (104 loc) · 3.78 KB
/
Div.jl
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
function Div(Px, Py, bound = "sym", order = 1)
"""
div - divergence operator
fd = div(Px,Py, options);
fd = div(P, options);
options.bound = 'per' or 'sym'
options.order = 1 (backward differences)
= 2 (centered differences)
Note that the -div and grad operator are adjoint
of each other such that
<grad(f),g>=<f,-div(g)>
See also: grad.
Copyright (c) 2007 Gabriel Peyre
"""
# retrieve number of dimensions
nbdims = ndims(Px)
if nbdims >= 3
if nbdims == 3
Py = P[:, :, 2]
Px = P[:, :, 1]
nbdims = 2
else
Pz = P[:, :, :, 3]
Py = P[:, :, :, 2]
Px = P[:, :, :, 1]
nbdims = 3
end
end
if bound == "sym"
nx = size(Px)[1]
if order == 1
fx = Px .- Px[vcat(1, collect(1, nx - 1)), :]
fx[1, :] = Px[1, :] # boundary
fx[nx, :] = - Px[nx - 1, :]
if nbdims >= 2
ny = size(Py)[2]
fy = Py - Py[:, vcat(1, collect(1, ny - 1))]
fy[:, 1] = Py[:, 1] # boundary
fy[:, ny] = - Py[:, ny - 1]
end
if nbdims >= 3
nz = size(Pz)[3]
fz = Pz - Pz[:, :, vcat(1, collect(1 : nz - 1))]
fz[:, :, 1] = Pz[:, :, 1] # boundary
fz[:, :, nz] = - Pz[:, :, nz - 1]
end
else
fx = (Px[vcat(collect(2 : nx), nx), :] - Px[vcat([1], collect(1 : nx-1)), :])./2.
fx[1, :] = + Px[2, :]./2. + Px[1, :] # boundary
fx[2, :] = + Px[3, :]./2. - Px[1, :]
fx[nx, :] = - Px[nx, :] - Px[nx - 1, :]./2.
fx[nx - 1, :] = + Px[nx, :] - Px[nx - 2, :]./2.
if nbdims >= 2
ny = size(Py)[2]
fy = (Py[:, vcat(collect(2 : ny), ny)] - Py[:, vcat(1, collect(1 : ny-1))])./2.
fy[:, 1] = + Py[:, 2]./2. + Py[:, 1] # boundary
fy[:, 2] = + Py[:, 3]./2. - Py[:,1]
fy[:, ny] = - Py[:, ny] - Py[:, ny-1]./2.
fy[:, ny - 1] = + Py[:, ny] - Py[:, ny - 2]./2.
end
if nbdims >= 3
nz = size(Pz)[3]
fz = (Pz[:, :, vcat(collect(2 : nz), nz)] - Pz[:, :, vcat(1, collect(1 : nz - 1))])./2.
fz[:, :, 1] = + Pz[:, :, 2]./2. + Pz[:, :, 1] # boundary
fz[:, :, 2] = + Pz[:, :, 3]./2. - Pz[:, :, 1]
fz[:, :, ny] = - Pz[:, :, nz] - Pz[:, :, nz - 1]./2.
fz[:, :, ny - 1] = + Pz[:, :, nz] - Pz[:, :, nz - 2]./2.
end
end
else
if order == 1
nx = size(Px)[1]
fx = Px - Px[vcat(nx, collect(1 : nx-1)), :]
if nbdims >= 2
ny = size(Py)[2]
fy = Py - Py[:, vcat(ny, collect(1 : ny-1))]
end
if nbdims>=3
nz = size(Pz)[3]
fz = Pz - Pz[:, :, vcat(nz, collect(1 : nz-1))]
end
else
nx = size(Px)[1]
fx = (Px[vcat(collect(2 : nx), 1), :]) - (Px[vcat(nx, collect(1 : nx - 1)), :])
if nbdims >= 2
ny = size(Py)[2]
fy = (Py[:, vcat(collect(2 : ny), 1)]) - (Py[:, vcat(ny, collect(1 : ny - 1))])
end
if nbdims >= 3
nz = size(Pz)[3]
fz = (Pz[:, :, vcat(collect(2 : nz), 1)]) - (Pz[:, :, vcat(nz, collect(1 : nz - 1))])
end
end
end
# gather result
if nbdims == 3
fd = fx + fy + fz
elseif nbdims == 2
fd = fx + fy
else
fd = fx
end
return fd
end