code: split off FD project

FD/Manifest.toml

+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distributed]]
+deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[InteractiveUtils]]
+deps = ["LinearAlgebra", "Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[OffsetArrays]]
+deps = ["DelimitedFiles", "Test"]
+git-tree-sha1 = "f8cd140824f74f2948ff8660bc2292e16d29c1b7"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.8.1"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[StaticArrays]]
+deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"]
+git-tree-sha1 = "d432c79bef174a830304f8601427a4357dfdbfb7"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "0.8.3"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

FD/Project.toml

+name = "FD"
+uuid = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
+authors = ["Stephan Hilb <stephan.hilb@mathematik.uni-stuttgart.de>"]
+version = "0.1.0"
+
+[deps]
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

FD/src/FD.jl

+module FD
+
+using StaticArrays
+#using SparseArrays
+using OffsetArrays
+#using LinearMaps
+#using IterativeSolvers
+
+include("filtering.jl")
+
+"Finite Differences"
+module FDiff
+
+using OffsetArrays
+
+const forward = OffsetVector([-1, 1], -1)
+const backward = OffsetVector([-1, 1], -2)
+const central = OffsetVector([-1, 0, 1], -2)
+function dirfd(fd, ndims, dir)
+  dir <= ndims || throw(ArgumentError("need dir <= ndims!"))
+  dfd = OffsetArray(zeros(ntuple(i -> i == dir ? length(fd) : 1, ndims)), ntuple(i -> i == dir ? fd.offsets[1] : -1, ndims))
+  # Workaround: LinearIndices(::OffsetVector) are not 1-based, so copyto! chockes
+  copyto!(dfd, fd.parent)
+end
+
+end
+
+
+using .Filtering
+using .FDiff
+
+"""
+    _fdkern(fd, ndims, dir)
+
+Create from `fd` a finite difference kernel usable in `ndims` dimension
+oriented in dimension `dir`.
+"""
+# TODO: make it typestable on ndims
+function _fdkern(fd::OffsetArray{T,1,Array{T,1}}, ndims::Integer, dir::Integer) where T
+  1 <= dir <= ndims || throw(ArgumentError("dimesion $(dir) out of range 1:$(ndims)"))
+  dfd = OffsetArray{T}(undef, Tuple(i == dir ? axes(fd, 1) : (0:0) for i in 1:ndims))
+  # copyto! from OffsetArray fails, broadcasting for differing axes fails too
+  # workaround: copyto! from parent
+  copyto!(dfd, fd.parent)
+end
+
+
+# Grid
+
+
+abstract type AbstractGrid{d} <: AbstractArray{Float64, d} end
+
+Base.ndims(grid::AbstractGrid{d}) where d = d
+Base.axes(grid::AbstractGrid) = Base.OneTo.(size(grid))
+Base.length(grid::AbstractGrid) = prod(size(grid))
+Base.IndexStyle(::Type{<:AbstractGrid}) = IndexCartesian()
+
+# aka `StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}`
+const FDStepRange = typeof(0:0.1:1)
+
+struct Grid{d} <: AbstractGrid{d}
+  domain::NTuple{d, FDStepRange}
+end
+Grid(domain::FDStepRange...) = Grid(domain)
+Base.show(io::IO, grid::Grid) =
+    print("Grid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
+          "with spacing $(step.(grid.domain))")
+
+Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
+
+Base.iterate(grid::Grid) = iterate(Iterators.product(grid.domain...))
+Base.iterate(grid::Grid, state) = iterate(Iterators.product(grid.domain...), state)
+
+Base.getindex(grid::Grid, idx...) = SVector(getindex.(grid.domain, idx))
+Base.getindex(grid::Grid, I::CartesianIndex) = SVector(getindex.(grid.domain, Tuple(I)))
+
+
+# GridView
+
+struct GridView{d} <: AbstractGrid{d}
+  grid::AbstractGrid{d}
+  indices::NTuple{d, UnitRange{Int}}
+end
+GridView(grid, indices::UnitRange{Int}...) = GridView(grid, indices)
+
+Base.size(grid::GridView) = length.(grid.indices)
+
+# FIXME: this is terrible!
+Base.iterate(grid::GridView) = iterate(grid.grid[I] for I in CartesianIndices(grid.indices))
+Base.iterate(grid::GridView, state) = iterate([grid.grid[I] for I in CartesianIndices(grid.indices)], state)
+
+Base.getindex(grid::GridView, idx...) = getindex(grid.grid, (first.(grid.indices) .- 1 .+ Base.to_indices(grid, idx))...)
+
+# MetaGrid
+
+
+struct MetaGrid{d} <: AbstractGrid{d}
+  grid::Grid{d}
+  indices::Array{NTuple{d, UnitRange{Int}}, d}
+  overlap::Int
+end
+
+"""
+    MetaGrid(domain::NTuple{d, FDStepRange}, pnum::NTuple{d, Int}, overlap::Int)
+
+Create a grid on `domain`, partitioned along dimensions according to `pnum`
+respecting `overlap`.
+"""
+function MetaGrid(domain::NTuple{d, FDStepRange}, pnum::NTuple{d, Int}, overlap::Int) where d
+  grid = Grid(domain)
+  tsize = size(grid) .+ (pnum .- 1) .* overlap
+
+  psize = tsize .÷ pnum
+  osize = tsize .- pnum .* psize
+
+  overhang(I, j) = I[j] == pnum[j] ? osize[j] : 0
+
+  indices = Array{NTuple{d, UnitRange{Int}}, d}(undef, pnum)
+  for I in CartesianIndices(pnum)
+    indices[I] = ntuple(j -> ((I[j] - 1) * psize[j] - (I[j] - 1) * overlap + 1) :
+                             (      I[j] * psize[j] - (I[j] - 1) * overlap + overhang(I, j)), d)
+  end
+  MetaGrid{d}(grid, indices, overlap)
+end
+
+Base.size(grid::MetaGrid) = size(grid.grid)
+
+Base.iterate(grid::MetaGrid) = iterate(grid.grid)
+Base.iterate(grid::MetaGrid, state) = iterate(grid.grid, state)
+
+Base.getindex(grid::MetaGrid, idx...) = getindex(grid.grid, idx...)
+
+indices(grid::MetaGrid, idx...) = getindex(grid.indices, idx...)
+
+Base.view(grid::MetaGrid, idx...) = GridView(grid, indices(grid, idx...))
+
+parts(grid::MetaGrid) = CartesianIndices(grid.indices)
+
+
+## GridFunction
+
+abstract type AbstractGridFunction{G,T,d,n} <: AbstractArray{T,1} end
+
+rangedim(f::AbstractGridFunction{G,T,d,n}) where {G,T,d,n} = n
+domaindim(f::AbstractGridFunction{G,T,d}) where {G,T,d} = d
+
+Base.similar(f::F, ::Type{S}, dims::Int...) where {F<:AbstractGridFunction, S} =
+    GridFunction(f.grid, Array{S}(undef, length(f)))
+
+
+"""
+    GridFunction{G<:AbstractGrid,T,d,n,D}
+
+`n`-dimensional function values of type `T` on a Grid `G` of dimension `d`.
+
+`D == d + 1`
+"""
+struct GridFunction{G<:AbstractGrid,T,d,n,A,D} <: AbstractGridFunction{G,T,d,n}
+  grid::G
+  data::A
+  function GridFunction{G,T,d,n,A,D}(grid::G, data::A) where {G,T,d,n,D,A<:AbstractArray{T,D}}
+    D == d + 1 || throw(ArgumentError)
+    d == ndims(grid) || throw(ArgumentError)
+    (size(grid)..., n) == size(data) || throw(ArgumentError("data does not match given size"))
+    new(grid, data)
+  end
+end
+# TODO: check size(data)
+#GridFunction(grid::G, data::AbstractArray{T,D}) where {d,G<:AbstractGrid{d},T,D} =
+#    GridFunction{G,T,d,size(data, d+1),typeof(data),d+1}(grid, data)
+function GridFunction(grid::G, data::AbstractArray, ::Val{n}) where {d,G<:AbstractGrid{d}, n}
+  if length(data) ÷ length(grid) == 1 && n > 1
+    data = repeat(data, inner=ntuple(i -> i == d+1 ? n : 1, d+1))
+  else
+    # TODO: check size(data)
+    data = reshape(data, size(grid)..., n)
+  end
+  length(data) ÷ length(grid) == n || throw(ArgumentError("dimensions don't match"))
+  GridFunction{G,eltype(data),d,n,typeof(data),d+1}(grid, data)
+end
+function GridFunction(grid::G, data::AbstractArray) where {d,G<:AbstractGrid{d}}
+  GridFunction(grid, data, Val(length(data) ÷ length(grid)))
+end
+function GridFunction(f::Function, n::Int, grid::G) where {G<:AbstractGrid}
+  data = Array{Float64,ndims(grid)+1}(undef, size(grid)..., n)
+  for I in CartesianIndices(grid)
+    data[Tuple(I)..., :] .= f(grid[I])
+  end
+  GridFunction(grid, data)
+end
+function GridFunction(grid::AbstractGrid, ::UndefInitializer, ::Val{n} = Val(1)) where n
+  data = Array{Float64}(undef, size(grid)..., 1)
+  GridFunction{typeof(grid), eltype(data), ndims(grid), n, typeof(data), ndims(grid) + 1}(grid, data)
+end
+
+#Base.show(io::IO, grid::GridFunction) =
+#    print("GridFunction on $(length(grid)) grid points")
+
+Base.size(f::GridFunction) = (length(f.data),)
+#Base.axes(f::GridFunction) = (Base.OneTo(length(f.data)),)
+Base.getindex(f::GridFunction, idx::Int) = getindex(f.data, idx)
+Base.setindex!(f::GridFunction, v, idx::Int) = setindex!(f.data, v, idx)
+
+## TODO: using this prints #undef for display(f)
+#Base.getindex(f::GridFunction, idx...) = getindex(f.data, idx...)
+#Base.setindex!(f::GridFunction, v, idx...) = setindex!(f.data, v, idx...)
+
+#Base.IndexStyle(::Type{<:GridFunction}) = IndexLinear()
+
+#component(f::GridFunction, k::Integer) = view(reshape(f.data,
+
+Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
+    GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
+
+grid(f::GridFunction) = f.grid
+
+
+Base.:*(a::Number, b::GridFunction) = GridFunction(b.grid, a * b.data, Val(rangedim(b)))
+
+
+
+
+## Operations
+
+function fdfilter(f::GridFunction{G,T,d,n}, kern) where {G,T,d,n}
+  A = zeros(axes(f.data))
+  for k = 1:n
+    vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
+    imfilter!(view(A, vidx...), f.data[vidx...], kern)
+  end
+  GridFunction(f.grid, A)
+end
+
+derivate(f::GridFunction, dir::Integer; fd = FDiff.forward) = fdfilter(f, _fdkern(fd, ndims(f.grid), dir))
+
+
+# currently gradient and divergence are specific implementations (as needed for
+# chambolles algorithm), other choices may be reasonable
+
+function gradient(f::GridFunction{G,T,d,1}; fd=FDiff.forward) where {G,T,d}
+  #A = Array{T,d+1}(undef, (size(f.grid)..., d))
+  A = zeros(size(f.grid)..., d)
+  for k = 1:d
+    vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
+    vidx2 = ntuple(i -> i <= d ? Colon() : 1, d + 1)
+    imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
+    # FIXME: hackish, only for forward fd
+    selectdim(selectdim(A, d+1, k), k, size(f.grid, k)) .= 0
+  end
+  GridFunction(f.grid, A)
+end
+
+function divergence(f::GridFunction{G,T,d}; fd=FDiff.backward) where {d,G<:AbstractGrid{d},T}
+  A = zeros(size(f.grid)..., 1)
+  for k = 1:d
+    vidx = ntuple(i -> i <= d ? Colon() : 1, d + 1)
+    vidx2 = ntuple(i -> i <= d ? Colon() : k, d + 1)
+    imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
+  end
+  GridFunction(f.grid, A)
+end
+
+function pw_norm(f::GridFunction)
+  d = domaindim(f)
+  GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(1))
+end
+
+function my_norm(f::GridFunction)
+  d = domaindim(f)
+  GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(rangedim(f)))
+end
+
+function pwmul(f::GridFunction{G,T,d,1}, g::GridFunction{G,T,d}) where {G,T,d}
+  # TODO: maybe comprehension?
+  data = similar(g.data)
+  for I in CartesianIndices(f.grid)
+    data[Tuple(I)..., :] = f.data[Tuple(I)..., 1] * g.data[Tuple(I)..., :]
+  end
+  GridFunction(f.grid, data, Val(rangedim(g)))
+end
+
+
+
+
+"""
+    neuman_kernel(kernel, bidx)
+
+Modify `kernel` to be usable on a boundary `bixd` in a neumannn boundary setting.
+
+Ghost nodes are eliminated through central finite differences. Currently only
+the basic 5-star Laplace is supported.
+"""
+function neumann_kernel(kernel, bidx)
+  out = copy(kernel)
+  for (i, bi) in enumerate(bidx)
+    colrange = (i == j ? Colon() : 0 for j in 1:ndims(kernel))
+    out[colrange...] .+= bi * FDiff.forward
+  end
+  out /= 2 ^ count(!iszero, bidx)
+  trim_kernel(out, bidx)
+end
+
+function divergence(grid::Grid{d}, field) where d
+  field = reshape(field, (size(grid)..., ndims(grid)))
+  size(field, d+1) == d || throw(ArgumentError("need a d-vectorfield on a d-grid"))
+  dst = zeros(size(field))
+  for k in 1:d
+    kern = FDiff.dirfd(FDiff.backward, d, k)
+    imfilter!(selectdim(dst, d+1, k), selectdim(field, d+1, k), bidx -> trim_kernel(kern, bidx))
+  end
+  out = dropdims(sum(dst, dims=d+1), dims=d+1)
+  out[:]
+end
+
+function gradient(grid::Grid, img::AbstractArray)
+  img = reshape(img, size(grid))
+  d = ndims(img)
+  dst = zeros(size(img)..., d)
+  for k in 1:d
+    kern = FDiff.dirfd(FDiff.forward, d, k)
+    imfilter!(selectdim(dst, d+1, k), img, bidx -> trim_kernel(kern, bidx))
+    # FIXME: hack
+    selectdim(selectdim(dst, d+1, k), k, size(grid, k)) .= 0
+  end
+  reshape(dst, length(grid) * ndims(grid))
+end
+
+function my_norm(grid::Grid, img::AbstractArray)
+  img = reshape(img, (size(grid)..., ndims(grid)))
+  d = ndims(grid)
+
+  out = sqrt.(sum(img.^2, dims=d+1))
+  out = repeat(out, inner=ntuple(i -> i == d+1 ? size(img, d+1) : 1, d+1))
+
+  out[:]
+end
+
+
+function apply(kern, flatimg, size)
+  img = reshape(flatimg, size)
+  out = imfilter(kern, img)
+  reshape(out, prod(size))
+end
+
+function solve(pde, shape::NTuple{d, Integer}) where d
+  ndims(pde.stencil) == d || throw("dimension mismatch")
+  length(pde.f) == prod(shape) || throw("bounds mismatch")
+
+  A = LinearMap{Float64}(u -> apply(pde.stencil, u, shape), prod(shape), issymmetric=true, isposdef=true)
+
+  u = zeros(prod(shape))
+  cg!(u, A, pde.f)
+  return u
+end
+
+end
+
+
+#module FD
+#
+#using OffsetArrays
+#include("filtering.jl")
+#
+#
+#"FD-Kernels"
+#module Kernels
+#
+#using OffsetArrays
+#
+#const forward = OffsetVector([-1, 1], 0:1)
+#const backward = OffsetVector([-1, 1], -1:0)
+#const central = OffsetVector([-1, 0, 1], -1:1)
+#
+#end
+#
+#using .Filtering
+#using .Kernels
+#
+#
+#
+### Implementation
+#
+#
+### Grid
+#
+## aka `StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}`
+#const FDStepRange = typeof(0:0.1:1)
+#
+#struct Grid{d}
+#  domain::NTuple{d, FDStepRange}
+#end
+#
+#Grid(domain::FDStepRange...) = Grid(domain)
+#Grid(ndims::Integer, range::FDStepRange) = Grid(ntuple(i -> range, ndims))
+#
+#Base.show(io::IO, grid::Grid) =
+#    print("FDGrid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
+#          "with spacing $(step.(grid.domain))")
+#
+## this is not type-stable?
+##Base.size(grid::Grid) = getproperty.(grid.domain, :len)
+#Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
+#Base.getindex(grid::Grid, idx...) = SVector(getindex.(grid.domain, idx))
+#
+#Base.ndims(grid::Grid{d}) where d = d
+#Base.axes(grid::Grid) = Base.OneTo.(size(grid))
+#Base.length(grid::Grid) = prod(size(grid))
+#
+## TODO: there might be a better way to wrap that product iterator
+#Base.iterate(grid::Grid) = iterate(Iterators.product(grid.domain...))
+#Base.iterate(grid::Grid, state) = iterate(Iterators.product(grid.domain...), state)
+#
+#
+### GridFunction
+#
+#struct GridFunction{T,d} <: AbstractArray{T,1}
+#  grid::Grid{d}
+#  data::Array{T,1}
+#end
+#
+#Base.show(io::IO, grid::GridFunction) =
+#    print("GridFunction on $(length(grid)) grid points")
+#
+#Base.size(f::GridFunction) = size(f.data)
+#Base.getindex(f::GridFunction, idx...) = getindex(f.data, idx...)
+#Base.setindex!(f::GridFunction, v, idx...) = setindex!(f.data, v, idx...)
+#
+#Base.IndexStyle(::Type{<:GridFunction}) = IndexLinear()
+#Base.similar(f::GridFunction, ::Type{S}, size::Int...) where S =
+#    GridFunction(f.grid, Array{S}(undef, size))
+#
+### Operations
+#
+#function fdfilter(f::GridFunction, kern)
+#  A = imfilter(reshape(f.data, size(f.grid)), kern)
+#  GridFunction(f.grid, A[:])
+#end
+#
+#derivate(f::GridFunction, dir::Integer; fd = FDiff.forward) = fdfilter(f, _fdkern(fd, ndims(f.grid), dir))
+#
+##gradient(f::GridFunction) = GridFunction(f.grid,
+#
+#
+#"""
+#    neuman_kernel(kernel, bidx)
+#
+#Modify `kernel` to be usable on a boundary `bixd` in a neumannn boundary setting.
+#
+#Ghost nodes are eliminated through central finite differences. Currently only
+#the basic 5-star Laplace is supported.
+#"""
+#function neumann_kernel(kernel, bidx)
+#  out = copy(kernel)
+#  for (i, bi) in enumerate(bidx)
+#    colrange = (i == j ? Colon() : 0 for j in 1:ndims(kernel))
+#    out[colrange...] .+= bi * FDiff.forward
+#  end
+#  out /= 2 ^ count(!iszero, bidx)
+#  trim_kernel(out, bidx)
+#end
+#
+#function divergence(grid::Grid{d}, field) where d
+#  field = reshape(field, (size(grid)..., ndims(grid)))
+#  size(field, d+1) == d || throw(ArgumentError("need a d-vectorfield on a d-grid"))
+#  dst = zeros(size(field))
+#  for k in 1:d
+#    kern = FDiff.dirfd(FDiff.backward, d, k)
+#    imfilter!(selectdim(dst, d+1, k), selectdim(field, d+1, k), kern)
+#  end
+#  out = dropdims(sum(dst, dims=d+1), dims=d+1)
+#  reshape(out, length(grid))
+#end
+#
+#function gradient(grid::Grid, img::AbstractArray)
+#  img = reshape(img, size(grid))
+#  d = ndims(img)
+#  dst = zeros(size(img)..., d)
+#  for k in 1:d
+#    kern = FDiff.dirfd(FDiff.forward, d, k)
+#    imfilter!(selectdim(dst, d+1, k), img, kern)
+#  end
+#  reshape(dst, length(grid) * ndims(grid))
+#end
+#
+#end

