Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 93:93df72e2b135 stencil_index

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement apply for 2D-Laplace which takes an StencilIndex as input
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 04 Feb 2019 09:13:48 +0100
parents 8d505e9bc715
children 84b1ad5a3755
comparison
equal deleted inserted replaced
92:b8c9e2db126f 93:93df72e2b135
35 # DiffOp must have a grid!!! 35 # DiffOp must have a grid!!!
36 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector) 36 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
37 for i ∈ 1:Grid.numberOfPoints(D.grid) 37 for i ∈ 1:Grid.numberOfPoints(D.grid)
38 u[i] = apply(D, v, i) 38 u[i] = apply(D, v, i)
39 end 39 end
40
41 return nothing 40 return nothing
42 end 41 end
43 42
44 function apply(D::DiffOp, v::AbstractVector)::AbstractVector 43 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
45 u = zeros(eltype(v), size(v)) 44 u = zeros(eltype(v), size(v))
61 end 60 end
62 61
63 using UnsafeArrays 62 using UnsafeArrays
64 63
65 # u = L*v 64 # u = L*v
65 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
66 h = Grid.spacings(L.grid)[1]
67 uᵢ = L.a * apply(L.op, h, v, i)
68 return uᵢ
69 end
70
71 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
72 lowerind, upperind, interiorind = stencilindices(L.grid)
73 for si ∈ lowerind
74 u[si.globalindex] = apply(L, v, si)
75 end
76 for si ∈ upperind
77 u[si.globalindex] = apply(L, v, si)
78 end
79 for si ∈ interiorind
80 u[si.globalindex] = apply(L, v, si)
81 end
82 return nothing
83 end
84
85 @inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
86 h = Grid.spacings(L.grid)
87 li = LinearIndices(L.grid.numberOfPointsPerDim)
88 # 2nd x-derivative
89 @inbounds vx = uview(v, uview(li,:,si.gridindex[2]))
90 uᵢ = apply(L.op, h[1], vx , si.gridindex[1], si)
91 # 2nd y-derivative
92 @inbounds vy = uview(v, uview(li,si.gridindex[1],:))
93 uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si)
94 return uᵢ
95 end
96
66 function apply(L::Laplace{2}, v::AbstractVector, i::Int) 97 function apply(L::Laplace{2}, v::AbstractVector, i::Int)
67 h = Grid.spacings(L.grid) 98 h = Grid.spacings(L.grid)
68 99
69 li = LinearIndices(L.grid.numberOfPointsPerDim) 100 li = LinearIndices(L.grid.numberOfPointsPerDim)
70 ci = CartesianIndices(L.grid.numberOfPointsPerDim) 101 ci = CartesianIndices(L.grid.numberOfPointsPerDim)