Mercurial > repos > public > sbplib_julia
view sbpD2.jl @ 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 | 55fea1ceb6aa |
children | 8b04efde1a46 |
line wrap: on
line source
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 return length(quadratureClosure) end function readOperator(D2fn, Hfn) d = readSectionedFile(D2fn) h = readSectionedFile(Hfn) # Create inner stencil innerStencilWeights = stringToVector(Float64, d["inner_stencil"]) width = length(innerStencilWeights) r = (-width//2, width//2) innerStencil = Stencil(r, innerStencilWeights) # Create boundary stencils boundarySize = length(d["boundary_stencils"]) closureStencils = Vector{Stencil}() for i ∈ 1:boundarySize stencilWeights = stringToVector(Float64, d["boundary_stencils"][i]) end d2 = D2( stringToVector(Float64, h["closure"]), innerStencil, closureStencils, stringToVector(Float64, d["e"]), stringToVector(Float64, d["d1"]), even ) # Return d2! return nothing end function readSectionedFile(filename)::Dict{String, Vector{String}} f = open(filename) sections = Dict{String, Vector{String}}() currentKey = "" for ln ∈ eachline(f) if ln == "" || ln[1] == '#' # Skip comments and empty lines continue end if isletter(ln[1]) # Found start of new section if ~haskey(sections, ln) sections[ln] = Vector{String}() end currentKey = ln continue end push!(sections[currentKey], ln) end return sections end function stringToVector(T::DataType, s::String; delimiter = " ") return parse(T, split(s, delimiter)) end