Skip to content

Commit

Permalink
move binnedpointlist to ExtendableGrids
Browse files Browse the repository at this point in the history
  • Loading branch information
j-fu committed Mar 21, 2024
1 parent bbddcd7 commit 55a1c66
Show file tree
Hide file tree
Showing 9 changed files with 328 additions and 28 deletions.
14 changes: 7 additions & 7 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "ExtendableGrids"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
authors = ["Juergen Fuhrmann <[email protected]>", "Christian Merdon <[email protected]>", "Johannes Taraz <[email protected]>"]
version = "1.3.2"
version = "1.4.0"

[deps]
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
Expand All @@ -20,6 +20,12 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"

[weakdeps]
Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"

[extensions]
ExtendableGridsGmshExt = "Gmsh"

[compat]
AbstractTrees = "0.3,0.4"
Bijections = "0.1.4"
Expand All @@ -39,11 +45,5 @@ UUIDs = "1.6"
WriteVTK = "1.14"
julia = "1.6"

[extensions]
ExtendableGridsGmshExt = "Gmsh"

[extras]
Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"

[weakdeps]
Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ function mkdocs()
"refinement.md",
"regionedit.md",
"tokenstream.md",
"binnedpointlist.md",
"gmsh.md",
"allindex.md",
"Examples" => generated_examples,
Expand Down
8 changes: 8 additions & 0 deletions docs/src/binnedpointlist.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
# BinnedPointList
Used to find and identify points in space

## API
```@autodocs
Modules = [ExtendableGrids]
Pages = ["binnedpointlist.jl"]
```
3 changes: 3 additions & 0 deletions src/ExtendableGrids.jl
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,9 @@ export cellmask!, bfacemask!, bedgemask!, rect!
include("arraytools.jl")
export glue, geomspace, linspace

include("binnedpointlist.jl")
export BinnedPointList

include("simplexgrid.jl")
export simplexgrid, geomspace, glue
export XCoordinates, YCoordinates, ZCoordinates
Expand Down
236 changes: 236 additions & 0 deletions src/binnedpointlist.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,236 @@
"""
$(TYPEDEF)
Binned point list structure allowing for fast
check for already existing points.
This provides better performance for indendifying
already inserted points than the naive linear search.
OTOH the implementation is still quite naive - it dynamically maintains
a cuboid bining region with a fixed number of bins.
Probably tree based adaptive methods (a la octree) will be more efficient,
however they will be harder to implement.
In an ideal world, we would maintain a dynamic Delaunay triangulation, which
at once could be the starting point of mesh generation which will follow here
anyway.
"""
mutable struct BinnedPointList{T}

# Space dimension
dim::Int32

# Point distance tolerance. Points closer than tol
# (in Euclidean distance) will be identified, i.e.
# are collapsed to the first inserted.
tol::T

# The union of all bins is the binning region -
# a cuboid given by two of its corners. It is calculated
# dynamically depending on the inserted points.
binning_region_min::Vector{T}
binning_region_max::Vector{T}

# Increase factor of binning region (with respect to
# the cuboid defined by the coordinates of the binned points)
binning_region_increase_factor::T

# The actual point list
points::ElasticArray{T, 2}

# The bins are vectors of indices of points in the point list
# We store them in a dim-dimensional array of length number_of_directional_bins^dim
bins::Array{Vector{Int32}}

# Number of bins in each space dimension
number_of_directional_bins::Int32

# Some points will fall outside of the binning region.
# We collect them in vector of ubinned point indices
unbinned::Vector{Int32}

# Number of unbinned points tolerated without rebinning
num_allowed_unbinned_points::Int32

# Maximum ratio of unbinned points in point list
max_unbinned_ratio::T

# Storage of current point bin
current_bin::Vector{Int32}

BinnedPointList{T}(::Nothing) where {T} = new()
end

"""
$(SIGNATURES)
Create and initialize binned point list
"""
function BinnedPointList(::Type{T}, dim;
tol = 1.0e-12,
number_of_directional_bins = 10,
binning_region_increase_factor = 0.01,
num_allowed_unbinned_points = 5,
max_unbinned_ratio = 0.05) where {T}
bpl = BinnedPointList{T}(nothing)

bpl.dim = dim
bpl.tol = tol
bpl.number_of_directional_bins = number_of_directional_bins
bpl.binning_region_increase_factor = binning_region_increase_factor
bpl.num_allowed_unbinned_points = num_allowed_unbinned_points
bpl.max_unbinned_ratio = max_unbinned_ratio

bpl.points = ElasticArray{Cdouble}(undef, dim, 0)

bpl.binning_region_min = fill(floatmax(T), dim)
bpl.binning_region_max = fill(-floatmax(T), dim)

