Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 56:27a8d3021a1c cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Convert apply functions to cell-based
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 15 Jan 2019 12:26:06 +0100
parents c62ea0112d4d
children 178a203f3e6d
comparison
equal deleted inserted replaced
55:c62ea0112d4d 56:27a8d3021a1c
1 abstract type DiffOp end 1 abstract type DiffOp end
2 2
3 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector) 3 function apply(D::DiffOp, v::AbstractVector, i::Int)
4 error("not implemented") 4 error("not implemented")
5 end 5 end
6 6
7 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real 7 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
8 error("not implemented") 8 error("not implemented")
30 30
31 function apply(c::Penalty, g, i::Int) 31 function apply(c::Penalty, g, i::Int)
32 error("not implemented") 32 error("not implemented")
33 end 33 end
34 34
35 # DiffOp must have a grid!!!
36 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
37 for i ∈ 1:Grid.numberOfPoints(D.grid)
38 u[i] = apply(D, v, i)
39 end
40
41 return nothing
42 end
43
44 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
45 u = zeros(eltype(v), size(v))
46 apply!(D,v,u)
47 return u
48 end
49
35 # Differential operator for a*d^2/dx^2 50 # Differential operator for a*d^2/dx^2
36 struct Laplace1D <: DiffOp 51 struct Laplace1D <: DiffOp
37 grid 52 grid
38 a 53 a
39 op 54 op
40 end 55 end
41 56
42 # u = L*v 57 # u = L*v
43 function apply!(L::Laplace1D, u::AbstractVector, v::AbstractVector) 58 function apply(L::Laplace1D, v::AbstractVector, i::Int)
44 h = Grid.spacings(L.grid)[1] 59 h = Grid.spacings(L.grid)[1]
45 apply!(L.op, u, v, h) 60 uᵢ = L.a * apply(L.op, h, v, i)
46 u .= L.a * u 61 return uᵢ
47 return nothing
48 end 62 end
49 63
50 64
51 # Differential operator for a*d^2/dx^2 + a*d^2/dy^2 65 # Differential operator for a*d^2/dx^2 + a*d^2/dy^2
52 struct Laplace2D <: DiffOp 66 struct Laplace2D <: DiffOp
54 a 68 a
55 op 69 op
56 end 70 end
57 71
58 # u = L*v 72 # u = L*v
59 function apply!(L::Laplace2D, u::AbstractVector, v::AbstractVector) 73 function apply(L::Laplace2D, v::AbstractVector, i::Int)
60 u .= 0*u
61 h = Grid.spacings(L.grid) 74 h = Grid.spacings(L.grid)
62 75
63 li = LinearIndices(L.grid.numberOfPointsPerDim) 76 li = LinearIndices(L.grid.numberOfPointsPerDim)
64 n_x, n_y = L.grid.numberOfPointsPerDim 77 ci = CartesianIndices(L.grid.numberOfPointsPerDim)
78 I = ci[i]
65 79
80 # 2nd x-derivative
81 uᵢ = apply(L.op, h[1], view(v, li[:,I[2]]), I[1])
82 # 2nd y-derivative
83 uᵢ += apply(L.op, h[2], view(v, li[I[1],:]), I[2])
66 84
67 # For each x 85 return uᵢ
68 temp = zeros(eltype(u), n_y)
69 for i ∈ 1:n_x
70
71 v_i = view(v, li[i,:])
72 apply!(L.op, temp, v_i, h[2])
73
74 u[li[i,:]] += temp
75 end
76
77 # For each y
78 temp = zeros(eltype(u), n_x)
79 for i ∈ 1:n_y
80 v_i = view(v, li[:,i])
81 apply!(L.op, temp, v_i, h[1])
82
83 u[li[:,i]] += temp
84 end
85
86 u .= L.a*u
87
88 return nothing
89 end 86 end