forked from JuliaLang/julia
-
Notifications
You must be signed in to change notification settings - Fork 0
/
getdivgrad.jl
48 lines (36 loc) · 1.16 KB
/
getdivgrad.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
# This file is a part of Julia. License is MIT: http://julialang.org/license
# https://github.com/JuliaLang/julia/issues/4707
#----------------- Get the A matrix
function getDivGrad(n1,n2,n3)
# the Divergence
D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))
# DIV from faces to cell-centers
Div = [D1 D2 D3]
return Div*Div';
end
#----------------- 1D finite difference on staggered grid
function ddx(n)
# generate 1D derivatives
return d = spdiags(ones(n)*[-1 1],[0,1],n,n+1)
end
#------------- Build a diagonal matrix
function spdiags(B,d,m,n)
# spdiags(B,d,m,n)
# creates a sparse matrix from its diagonals
d = d[:]
p = length(d)
len = zeros(Int, p+1, 1)
for k = 1:p
len[k+1] = len[k] + length(max(1,1-d[k]):min(m,n-d[k]))
end
a = zeros(Int, len[p+1], 3)
for k = 1:p
# Append new d[k]-th diagonal to compact form
i = max(1,1-d[k]):min(m,n-d[k])
a[(len[k]+1):len[k+1],:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
end
A = sparse(a[:,1],a[:,2],a[:,3],m,n)
return A
end