Mercurial > repos > public > sbplib_julia
changeset 34:bb841977d198
Move stencil operator application to its own function
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 10 Jan 2019 15:49:44 +0100 |
| parents | e8d7137b3f07 |
| children | 8b04efde1a46 |
| files | diffOp.jl sbpD2.jl |
| diffstat | 2 files changed, 27 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
diff -r e8d7137b3f07 -r bb841977d198 diffOp.jl --- a/diffOp.jl Thu Jan 10 14:46:15 2019 +0100 +++ b/diffOp.jl Thu Jan 10 15:49:44 2019 +0100 @@ -35,17 +35,6 @@ h = scaling(L.grid) - for i ∈ 1:N - u[i] = apply(L.op.closureStencils[i], v, i)/h^2 - end - - for i ∈ N+1:M-N - u[i] = apply(L.op.innerStencil, i)/h^2 - end - - for i ∈ M:-1:M-N+1 - u[i] = apply(flip(L.op.closureStencils[M-i+1]), v, i)/h^2 - end - + apply!(L.op, u, v, grid.spacings(L.grid)[1], 1, L.grid.numberOfPointsPerDim, stride=1) return nothing end
diff -r e8d7137b3f07 -r bb841977d198 sbpD2.jl --- a/sbpD2.jl Thu Jan 10 14:46:15 2019 +0100 +++ b/sbpD2.jl Thu Jan 10 15:49:44 2019 +0100 @@ -1,9 +1,33 @@ -struct D2{T} +abstract type ConstantStencilOperator end + +function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real, start::Int, N::Int; stride::Int=1) + cSize = closureSize(op) + + for i ∈ range(start; length=cSize, step=stride) + u[i] = apply!(op.closureStencils[i], v, i; stride=stride)/h^2 + end + + innerStart = start + cSize*stride + innerEnd = N - cSize*stride-1 + for i ∈ range(innerStart, stop=innerEnd, step=stride) + u[i] = apply(op.innerStencil, v, i; stride=stride)/h^2 + end + + for i ∈ range(innerEnd+1, length=cSize, step=cSize) + u[i] = op.parity*apply(flip(op.closureStencils[M-i+1]), v, i; stride=stride)/h^2 + end +end + +odd = -1 +even = 1 + +struct D2{T} <: ConstantStencilOperator quadratureClosure::Vector{T} innerStencil::Stencil closureStencils::Vector{Stencil} # TBD: Should this be a tuple? eClosure::Vector{T} dClosure::Vector{T} + parity::Int end function closureSize(D::D2)::Int @@ -34,6 +58,7 @@ closureStencils, stringToVector(Float64, d["e"]), stringToVector(Float64, d["d1"]), + even ) # Return d2!