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