bpl.unbinned = Vector{Int32}(undef, 0)
bpl.bins = Array{Vector{Int32}, dim}(undef, zeros(Int32, dim)...)

bpl.current_bin = zeros(Int32, bpl.dim)
bpl
end

"""
$(SIGNATURES)
Create and initialize binned point list
"""
BinnedPointList(dim; kwargs...) = BinnedPointList(Cdouble, dim; kwargs...)

#
# Find point in index list (by linear search)
# Return its index, or zero if not found
#
function _findpoint(bpl, index, p)
for i = 1:length(index)
@views if norm(bpl.points[:, index[i]] - p) < bpl.tol
return index[i]
end
end
return 0
end

#
# Calculate the bin of the point. Result is stored
# in bpl.current_bin
#
function _bin_of_point!(bpl, p)
for idim = 1:(bpl.dim)
# scaled value for particular dimension
s = (p[idim] - bpl.binning_region_min[idim]) / (bpl.binning_region_max[idim] - bpl.binning_region_min[idim])
if s > 1 || s < 0
# point lies outside of binning area
bpl.current_bin[idim] = 0
else
# calculate bin in particular dimension
bpl.current_bin[idim] = ceil(Int32, bpl.number_of_directional_bins * s)
end
end
end

#
# Re-calculate binning if there are too many unbinned points
# This amounts to two steps:
# - Enlarge binning area in order to include all points
# - Re-calculate all point bins
function _rebin_all_points!(bpl)
if length(bpl.unbinned) > max(bpl.num_allowed_unbinned_points,
bpl.max_unbinned_ratio * size(bpl.points, 2))

# Calculate extrema of unbinned points
@views e = extrema(bpl.points[:, bpl.unbinned]; dims = 2)

for i = 1:(bpl.dim)
# Increase binning region according to unbinned extrema
bpl.binning_region_min[i] = min(bpl.binning_region_min[i], e[i][1])
bpl.binning_region_max[i] = max(bpl.binning_region_max[i], e[i][2])

# Slightly increase binning region further in order to
# include all existing points with tolerance
delta = max(bpl.binning_region_max[i] - bpl.binning_region_min[i], bpl.tol)
bpl.binning_region_min[i] -= bpl.binning_region_increase_factor * delta
bpl.binning_region_max[i] += bpl.binning_region_increase_factor * delta
end

# Re-allocate all bins
if bpl.dim == 1
bpl.bins = [zeros(Int32, 0) for i = 1:(bpl.number_of_directional_bins)]
elseif bpl.dim == 2
bpl.bins = [zeros(Int32, 0) for i = 1:(bpl.number_of_directional_bins), j = 1:(bpl.number_of_directional_bins)]
elseif bpl.dim == 3
bpl.bins = [zeros(Int32, 0)
for i = 1:(bpl.number_of_directional_bins), j = 1:(bpl.number_of_directional_bins),
k = 1:(bpl.number_of_directional_bins)]
end

# Register all points in their respctive bins
for i = 1:size(bpl.points, 2)
@views _bin_of_point!(bpl, bpl.points[:, i])
push!(bpl.bins[bpl.current_bin...], i)
end

# Re-allocate unbinned index
bpl.unbinned = zeros(Int32, 0)
end
end

"""
$(SIGNATURES)
If another point with distance less the tol from p is
in pointlist, return its index. Otherwise, insert point into pointlist.
"""
function Base.insert!(bpl::BinnedPointList{T}, p) where {T}
_rebin_all_points!(bpl)
_bin_of_point!(bpl, p)
if reduce(*, bpl.current_bin) > 0
i = _findpoint(bpl, bpl.bins[bpl.current_bin...], p)
if i > 0
return i
else
append!(bpl.points, p)
i = size(bpl.points, 2)
push!(bpl.bins[bpl.current_bin...], i)
return i
end
else
i = _findpoint(bpl, bpl.unbinned, p)
if i > 0
return i
else
append!(bpl.points, p)
i = size(bpl.points, 2)
push!(bpl.unbinned, i)
return i
end
end
end

Base.insert!(bpl::BinnedPointList{T}, x::Number) where {T} = insert!(bpl, (x))
Base.insert!(bpl::BinnedPointList{T}, x::Number, y::Number) where {T} = insert!(bpl, (x, y))
Base.insert!(bpl::BinnedPointList{T}, x::Number, y::Number, z::Number) where {T} = insert!(bpl, (x, y, z))

"""
$(SIGNATURES)
Return the array of points in the point list.
"""
points(bpl::BinnedPointList{T}) where {T} = bpl.points

# Just for being able to check of all of the above was worth the effort...
function naiveinsert!(bpl::BinnedPointList{T}, p) where {T}
for i = 1:size(bpl.points, 2)
@views if norm(bpl.points[:, i] - p) < bpl.tol
return i
end
end
append!(bpl.points, p)
size(bpl.points, 2)
end
1 change: 1 addition & 0 deletions test/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"

