Mercurial > repos > public > sbplib_julia

comparison StencilIndex.jl @ 94:84b1ad5a3755 stencil_index

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Made everything work(?) but also go really slow. And also not type-stable.
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 04 Feb 2019 16:09:07 +0100
parents 93df72e2b135
children
comparison
equal deleted inserted replaced
93:93df72e2b135 94:84b1ad5a3755
1 abstract type StencilIndex end 1 abstract type Region end
2 struct Interior <: Region end
3 struct Lower <: Region end
4 struct Upper <: Region end
5
6 struct StencilIndex{R<:Region, T<:Integer}
7 localindex::CartesianIndex
8 globalindex::T
9
10 StencilIndex{R}(li::CartesianIndex, gi::T) where {R<:Region,T<:Integer} = new{R, T}(li, gi)
11 StencilIndex(li::CartesianIndex, gi::T, ::Type{R}) where {R<:Region,T<:Integer} = StencilIndex{R}(li, gi)
12
13 # Index(t::Tuple{T, Type{R}}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
14 # Above doesn't work, below does but is less type strict
15 #Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
16 end
2 17
3 function Base.getindex(si::StencilIndex, i::Int) 18 function Base.getindex(si::StencilIndex, i::Int)
4 return si.gridindex[i] 19 return si.localindex[i]
5 end 20 end
6 21
7 struct LowerClosureIndex <: StencilIndex 22 #Index(t::Vararg{Tuple{T, DataType}}) where T = Index.(t)
8 globalindex::Integer 23 # TODO: Where to place this function?
9 gridindex::CartesianIndex
10 end
11 24
12 struct UpperClosureIndex <: StencilIndex 25 function stencilindices(diffOp)
13 globalindex::Integer 26 N = diffOp.grid.numberOfPointsPerDim
14 gridindex::CartesianIndex
15 end
16 27
17 struct InteriorIndex <: StencilIndex 28 lowerclosure = Vector{Vector{StencilIndex{Lower, Int64}}}(undef,0)
18 globalindex::Integer 29 upperclosure = Vector{Vector{StencilIndex{Upper, Int64}}}(undef,0)
19 gridindex::CartesianIndex 30 interior = Vector{Vector{StencilIndex{Interior, Int64}}}(undef,0)
20 end 31 cSize = closureSize(diffOp.op)
32 ci = CartesianIndices(diffOp.grid.numberOfPointsPerDim)
21 33
22 # TODO: The design of StencilIndex is wrong. Use Jonatans design instead. 34 # TODO: Loop over all points or one loop for each region?
23 # TODO: This should take a Stencil or DiffOp so that we can extract all the 35 for j = 1:Grid.numberOfDimensions(diffOp.grid)
24 # indices in the closures. 36 templ = Vector{StencilIndex{Lower,Int64}}(undef, 0)
25 # TODO: Where to place this function? 37 tempu = Vector{StencilIndex{Upper,Int64}}(undef, 0)
26 function stencilindices(grid::Grid.EquidistantGrid) 38 tempi = Vector{StencilIndex{Interior,Int64}}(undef, 0)
27 lowerclosure = Vector{LowerClosureIndex}(undef, 0) 39 for i ∈ 1:Grid.numberOfPoints(diffOp.grid)
28 upperclosure = Vector{UpperClosureIndex}(undef, 0) 40 val = ci[i][j]
29 interior = Vector{InteriorIndex}(undef, 0) 41 if val ∈ range(1; length=cSize)
30 # TODO: Fix such that the indices of the entire closure width is included. 42 push!(templ, StencilIndex{Lower}(ci[i],i))
31 islower = x -> (x == 1) 43 elseif val ∈ range(N[j] - cSize+1, length=cSize)
32 isupper = x -> (x in grid.numberOfPointsPerDim) 44 push!(tempu, StencilIndex{Upper}(ci[i],i))
33 ci = CartesianIndices(grid.numberOfPointsPerDim) 45 else
34 for i ∈ 1:Grid.numberOfPoints(grid) 46 push!(tempi, StencilIndex{Interior}(ci[i],i))
35 I = Tuple(ci[i]) 47 end
36 if any(islower, I)
37 push!(lowerclosure, LowerClosureIndex(i,ci[i]))
38 elseif any(isupper, I)
39 push!(upperclosure, UpperClosureIndex(i,ci[i]))
40 else
41 push!(interior, InteriorIndex(i,ci[i]))
42 end 48 end
49 push!(lowerclosure,templ)
50 push!(upperclosure,tempu)
51 push!(interior,tempi)
43 end 52 end
44 return lowerclosure, upperclosure, interior 53 return lowerclosure, upperclosure, interior
45 end 54 end
55