FD/src/filtering.jl

+module Filtering
+
+export imfilter, imfilter!, trim_kernel
+
+@inline function _shift_range(range, shift::Integer)
+  Base.UnitRange(first(range)+shift : last(range)+shift)
+end
+
+## TODO: allow border to be tuple
+#function _innershift_indices(a::Array{T, d}, shift::NTuple{d, Integer}, border = abs(maximum(shift))) where {T, d}
+#  abs(maximum(shift)) <= border || throw("cannot shift by $shift with border $border")
+#  range(i) = (1 + shift[i] + border):(lastindex(a, i) + shift[i] - border)
+#  ntuple(range, Val(ndims(a)))
+#end
+#
+#function innershift_view(a::Array{T, d}, shift::NTuple{d, Integer} = ntuple(i -> 0, d), border = abs(maximum(shift))) where {T, d}
+#  view(a, _innershift_indices(a, shift, border)...)
+#end
+
+@inline function _bidx_to_range(range::AbstractUnitRange, bidx::Integer)
+  if bidx > 0
+    return last(range):last(range)
+  elseif bidx < 0
+    return first(range):first(range)
+  else
+    return first(range)+1:last(range)-1
+  end
+end
+
+
+# we return Base.Slice since that is what `OffsetArray` uses for `axes()`
+function _trim_kernel_axes(kaxes, bidx)
+  ntuple(Val(length(kaxes))) do i
+    if bidx[i] < 0
+      return Base.Slice(0:last(kaxes[i]))
+    elseif bidx[i] > 0
+      return Base.Slice(first(kaxes[i]):0)
+    else
+      return Base.Slice(kaxes[i])
+    end
+  end
+end
+
+"""
+    trim_kernel(kernel, bidx)
+
+Trim `kernel` so that the resulting kernel is usable on the boundary
+indicated by `bidx`
+"""
+function trim_kernel(kernel, bidx)
+  dst_axes = _trim_kernel_axes(axes(kernel), bidx)
+  out = similar(kernel, dst_axes)
+  range = CartesianIndices(out)
+  copyto!(out, range, kernel, range)
+end
+
+
+"""
+    imfilter!(dst, img, kern, slice=:inner)
+
+Filters `img` applying `kern` and write to `dst`. `dst` is added to, so you
+should initialize with zeros.
+"""
+# TODO: weirdly Base.Slice instead of Base.UnitRange leads to axes that are not
+# compatible when broadcasting, figure out why!
+function imfilter!(dst, img, kern, slice=_bidx_to_range.(axes(img), ntuple(i->0, ndims(img))))
+  ndims(dst) == ndims(img) == ndims(kern) == length(slice) ||
+      throw(ArgumentError("dst, image, kernel or slice dimensions disagree ($(ndims(dst)) vs $(ndims(img)) vs $(ndims(kern)) vs $(length(slice)))"))
+  axes(dst) == axes(img) || throw(ArgumentError("output and input image axes disagree"))
+  all(size(kern) .<= size(img)) || throw(ArgumentError("kernel bigger than image"))
+  last.(slice) .+ last.(axes(kern)) <= last.(axes(img)) || first.(slice) .+ first.(axes(kern)) >= first.(axes(img)) ||
+      throw(ArgumentError("kern and/or slice out of bounds for image"))
+
+  dstslice = view(dst, UnitRange.(slice)...)
+  for i in CartesianIndices(kern)
+    iszero(kern[i]) && continue
+    dstslice .+= kern[i] .* view(img, (_shift_range.(slice, Tuple(i)))...)
+  end
+  dst
+end
+
+function imfilter!(dst, img, kernelf::Function)
+  for bidx in Iterators.product(ntuple(i -> -1:1, ndims(img))...)
+    range = _bidx_to_range.(axes(img), bidx)
+    imfilter!(dst, img, kernelf(bidx), range)
+  end
+  dst
+end
+
+"""
+    imfilter(img, kernelf)
+
+Filter `img` using a kernel function `kernelf(bidx) = kernel` specifying a
+kernel for each type of node.
+"""
+# need zeros here, because `imfilter!` acts additively
+imfilter(img, kernelf::Function) = imfilter!(zeros(size(img)), img, kernelf)
+
+
+"""
+    imfilter(img, kern)
+
+Filter `img` by `kern` using correlation.
+Boundary is zeroed out.
+"""
+imfilter(img, kern) = imfilter(img, bidx -> trim_kernel(kern, bidx))
+# TODO: use keyword arguments
+imfilter_inner(img, kern) = imfilter(img, bidx -> all(bidx .== 0) ? kern : zeros(ntuple(i -> 0, ndims(img))))
+
+end