[compat]
Expand Down
43 changes: 22 additions & 21 deletions test/gmsh.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,48 +43,48 @@ using Gmsh: gmsh
end

@testset "Read/write simplex gmsh 2d / 3d" begin
path = ""
path = joinpath(pkgdir(ExtendableGrids),"test")

X = collect(0:0.02:2)
Y = collect(0:2:4)
grid1 = simplexgrid(X, Y) #ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh")
ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh")
ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = joinpath(path,"testfile.msh"))

grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float32, Ti = Int64)
grid2 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"testfile.msh"); Tc = Float32, Ti = Int64)

@test seemingly_equal(grid2, grid1; sort = true, confidence = :low)
@test seemingly_equal(grid2, grid1; sort = true, confidence = :full)

grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_2d.msh"; Tc = Float64, Ti = Int64)
grid1 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"sto_2d.msh"); Tc = Float64, Ti = Int64)

ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh")
grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float64, Ti = Int64)
ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = joinpath(path,"testfile.msh"))
grid2 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"testfile.msh"); Tc = Float64, Ti = Int64)

@test seemingly_equal(grid1, grid2; sort = true, confidence = :low)
@test seemingly_equal(grid1, grid2; sort = true, confidence = :full)

grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_3d.msh"; Tc = Float32, Ti = Int64)
grid1 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"sto_3d.msh"); Tc = Float32, Ti = Int64)

ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh")
grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float64, Ti = Int32)
ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = joinpath(path,"testfile.msh"))
grid2 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"testfile.msh"); Tc = Float64, Ti = Int32)

@test seemingly_equal(grid1, grid2; sort = true, confidence = :low)
@test seemingly_equal(grid1, grid2; sort = true, confidence = :full)

grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_2d.msh")
grid1 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"sto_2d.msh"))

grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_3d.msh"; Tc = Float32, Ti = Int32)
grid2 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"sto_3d.msh"); Tc = Float32, Ti = Int32)

@test !seemingly_equal(grid1, grid2; sort = true, confidence = :low)
@test !seemingly_equal(grid1, grid2; sort = true, confidence = :full)

grid1 = ExtendableGrids.simplexgrid_from_gmsh("testmesh.gmsh"; incomplete = true)
grid1 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"testmesh.gmsh"); incomplete = true)
ExtendableGrids.seal!(grid1; encode = false)

ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = "completed_testfile.msh")
grid2 = ExtendableGrids.simplexgrid_from_gmsh("completed_testfile.msh")
ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = joinpath(path,"completed_testfile.msh"))
grid2 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"completed_testfile.msh"))

grid3 = ExtendableGrids.simplexgrid_from_gmsh("testmesh.gmsh"; incomplete = true)
grid3 = ExtendableGrids.simplexgrid_from_gmsh(joinpath(path,"testmesh.gmsh"); incomplete = true)
ExtendableGrids.seal!(grid3; encode = true)

@test seemingly_equal(grid1, grid2; sort = true, confidence = :low)
Expand All @@ -107,19 +107,20 @@ end
end

@testset "Read/write mixed gmsh 2d" begin
path = ""
grid1 = ExtendableGrids.mixedgrid_from_gmsh(path * "mixedgrid_2d.msh"; Tc = Float64, Ti = Int64)
path = joinpath(pkgdir(ExtendableGrids),"test")
grid1 = ExtendableGrids.mixedgrid_from_gmsh(joinpath(path,"mixedgrid_2d.msh"); Tc = Float64, Ti = Int64)

ExtendableGrids.mixedgrid_to_gmsh(grid1; filename = path * "testfile.msh")
grid2 = ExtendableGrids.mixedgrid_from_gmsh(path * "testfile.msh"; Tc = Float32, Ti = UInt64)
ExtendableGrids.mixedgrid_to_gmsh(grid1; filename = joinpath(path, "testfile.msh"))
grid2 = ExtendableGrids.mixedgrid_from_gmsh(joinpath(path,"testfile.msh"); Tc = Float32, Ti = UInt64)

@test seemingly_equal(grid1, grid2; sort = true, confidence = :low)
@test seemingly_equal(grid1, grid2; sort = true, confidence = :full)
end

@testset "Read .geo files" begin
grid = simplexgrid("disk1hole.geo")
path = joinpath(pkgdir(ExtendableGrids),"test")
grid = simplexgrid(joinpath(path,"disk1hole.geo"))
@test num_cells(grid) > 0
grid = simplexgrid("cube6.geo")
grid = simplexgrid(joinpath(path,"cube6.geo"))
@test num_cells(grid) > 0
end
Loading

0 comments on commit 55a1c66

Please sign in to comment.