Mercurial > repos > public > sbplib_julia
changeset 91:c0f33eccd527 cell_based_test
Create types StencilIndex, LowerClosureIndex, UpperClosureIndex and InteriorIndex. First attempt at seperating out interior and closure indices from. Not fully implemented.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Tue, 29 Jan 2019 14:32:28 +0100 |
| parents | 8d505e9bc715 |
| children | b8c9e2db126f 98c699eb1d3e |
| files | EquidistantGrid.jl StencilIndex.jl sbp.jl sbpD2.jl |
| diffstat | 4 files changed, 59 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/EquidistantGrid.jl Fri Jan 25 15:26:47 2019 +0100 +++ b/EquidistantGrid.jl Tue Jan 29 14:32:28 2019 +0100 @@ -56,11 +56,11 @@ # @Input: grid - an EquidistantGrid # @Return: points - the points of the grid. function points(grid::EquidistantGrid) + # Signed grid spacing dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1) - - points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid)) # Compute the points based on their Cartesian indices and the signed # grid spacings + points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid)) cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim) for i ∈ 1:numberOfPoints(grid) ci = Tuple(cartesianIndices[i]) .-1 @@ -80,6 +80,7 @@ points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim]) end +# TODO: Move to own plotting module. using PyPlot, PyCall function plotgridfunction(grid::EquidistantGrid, gridfunction)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/StencilIndex.jl Tue Jan 29 14:32:28 2019 +0100 @@ -0,0 +1,43 @@ +abstract type StencilIndex end + +function Base.getindex(si::StencilIndex, i::Int) + return si.index[i] +end + +struct LowerClosureIndex <: StencilIndex + index::CartesianIndex +end + +struct UpperClosureIndex <: StencilIndex + index::CartesianIndex +end + +struct InteriorIndex <: StencilIndex + index::CartesianIndex +end + +# TODO: This should take a Stencil or DiffOp so that we can extract all the +# indices in the closures. +# TODO: Where to place this function? +function stencilindices(grid::Grid.EquidistantGrid) + lowerclosure = Vector{LowerClosureIndex}(undef, 0) + upperclosure = Vector{UpperClosureIndex}(undef, 0) + interior = Vector{InteriorIndex}(undef, 0) + # TODO: Fix such that the indices of the entire closure width is included. + islower = x -> (x == 1) + isupper = x -> (x in grid.numberOfPointsPerDim) + ci = CartesianIndices(grid.numberOfPointsPerDim) + for i ∈ ci + I = Tuple(i) + if any(islower, I) + push!(lowerclosure, LowerClosureIndex(i)) + # TODO: Corner points should be in both Lower and Upper? + # Should they have a separate type? + elseif any(isupper, I) + push!(upperclosure, UpperClosureIndex(i)) + else + push!(interior, InteriorIndex(i)) + end + end + return lowerclosure, upperclosure, interior +end
--- a/sbp.jl Fri Jan 25 15:26:47 2019 +0100 +++ b/sbp.jl Tue Jan 29 14:32:28 2019 +0100 @@ -1,5 +1,6 @@ module sbp include("grid.jl") +include("StencilIndex.jl") include("stencil.jl") include("sbpD2.jl") include("diffOp.jl")
--- a/sbpD2.jl Fri Jan 25 15:26:47 2019 +0100 +++ b/sbpD2.jl Tue Jan 29 14:32:28 2019 +0100 @@ -15,6 +15,18 @@ return uᵢ end +Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::InteriorIndex) + return apply(op.innerStencil, v, i)/h^2 +end + +Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::LowerClosureIndex) + return apply(op.closureStencils[i], v, i)/h^2 +end + +Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::UpperClosureIndex) + return Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2 #TODO: Write an applybackwards instead? +end + @enum Parity begin odd = -1 even = 1