Manifest.toml

+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[ColorTypes]]
+deps = ["FixedPointNumbers", "Random", "Test"]
+git-tree-sha1 = "0e3209ba7418aed732e5c3818076b4400ee36c08"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.7.4"
+
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Pkg", "Printf", "Reexport", "Test"]
+git-tree-sha1 = "8c89e0a9a583954eae3efcf6a531e51c02b38cee"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.9.4"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "ae262fa91da6a74e8937add6b613f58cd56cdad4"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "1.1.0"
+
+[[Contour]]
+deps = ["LinearAlgebra", "StaticArrays", "Test"]
+git-tree-sha1 = "b974e164358fea753ef853ce7bad97afec15bb80"
+uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
+version = "0.5.1"
+
+[[DataStructures]]
+deps = ["InteractiveUtils", "REPL", "Random", "Serialization", "Test"]
+git-tree-sha1 = "5162c360c52e9265e0bfafb71415c31f746f7e04"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.12.0"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distributed]]
+deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[FD]]
+deps = ["OffsetArrays", "StaticArrays"]
+path = "FD"
+uuid = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
+version = "0.1.0"
+
+[[FixedPointNumbers]]
+deps = ["Pkg", "Test"]
+git-tree-sha1 = "b8045033701c3b10bf2324d7203404be7aef88ba"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.5.3"
+
+[[GR]]
+deps = ["Base64", "DelimitedFiles", "Pkg", "Random", "Serialization", "Sockets", "Test"]
+git-tree-sha1 = "0437921cbf8ee555ab7a2a3115c1d0907289e816"
+uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
+version = "0.34.1"
+
+[[InteractiveUtils]]
+deps = ["LinearAlgebra", "Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JSON]]
+deps = ["Dates", "Distributed", "Mmap", "Pkg", "Sockets", "Test", "Unicode"]
+git-tree-sha1 = "fec8e4d433072731466d37ed0061b3ba7f70eeb9"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.19.0"
+
+[[LibGit2]]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Measures]]
+deps = ["Pkg", "Test"]
+git-tree-sha1 = "ddfd6d13e330beacdde2c80de27c1c671945e7d9"
+uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
+version = "0.3.0"
+
+[[Missings]]
+deps = ["Dates", "InteractiveUtils", "SparseArrays", "Test"]
+git-tree-sha1 = "196528fa1df03df435025f52f9c1ff8356f94738"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.3.0"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[NaNMath]]
+deps = ["Compat"]
+git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.2"
+
+[[OffsetArrays]]
+deps = ["DelimitedFiles", "Test"]
+git-tree-sha1 = "f8cd140824f74f2948ff8660bc2292e16d29c1b7"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.8.1"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[PlotThemes]]
+deps = ["PlotUtils", "Requires", "Test"]
+git-tree-sha1 = "f3afd2d58e1f6ac9be2cea46e4a9083ccc1d990b"
+uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
+version = "0.3.0"
+
+[[PlotUtils]]
+deps = ["Colors", "Dates", "Printf", "Random", "Reexport", "Test"]
+git-tree-sha1 = "78553a920c4869d20742ba66193e3692369086ec"
+uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
+version = "0.5.4"
+
+[[Plots]]
+deps = ["Base64", "Contour", "Dates", "FixedPointNumbers", "GR", "JSON", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "Random", "RecipesBase", "Reexport", "Requires", "Showoff", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "Test", "UUIDs"]
+git-tree-sha1 = "fa4f80b7e2173b8605353bdadf10d3ee8569efc6"
+uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+version = "0.20.2"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[RecipesBase]]
+deps = ["Pkg", "Random", "Test"]
+git-tree-sha1 = "0b3cb370ee4dc00f47f1193101600949f3dcf884"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "0.6.0"
+
+[[Reexport]]
+deps = ["Pkg"]
+git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "0.2.0"
+
+[[Requires]]
+deps = ["Test"]
+git-tree-sha1 = "f6fbf4ba64d295e146e49e021207993b6b48c7d1"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "0.5.2"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[Showoff]]
+deps = ["Compat", "Pkg"]
+git-tree-sha1 = "276b24f3ace98bec911be7ff2928d497dc759085"
+uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
+version = "0.2.1"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures", "Random", "Test"]
+git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "0.3.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[StaticArrays]]
+deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"]
+git-tree-sha1 = "d432c79bef174a830304f8601427a4357dfdbfb7"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "0.8.3"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsBase]]
+deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "Test"]
+git-tree-sha1 = "723193a13e8078cec6dcd0b8fe245c8bfd81690e"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.25.0"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[UUIDs]]
+deps = ["Random"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

Project.toml

+[deps]
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+FD = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

chambollefp.jl

-include("fd.jl")
-
 using StaticArrays
 using SparseArrays
 using OffsetArrays
@@ -7,16 +5,18 @@ using LinearAlgebra
 using Colors
 using Plots
 
+using FD
+
 
 imf(f, grid) = [f(x) for x in grid]
 
 imshow(grid, img) = plot(Gray.((reshape(float.(img), size(grid)) .- minimum(float.(img))) / (maximum(float.(img)) - minimum(float.(img))) ), show=true)
 
 function _fpiter(grid, p, f, A, α, β, τ)
-  x = (A'*A + β*I) \ (FDM.divergence(grid, p) .+ A'*f ./ α)
-  q = FDM.gradient(grid, x)
+  x = (A'*A + β*I) \ (FD.divergence(grid, p) .+ A'*f ./ α)
+  q = FD.gradient(grid, x)
 
-  p_new = (p .+ τ .* q) ./ (1 .+ τ .* FDM.my_norm(grid, q))
+  p_new = (p .+ τ .* q) ./ (1 .+ τ .* FD.my_norm(grid, q))
 
   reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
   (p_new, reserr)
@@ -32,7 +32,7 @@ function fpalg(grid, f, A; α, β, τ = 1/8, ϵ = 1e-4)
     #sleep(0.1 / k)
     (p, reserr) = _fpiter(grid, p, f, A, α, β, τ)
 
-    uimg = (A'*A + β*I) \ (α * FDM.divergence(grid, p) .+ A'*f)
+    uimg = (A'*A + β*I) \ (α * FD.divergence(grid, p) .+ A'*f)
 
     imshow(grid, uimg)
     println("[$(lpad(k, 4))] res = $(round(reserr, digits=5))")
@@ -43,8 +43,8 @@ end
 
 #st = OffsetArray([0 -1 0; -1 4 -1; 0 -1 0], (-2, -2))
 
-grid = FDM.Grid(0:0.01:1, 0:0.01:1)
-#grid = FDM.Grid(0:0.5:1, 0:0.5:1)
+grid = FD.Grid(0:0.01:1, 0:0.01:1)
+#grid = FD.Grid(0:0.5:1, 0:0.5:1)
 
 n = 1
 #A = spdiagm(collect(k => fill(1/(2n+1), length(grid)-abs(k)) for k in -n:n)...)
@@ -68,9 +68,9 @@ img = fpalg(grid, f, A, α=1, β=1, τ=1/100, ϵ=1e-7)
 #A = Array{Float64}(undef, (length(grid), length(grid)))
 #for i in 1:length(grid)
 #  ei = float.(1:length(grid) .== i)
-#  #A[:, i] = FDM.divergence(grid, FDM.gradient(grid, ei))
-#  #A[:, i] = FDM.divergence(grid, append!(ei, zeros(length(grid))))
-#  A[:, i] = reshape(reshape(FDM.gradient(grid, ei), size(grid)..., 2)[:,:,2], length(grid))
+#  #A[:, i] = FD.divergence(grid, FD.gradient(grid, ei))
+#  #A[:, i] = FD.divergence(grid, append!(ei, zeros(length(grid))))
+#  A[:, i] = reshape(reshape(FD.gradient(grid, ei), size(grid)..., 2)[:,:,2], length(grid))
 #end
 #A
 

chambollefp_dd.jl

-include("fd.jl")
+#include("fd.jl")
 
 using LinearAlgebra
 using SparseArrays
-import .FDM
+import FD
+
+
+import Colors: Gray
+import Plots: plot
+
+# TODO: add clamp option
+my_plot(img::FD.AbstractGridFunction{G,T,d,1}) where {G,T,d} =
+    plot(reshape(Gray.((img .- minimum(img)) ./ (maximum(img) - minimum(img))), size(img.grid)), show=true)
+
 
 "partition of unity function"
 function theta(grid, idx...)
-  vidx = FDM.indices(grid, idx...)
-
-  slices = Vector{Vector{Float64}}(undef, length(idx))
-  for (k, kglobrange) in enumerate(vidx)
-    krange = Base.OneTo(length(kglobrange))
-    slices[k] = min.(krange, reverse(krange), grid.overlap) ./ grid.overlap
-    # FIXME: dirty hack
-    if first(kglobrange) == 1
-      slices[k][1:grid.overlap] .= 1.0
-    end
-    if last(kglobrange) == size(grid, k)
-      slices[k][end-grid.overlap:end] .= 1.0
-    end
+  function h(r, rp, i)
+    i <= first(r) + grid.overlap && return 1.0
+    i >= last(r) - grid.overlap && return 1.0
+    (min(i - first(rp), last(rp) - i, grid.overlap) + 1) / (grid.overlap + 1)
   end
 
-  #dst = Array{Float64}(undef, size(grid)...)
+  vidx = FD.indices(grid, idx...)
+  ctheta = map(x -> [h(x[1], x[2], i) for i in x[2]], zip(axes(grid), vidx))
+
   dst = zeros(size(grid))
   for (I,J) in zip(CartesianIndices(vidx), CartesianIndices(length.(vidx)))
-    dst[I] = prod(getindex.(slices, Tuple(J)))
+    dst[I] = prod(getindex.(ctheta, Tuple(J)))
   end
-  dst
-
-  FDM.GridFunction(grid, dst)
+  FD.GridFunction(grid, dst)
 end
 
 function p2u(grid, p, A, f; α, β)
-  FDM.GridFunction(grid, (A'*A + β*I) \ (α * FDM.divergence(p).data[:] .+ A'*f.data[:]))
+  FD.GridFunction(grid, (A'*A + β*I) \ (α .* FD.divergence(p).data[:] .+ A'*f.data[:]))
 end
 
 # TODO: investigate matrix multiplication performance with GridFunctions
-function fpiter(grid, p, f, A, α, β, τ)
-  x = FDM.GridFunction(grid, (A'*A + β*I) \ (FDM.divergence(p).data[:] .+ A'*f.data[:] ./ α))
-  q = FDM.gradient(x)
+function fpiter(grid, p, w, A, α, β, τ)
+  x = FD.GridFunction(grid, (A'*A + β*I) \ (FD.divergence(p).data[:] .+ w ./ α))
+  q = FD.gradient(x)
 
   #display(x.data)
-  #FDM.my_plot(x)
+  #FD.my_plot(x)
   #checknan(q)
-  p_new = FDM.GridFunction(grid, (p .+ τ .* q) ./ (1 .+ τ .* FDM.my_norm(q)))
+  p_new = FD.GridFunction(grid, FD.pwmul(FD.GridFunction(grid, 1 ./ (1 .+ τ .* FD.pw_norm(q))), FD.GridFunction(grid, p .+ τ .* q)))
 
   reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
   (p_new, reserr)
 end
 
-function fpalg(grid, f, A; α, β, τ = 1/8, ϵ = 1e-4)
+function fpalg(grid, w, A; α, β, τ = 1/8, ϵ = 1e-4)
   local uimg
-  p = FDM.GridFunction(x -> 0, ndims(grid), grid)
+  p = FD.GridFunction(x -> 0, ndims(grid), grid)
   reserr = 1
   k = 0
   while reserr > ϵ && k < 30e1
     k += 1
-    #sleep(0.1 / k)
-    (p, reserr) = fpiter(grid, p, f, A, α, β, τ)
+    (p, reserr) = fpiter(grid, p, w, A, α, β, τ)
 
-    u = p2u(grid, p, A, f, α=α, β=β)
-
-    #FDM.my_plot(uimg)
-    println("[$(lpad(k, 4))] res = $(round(reserr, digits=5))")
-    println("min/max: $(minimum(uimg)) / $(maximum(uimg))")
-    println("min/max: $(minimum(p)) / $(maximum(p))")
+    #FD.my_plot(uimg)
+    #println("[$(lpad(k, 4))] res = $(round(reserr, digits=5))")
+    #println("min/max: $(minimum(uimg)) / $(maximum(uimg))")
+    #println("min/max: $(minimum(p)) / $(maximum(p))")
   end
   p
   #uimg
 end
 
-grid = FDM.MetaGrid((0:0.01:1, 0:0.01:1), (3, 3), 5)
+grid = FD.MetaGrid((0:0.01:1, 0:0.01:1), (3, 3), 5)
 
 A = I
 #A = spdiagm(0 => ones(length(grid)), 1 => ones(length(grid)-1) )
 #f = A * f
 
 
-f = FDM.GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
-g = FDM.GridFunction(x -> norm(x), 2, grid)
+f = FD.GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
+g = FD.GridFunction(x -> norm(x), 2, grid)
+
+α=1
+β=1
+τ=1/100
+ϵ=1e-7
+ρ=0.25
+
+p = FD.GridFunction(x -> 0, 2, grid)
+
+function dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
+  p .= (1-ρ) .* p
 
+  for part in FD.parts(grid)
+    gridv = view(grid, Tuple(part)...)
+    fv = view(f, Tuple(part)...)
+    #p.data .= 0
+    pv = view(p, Tuple(part)...)
 
-for part in parts(grid)
-  gridv = view(grid, part...)
-  fv = view(f, part...)
-  fpalg(gridv, fv, A, α=1, β=1, τ=1/100, ϵ=1e-7)
+    θ = theta(grid, Tuple(part)...)
+    θv = view(θ, Tuple(part)...)
+    w = FD.GridFunction(grid, A' * f.data[:] - FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p)))
+    wv = view(w, Tuple(part)...)
+
+    #wv = A' * fv
+
+    pv .+= ρ .* fpalg(gridv, wv, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
+
+    #u = p2u(grid, p, A, f, α=α*θ, β=β)
+    u = p2u(grid, p, A, f, α=α, β=β)
+    uv = view(u, Tuple(part)...)
+
+    my_plot(θv)
+    #sleep(1)
+    #my_plot(FD.pw_norm(pv))
+    @show maximum(pv), minimum(pv)
+    #sleep(1)
+    #my_plot(θv)
+    return
+  end
 end
 
+dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
+
+#for n = 1:10
+#  println("iteration $n")
+#  dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
+#end
+
+
+
+
 #gridv = view(grid, 1, 1)
 #fv = view(f, 1, 1)
 #g = view(g, 1, 1)
 
 
-FDM.my_plot(f)
+#FD.my_plot(f)
 
 
 
 
-#FDM.plot(f)
+#FD.plot(f)
 
 
 

fd.jl

@@ -210,10 +210,13 @@ Base.setindex!(f::GridFunction, v, idx::Int) = setindex!(f.data, v, idx)
 #component(f::GridFunction, k::Integer) = view(reshape(f.data,
 
 Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
-    GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)...))
+    GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
 
+grid(f::GridFunction) = f.grid
 
 
+Base.:*(a::Number, b::GridFunction) = GridFunction(b.grid, a * b.data, Val(rangedim(b)))
+
 
 
 import Colors: Gray
@@ -263,11 +266,25 @@ function divergence(f::GridFunction{G,T,d}; fd=FD.backward) where {d,G<:Abstract
   GridFunction(f.grid, A)
 end
 
+function pw_norm(f::GridFunction)
+  d = domaindim(f)
+  GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(1))
+end
+
 function my_norm(f::GridFunction)
   d = domaindim(f)
   GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(rangedim(f)))
 end
 
+function pwmul(f::GridFunction{G,T,d,1}, g::GridFunction{G,T,d}) where {G,T,d}
+  # TODO: maybe comprehension?
+  data = similar(g.data)
+  for I in CartesianIndices(f.grid)
+    data[Tuple(I)..., :] = f.data[Tuple(I)..., 1] * g.data[Tuple(I)..., :]
+  end
+  GridFunction(f.grid, data, Val(rangedim(g)))
+end
+