Mercurial > repos > public > sbplib_julia
changeset 11:7f075bacbd68
merge with default
author | Ylva Rydin <ylva.rydin@telia.com> |
---|---|
date | Mon, 17 Dec 2018 14:42:38 +0100 |
parents | bed51234616b (current diff) 433008d3b7d3 (diff) |
children | c4bc1165fdf7 |
files | sbp.jl |
diffstat | 5 files changed, 93 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/diffOp.jl Mon Dec 17 14:42:11 2018 +0100 +++ b/diffOp.jl Mon Dec 17 14:42:38 2018 +0100 @@ -0,0 +1,47 @@ +abstract type DiffOp end + +function apply(D::DiffOp, v::AbstractVector) + error("not implemented") +end + +function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real + error("not implemented") +end + +function matrixRepresentation(D::DiffOp) + error("not implemented") +end + +function boundaryCondition(D::DiffOp) + error("not implemented") +end + +function interface(Du::DiffOp, Dv::DiffOp, b::BoundaryID; type) + error("not implemented") +end + + +# Differential operator for a*d^2/dx^2 +struct Laplace1D <: DiffOp + grid + a + op +end + +# u = L*v +function apply(L::Laplace1D, u::AbstractVector, v::AbstractVector)::AbstractVector + N = closureSize(L.op) + M = length(v) + + for i ∈ 1:N + u[i] = apply(L.op.closureStencils[i], v, i) + end + + for i ∈ N+1:M-N + u[i] = apply(L.op.innerStencil, i); + end + + for i ∈ M:-1:M-N+1 + u[i] = apply(flip(L.op.closureStencils[M-i+1]), v, i) + end +end
--- a/grid.jl Mon Dec 17 14:42:11 2018 +0100 +++ b/grid.jl Mon Dec 17 14:42:38 2018 +0100 @@ -0,0 +1,1 @@ +abstract type BoundaryID end
--- a/sbp.jl Mon Dec 17 14:42:11 2018 +0100 +++ b/sbp.jl Mon Dec 17 14:42:38 2018 +0100 @@ -1,6 +1,7 @@ module sbp include("grid.jl") include("diffOp.jl") +include("stencil.jl") include("sbpD2.jl") include("TimeStepper.jl") end # module
--- a/sbpD2.jl Mon Dec 17 14:42:11 2018 +0100 +++ b/sbpD2.jl Mon Dec 17 14:42:38 2018 +0100 @@ -1,1 +1,12 @@ +struct D2{T} + quadratureClosure::Vector{T} + innerStencil::Stencil + closureStencils::Vector{Stencil} # TBD: Should this be a tuple? + eClosure::Vector{T} + dClosure::Vector{T} +end + +function closureSize(D::D2)::Int + return length(quadratureClosure) +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stencil.jl Mon Dec 17 14:42:38 2018 +0100 @@ -0,0 +1,33 @@ +struct Stencil + range::NTuple{2,Int} + weights::Vector # TBD: Should this be a tuple? + function Stencil(range, weights) + width = range[2]-range[1]+1 + if width != length(weights) + error("The width and the number of weights must be the same") + end + new(range, weights) + end +end + +function flip(s::Stencil) + range = (-s.range[2], -s.range[1]) + s = Stencil(range, s.weights[end:-1:1]) +end + +# Provides index into the Stencil based on offset for the root element +function Base.getindex(s::Stencil, i::Int) + if s.range[1] <= i <= s.range[2] + return s.weights[1 + i - s.range[1]] + else + return 0 + end +end + +function apply(s::Stencil, v::AbstractVector, i::Int) + w = zero(v[0]) + for j ∈ i+(s.range[1]:s.range[2]) + w += v[j] + end + return w +end