From 55a1c66cfc16a910e05f6bb4f98e8af4c86502a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Thu, 21 Mar 2024 17:12:02 +0100 Subject: [PATCH 1/3] move binnedpointlist to ExtendableGrids --- Project.toml | 14 +-- docs/make.jl | 1 + docs/src/binnedpointlist.md | 8 ++ src/ExtendableGrids.jl | 3 + src/binnedpointlist.jl | 236 ++++++++++++++++++++++++++++++++++++ test/Project.toml | 1 + test/gmsh.jl | 43 +++---- test/runtests.jl | 7 ++ test/testbinned.jl | 43 +++++++ 9 files changed, 328 insertions(+), 28 deletions(-) create mode 100644 docs/src/binnedpointlist.md create mode 100644 src/binnedpointlist.jl create mode 100644 test/testbinned.jl diff --git a/Project.toml b/Project.toml index 0dc45324..e9ccc241 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableGrids" uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" authors = ["Juergen Fuhrmann ", "Christian Merdon ", "Johannes Taraz "] -version = "1.3.2" +version = "1.4.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -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" @@ -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" diff --git a/docs/make.jl b/docs/make.jl index 64246621..02eadcad 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,6 +34,7 @@ function mkdocs() "refinement.md", "regionedit.md", "tokenstream.md", + "binnedpointlist.md", "gmsh.md", "allindex.md", "Examples" => generated_examples, diff --git a/docs/src/binnedpointlist.md b/docs/src/binnedpointlist.md new file mode 100644 index 00000000..2b4b1ac9 --- /dev/null +++ b/docs/src/binnedpointlist.md @@ -0,0 +1,8 @@ +# BinnedPointList +Used to find and identify points in space + +## API +```@autodocs +Modules = [ExtendableGrids] +Pages = ["binnedpointlist.jl"] +``` diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl index 1cb082d1..f7edfb79 100644 --- a/src/ExtendableGrids.jl +++ b/src/ExtendableGrids.jl @@ -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 diff --git a/src/binnedpointlist.jl b/src/binnedpointlist.jl new file mode 100644 index 00000000..7ce702ec --- /dev/null +++ b/src/binnedpointlist.jl @@ -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 diff --git a/test/Project.toml b/test/Project.toml index 5023a824..956a752d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -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] diff --git a/test/gmsh.jl b/test/gmsh.jl index e10825d5..6e565dfd 100644 --- a/test/gmsh.jl +++ b/test/gmsh.jl @@ -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) @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 11786b49..2047a5ec 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,8 @@ using ExtendableGrids, SHA using AbstractTrees, StatsBase +using TetGen, SimplexGridFactory + @testset "Aqua" begin Aqua.test_ambiguities([ExtendableGrids, Base, Core], exclude=[view, ==, StatsBase.TestStat, copyto!]) Aqua.test_unbound_args(ExtendableGrids) @@ -239,6 +241,11 @@ end end + +@testset " BinnedPointList" begin + include("testbinned.jl") +end + function tglue(; dim = 2, breg = 0) X1 = linspace(0, 1, 5) Y1 = linspace(0, 1, 5) diff --git a/test/testbinned.jl b/test/testbinned.jl new file mode 100644 index 00000000..6163a333 --- /dev/null +++ b/test/testbinned.jl @@ -0,0 +1,43 @@ + function testbinning(a) + dim = size(a, 1) + n = size(a, 2) + idx = rand(1:n, n ÷ 4) + a1 = a[:, idx] + + bpl = BinnedPointList(dim) + for i = 1:size(a, 2) + insert!(bpl, a[:, i]) + end + + for i = 1:size(a1, 2) + ix = insert!(bpl, a1[:, i]) + if ix != idx[i] + return false + end + end + true + end + function testbinning2() + A = 50 + B = 100 + Z = 20 + d = 5 + builder = SimplexGridBuilder(; Generator = TetGen) + p1 = point!(builder, 0, 0, 0) + p2 = point!(builder, B, 0, 0) + p3 = point!(builder, B, A, 0) + p4 = point!(builder, B, A + d, 0) + p5 = point!(builder, B, 2 * A + d, 0) + p6 = point!(builder, 0, 2 * A + d, 0) + p7 = point!(builder, 0, A + d, 0) + p8 = point!(builder, 0, A, 0) + end + + @test testbinning(rand(1, 10)) + @test testbinning(rand(1, 10000)) + @test testbinning(rand(2, 10)) + @test testbinning(rand(2, 10000)) + @test testbinning(rand(3, 10)) + @test testbinning(rand(3, 10000)) + + @test testbinning2() == 8 From 5825ca6e8ce6dc570da0dcc191008b6c8c9e9092 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Thu, 21 Mar 2024 20:52:00 +0100 Subject: [PATCH 2/3] new grid glueing algorithm using BinnedPointList --- src/binnedpointlist.jl | 10 ++++++ src/simplexgrid.jl | 76 ++++++++++++++++++++++++++++++++++++------ 2 files changed, 76 insertions(+), 10 deletions(-) diff --git a/src/binnedpointlist.jl b/src/binnedpointlist.jl index 7ce702ec..50f99e10 100644 --- a/src/binnedpointlist.jl +++ b/src/binnedpointlist.jl @@ -213,6 +213,16 @@ function Base.insert!(bpl::BinnedPointList{T}, p) where {T} end end +function findpoint(bpl::BinnedPointList{T}, p) where {T} + _rebin_all_points!(bpl) + _bin_of_point!(bpl, p) + if reduce(*, bpl.current_bin) > 0 + return _findpoint(bpl, bpl.bins[bpl.current_bin...], p) + else + return _findpoint(bpl, bpl.unbinned, p) + 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)) diff --git a/src/simplexgrid.jl b/src/simplexgrid.jl index 5092923e..571e7f94 100644 --- a/src/simplexgrid.jl +++ b/src/simplexgrid.jl @@ -10,7 +10,7 @@ function simplexgrid(coord::Array{Tc,2}, Create d-dimensional simplex grid from five arrays. -- coord: d ``\\times`` n_points matrix of coordinates +- coord: d ``\\times`` n_points matrix of coordinates - cellnodes: d+1 ``\\times`` n_tri matrix of triangle - point incidence - cellregions: n_tri vector of cell region markers - bfacenodes: d ``\\times`` n_bf matrix of boundary facet - point incidences @@ -749,13 +749,72 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid; end end + function barycenter!(bc,bfnodes,coord,iface) + for idim=1:dim + bc[idim]=0.0 + for jdim=1:dim + bc[idim]+=coord[idim,bfnodes[jdim,iface]]/dim + end + end + bc + end + + nbfreg1=num_bfaceregions(g1) + nbfreg2=num_bfaceregions(g2) + breg1used=zeros(Bool,nbfreg1) + breg2used=zeros(Bool,nbfreg2) + for reg in g1regions + breg1used[reg]=true + end + for reg in g2regions + breg2used[reg]=true + end - use_region(ireg,regionlist) = findfirst(i->i==ireg,regionlist)!=nothing - # Run over all pairs of boundary faces and try to match them - for ibf1=1:nbf1 - if use_region(bfreg1[ibf1],g1regions) - for ibf2=1:nbf2 - if use_region(bfreg2[ibf2],g2regions) + if false + # Run over all pairs of boundary faces and try to match them + # This can scale catastrophically... + for ibf1=1:nbf1 + if breg1used[bfreg1[ibf1]] + for ibf2=1:nbf2 + if breg2used[bfreg2[ibf2]] + match_faces(ibf1,ibf2) + end + end + end + end + else + # Use a binned point list for speed up. + # Facets are found through their barycenters. + bcenter=zeros(dim) + bpl=BinnedPointList(eltype(coord1),dim) + nbf1used=0 + + for ibf1=1:nbf1 + if breg1used[bfreg1[ibf1]] + nbf1used+=1 + end + end + + # marker for facet numbers + bf1used=zeros(Int,nbf1used) + + # insert all barycenters of used bfaces from grid 1 into the binned pointlist + for ibf1=1:nbf1 + if breg1used[bfreg1[ibf1]] + i=insert!(bpl,barycenter!(bcenter,bfn1,coord1,ibf1)) + bf1used[i]=ibf1 + end + end + + # loop over used grid2 bfaces and find if barycenters match with a point in + # the binned pointlist + for ibf2=1:nbf2 + if breg2used[bfreg2[ibf2]] + i=findpoint(bpl,barycenter!(bcenter,bfn2,coord2,ibf2)) + # if a match has been found, identify facets + if i>0 + ibf1=bf1used[i] + # use "old" matching, double check at once match_faces(ibf1,ibf2) end end @@ -763,9 +822,6 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid; end @info "glue: matches found: nodes $(n_matching_nodes) bfaces: $(n_matching_faces)" - - - creg1=g1[CellRegions] creg2=g2[CellRegions] From acf1c1bb9f4de07bb1f2070bce3dc0236f181800 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Thu, 21 Mar 2024 22:37:36 +0100 Subject: [PATCH 3/3] fix docs, test --- docs/src/binnedpointlist.md | 17 +++- src/ExtendableGrids.jl | 2 +- src/binnedpointlist.jl | 182 ++++++++++++++++++++++++------------ src/simplexgrid.jl | 14 +-- test/Project.toml | 1 - test/runtests.jl | 2 - test/testbinned.jl | 62 +++++------- 7 files changed, 167 insertions(+), 113 deletions(-) diff --git a/docs/src/binnedpointlist.md b/docs/src/binnedpointlist.md index 2b4b1ac9..f28eb97e 100644 --- a/docs/src/binnedpointlist.md +++ b/docs/src/binnedpointlist.md @@ -2,7 +2,18 @@ Used to find and identify points in space ## API -```@autodocs -Modules = [ExtendableGrids] -Pages = ["binnedpointlist.jl"] + +```@docs +BinnedPointList +findpoint +Base.insert! +``` + +## Internal + +```@docs +ExtendableGrids._findpoint +ExtendableGrids._bin_of_point! +ExtendableGrids._rebin_all_points! +ExtendableGrids.naiveinsert! ``` diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl index f7edfb79..8949c0d5 100644 --- a/src/ExtendableGrids.jl +++ b/src/ExtendableGrids.jl @@ -164,7 +164,7 @@ include("arraytools.jl") export glue, geomspace, linspace include("binnedpointlist.jl") -export BinnedPointList +export BinnedPointList, findpoint include("simplexgrid.jl") export simplexgrid, geomspace, glue diff --git a/src/binnedpointlist.jl b/src/binnedpointlist.jl index 50f99e10..ddaa78de 100644 --- a/src/binnedpointlist.jl +++ b/src/binnedpointlist.jl @@ -2,70 +2,99 @@ $(TYPEDEF) -Binned point list structure allowing for fast -check for already existing points. +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. +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. +OTOH the implementation is still quite naive - it dynamically +maintains a cuboid binning 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. +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. +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. + +$(TYPEDFIELDS) """ mutable struct BinnedPointList{T} - # Space dimension + """ + 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. + """ + 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. + """" + 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) + """ + 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 + """ + 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 + """ + 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 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 + """ + 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 + """ + Number of unbinned points tolerated without rebinning + """ num_allowed_unbinned_points::Int32 - # Maximum ratio of unbinned points in point list + """ + Maximum ratio of unbinned points in point list + """ max_unbinned_ratio::T - # Storage of current point bin + """ + Storage of current point bin + """ current_bin::Vector{Int32} BinnedPointList{T}(::Nothing) where {T} = new() end """ - $(SIGNATURES) + 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} Create and initialize binned point list """ @@ -97,17 +126,19 @@ function BinnedPointList(::Type{T}, dim; end """ -$(SIGNATURES) + BinnedPointList(dim; kwargs...) Create and initialize binned point list """ -BinnedPointList(dim; kwargs...) = BinnedPointList(Cdouble, dim; kwargs...) +BinnedPointList(dim; kwargs...) = BinnedPointList(Float64, dim; kwargs...) + +""" + _findpoint(binnedpointlist, index, p) -# -# Find point in index list (by linear search) -# Return its index, or zero if not found -# +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 @@ -117,10 +148,12 @@ function _findpoint(bpl, index, p) return 0 end -# -# Calculate the bin of the point. Result is stored -# in bpl.current_bin -# +""" + _bin_of_point!(binnedpointlist, p) + +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 @@ -135,11 +168,14 @@ function _bin_of_point!(bpl, p) 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 +""" + _rebin_all_points!(bpl) + +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)) @@ -181,11 +217,29 @@ function _rebin_all_points!(bpl) end end + +""" + findpoint(binnedpointlist, p) + +Find point in binned point list. Return its index in the point list if found, +otherwise return 0. +""" +function findpoint(bpl::BinnedPointList{T}, p) where {T} + _rebin_all_points!(bpl) + _bin_of_point!(bpl, p) + if reduce(*, bpl.current_bin) > 0 + return _findpoint(bpl, bpl.bins[bpl.current_bin...], p) + else + return _findpoint(bpl, bpl.unbinned, p) + end +end + """ -$(SIGNATURES) + Base.insert!(binnedpointlist,p) If another point with distance less the tol from p is in pointlist, return its index. Otherwise, insert point into pointlist. +`p` may be a vector or a tuple. """ function Base.insert!(bpl::BinnedPointList{T}, p) where {T} _rebin_all_points!(bpl) @@ -213,28 +267,34 @@ function Base.insert!(bpl::BinnedPointList{T}, p) where {T} end end -function findpoint(bpl::BinnedPointList{T}, p) where {T} - _rebin_all_points!(bpl) - _bin_of_point!(bpl, p) - if reduce(*, bpl.current_bin) > 0 - return _findpoint(bpl, bpl.bins[bpl.current_bin...], p) - else - return _findpoint(bpl, bpl.unbinned, p) - end -end - +""" + Base.insert!(binnedpointlist,x) + +Insert 1D point via coordinate. +""" Base.insert!(bpl::BinnedPointList{T}, x::Number) where {T} = insert!(bpl, (x)) + +""" + Base.insert!(binnedpointlist,x,y) + +Insert 2D point via coordinates. +""" + Base.insert!(bpl::BinnedPointList{T}, x::Number, y::Number) where {T} = insert!(bpl, (x, y)) +""" + Base.insert!(binnedpointlist,x,y,z) + +Insert 3D point via coordinates. +""" 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 + naiveinsert(binnedpointlist, p) -# Just for being able to check of all of the above was worth the effort... +Insert via linear search, without any binning. +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 diff --git a/src/simplexgrid.jl b/src/simplexgrid.jl index 571e7f94..79c84378 100644 --- a/src/simplexgrid.jl +++ b/src/simplexgrid.jl @@ -651,7 +651,8 @@ end g1regions=1:num_bfaceregions(g1), g2regions=1:num_bfaceregions(g2), interface=0, - tol=1.0e-10) + tol=1.0e-10, + naive=false) Merge two grids along their common boundary facets. @@ -661,7 +662,7 @@ Merge two grids along their common boundary facets. - g2regions: boundary regions to be used from grid2. Default: all. - interface: if nonzero, create interface region in new grid, otherwise, ignore - tol: Distance below which two points are seen as identical. Default: 1.0e-10 - +- naive: use naive quadratic complexity matching (for checking backward compatibility). Default: false Deprecated: - breg: old notation for interface @@ -671,8 +672,9 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid; g2regions=1:num_bfaceregions(g2), breg=nothing, interface=0, - tol=1.0e-10) - Ti=Cint + tol=1.0e-10, + naive = false) + Ti=eltype(g1[CellNodes]) dim=dim_space(g1) if breg!=nothing @@ -770,7 +772,7 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid; breg2used[reg]=true end - if false + if naive # Run over all pairs of boundary faces and try to match them # This can scale catastrophically... for ibf1=1:nbf1 @@ -821,7 +823,7 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid; end end - @info "glue: matches found: nodes $(n_matching_nodes) bfaces: $(n_matching_faces)" + @info "glue: $(n_matching_faces) matching bfaces found" creg1=g1[CellRegions] creg2=g2[CellRegions] diff --git a/test/Project.toml b/test/Project.toml index 956a752d..5023a824 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,7 +8,6 @@ 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] diff --git a/test/runtests.jl b/test/runtests.jl index 2047a5ec..32942497 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,8 +5,6 @@ using ExtendableGrids, SHA using AbstractTrees, StatsBase -using TetGen, SimplexGridFactory - @testset "Aqua" begin Aqua.test_ambiguities([ExtendableGrids, Base, Core], exclude=[view, ==, StatsBase.TestStat, copyto!]) Aqua.test_unbound_args(ExtendableGrids) diff --git a/test/testbinned.jl b/test/testbinned.jl index 6163a333..34f79b66 100644 --- a/test/testbinned.jl +++ b/test/testbinned.jl @@ -1,43 +1,27 @@ - function testbinning(a) - dim = size(a, 1) - n = size(a, 2) - idx = rand(1:n, n ÷ 4) - a1 = a[:, idx] - - bpl = BinnedPointList(dim) - for i = 1:size(a, 2) - insert!(bpl, a[:, i]) - end - - for i = 1:size(a1, 2) - ix = insert!(bpl, a1[:, i]) - if ix != idx[i] - return false - end - end - true +function testbinning(a) + dim = size(a, 1) + n = size(a, 2) + idx = rand(1:n, n ÷ 4) + a1 = a[:, idx] + + bpl = BinnedPointList(dim) + for i = 1:size(a, 2) + insert!(bpl, a[:, i]) end - function testbinning2() - A = 50 - B = 100 - Z = 20 - d = 5 - builder = SimplexGridBuilder(; Generator = TetGen) - p1 = point!(builder, 0, 0, 0) - p2 = point!(builder, B, 0, 0) - p3 = point!(builder, B, A, 0) - p4 = point!(builder, B, A + d, 0) - p5 = point!(builder, B, 2 * A + d, 0) - p6 = point!(builder, 0, 2 * A + d, 0) - p7 = point!(builder, 0, A + d, 0) - p8 = point!(builder, 0, A, 0) + + for i = 1:size(a1, 2) + ix = insert!(bpl, a1[:, i]) + if ix != idx[i] + return false + end end + true +end - @test testbinning(rand(1, 10)) - @test testbinning(rand(1, 10000)) - @test testbinning(rand(2, 10)) - @test testbinning(rand(2, 10000)) - @test testbinning(rand(3, 10)) - @test testbinning(rand(3, 10000)) +@test testbinning(rand(1, 10)) +@test testbinning(rand(1, 10000)) +@test testbinning(rand(2, 10)) +@test testbinning(rand(2, 10000)) +@test testbinning(rand(3, 10)) +@test testbinning(rand(3, 10000)) - @test testbinning2() == 8