Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 49:947f7579ba9c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Laplace2d
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 11 Jan 2019 11:16:37 +0100
parents be53b58581a6
children 4300a3fbd818
comparison
equal deleted inserted replaced
48:be53b58581a6 49:947f7579ba9c
33 h = grid.spacings(L.grid)[1] 33 h = grid.spacings(L.grid)[1]
34 apply!(L.op, u, v, h) 34 apply!(L.op, u, v, h)
35 u .= L.a * u 35 u .= L.a * u
36 return nothing 36 return nothing
37 end 37 end
38
39
40 # Differential operator for a*d^2/dx^2 + a*d^2/dy^2
41 struct Laplace2D <: DiffOp
42 grid
43 a
44 op
45 end
46
47 # u = L*v
48 function apply!(L::Laplace2D, u::AbstractVector, v::AbstractVector)
49 u .= 0*u
50 h = grid.spacings(L.grid)
51
52 li = LinearIndices(L.grid.numberOfPointsPerDim)
53 n_x, n_y = L.grid.numberOfPointsPerDim
54
55
56 # For each x
57 temp = zeros(eltype(u), n_y)
58 for i ∈ 1:n_x
59
60 v_i = view(v, li[i,:])
61 apply!(L.op, temp, v_i, h[2])
62
63 u[li[i,:]] += temp
64 end
65
66 # For each y
67 temp = zeros(eltype(u), n_x)
68 for i ∈ 1:n_y
69 v_i = view(v, li[:,i])
70 apply!(L.op, temp, v_i, h[1])
71
72 u[li[:,i]] += temp
73 end
74
75 u .= L.a*u
76
77 return nothing
78 end