Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 52:0236f8e90567

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 11 Jan 2019 11:55:48 +0100
parents 947f7579ba9c
children 4300a3fbd818
comparison
equal deleted inserted replaced
51:614b56a017b9 52:0236f8e90567
30 30
31 # u = L*v 31 # u = L*v
32 function apply!(L::Laplace1D, u::AbstractVector, v::AbstractVector) 32 function apply!(L::Laplace1D, u::AbstractVector, v::AbstractVector)
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 return nothing 36 return nothing
36 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