From 93a425111445f8bdd8680a6ecf84659f70a50061 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 08:55:07 +0100 Subject: [PATCH 01/19] Update make.jl include new documentary parts --- docs/make.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 8d79ed3..5937cd2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,10 +19,12 @@ makedocs(; ), pages=[ "Home" => "index.md", + "Convex Hulls" => "convexhull.md", "Manual Voronoi" => Any[ "Examples Voronoi Generation" => "man/short.md", "man/workflowmesh.md", "man/geometry.md", + "man/sphere.md", "man/raycast.md", "man/multithread.md", "man/database.md", @@ -50,6 +52,6 @@ makedocs(; deploydocs(; repo="github.com/martinheida/HighVoronoi.jl", devbranch="main", - devurl = "v1.3.0",#"v1.0.2", + devurl = "v1.4.0",#"v1.0.2", # versions = ["stable" => "v^", "v#.#.#",], ) From d39c21787657bacb415c23871995c2f7545b7209 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 08:58:09 +0100 Subject: [PATCH 02/19] new index file --- docs/src/index.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index f03cd4c..fbe5f3c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,7 +2,7 @@ CurrentModule = HighVoronoi ``` -# HighVoronoi 1.3.0: $N\log N$ complexity Parallel Computed Voronoi Grids in $\mathbb{R}^d$ +# HighVoronoi 1.4.0: $N\log N$ complexity Parallel Computed Voronoi Grids in $\mathbb{R}^d$ and on $\mathbb{S}^d$ Documentation for [HighVoronoi](https://github.com/martinheida/HighVoronoi.jl). Voronoi mesh generation in arbitrary dimensions + Finite Volume setup, also for vertices with $d+k$, $k>1$ generators. @@ -10,6 +10,11 @@ Documentation for [HighVoronoi](https://github.com/martinheida/HighVoronoi.jl). - [QUICK START on FINITE VOLUME methods: Click here](@ref QuickFV) / [The ABSTRACT WORKFLOW is here](@ref workflowfv) - [Toy file for testing numerical solver](@ref toyfile) +### News to version 1.4.0: + +- [Parallelized computation of convex hulls](@ref convexhull) based on the raycast algorithm (no volumes / integrals) +- [Parallelized computation of Voronoi diagrams on spheres](@ref voronoisphere), including volumes and integrals (spheres are the only non-euclidian geometry where the raycast algorithm can be applied) + ### News to version 1.3.0: - Parallelized compuation of Voronoi Diagrams, volumes and integrals From d94fffc293685110ba1e9610a466605f83a248c9 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 08:59:00 +0100 Subject: [PATCH 03/19] sphere and convexhull docs --- docs/src/man/convexhull.md | 13 +++++++ docs/src/man/sphere.md | 70 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) create mode 100644 docs/src/man/convexhull.md create mode 100644 docs/src/man/sphere.md diff --git a/docs/src/man/convexhull.md b/docs/src/man/convexhull.md new file mode 100644 index 0000000..6674c50 --- /dev/null +++ b/docs/src/man/convexhull.md @@ -0,0 +1,13 @@ + +# [Convex Hull](@id convexhull) + +Use the following to compute the convex hull of the `VoronoiNodes` `xs`: + +```@julia +cv = ConvexHull(xs,"computing convex hull: "; nthreads=Threads.nthreads(), method = RCOriginal) +``` + +- Note that only the following methods are available: `RCOriginal`, `RCNonGeneral`. Again, `RCOriginal` is faster and for nodes in general position (that means there are only $d$ points in $\mathbb{R}^d$ that generate one surface element) while `RCNonGeneral` can be used for nodes in non-general position. +- `nthreads` will be initialized with `Threads.nthreads()` by default. If you want less threads, provide a positive `Int64`. You can set the string to ` "" ` in order to have only a progress bar. +- `length(cv)` will give you the number of surface elements +- `cv[i]` will then give you `(sig,r,u)` where `sig` is a list of generating nodes, `r` is a point on the plane equidistant to all `sig` and `u` is the outer normal of that plane. diff --git a/docs/src/man/sphere.md b/docs/src/man/sphere.md new file mode 100644 index 0000000..9c6e747 --- /dev/null +++ b/docs/src/man/sphere.md @@ -0,0 +1,70 @@ + +# [Voronoi diagrams on spheres](@id voronoisphere) + +## Simple example and explanation + +You can generate a partition of the full sphere as follows: + +```@julia +A=randn(3,400) +for k in 1:400 + A[:,k] .= normalize(A[:,k]) +end +vs = VoronoiSphere(VoronoiNodes(A), integrate=true, integrator=VI_FAST_POLYGON) +vd = VoronoiData(vs) +println(sum(vd.volume)) +``` +In the above example, the output should be close to $4\pi$, but will be lower by a systemic error (see following comments for a fix). + + +Some comments: +- `VI_MONTECARLO` may not be very accurate. You either need a high number of rays or you may be better of with `VI_FAST_POLYGON` and alike +- `VoronoiSphere` supports all commands of `VoronoiGeometry` except for the following: `periodic_grid`, `improving` +- `VoronoiSphere` comes with the additional commands: + - `total_area`: provide a `Float64` if you know the total area. This is helpful as it will correct the surface integral of functions according to a scaling factor that arises from computed area vs. provided `total_area`. + - `transformations`: provide a list of maps that will copy the given nodes and identify the copies with the originals at the stage of `VoronoiData` computation. Application: In 4 dimensions, the upper half of the $\mathbb{S}^3$ with "periodic boundary conditions" is equivalent to the $SO(3)$. See the following example. + - `center`: if you know the center of the sphere, you can provide it. Otherwise, `HighVoronoi` will try to figure out the center from the data. + - `systematic_error=0.0001`: This parameter adjusts an unavoidable error at the interface integrals. a value of `0.0001` means a systematic error of `0.01`%. If you set it to zero, it will crash. + +## Advanced examples + +The area of the upper half of the sphere using periodic boundary conditions: + +```@julia + + A=randn(3,400) + for k in 1:400 + A[:,k] .= normalize(A[:,k]) + if A[1,k] < 0.0 + A[1,k] = -A[1,k] + end + end + vs = VoronoiSphere(VoronoiNodes(A),transformations=(x->-x,), integrate=true, integrator=VI_FAST_POLYGON) + vd = VoronoiData(vs) + println(sum(vd.volume)) + +``` + +The area of the upper half of the sphere using a mirror of the original points. In this case, if cell 1 has neighbor 1 ( use `vd.neighbors[1]` ), this segment will be part of the 2-dimensional boundary of the upper hull. + +```@julia + using LinearAlgebra + + A = randn(3,400) + for k in 1:400 + A[:,k] .= normalize(A[:,k]) + if A[1,k] < 0.0 + A[1,k] = -A[1,k] + end + end + + function mirror(x) + return diagm( [-1.0, 1.0, 1.0] ) * x + end + + vs = VoronoiSphere(VoronoiNodes(A),transformations=(x->mirror(x),), integrate=true, integrator=VI_FAST_POLYGON) + vd = VoronoiData(vs) + println(sum(vd.volume)) + +``` + From 75832bf0803ea7dd7089c11776b5573cdde5da62 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:14:28 +0100 Subject: [PATCH 04/19] Update Project.toml --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index 19a8fa8..b52c443 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" + [compat] julia = "1.7" @@ -34,6 +36,7 @@ StaticArrays = "1.5" Plots = "1.30" Distances = "0.10.11" ProgressMeter = "1.10" +DoubleFloats = "1.0" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 2f71e0403a58195066d7fd05e0a0f2a5a4510dd5 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:14:56 +0100 Subject: [PATCH 05/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b52c443..851a900 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "HighVoronoi" uuid = "1d30c219-805f-47bd-9cb8-33f291c94a4c" authors = ["Martin Heida"] -version = "1.3.1" +version = "1.4.0" [deps] Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" From a7c986e6c2eb9f8f6b072e142bf5c5d4fc074300 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:15:56 +0100 Subject: [PATCH 06/19] update tests --- test/convexhull.jl | 9 +++++++++ test/jld.jl | 3 +++ test/runtests.jl | 4 +++- test/sphere.jl | 16 ++++++++++++++++ 4 files changed, 31 insertions(+), 1 deletion(-) create mode 100644 test/convexhull.jl create mode 100644 test/sphere.jl diff --git a/test/convexhull.jl b/test/convexhull.jl new file mode 100644 index 0000000..b3678c7 --- /dev/null +++ b/test/convexhull.jl @@ -0,0 +1,9 @@ + +@testset "convex hull" begin + xs = VoronoiNodes(rand(6,1000)) + + ch = HighVoronoi.ConvexHull(xs)#,method=RCNonGeneral) + println(ch[1]) + @test true +end + diff --git a/test/jld.jl b/test/jld.jl index 6085a15..d302c2c 100644 --- a/test/jld.jl +++ b/test/jld.jl @@ -15,6 +15,9 @@ function test_jld(db) xs = VoronoiNodes(rand(5,100)) vg = VoronoiGeometry(xs,cuboid(5,periodic=[1]),vertex_storage=db,search_settings=(method=RCOriginal,),integrate=true,integrator=VI_FAST_POLYGON,integrand=x->[x[1]],silence=false) + vd = VoronoiData(vg) + vol = sum(vd.volume) + abs(vol-1.0)>0.01 && error("") println("Step 1") jldopen("geometry5d.jld2","w") do file file["geo"] = vg diff --git a/test/runtests.jl b/test/runtests.jl index 6a82c78..7e0afc4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ using Test -#using Revise +using Revise using HighVoronoi using SpecialFunctions @@ -21,6 +21,8 @@ const global_silence = false end include("tools.jl") include("basics.jl") + include("convexhull.jl") + include("sphere.jl") include("voronoidata.jl") include("statistics.jl") include("fraud.jl") diff --git a/test/sphere.jl b/test/sphere.jl new file mode 100644 index 0000000..824c0ac --- /dev/null +++ b/test/sphere.jl @@ -0,0 +1,16 @@ + +@testset "sphere" begin + + A=randn(3,400) + for k in 1:400 + A[:,k] .= normalize(A[:,k]) + if A[1,k] < 0.0 + A[1,k] = -A[1,k] + end + end + vs = VoronoiSphere(VoronoiNodes(A),transformations=(x->-x,), integrate=true, integrator=VI_FAST_POLYGON) + #error() + vd = VoronoiData(vs) + @test sum(vd.volume)>6.0 +end + From 9fba1ac7c5b17c8406b630697ccfc5a219d145e7 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:20:12 +0100 Subject: [PATCH 07/19] Update README.md --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 880d862..dc20b72 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ -# HighVoronoi.jl -A Julia Package for parallelized computing of high dimensional (i.e. any dimension >= 2) Voronoi Meshes and setting up Finite Volume computations on these Meshes +# HighVoronoi.jl new version 1.4.0 +A Julia Package for parallelized computing of high dimensional (i.e. any dimension >= 2) Euclidean and Spherical Voronoi Meshes and setting up Finite Volume computations on these Meshes. Also allows parallelized convex hull computation. [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://martinheida.github.io/HighVoronoi.jl/stable/) [![Dev(broken)](https://img.shields.io/badge/docs-dev-blue.svg)](https://martinheida.github.io/HighVoronoi.jl/dev/) @@ -9,4 +9,5 @@ A Julia Package for parallelized computing of high dimensional (i.e. any dimensi Refer to the manual for detailed information on how to use. +The new version 1.4.0 provides spherical meshes and convex hull computation. The new version 1.3.0 has improved performance due to new internal data structure and parallelization. From 443a64ac246e7ac95960e9f96cfafb6169496505 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:37:31 +0100 Subject: [PATCH 08/19] modifications for version 1.4 --- src/FVmesh.jl | 2 +- src/HighVoronoi.jl | 12 +- src/chull.jl | 1224 ++++++++++++++++++++++++++++++++++++ src/domain.jl | 2 +- src/edgehashing.jl | 31 + src/edgeiterate.jl | 8 +- src/extended.jl | 82 +++ src/fastpolyintegrator.jl | 8 +- src/geometry.jl | 44 +- src/heuristic.jl | 4 +- src/heuristic_mc.jl | 4 +- src/hvdatabase.jl | 72 ++- src/hvlocks.jl | 62 +- src/improving.jl | 2 +- src/integrate.jl | 142 ++++- src/mcintegrator.jl | 12 +- src/nodes.jl | 2 +- src/periodicmesh.jl | 2 +- src/polyintegrator.jl | 4 +- src/progress.jl | 5 + src/queuehashing.jl | 12 +- src/raycast-types.jl | 156 ++++- src/raycast.jl | 573 ++++++++++++++++- src/searchtrees.jl | 4 + src/serialdomain.jl | 27 +- src/sphere.jl | 165 +++++ src/sphericalmeshview.jl | 275 ++++++++ src/sphericalpublicview.jl | 358 +++++++++++ src/statistics.jl | 2 +- src/sysvoronoi.jl | 25 +- src/tools.jl | 31 + src/voronoi_mesh.jl | 5 +- src/voronoidata.jl | 36 +- src/voronoidomain.jl | 2 + 34 files changed, 3226 insertions(+), 169 deletions(-) create mode 100644 src/chull.jl create mode 100644 src/sphere.jl create mode 100644 src/sphericalmeshview.jl create mode 100644 src/sphericalpublicview.jl diff --git a/src/FVmesh.jl b/src/FVmesh.jl index 750c2ec..f602a81 100644 --- a/src/FVmesh.jl +++ b/src/FVmesh.jl @@ -176,7 +176,7 @@ end @inline Base.getproperty(cd::SerialMesh, prop::Symbol) = dyncast_get(cd,Val(prop)) @inline @generated dyncast_get(cd::SerialMesh, ::Val{:length}) = :(getfield(cd,:data)[1]) @inline @generated dyncast_get(cd::SerialMesh, d::Val{S}) where S = :( getfield(cd, S)) -@inline @generated dyncast_get(cd::SerialMesh, d::Val{:boundary_Vertices}) where S = :( getfield(cd, :meshes)[1].mesh.boundary_Vertices) +@inline @generated dyncast_get(cd::SerialMesh, d::Val{:boundary_Vertices}) = :( getfield(cd, :meshes)[1].mesh.boundary_Vertices) @inline Base.setproperty!(cd::SerialMesh, prop::Symbol, val) = dyncast_set(cd,Val(prop),val) diff --git a/src/HighVoronoi.jl b/src/HighVoronoi.jl index 1a22686..42b0a81 100644 --- a/src/HighVoronoi.jl +++ b/src/HighVoronoi.jl @@ -14,9 +14,11 @@ using GLPK using Plots using Distances using ProgressMeter -#using Cthulhu +using Cthulhu using Base.Threads using Base.Threads: Atomic, atomic_cas! +using DoubleFloats +#using LoggingExtras #using Traceur @@ -183,6 +185,9 @@ include("parallelintegral.jl") include("voronoidomain.jl") include("serialdomain.jl") include("geometry.jl") # VoronoiGeometry +include("sphericalmeshview.jl") +include("sphere.jl") +include("sphericalpublicview.jl") include("voronoidata.jl") include("discretefunctions.jl") include("substitute.jl") # refinement by substitution @@ -228,7 +233,7 @@ include("l1projection_new.jl") include("R-projection.jl") include("finitevolume.jl") include("statistics.jl") -#include("chull.jl") +include("chull.jl") # What will be exported: export DensityRange export VoronoiNodes @@ -241,6 +246,8 @@ export Voronoi_Integral export VoronoiData export CompactVoronoiData export VoronoiFV +export ConvexHull +export VoronoiSphere export refine! export refine @@ -292,6 +299,7 @@ export RaycastParameter export RCOriginal export RCCombined export RCNonGeneral +export RCNonGeneralHP export MultiThread export SingleThread export AutoThread diff --git a/src/chull.jl b/src/chull.jl new file mode 100644 index 0000000..bc9f875 --- /dev/null +++ b/src/chull.jl @@ -0,0 +1,1224 @@ +################################################################################################################## +################################################################################################################## + +## PrependedVector + +################################################################################################################## +################################################################################################################## + +mutable struct PrependedVector{T, AV <: AbstractVector{T}} <: AbstractVector{T} + first::T + rest::AV +end + +# Implementierung der Länge (1 für first + Länge des restlichen Vektors) +Base.length(pv::PrependedVector) = 1 + length(pv.rest) +Base.size(pv::PrependedVector) = (1 + size(pv.rest)[1],) + +# Zugriff auf Elemente mit `getindex` +function Base.getindex(pv::PV, i::Int) where {PV<:PrependedVector} + if i == 1 + return pv.first + elseif i > 1 && i <= length(pv) + @inbounds return pv.rest[i - 1] + else + throw(BoundsError(pv, i)) + end +end + +function Base.setindex!(pv::PV, value, i::Int) where {PV<:PrependedVector} + if i == 1 + pv.first = value + elseif i > 1 && i <= length(pv) + @inbounds pv.rest[i - 1] = value + else + throw(BoundsError(pv, i)) + end + return pv +end + +function Base.getindex(pv::PV, i::Int) where {T<:Int,PV<:PrependedVector{T}} + if i == 1 + return pv.first + 1 + elseif i > 1 && i <= length(pv) + return @inbounds pv.rest[i - 1] + 1 + else + throw(BoundsError(pv, i)) + end +end + +function Base.setindex!(pv::PV, value, i::Int) where {T<:Int,PV<:PrependedVector{T}} + if i == 1 + pv.first = value - 1 + elseif i > 1 && i <= length(pv) + #@inbounds + pv.rest[i - 1] = value - 1 + else + throw(BoundsError(pv, i)) + end + return pv +end + +# Unterstützung für Iteration +Base.iterate(pv::PrependedVector) = (pv.first, 2) + +function Base.iterate(pv::PrependedVector, state) + if state<=length(pv) + return pv[state],state+1 + else + return nothing + end +end + +# Optional: print für eine anschauliche Ausgabe +Base.show(io::IO, pv::PV) where {PV<:PrependedVector} = begin + print(io, "[$(pv.first)") + for x in pv.rest + print(io, ", $(x)") + end + print(io, "]") +end + + +################################################################################################################## +################################################################################################################## + +## EasyExtendableVector + +################################################################################################################## +################################################################################################################## + + +struct EasyExtendableVector{T} <: AbstractVector{T} + data::Vector{Vector{T}} # Speicher für die Blöcke + length::Atomic{Int64} # Länge des gesamten Vektors, gespeichert in einer MVector + blocklength::Int64 # Länge eines Blocks + lock::ReadWriteLock + # Konstruktor für die Initialisierung + function EasyExtendableVector{T}(blocklength::Int64) where T + new{T}(Vector{Vector{T}}[], Atomic{Int64}(0), blocklength,ReadWriteLock()) + end +end + +# Funktion für die Länge des Vektors +@inline function Base.length(vec::EasyExtendableVector) + return vec.length[] +end + +# Funktion für getindex +function Base.getindex(vec::EasyExtendableVector{T}, i::Int) where T + @assert 1 <= i <= length(vec) "Index out of bounds" + block_idx = (i - 1) ÷ vec.blocklength + 1 + elem_idx = ((i - 1) % vec.blocklength) + 1 + readlock(vec.lock) + ret = vec.data[block_idx][elem_idx] + readunlock(vec.lock) + return ret +end + +# Funktion für push! +function Base.push!(vec::EasyExtendableVector{T}, x::T) where T + writelock(vec.lock) + l_vec = vec.length[] + if l_vec >= length(vec.data) * vec.blocklength + # Neuen Block hinzufügen, wenn der aktuelle Speicher voll ist + push!(vec.data, Vector{T}(undef, vec.blocklength)) + end + # Position des neuen Elements bestimmen + block_idx = (l_vec) ÷ vec.blocklength + 1 + elem_idx = (l_vec % vec.blocklength) + 1 + vec.data[block_idx][elem_idx] = x + atomic_add!(vec.length, 1) # Länge erhöhen + writeunlock(vec.lock) +end + + + +################################################################################################################## +################################################################################################################## + +## routines + +################################################################################################################## +################################################################################################################## + + + +function search_max(xs) + max_ = typemin(Float64) + idx = 0 + for i in 1:length(xs) + x = xs[i] + if x[1]>max_ + max_ = x[1] + idx = i + end + end + return idx +end + +function swap_with_first!(sig, s) + # Finde den Index i, so dass sig[i] == s + i = findfirst(x -> x == s, sig) + + if i !== nothing # Nur, wenn s in sig gefunden wurde + # Tausche den Eintrag an sig[i] mit dem ersten Eintrag sig[1] + sig[i], sig[1] = sig[1], sig[i] + else + println("Wert s nicht im Vektor gefunden.") + end + + return sig +end + +################################################################################################################## +################################################################################################################## + +## Status + +################################################################################################################## +################################################################################################################## + +const STATUS = [Threads.Atomic{Int64}(0) for i in 1:8] +const PRINT_STATUS = Threads.Atomic{Bool}(false) + +#status(i) = @info("($(Threads.threadid()),$i),") #nothing#Threads.atomic_xchg!(STATUS[Threads.threadid()],i) +@inline function status(i) + id = Threads.threadid() + Threads.atomic_xchg!(STATUS[id],i) + PRINT_STATUS[] && print("($id,$i) ") +end + +@inline function status_2(i) + id = Threads.threadid() + Threads.atomic_xchg!(STATUS[id],i) + PRINT_STATUS[] && print("($id,$i) ") +end + +#const demux_logger = TeeLogger( +# MinLevelLogger(FileLogger("info.log"), Logging.Info), +# MinLevelLogger(FileLogger("warn.log"), Logging.Warn), +#); +################################################################################################################## +################################################################################################################## + +## Main Part + +################################################################################################################## +################################################################################################################## + +struct ConvexHull{DB, II, XX} + data::DB + indices::II + xs::XX +end + +# Define length for the structure +@inline Base.length(c::ConvexHull) = length(c.indices) + +# Define getindex to access elements by index +function Base.getindex(c::ConvexHull, i::Int) + f = get_facet(c.data, c.indices[i]) + sig = copy(f[1]) + r = f[2] + u = f[3] + r2 = r + u * dot(u,c.xs[sig[1]]-r) + r3 = r2 + u * dot(u,c.xs[sig[1]]-r2) + return (sig,r3,u) +end + +get_internal_data(c::ConvexHull, i::Int) = get_facet(c.data, c.indices[i]) + +# Define an iterator +@inline Base.iterate(c::ConvexHull, state=1) = state > length(c) ? nothing : (getindex(c,state), state + 1) + +ConvexHull(xs::HV1,intro="computing convex hull: ";nthreads=Threads.nthreads(), method = RCOriginal) where {P,HV1<:HVNodes{P}} = systematic_chull(xs,intro; nthreads=nthreads, method = method) + + +function systematic_chull(xs::HV1,intro="computing convex hull: ";nthreads=Threads.nthreads(), method = RCOriginal) where {P,HV1<:HVNodes{P}}#,HV2<:HVNodes{P}} + prog = ThreadsafeProgressMeter(intro) + dimension = size(P)[1] + + start = search_max(xs) + #println(typeof(RaycastHP(method))," vs. ",typeof(method)) + search = RaycastParameter(Float64,method=RaycastHP(method)) + searcher = Raycast(xs,options=search) + low, hi = HVNearestNeighbors.reduction!(searcher.tree.tree.tree) + RADIUS = norm(hi-low) + #println("$start: $(xs[start])") + + db = HeapDataBase(P,8*1024) + first_result = descent_chull(xs,searcher,start) + #println(first_result) + indices = EasyExtendableVector{Int64}(1024) + + threading = MultiThread(min(nthreads,Threads.nthreads()),1) + + _threads = create_multithreads(threading) + node_threads = length(_threads) + searchers = Vector{typeof(searcher)}(undef,node_threads) + searchers[1] = searcher #chull_Raycasters(searcher,node_threads) + #searchers = [searcher] #chull_Raycasters(searcher,node_threads) + for i in 2:node_threads + searchers[i] = Raycast(xs,options=search) + HVNearestNeighbors.reduction!(searchers[i].tree.tree.tree) + end + queue_position = Atomic{Int64}(1) + global_end = Atomic{Bool}(false) + next!(prog) #first element calculated + lower_b = round(Int,lowerbound(dimension,dimension)) + #edgecount_global = EdgeHashTable(8*lower_b*dimension^2,threading) + edgecount_global = HashEdgeContainer(8*lower_b*dimension^2,4*node_threads,threading) + + zero = zeros(P) + pxs = PrependedVector(zero,xs) + queue_facet(first_result,db,indices,searcher,pxs,edgecount_global) + confirmed_hull = falses(length(xs)) + #with_logger(HighVoronoi.demux_logger) do + if node_threads>1 + Threads.@threads for i in 1:node_threads + __chull(db,indices,queue_position,prog,searchers[i],edgecount_global,node_threads,global_end,PrependedVector(zeros(P),xs),RADIUS,confirmed_hull) + end + else + for i in 1:node_threads + __chull(db,indices,queue_position,prog,searchers[i],edgecount_global,node_threads,global_end,pxs,RADIUS,confirmed_hull) + end + end + #end + return ConvexHull(db,indices,xs) +end + +function queue_facet(first_result, db, indices, searcher, xs, edgecount_global) + status(31) + first_index = push_facet!(db,first_result) + first_index == 0 && return 0 + status(32) + push!(indices,first_index) + status(33) + queue_convex_edges_on_find(first_result,searcher,xs,edgecount_global) # queue edges of very first vertex + return first_index +end + +function __chull(db,indices,queue_position,prog,searcher,edgecount,node_threads,global_end,xs::HVN,RADIUS,confirmed_hull) where {P,HVN<:PrependedVector{P}} + onb = [MVector(zeros(P)) for i in 1:size(P)[1]] + final_edge = [0] + start_time = time_ns() + atomic_xchg!(PRINT_STATUS,false) + this_sig = [0] + ps_data = PlaneSearchData(P,length(xs.rest),xs.rest) + try + while !global_end[] + status(1) + this_entry = atomic_add!(queue_position,1) + if this_entry == length(indices)+node_threads + atomic_xchg!(global_end,true) + end + status(2) + ii = 0 + while this_entry>length(indices) && !global_end[] + active_wait(100) + ii += 1 + if mod(ii,100)==0 + #print("$(Threads.threadid())[$(HighVoronoi.STATUS[1][]),$(HighVoronoi.STATUS[2][]),$(HighVoronoi.STATUS[3][]),$(HighVoronoi.STATUS[4][])]") + myid = Threads.threadid() + #print("-$(Threads.threadid())") + #yield() + #print("+$(Threads.threadid())($myid)$(HighVoronoi.STATUS) ") + end + my_time = time_ns() + if my_time-start_time>1E10 + + println("e:$(Threads.threadid())[$(HighVoronoi.STATUS[1][]),$(HighVoronoi.STATUS[2][]),$(HighVoronoi.STATUS[3][]),$(HighVoronoi.STATUS[4][])]") + println(Threads.threadid(),": ",STATUS) + #atomic_xchg!(PRINT_STATUS,true) + atomic_xchg!(global_end,true) + error("Ende: $(Threads.threadid())") + end + end + status_2(3) + global_end[] && continue + this_index = 0 + iii = 0 + while this_index==0 && iii<100 + this_index = indices[this_entry] + if this_index==0 + iii += 1 + active_wait(100) + end + end + this_index==0 && println("Mist") + facet_data2 = get_facet(db,this_index) + resize!(this_sig,length(facet_data2[1])) + this_sig .= facet_data2[1] + facet_data = (this_sig,facet_data2[2],facet_data2[3]) + #print("--- $this_entry ")#--- ",facet_data) + #this_entry>500 && break + sig,r = preparate_convex_data!(xs,facet_data) # r is the vertex of the outer extended tetrahedron + status(4) + +# print("$(Threads.threadid())[$(HighVoronoi.STATUS[1][]),$(HighVoronoi.STATUS[2][]),$(HighVoronoi.STATUS[3][]),$(HighVoronoi.STATUS[4][])]") + #println(sig,r,facet_data[1]) + if length(sig)>0 + ei = get_EdgeIterator(sig,r,searcher,1,xs,OnSysVoronoi()) + status(5) + #@descend explore_chull_vertex(ei,facet_data,xs,onb,sig,r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) + #error("") + nv = explore_chull_vertex(ei,facet_data,xs,onb,sig,r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) + end + status_2(6) + #facet_data[1]==[18, 174, 243, 333, 377, 418] && error("KLAPPT!") + + #println() + end +catch e + println(Threads.threadid(),": ",STATUS) + #atomic_xchg!(PRINT_STATUS,true) + open("error_log_chull_$(Threads.threadid()).txt", "w") do f + # Stacktrace speichern + Base.showerror(f, e, catch_backtrace()) + end + atomic_xchg!(global_end,true) + println("hier $(Threads.threadid()) fehler!") + rethrow() +end + #println("ending: $(Threads.threadid())") +end + +function verify(sig,r,u,xs) + ret = true + r += u * dot(u,xs[sig[2]]-r) + c = maximum(dot(xs[g], u) for g in sig) + c = max(c,dot(r,u))*1.00000000001 + for i in 1:length(xs) + i in sig && continue + if dot(xs[i],u)>c + print(dot(xs[i],u)-c,", ",sig,", ",i," ") + ret = false + end + #print(" , ") + end + dim = length(r) + for k in 1:(length(r)) + edge = copy(sig) + skip = edge[k] + deleteat!(edge,k) + onb =[MVector(r) for _ in 1:length(r)] + x0 = xs[edge[1]] + for i in 1:dim + if i<=dim-2 + x = xs[edge[i+1]] + onb[i] .= x .- x0 + elseif i==dim-1 + onb[i] .= u + continue + else + onb[i] .= x0 .- xs[skip] + end + for j in 1:(i-1) + onb[i] .-= onb[j] .* dot(onb[j],onb[i]) + onb[i] .-= onb[j] .* dot(onb[j],onb[i]) + end + normalize!(onb[i]) + end + + end + return ret +end + +@inline function queue_convex_edges_on_find(raw_data,searcher,xs::PV_P,edgecount) where {P,PV_P<:PrependedVector{P}} + sig,r = preparate_convex_data!(xs,raw_data) + queue_edges_OnFind(sig,r,searcher,1,xs,edgecount) +end + +function preparate_convex_data!(xs,data) + sig = copy(PrependedVector(1,data[1])) + r = data[2]+data[3] + dist = norm(r-xs[sig[2]]) + xs[1] = r + dist*data[3] + return sig,r +end + +function get_r_u_o(sig,original_r,edge,edgeIterator,xs,_normal,searcher,RADIUS) + full_edge, u = get_full_edge(sig,original_r,edge,edgeIterator,xs) + deleteat!(full_edge,1) + c1 = maximum(dot(xs[g], u) for g in full_edge) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + r = original_r + + bbb = HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,u,c,true)==0 + if bbb + para = norm(u - _normal * dot(_normal,u)) + r2 = r + (RADIUS/para) * u # move outwards in direction u + print(" !!!! ") + return full_edge,r2, -_normal, r # why not r2 ?? + end + return full_edge, original_r, u, original_r +end + +@inline function delta_u(onb,dim) + s = 0.0 + for i in 2:dim + for j in 1:(i-1) + s += abs(dot(onb[i],onb[j])) + end + end + return s +# return @inbounds sum(i->abs(dot(onb[i],onb[dim])),1:(dim-1)) +end + +@inline function samecount(big,small,expected) + _samecount = 0 + for ed in big + _samecount += ed in small ? 1 : 0 + end + return _samecount==expected +end + + +@inline raycast_des_chull(a,b,c,d,e,f,g,h,i,j::Raycast_Original_HP,k,l,m) = raycast_des2(a,b,c,d,e,f,g,h,i,j,k) +@inline raycast_des_chull(a,b,c,d,e,f,g,h,i,j::Raycast_Non_General_HP,k,l,m) = raycast_des2(a,b,c,d,e,f,g,h,i,j,l) +function raycast_des3(full_edge, r, u, xs, searcher,old,fe,facet_data,cast,method,safety=2^length(r),debug=false,data=nothing) + + #_generator_1, _t, _r2 = raycast_des3(copy(full_edge), r, u, xs, searcher,copy(old),copy(fe),copy(facet_data),cast,RCOriginalSafety,safety,debug,data) + generator_1, t, r2 = raycast_des_chull(full_edge, r, u, xs, searcher,old,fe,facet_data,cast,method,safety,debug,data) + #_generator_1!=generator_1 && error("$generator_1 vs. $_generator_1") + #print("+") + #generator_1, t, r2 = raycast_des2(full_edge, r, u, xs, searcher,old,fe,facet_data,cast,RCCombined) + debug && println(full_edge,",",generator_1,",",t) + if t==Inf && generator_1!=0 + filter!(i->i!=generator_1,full_edge) + return raycast_des2(full_edge, r, u, xs, searcher,old,fe,facet_data,cast,RCNonGeneralHP,debug) + else + return generator_1, t, r2 + end +end + +#=function raycast_des3(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,method::Raycast_Original_Safety,safety=2^length(r),debug=false, data = 3) + x0 = xs[edge[1]] + maxiter = safety + + reset(data,x0,r,u,xs,origin,searcher) + + data.r = r + u * dot(u , (x0-r)) + i, t = search_vertex_plane(searcher.tree,data) + println(i," ",dot(u,xs[i]),", ",data.c0) + t == Inf && return 0, Inf, r + i2 = i + #if i2 in origin + # HVNearestNeighbors.global_i2[1] = i2 + # i, t = _nn(searcher.tree, vvv, skip) + #end + #i2==195 && print("!") + #sk = skip(i2) + #HVNearestNeighbors.global_i2[1] = 0 + #(i2 in origin) && error("????????? $i2, $(dot(xs[i2],u)) $c $(sk) $(typeof(searcher.tree.tree))") + count = 0 + while i!=0 + count==maxiter && (break) + count += 1 + x = xs[i] + t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + vvv = r+t*u + data.r = vvv + i, t = search_vertex_plane(searcher.tree,data) + if (i==i2 || (i in origin)) + i=0 + end + if i!=0 + i2 = i + end + end + x2 = xs[i2] + t = get_t(r,u,x0,x2) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + _r = r+t*u + + append!(sig,i2) + sort!(sig) + + return i2, (count==maxiter && i!=0) ? Inf : t, _r + +end=# + +function explore_chull_vertex(edgeIterator,facet_data,xs,onb,sig,original_r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) + dim = size(eltype(xs))[1] + debug = false#facet_data[1]==[477, 1457, 1540, 1605, 1660, 1980] #[33, 104, 210, 285, 392, 458] + #println(facet_data) + sss = copy(edgeIterator.sig) + sss .-= 1 + fd1 = facet_data[1] + #HighVoronoi.samecount(fd1,sss,dim)==false && error("also hier geht's los... $(facet_data[1]), $edgeIterator") + for (edge,skip) in edgeIterator + status_2(11) + _normal = facet_data[3] + b = pushedge!(edgecount,edge,1) + (edge[1]!=1 || b ) && continue + ee2 = Vector{Int64}(edge) + ee2 .-=1 + #HighVoronoi.samecount(fd1,ee2,dim-1)==false && error("hier schon Mist: $(typeof(edgeIterator)) $ee2, $(fd1),$edgeIterator") + le = length(edge) + edge2 = zeros(Int64,le+1)#Vector{Int64}(undef, le+1) + view(edge2,1:(le-1)) .= view(edge,2:le) + edge2 .-= 1 + + full_edge, r,u, orientation = get_r_u_o(sig,original_r,edge,edgeIterator,xs,_normal,searcher,RADIUS) + x0 = xs[edge[2]] + debug && println("klassicher fall: $(r==original_r) ") + + + #full_edge, u_1, onb = get_full_edge_basis(sig,r,edge,edgeIterator,xs) + full_edge .-= 1 + start_edge = full_edge[1] + debug && println("###################################################") + debug && println(full_edge) + #@descend raycast_des2(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginalSafety,20,debug) + #error("") + #generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginal,20,debug,ps_data) + generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),searcher.parameters.method,20,debug,ps_data) + #generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginalSafety,20,debug,ps_data) + + generator_1==0 && continue + + #generator_1, t, r2 = raycast_des2(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginal) + #generator_1, t, r2 = raycast_des2(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCCombined) + t2 = 0.0 + #try + t2 = get_t(r2,u,x0,xs.rest[generator_1])# (sum(abs2, r - xs[generator]) - sum(abs2, r - x0)) / (2 * u' * (xs[generator]-x0)) + #catch + # println("$fd1") + # rethrow() + #end + debug && println("First step: ",full_edge,", ir=",sort!(inrange(searcher.tree.tree,r2,norm(r2-xs.rest[full_edge[1]])*1.0000000001)),", r2=",r2," from r=",r) + orientation==original_r && (orientation = r2) + + HighVoronoi.samecount(full_edge,edge2,dim-1)==false && continue + + edge2[le] = generator_1 + for i in 1:dim + if i1E-13 && print("!") + if du>1E-13 + u2 = u + u1 = u_qr_onb(onb,x0)#u_qr(sig,xs,1) + onb[dim] .= dot(u2,u1)>0 ? u1 : -1.0*u1 + du2 = delta_u(onb,dim) + #du2>10*du && error("$du, $du2") + du = du2 + end + u = v = SVector(onb[dim]) + + c1 = maximum(dot(xs.rest[g], v) for g in full_edge) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + _valid = HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,v,c,true)!=0 + + debug && println("valid=$_valid") + generator, t, r3 = _valid ? raycast_des2(full_edge, r2, v, xs.rest, searcher,0,full_edge,full_edge,Raycast_By_Walkray(),RCNonGeneralHP,false,du) : (0, Inf64, r2) + + HighVoronoi.samecount(full_edge,edge2,dim)==false && continue #error("what the fuck??? $samecount $final_edge, $(facet_data[1])") + #HighVoronoi.samecount(facet_data[1],edge2,dim-1)==false && error("hier schon Mist: $edge2, $(facet_data[1]), $fd1") + + edge2[le+1] = generator + if t==Inf + #=if !_valid + __t = false# (full_edge == [165, 199, 225, 311, 370, 416]) + __t && println("----------------------------------------------") + my_r1 = r2 + v * dot(v,x0-r2) + my_r2 = my_r1 + v * dot(v,x0-my_r1) + mymax = maximum(norm(my_r2-xs.rest[g]) for g in full_edge) + closest = HVNearestNeighbors.inrange(searcher.tree.tree.tree,my_r2,mymax) + lclos = length(closest) + myrad = maximum(norm(r2-xs.rest[g]) for g in full_edge) + for __i in 1:length(closest) + ind = closest[__i] + if ind in full_edge + closest[__i] = 0 + continue + end + xx = xs.rest[ind] + __t && print("$ind: $(norm(xx-r2)) ") + if abs(myrad-norm(xx-r2))/myrad>1E-15 + __t && print("true ") + closest[__i] = 0 + end + __t && println(",") + end + filter!(i->i!=0,closest) + if length(closest)!=0 + println(m,", ",lclos) + for __i in 1:length(closest) + ind = closest[__i] + println("$ind: $(norm(xs.rest[ind]-r2)), $(norm(r2-xs.rest[ind]))") + end + error("$full_edge, $closest") + end + append!(full_edge,closest) + sort!(full_edge) + #abs(dot(v,my_r2-x0))>1E-10 && println("fehler") + end=# + if norm(r3-xs.rest[full_edge[1]])>10000 + debug && println("A0: ",full_edge) + generator, t, r3 = raycast_des2(full_edge, r2, -v, xs.rest, searcher,0,full_edge,full_edge,Raycast_By_Walkray(),RCNonGeneralHP,debug,du) + generator==0 && continue + #generator, t, r3 = raycast_des2(full_edge, r2, -v, xs.rest, searcher,0,full_edge,full_edge,Raycast_By_Walkray(),RCOriginal) + u = v + debug && print("A1: $edge2, $(vertex_variance(full_edge,r3,xs.rest))") + #=if debug + for s in full_edge + print(" $s: $(norm(r3-xs.rest[s])),") + end + println() + println("--------------") + end=# + r,_ = walkray_correct_vertex(r3, full_edge, searcher, edge2, generator) + r3 = r + debug && print("A2: $edge2, $(vertex_variance(full_edge,r3,xs.rest))") + HighVoronoi.samecount(final_edge,edge2,dim)==false && continue + #=if debug + for s in full_edge + print(" $s: $(norm(r3-xs.rest[s])),") + end + end=# + end + debug && println("hier") + else + debug && println("rotate: $final_edge,$full_edge") + debug && print("B: ")# $generator_1,$generator, $final_edge, $edge2 ||| ") + r,_ = walkray_correct_vertex(r3, full_edge, searcher, edge2, generator) + #rotate_data = [original_r,original_r] + r2, u22 = rotate_outwards(onb,searcher,edge2,r,xs.rest,x0, orientation-xs.rest[edge2[1]], full_edge,final_edge,debug) + c1 = maximum(dot(xs.rest[g], u22) for g in full_edge) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + _valid = HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,u22,c,true)!=0 + HighVoronoi.samecount(final_edge,edge2,dim)==false && continue + r3 = r2 + u = u22 + end + #=ir = sort!(copy(inrange(searcher.tree.tree,r3,norm(r3-xs.rest[final_edge[1]])*1.000000000001))) + if sort!(full_edge) != ir && false + m = 0.0 + for ii in full_edge + m = max(m, norm(xs.rest[ii]-r3)) + end + bb = false + for ii in ir + ii in full_edge && continue + bb |= norm(xs.rest[ii]-r3)6 + println() + for ii in ir + print("$ii: $(norm(xs.rest[ii]-r3)), ") + end + #println(" vv=$(vertex_variance(ir,r3,searcher.tree.extended_xs))") + for ii in full_edge + print("$ii: $(norm(xs.rest[ii]-r3)), ") + end + println() + for ii in final_edge + print("$ii: $(norm(xs.rest[ii]-r3)), ") + end + status(22) + error("$(facet_data[1]), $edge2, $final_edge")# $(vertex_variance(edge2,r3,xs.rest)), $(vertex_variance(final_edge,r3,xs.rest)),") + end + #println() + #nns = sort!(inrange(searcher.tree.tree,r2,norm(r2-x0)*1.000000001)) + #(!vvvv ) && error("$edge2, $(nns), $r2") + status(23) + (!vvvv) && error("$(facet_data[1]), $final_edge, $edge2, $r2, $(vertex_variance(final_edge,r2,xs.rest,length(final_edge)-1)), $(dot(r2,r2))") + status(17)=# + #if final_edge==[683, 1457, 1540, 1605, 1660, 1916, 1980] + # println("Fehler: ",facet_data) + #end + queue_facet((final_edge,r3,u), db, indices, searcher, xs, edgecount) + status_2(18) + #view(confirmed_hull,final_edge) .= true + #= + + =# + end +end + +function rotate_outwards(onb,searcher,sig,r_::P,xs,x0,orientation,final_sig,final_edge,debug=false) where {P} + #println("variance: ",vertex_variance(final_sig,r_,xs),final_sig," ",r_) + i = 0 + i_ = 0 + r = r_ + u = Ref(r_) + dim = length(r_) + last_sig = copy(final_sig) + kkk = 0 + #debug = false + while kkk<30 + debug && println(last_sig) + kkk += 1 + #b = false + for kk in 0:1 + i = dim+(1-kk) + i_ = dim + kk + #println("$i,$i_, $sig, $last_sig - ") + onb[dim] .= x0 .- xs[sig[i_]] + onb[dim-1] .= xs[sig[i]] .- x0 + for i in 1:(dim-2) + onb[dim-1] .-= dot(onb[dim-1],onb[i]) .* onb[i] + onb[dim-1] .-= dot(onb[dim-1],onb[i]) .* onb[i] + end + normalize!(onb[dim-1]) + for i in 1:(dim-1) + onb[dim] .-= dot(onb[dim],onb[i]) .* onb[i] + onb[dim] .-= dot(onb[dim],onb[i]) .* onb[i] + onb[dim] .-= dot(onb[dim],onb[i]) .* onb[i] + end + normalize!(onb[dim]) + u[] = P(onb[dim]) + #println(dot(u,orientation),u) + dot(onb[dim],orientation)<=0 && continue + #dot(u,orientation)<=0 && continue + + full_sig = copy(sig) + deleteat!(full_sig,i_) + + if debug + println("+",sig,inrange(searcher.tree.tree,r,norm(r-x0)*1.00000001),", ",full_sig) + for iii in 1:dim + dd = dot(onb[dim],onb[iii]) + if abs(dd)>1E-8 + print("($iii,$dim: $dd), ") + end + end + for iii in 1:dim + dd = dot(onb[dim-1],onb[iii]) + if abs(dd)>1E-8 + print("($iii,$(dim-1): $dd), ") + end + end + println() + end + #if sig==[248, 294, 379, 397, 418, 21, 406] + # raycast_hull(full_sig, r, u, xs, searcher,0,full_sig,last_sig,Raycast_By_Walkray()) + # println(last_sig) + # bb = true +# # error("hier") + #end + du = delta_u(onb,dim) + if du>1E-13 + if kk==1 + sig[dim+1], sig[dim] = sig[dim], sig[dim+1] + end + u1 = u_qr_onb(onb,x0)#u_qr(sig,xs,1) + #u2 = dot(u,u1)>0 ? u1 : -1.0*u1 + #u = u2 + #=println(dot(u,u2)) + for iii in 1:dim + #s = sig[i] + print("$iii: $(dot(onb[iii],u2)), ") + end + println() + for iii in 1:dim + #s = sig[i] + print("$iii: $(dot(onb[iii],u)), ") + end + println() + for iii in 2:dim+1 + #s = sig[i] + print("$iii: $(dot(normalize(xs[sig[iii]]-x0),u2)), ") + end + println() + for iii in 2:dim+1 + #s = sig[i] + print("$iii: $(dot(normalize(xs[sig[iii]]-x0),u)), ") + end + println() + print("kk=$kk, du=$du -> ")=# + onb[dim] .= dot(onb[dim],u1)>0 ? u1 : -1.0*u1 + du2 = delta_u(onb,dim) + if kk==1 + sig[dim+1], sig[dim] = sig[dim], sig[dim+1] + end + du2>du && error("") + du = du2 + end + u[] = SVector{size(P)[1],Float64}(onb[dim]) + + __c1 = maximum(dot(xs[g], u[]) for g in last_sig) + __c2 = abs(__c1) + __c = __c1 + __c2*searcher.plane_tolerance + _valid = HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,u[],__c,true)!=0 + +# debug && println("valid=$_valid") +# generator, t, r3 = _valid ? raycast_des2(full_edge, r2, v, xs.rest, searcher,0,full_edge,full_edge,Raycast_By_Walkray(),searcher.parameters.method,false,du) : (0, Inf64, r2) + + debug && println("***************************************************") + if _valid + generator, t, r2 = raycast_des2(full_sig, P(r), u[], xs, searcher,0,full_sig,last_sig,Raycast_By_Walkray(),RCNonGeneralHP,debug,du) + last_sig = full_sig + sig[i_] = generator + r = r2 + r,_ = walkray_correct_vertex(r, last_sig, searcher, sig, generator) + break + else + resize!(final_sig,length(last_sig)) + final_sig .= last_sig + #println("Ende: $last_sig") + debug && println("* du = $du, full=$full_sig, last=$last_sig, r = $r ",inrange(searcher.tree.tree,r+10*u,norm(r+10*u-xs[final_sig[1]])*1.00000000001),inrange(searcher.tree.tree,r-0.5*u,norm(r-0.5*u-xs[final_sig[1]])*1.00000001)) + #print("F: $full_sig, ") + for ii in 1:length(final_sig) + aprod = abs(dot(normalize(xs[last_sig[ii]]-x0),u[])) + debug && print("$(last_sig[ii]): $aprod - ") + if xs[last_sig[ii]]!=x0 && aprod>max(1E-10,10*du) + debug && print("$ii: $(normalize(xs[last_sig[ii]]-x0)), $(dot(normalize(xs[last_sig[ii]]-x0),u)) -- ") + last_sig[ii] = 0 + end + debug && println() + end + debug && println(last_sig) + #println(length(full_sig)," vs. ",length(last_sig)) + filter!(f->f!=0, last_sig) + resize!(final_edge,length(last_sig)) + final_edge .= last_sig + sort!(final_edge) + #println(final_edge) + return r,u[] + end + end + end +end + +################################################################################################################## +################################################################################################################## + +## Descend + +################################################################################################################## +################################################################################################################## + + + + + +function descent_chull(xs::Points, searcher::RaycastIncircleSkip, start) + dim = searcher.dimension + sig = [start] + r = xs[start] + x0 = xs[start] + minimal_edge = zeros(Int64,dim+1) + nearest, t = _nn(searcher.tree, r, i->(i==start)) + direction = MVector(xs[nearest]-xs[start]) + direction[1] = 0.0 + normalize!(direction) + direction[1] = 1.0 + normalize!(direction) + MV = typeof(direction) + base = [zeros(MV) for i in 1:dim] + success = false + final_sig = [0] + + + base[end] .= direction + #minimal_edge .= 0 + minimal_edge[1] = start + my_vv = 1.0 + try + for k in 1:(dim-1) # find an additional generator for each dimension + u = SVector(base[end]) + + generator, t, r2 = raycast_des(sig, r, u, xs, searcher,0,sig,sig,Raycast_By_Descend(),RCNonGeneralHP) + b = false + if t == Inf + u = -u + generator, t, r2 = raycast_des(sig, r, u, xs, searcher,0,sig,sig,Raycast_By_Descend(),RCNonGeneralHP) + end + if t == Inf + error("Could not find a vertex in both directions of current point." * + "Consider increasing search range (tmax)") + end + base[k] .= xs[generator] .- x0 + normalize!(base[k]) + for i in 1:(k-1) + base[k] .-= dot(base[k],base[i]) .* base[i] + base[k] .-= dot(base[k],base[i]) .* base[i] + normalize!(base[k]) + end + base[end] .-= dot(base[k],base[end]) .* base[k] + normalize!(base[end]) + r = r2 + minimal_edge[k+1] = generator + my_vv = vertex_variance(view(minimal_edge,1:(k+1)),r,xs,k,view(searcher.ddd,1:(k+1))) + my_vv>searcher.variance_tol && error("$my_vv, $minimal_edge") + swap_with_first!(sig,start) + end + u = SVector(base[end]) + generator, t, r2 = raycast_des(sig, r, u, xs, searcher,0,sig,sig,Raycast_By_Descend(),RCNonGeneralHP) + b = false + if t == Inf # erfolgreich! + success = true + final_sig = copy(sig) + generator, t, r2 = raycast_des(sig, r, -u, xs, searcher,0,sig,sig,Raycast_By_Descend(),RCNonGeneralHP) + r = r2 + else # vertex, aber nicht äußere FLäche... + #println("hier 2, start = $start, sig = $sig") + minimal_edge[end] = generator + r = r2 + my_vv = vertex_variance(minimal_edge,r,xs, dim , view(searcher.ddd,1:(length(minimal_edge)))) + my_vv>searcher.variance_tol && error("$my_vv, $minimal_edge") + sort!(sig) + r,_ = walkray_correct_vertex(r, sig, searcher, minimal_edge,minimal_edge[dim+1]) + swap_with_first!(sig,start) + #println("sig = $sig, r = $r") + r, success = search_local_hull!(sig,start,searcher,r,base[end],SVector(x0+direction),xs,final_sig) + u = SVector(base[end]) + success = true + end + + catch + rethrow() + my_vv=1.0 + end + return (final_sig, r, SVector(base[end])) +end + +function search_local_hull!(sig,start,searcher,r,direction,y0,xs,final_sig) + stop = false + while !stop + if (length(sig)==length(r)+1) + #=print("HIER 1: $(sig[1]) vs $start: ") + ei = General_EdgeIterator(sig,r,start,searcher.general_edgeiterator) + for (edge, skip) in ei + print("($edge,$skip) - ") + end + println()=# + ei = General_EdgeIterator(sig,r,start,searcher.general_edgeiterator) + r2, stop = travel_right_direction(sig,r,ei,direction,y0,searcher,xs,final_sig) + r = r2 + else + #println("HIER 2") + ei = searcher.edgeiterator + fei = get!(searcher.FEIStorage_global,sig,FEIStorage(sig,r)) + HighVoronoi.reset(ei,sig,r,searcher.tree.extended_xs,_Cell,searcher,fei) + r2, stop = travel_right_direction(sig,r,ei,direction,y0,searcher,xs,final_sig) + r = r2 + end + swap_with_first!(sig,start) + #println("sig = $sig, stop = $stop, r = $r") + end + return r, true +end + +function travel_right_direction(sig,r,edgeIterator,direction,y0,searcher,xs,final_sig) + i = 0 + for (edge,skip) in edgeIterator + i += 1 + i==1 && continue + full_edge, u = get_full_edge(sig,r,edge,edgeIterator,xs) + + dot(y0,u)<=0 && continue + resize!(final_sig,length(full_edge)) + final_sig .= full_edge + #view(final_sig,2:(length(full_edge)+1)) .= full_edge + #final_sig[1] = skip[1] + return walkray_right_direction(sig,r,u,searcher,direction,xs,full_edge,edge) + end + #print("steps: $i , ") + return r, true +end + +function walkray_right_direction(sig,r,u,searcher,direction,xs,full_edge,edge) + c1 = maximum(dot(xs[g], u) for g in full_edge) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + if HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,u,c,true)==0 + direction .= u + (return r, true) + end + sig2, r2, _ = walkray(full_edge, r, xs, searcher, sig, u, edge ) # provide missing node "j" of new vertex and its coordinate "r" + resize!(sig,length(sig2)) + sig .= sig2 + direction .= u + return r2, false +end + +#= +function travel_right_direction(sig,r,edgeIterator,direction,y0,searcher,xs) + i = 0 + for (edge,_) in edgeIterator + i += 1 + i==1 && continue + full_edge, u = get_full_edge(sig,r,edge,edgeIterator,xs) + #=print("FE=$full_edge: ") + for s in sig + print("$s->$(dot(xs[s]-xs[sig[1]],u)), ") + end=# + dot(y0,u)<=0 && continue + c1 = maximum(dot(xs[g], u) for g in full_edge) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + if HVNearestNeighbors.peak_direction(searcher.tree.tree.tree,u,c,true)==0 + x0 = xs[sig[1]] + #println() + #println("H: u=$u , u*x0=$(dot(x0,u)) , xs[10]=$(xs[10]) ") + #HVNearestNeighbors.search_node_direction(searcher.tree.tree.tree,u,c,10) + #=for ii in 1:length(xs) + cc = dot(xs[ii], u) + if cc>c + print("$ii, ") + end + if (dot(xs[ii]-x0,u)>0) + print("f$ii=$(dot(xs[ii]-x0,u)), ") + end + end + + println("c = $c, sig = $sig, full_edge = $full_edge") + println(dot(xs[17],u))=# + direction .= u + (return r, true) + end + #print("h ") + sig2, r2, _ = walkray(full_edge, r, xs, searcher, sig, u, edge ) # provide missing node "j" of new vertex and its coordinate "r" + resize!(sig,length(sig2)) + sig .= sig2 + direction .= u + return r2, false + end + #print("steps: $i , ") + return r, true +end + +=# + +#= + +function systematic_chull(xs::HV1) where {P,HV1<:HVNodes{P},HV2<:HVNodes{P}} + @time search=RaycastParameter(eltype(P)) + DIM = size(P)[1] + queue = IndifferentQueue{Pair{Vector{Int64},P}} # store new verices in a queue + lower_b = round(Int,lowerbound(dimension,dimension)) + edgecount_global = EdgeHashTable(8*lower_b*dimension^2) # remember edges + vertexcount_global = VertexHashTable(8*lower_b*dimension^2) # remember vertices with more than dim+1 generators + + +end + + + +function chull(xs::HV1,hull::HV2) where {P,HV1<:HVNodes{P},HV2<:HVNodes{P}} + long_xs = DoubleVector(hull,xs) + mmm = cast_mesh(ExternalMemory(4),long_xs) + @time search=RaycastParameter(eltype(P)) + DIM = size(P)[1] + #@time voronoi(mmm,searcher=Raycast(long_xs;domain=cuboid(DIM,dimensions=20*ones(Float64,DIM),periodic=Int64[],offset=-10*ones(Float64,DIM)),options=search),intro="",Iter=1:12)#length(hull)) + @time voronoi(mmm,searcher=Raycast(copy(long_xs);options=search),intro="",Iter=1:14)#length(hull)) + b = falses(length(long_xs)) + for i in 1:14 + for (sig,r) in all_vertices_iterator(mmm,i) + #if sig[2]>6 + b[sig] .= true + #end + end + end + println("alle: ",sum(b)) +end + +function chull2(xs::HV1) where {P,HV1<:HVNodes{P},HV2<:HVNodes{P}} + mmm = cast_mesh(ExternalMemory(4),xs) + search=RaycastParameter(eltype(P)) + DIM = size(P)[1] + @time voronoi(mmm,searcher=Raycast(xs;options=search),intro="",Iter=1:1)#length(hull)) + +end + + + +function cube_vertices(dim::Int,stretch,offset,c) + # Anzahl der Eckpunkte eines DIM-dimensionalen Würfels ist 2^DIM + num_vertices = 2^dim + + # Initialisiere ein leeres Array für die Eckpunkte + vertices_b = zeros(Bool, num_vertices, dim) + + for i in 0:(num_vertices-1) + # Jede Zahl von 0 bis 2^DIM-1 entspricht einem Eckpunkt + # Konvertiere die Zahl in eine binäre Darstellung + bin_rep = digits(i, base=2, pad=dim) + # Speichere die binäre Darstellung als Koordinate im Array + vertices_b[i+1, :] .= bin_rep + end + + # Konvertiere die Bool-Werte in Float64 + vertices_ = convert(Array{Float64}, vertices_b) + vertices = copy(transpose(vertices_)) + + # Verschieben der Punkte um offset + vertices .+= offset + for i in eachindex(vertices) + vertices[i] += c*(2*rand()-1) + end + + v1 = VoronoiNodes(vertices) + + vectors = Vector{Float64}[] + + for i in 1:dim + # Erstelle einen Nullvektor + vec = zeros(Float64, dim) + # Setze die i-te Komponente auf 1 + vec[i] = 1.0 + for i in eachindex(vec) + vec[i] += 2*c*(2*rand()-1) + end + push!(vectors, vec) + # Erstelle den negativen Vektor + neg_vec = -vec + push!(vectors, neg_vec) + end + #append!(v1,VoronoiNodes(vectors)) + v1 = VoronoiNodes(vectors) + for i in eachindex(v1) + v1[i] = 10*normalize(v1[i]) + end + println(v1) + return v1 +end + +=# diff --git a/src/domain.jl b/src/domain.jl index 3970be3..a3d51e5 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -47,7 +47,7 @@ end=# ############################################################################################################################ -function add_virtual_points(domain::AD,new_xs::ReflectedNodes;search_settings::RP=RaycastParameter(Float64),do_refine=statictrue,kwargs...) where {AD<:AbstractDomain,RP<:RaycastParameter} +function add_virtual_points(domain::AD,new_xs::ReflectedNodes;search_settings::RP=RaycastParameter(Float64),do_refine=statictrue,kwargs...) where {AD<:AbstractDomain,RP} # <: RaycastParameter} length(new_xs)==0 && return Int64[] expand_internal_boundary(domain,new_xs) prepend!(domain,new_xs) diff --git a/src/edgehashing.jl b/src/edgehashing.jl index 82369a6..1f59a61 100644 --- a/src/edgehashing.jl +++ b/src/edgehashing.jl @@ -71,7 +71,9 @@ function pushedge!(ht::EdgeHashTable, key::K, cell::Int64, mode=true) where K i = UInt64(0) ret = false + #print("$(Threads.threadid())?") lock(ht.lock) + #print("#") #print("$value, $index2") while true idx = reinterpret(Int64,(index1 + i * index2) & ht.mylength[1] + 1) # try bitcast( ) instead @@ -103,6 +105,7 @@ function pushedge!(ht::EdgeHashTable, key::K, cell::Int64, mode=true) where K break end end + #print("$(Threads.threadid())!") unlock(ht.lock) return ret end @@ -146,6 +149,34 @@ function Base.empty!(ht::EdgeHashTable) fill!(ht.occupied,false) unlock(ht.lock) end +struct HashEdgeContainer{EHT} + container::Vector{EHT} + length::Int64 + function HashEdgeContainer(len::Int64,split::Int64,lock=SingleThread) + mysplit = reinterpret(Int64,next_power_of_two(split))-1 + mylength = div(len,mysplit) + eht = EdgeHashTable(mylength,lock) + cont = Vector{typeof(eht)}(undef,mysplit+1) + cont[1] = eht + for i in 2:(mysplit+1) + cont[i] = EdgeHashTable(mylength,lock) + end + return new{typeof(eht)}(cont,mysplit) + end +end + +function pushedge!(ht::HashEdgeContainer, key::K, cell::Int64, mode=true) where K + index = (sum(key) & ht.length) + 1 +# error("$(sum(key)), $(ht.length), $(sum(key) & ht.length)") + pushedge!(ht.container[index],key,cell,mode) +end +function Base.empty!(ht::HashEdgeContainer) + for eht in ht.container + empty!(eht) + end +end + + #= # Methode zum Prüfen, ob ein Schlüssel in der EdgeHashTable vorhanden ist function Base.haskey(ht::EdgeHashTable, key) diff --git a/src/edgeiterate.jl b/src/edgeiterate.jl index 96c0c26..2657d42 100644 --- a/src/edgeiterate.jl +++ b/src/edgeiterate.jl @@ -110,6 +110,12 @@ function get_full_edge(sig,r,edge,NF_::FastEdgeIterator,xs) return Vector{Int64}(view(NF.sig, fullview[1:count])), convert_SVector( -1 .* ray(NF_) ) end +#=function get_full_edge_basis(sig,r,edge,NF_::FastEdgeIterator,xs) + NF=NF_.iterators[1] + fullview,count = get_full_edge(NF.rays,NF.local_cone,NF.new_node,NF.active_nodes[1],length(NF.sig),NF.dim,NF_.ray_tol) + return Vector{Int64}(view(NF.sig, fullview[1:count])), convert_SVector( -1 .* ray(NF_) ), [copy(NF.rays[i]) for i in 1:length(NF)] +end=# + #=function get_full_edge_indexing(sig,r,edge,NF_::FastEdgeIterator,xs) NF=NF_.iterators[1] fullview, count = get_full_edge_indexing(NF.rays,NF.local_cone,NF.new_node,NF.active_nodes[1],length(NF.sig),NF.dim,NF_.ray_tol) @@ -535,7 +541,7 @@ function update_edge(NF_::FastEdgeIterator) end end # adjust orthogonal basis - nu[dim] .= my_cone[first_entry] + nu[dim] .= my_cone[first_entry] # take the random vector initially generated in this place b = true # print("$(my_minimal)") for i in 1:(start_index-1) diff --git a/src/extended.jl b/src/extended.jl index 6d57fb9..7322ee6 100644 --- a/src/extended.jl +++ b/src/extended.jl @@ -175,6 +175,88 @@ function search_vertex2(tree::ExtendedTree,point::MV,idx,dist) where {S,FLOAT<:R #error("") end +mutable struct PlaneSearchData{P,P2,HVN} + blocked::BitVector + c0::Float64 + direction::P + x0::P + r::P + vertex::P2 + num_nodes::Int64 + bestnode::MVector{1,Int64} + bestdist::MVector{1,Float64} + xs::HVN + PlaneSearchData(::Type{P},nn,xs) where P = new{P,MVector{size(P)[1],Float64},typeof(xs)}(falses(64+nn),0.0,zeros(P),zeros(P),zeros(P),zeros(MVector{size(P)[1],Float64}),nn,zeros(MVector{1,Int64}),zeros(MVector{1,Float64}),xs) +end +function block(data::PlaneSearchData,index) + len = length(data.blocked) + if index>len + append!(data.blocked,falses(index-len)) + end + data.blocked[index] = true +end + +function blocked(data::PlaneSearchData,index) + len = length(data.blocked) + if index>len + return false + else + return data.blocked[index] + end +end + +function reset(data,x0,r,u,xs,origin,searcher) + data.r = r + data.x0 = x0 + data.direction = u + data.vertex .= 0.0 + + tree = searcher.tree.tree.tree + + _min = tree.hyper_rec.mins + _max = tree.hyper_rec.maxes + m1 = data.vertex + l1 = length(m1) + for i in 1:l1 + if data.direction[i]>0 + m1[i] = _max[i] + end + end + + #resize!(data.blocked, bla ) + + c1 = maximum(dot(xs[g], u) for g in origin) + c2 = abs(c1) + data.c0 = c1 + c2*searcher.plane_tolerance + + len = length(searcher.tree.tree.tree.nodes) + #resize!(data.blocked,len+data.num_nodes) + fill!(data.blocked,false) +end + +function search_vertex_plane(tree::ExtendedTree,data) #where {S,FLOAT<:Real} # ,MV<:Union{MVector{S,FLOAT},SVector{S,FLOAT}} + #(a,b) = sum(i->(a[i]-b[i])^2,1:S) + exs = tree.extended_xs + s = tree.size + lm=tree.mirrors + data.bestnode[1]=0 + data.bestdist[1]=typemax(Float64) + for j in 1:lm + #b_skip(s+j) + x_new = exs[s+j] + mydist = sum(abs2, data.r - x_new) + if mydistdata.c0 + data.bestnode[1] = s+j + data.bestdist[1] = mydist + end + end + search_vertex_plane(tree.tree,data.r,data.bestnode,data.bestdist,data) + return data.bestnode[1], data.bestdist[1] + #error("") +end + + + function _nn(tree::ExtendedTree,x::Point,skip=(x->false))::Tuple{Int64,Float64} nn(tree,x,skip) diff --git a/src/fastpolyintegrator.jl b/src/fastpolyintegrator.jl index 034c7d9..2729dfe 100644 --- a/src/fastpolyintegrator.jl +++ b/src/fastpolyintegrator.jl @@ -24,7 +24,9 @@ struct Fast_Polygon_Integrator{T,TT,IDC<:IterativeDimensionChecker,D} function Fast_Polygon_Integrator(Integ::HVIntegral,integrand=nothing, bulk_integral=false) b_int=(typeof(integrand)!=Nothing) ? bulk_integral : false enable(Integ,volume=true,integral=b_int) - idc = IterativeDimensionChecker(dimension(mesh(Integ))) + idc = IterativeDimensionChecker(mesh(Integ)) + #@descend Fast_Polygon_Integrator{typeof(integrand),typeof(Integ),typeof(idc)}( integrand, b_int, Integ, idc ) + #error() PI=Fast_Polygon_Integrator{typeof(integrand),typeof(Integ),typeof(idc)}( integrand, b_int, Integ, idc ) return PI end @@ -42,8 +44,8 @@ function copy(I::Fast_Polygon_Integrator) return Fast_Polygon_Integrator{typeof(I._function),typeof(I.Integral)}(I._function,I.bulk,copy(I.Integral)) end -@inline function integrate(Integrator::Fast_Polygon_Integrator; domain, relevant, modified,progress) - _integrate(Integrator, domain=domain, calculate=modified, iterate=relevant,progress=progress) +@inline function integrate(Integrator::Fast_Polygon_Integrator, domain, relevant, modified, progress) + _integrate(Integrator, domain, modified, relevant, progress) end diff --git a/src/geometry.jl b/src/geometry.jl index 605fa85..0d83ac9 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -155,12 +155,36 @@ function VoronoiGeometry(xs::Points,b=Boundary(); vertex_storage=DatabaseVertexS mmm = cast_mesh(vertex_storage,copy(xs)) voronoi(mmm,searcher=Raycast(xs;domain=b,options=search),intro="",printsearcher=printevents, silence=silence) + #= + nod = nodes(mmm) + for i in 1:10 + print("$i: ") + for (sig,r) in vertices_iterator(mmm,i) + for s in sig + s>length(nod) && break + print("$(norm(r-nod[s])), ") + end + print(" || ") + end + println() + end + error() =# + d2 = Create_Discrete_Domain(mmm,b,intro="",search_settings=search) # periodized version including all boundary data improve_mesh(d2; improving..., printevents=printevents,search=search) lboundary = length(b) - integrate_geo(integrate,d2,myintegrator,integrand,mc_accurate,collect(1:public_length(d2)),collect(1:(length(mesh(d2))+lboundary)),silence) + relevant = collect(1:public_length(d2)) + modified = collect(1:(length(mesh(d2))+lboundary)) + #@descend integrate_geo(integrate,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence) + #error() + integrate_geo(integrate,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence) + #=m2,i2 = integrate_view(d2) + println(i2.neighbors) + for (sig,r) in vertices_iterator(m2,1) + print("$sig, ") + end=# result = VoronoiGeometry(myintegrator,d2,integrand,search,mc_accurate,nothing)#NoFile()) end catch err @@ -170,6 +194,7 @@ function VoronoiGeometry(xs::Points,b=Boundary(); vertex_storage=DatabaseVertexS redirect_stdout(oldstd) return result end +@inline internal_boundary(d2,myinte) = internal_boundary(d2) function integrate_geo(threading::MultiThread,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence=false) integral = integrate_view(d2).integral ref_mesh = mesh(integral) @@ -208,17 +233,28 @@ function integrate_geo(threading::MultiThread,d2,myintegrator,integrand,mc_accur progress = ThreadsafeProgressMeter(2*length(relevant),silence,intro) Threads.@threads for i in 1:length(parallel_integrals) #for i in 4:-1:1 - HighVoronoi.integrate(myintes[i],domain=internal_boundary(d2),relevant=relevants[i],modified=modifieds[i],progress=progress) + HighVoronoi.integrate(myintes[i],internal_boundary(d2,myintes[i]),relevants[i],modifieds[i],progress) end + return end parallelize_integrators(myintes2) = myintes2 @inline integrate_geo(integrate::SingleThread,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence=false) = integrate_geo(true,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence) @inline function integrate_geo(integrate::Bool,d2,myintegrator,integrand,mc_accurate,relevant,modified,silence=false) - !integrate && return + !integrate && return II2=Integrator(integrate_view(d2).integral,myintegrator,integrand=integrand,mc_accurate=mc_accurate) + + mmm = mesh(integrate_view(d2).integral) + nnn = nodes(mmm) + #println("START: Length=$(length(nnn)), $(nnn[81])") + #println(nnn) + #return + myinte=backup_Integrator(II2,true) intro="$(Integrator_Name(myinte))-integration over $(length(relevant)) cells:" - HighVoronoi.integrate(myinte,domain=internal_boundary(d2),relevant=relevant,modified=modified,progress=ThreadsafeProgressMeter(2*length(relevant),false,intro)) + #@descend HighVoronoi.integrate(myinte,internal_boundary(d2,myinte),relevant,modified,ThreadsafeProgressMeter(2*length(relevant),false,intro)) + #error() + HighVoronoi.integrate(myinte,internal_boundary(d2,myinte),relevant,modified,ThreadsafeProgressMeter(2*length(relevant),false,intro)) + return end function VoronoiGeometry(file::String,proto=nothing; _myopen=jldopen, offset="", search_settings=(__useless=0,), integrate=false, volume=true,area=true,bulk=false,interface=false,integrator=VI_GEOMETRY,integrand=nothing,mc_accurate=(1000,100,20),periodic_grid=nothing,silence=false) diff --git a/src/heuristic.jl b/src/heuristic.jl index 396b5d0..900ddef 100644 --- a/src/heuristic.jl +++ b/src/heuristic.jl @@ -31,8 +31,8 @@ function copy(I::Heuristic_Integrator) return Heuristic_Integrator{typeof(I._function),typeof(I.Integral)}(I._function,I.bulk,copy(I.Integral)) end -@inline function integrate(Integrator::Heuristic_Integrator; progress=ThreadsafeProgressMeter(0,true,""), domain=Boundary(), relevant=1:(length(Integrator.Integral)+length(domain)), modified=1:(length(Integrator.Integral))) - _integrate(Integrator; domain=domain, calculate=modified, iterate=relevant,progress=progress) +@inline function integrate(Integrator::Heuristic_Integrator, domain, relevant, modified,progress) + _integrate(Integrator, domain, modified, relevant, progress) end #function integrate(Integrator::Heuristic_Integrator; domain=Boundary(), relevant=1:(length(Integrator.Integral)+length(domain)), modified=1:(length(Integrator.Integral))) #println("PolyInt: ")#$(length(relevant)), $(length(modified))") diff --git a/src/heuristic_mc.jl b/src/heuristic_mc.jl index a04868a..508a4ee 100644 --- a/src/heuristic_mc.jl +++ b/src/heuristic_mc.jl @@ -31,8 +31,8 @@ function copy(I::HeuristicMCIntegrator) return HeuristicMCIntegrator{typeof(I.f)}(mc2.Integral,heu2,mc2,I.f) end -function integrate(Integrator::HeuristicMCIntegrator; progress=ThreadsafeProgressMeter(0,true,""), domain=Boundary(), relevant=1:(length(Integrator.Integral)+length(domain)), modified=1:(length(Integrator.Integral))) - _integrate(Integrator; domain=domain, calculate=modified, progress=progress, iterate=relevant) +function integrate(Integrator::HeuristicMCIntegrator, domain, relevant, modified, progress) + _integrate(Integrator, domain, modified, relevant, progress) end diff --git a/src/hvdatabase.jl b/src/hvdatabase.jl index 8bb0f38..d1dc5d1 100644 --- a/src/hvdatabase.jl +++ b/src/hvdatabase.jl @@ -41,6 +41,7 @@ mutable struct HeapDataBase{P, RI, Q<:QueueHashTable} <: HVDataBase{P, Vector{In buffer_r::Vector{Vector{Float64}} lock::ReadWriteLock nthreads::Int64 + set::Dict{Vector{Int64},Int64} # Konstruktor für HeapDataBase function HeapDataBase(::Type{P},_blocksize::Int64) where {P} @@ -50,7 +51,7 @@ mutable struct HeapDataBase{P, RI, Q<:QueueHashTable} <: HVDataBase{P, Vector{In queue = QueueHashTable(_blocksize,SingleThread()) nth = Threads.nthreads() #new{P, RI, typeof(queue)}(data, floatview, queue, _blocksize,_blocksize #=position=blocksize??=#,[Int64[] for i in 1:nth],[Float64[] for i in 1:nth],ReadWriteLock(),nth) - new{P, RI, typeof(queue)}(data, floatview, queue, _blocksize,_blocksize,[Int64[] for i in 1:nth],[Float64[] for i in 1:nth],ReadWriteLock(),nth) + new{P, RI, typeof(queue)}(data, floatview, queue, _blocksize,_blocksize,[Int64[] for i in 1:nth],[Float64[] for i in 1:nth],ReadWriteLock(),nth,Dict{Vector{Int64},Int64}()) end end indexvector(::HDB) where {HDB<:HeapDataBase} = Vector{Int64}() @@ -190,3 +191,72 @@ function Base.delete!(hdb::HDB,index,key) where {P,HDB<:HeapDataBase{P}} #print("/") end + +function push_facet!(hdb::HDB,data::P) where {HDB<:HeapDataBase,P<:Tuple} + sig = data[1] + r = data[2] + u = data[3] + sort!(sig) + pos = 0 + scs = copy(sig) + #if haskey(hdb.set,scs) + # println() + # println("----- KENN ICH SCHON !!! ") + #end + #push!(hdb.set,scs=>1) + status(41) + writelock(hdb.lock) + if !(pushqueue!(hdb.keys,sig)) # works like a haskey but also pushes entry if not known + status(42) + status(43) + push_entry!(hdb,length(sig)) + pos = hdb.position + hdb.blocksize*(length(hdb.data)-1) + push!(hdb,sig) + push!(hdb,r) + push!(hdb,u) + status(44) + else + status(45) + #println() + #println("KENN ICH SCHON !!! ") + end + writeunlock(hdb.lock) + return pos +end + + +function get_facet(hdb::HDB,index) where {P,HDB<:HeapDataBase{P}} + id = Threads.threadid() + index==0 && error("severe error: index shall never be '0' !") + id>hdb.nthreads && rescale_db(hdb) + pos = mod(index,hdb.blocksize) + block = div(index,hdb.blocksize)+1 + if pos==0 + pos = hdb.blocksize + block -= 1 + end + readlock(hdb.lock) + #if pos==0 || block==0 + # open("protokol.txt","a") do f + # write(f,"Pos: $pos, Block: $block, Index: $index \n") + # end + #end + #@inbounds + l = hdb.data[block][pos] + if l<0 + readunlock(hdb.lock) + return (view(hdb.buffer_sig[id],1:0),zeros(P)) + end + sig = read_from_db(hdb.data,hdb.buffer_sig[id],pos,block,l,hdb.blocksize) + pos2 = pos + l + if (pos2>=hdb.blocksize) + pos2 -= hdb.blocksize + block += 1 + end + data = read_from_db(hdb.floatview,hdb.buffer_r[id],pos2,block,2*length(P),hdb.blocksize) + r = P(view(data,1:length(P))) + u = P(view(data,(length(P)+1):(2*length(P)))) + readunlock(hdb.lock) + return (sig,r,u) +end + diff --git a/src/hvlocks.jl b/src/hvlocks.jl index 7eb8f17..1c8277f 100644 --- a/src/hvlocks.jl +++ b/src/hvlocks.jl @@ -126,13 +126,10 @@ struct ReadWriteLock head::Threads.Atomic{Int64} tail::Threads.Atomic{Int64} reads_count::Threads.Atomic{Int64} - all::Threads.Atomic{Int64} - index::Int64 - trace::Vector{RWLTrace} - lock::SpinLock - #counts_read::Vector{Threads.Atomic{Int64}} - #counts_write::Vector{Threads.Atomic{Int64}} - #lock::ReentrantLock + #all::Threads.Atomic{Int64} + #index::Int64 + #trace::Vector{RWLTrace} + #lock::SpinLock end @inline readlock(::Nothing) = nothing @@ -142,82 +139,41 @@ end @inline ReadWriteLock(::T) where T = nothing function ReadWriteLock() - #nn = Threads.nthreads() - #cr = [Threads.Atomic{Int64}(0) for i in 1:nn] - #cw = [Threads.Atomic{Int64}(0) for i in 1:nn] - return ReadWriteLock(Threads.Atomic{Int64}(0),Threads.Atomic{Int64}(0),Threads.Atomic{Int64}(0),Threads.Atomic{Int64}(0),atomic_add!(HighVoronoi.RWLCounter,1),Vector{RWLTrace}(undef,100),SpinLock())#,cr,cw,ReentrantLock()) + return ReadWriteLock(Threads.Atomic{Int64}(0),Threads.Atomic{Int64}(0),Threads.Atomic{Int64}(0))#,Threads.Atomic{Int64}(0),atomic_add!(HighVoronoi.RWLCounter,1),Vector{RWLTrace}(undef,100),SpinLock())#,cr,cw,ReentrantLock()) end @inline function readlock(rwl::ReadWriteLock) - #check(rwl) -# lock(rwl.lock) - ti = time_ns() this_tail = atomic_add!(rwl.tail,1) -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),1,a) -# unlock(rwl.lock) + #print("$(Threads.threadid())+") ii = 0 while atomic_add!(rwl.head,0)1000000000 -# lock(rwl.lock) -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),-1,a) -# s = "$(rwl.index): $(atomic_add!(rwl.head,0)), $(this_tail), $(atomic_add!(rwl.reads_count,0)), $ti \n", rwl.trace,"\n" -# unlock(rwl.lock) - error(s) - end end atomic_add!(rwl.reads_count,1) atomic_add!(rwl.head,1) end @inline function writelock(rwl::ReadWriteLock) - #check(rwl) -# lock(rwl.lock) - ti = time_ns() this_tail = atomic_add!(rwl.tail,1) - #print("$this_tail, $(rwl.head[]), $(rwl.reads_count[])") -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),3,a) -# unlock(rwl.lock) + #print("$(Threads.threadid())*") ii = 0 while atomic_add!(rwl.head,0)0 active_wait(100) ii += 1 mod(ii,100)==0 && yield() - if time_ns()-ti>1000000000 -# lock(rwl.lock) -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),-3,a) -# s = "$(rwl.index): $(atomic_add!(rwl.head,0)), $(this_tail), $(atomic_add!(rwl.reads_count,0)), $ti \n", rwl.trace,"\n" -# unlock(rwl.lock) - error(s) - end - #atomic_fence() end end @inline function readunlock(rwl::ReadWriteLock) -# lock(rwl.lock) + #print("$(Threads.threadid())-") atomic_sub!(rwl.reads_count,1) -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),2,a) -# unlock(rwl.lock) - - #atomic_sub!(rwl.counts_read[Threads.threadid()],1) end @inline function writeunlock(rwl::ReadWriteLock) -# lock(rwl.lock) + #print("$(Threads.threadid())_") atomic_add!(rwl.head,1) -# a = atomic_add!(rwl.all,1) -# rwl.trace[mod(a,100)+1] = RWLTrace(Threads.threadid(),4,a) -# unlock(rwl.lock) - - #atomic_sub!(rwl.counts_write[Threads.threadid()],1) end @inline Base.lock(rwl::ReadWriteLock) = writelock(rwl) diff --git a/src/improving.jl b/src/improving.jl index c180795..d71cbed 100644 --- a/src/improving.jl +++ b/src/improving.jl @@ -22,7 +22,7 @@ function improve_mesh(_domain; max_iterations=0, tolerance=1.0, printevents,sear modified = zeros(Bool,lmesh) modified_i = zeros(Bool,lmesh) buffer = zeros(Float64,dim) - integrate(Integrator,domain=b,relevant=1:plmesh,modified=1:plmesh) + integrate(Integrator,b,1:plmesh,1:plmesh) for i in 1:plmesh modified_i .= false buffer .= 0.0 diff --git a/src/integrate.jl b/src/integrate.jl index f7bbd09..73db336 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -72,9 +72,9 @@ end _NeighborFinder(dim,x) = NeighborFinder(dim,x)#dim>=UseNeighborFinderDimension ? NeighborFinder(dim,x) : nothing -struct IntegrateData{T,VP,TT} +struct IntegrateData{T,VP,TT,B} extended_xs::VP - domain::Boundary + domain::B size::Int64 active::BitVector float_vec_buffer::Vector{Float64} @@ -87,6 +87,8 @@ struct IntegrateData{T,VP,TT} deprecated::Vector{Bool} buffer_data::TT end +const SphereIntegrateData{T,VP,TT} = IntegrateData{T,VP,TT,SP} where {T,VP,TT,SP<:SphereBoundary} + function IntegrateData(xs::HV,dom,tt::TT) where {P,HV<:HVNodes{P},TT} return _IntegrateData(xs,dom,0) end @@ -94,13 +96,14 @@ end function _IntegrateData(xs::HV,dom,tt) where {P,HV<:HVNodes{P}} dim = size(P)[1]#length(xs[1]) l=length(dom) + #println("length: $(length(xs)) with $(xs[81])") m=append!(copy(xs),Vector{P}(undef,l)) a=falses(l) #BitVector(zeros(Int8,l)) nf = NeighborFinder(dim,zeros(P)) c = Vector{Int64}(undef,length(m)) a = Vector{Bool}(undef,length(m)) d = Vector{Bool}(undef,length(m)) - return IntegrateData{typeof(nf),typeof(m),typeof(tt)}(m,dom,length(xs),a,Float64[],(Vector{Float64})[],length(zeros(P)),nf,c,a,d,tt) + return IntegrateData{typeof(nf),typeof(m),typeof(tt),typeof(dom)}(m,dom,length(xs),a,Float64[],(Vector{Float64})[],length(zeros(P)),nf,c,a,d,tt) end function activate_data_cell(tree,_Cell,neigh) @@ -116,6 +119,102 @@ function activate_data_cell(tree,_Cell,neigh) end end +function activate_data_cell(tree::SphereIntegrateData,_Cell,neigh) + tree.active .= false + lxs=tree.size + extended_xs = tree.extended_xs + for n in neigh + if n>lxs + #plane=n-lxs + #tree.active[plane] && continue + #tree.active[plane]=true + sphere_boundary = tree.domain + mesh = sphere_boundary.mesh + + center = sphere_boundary.center + radius = sphere_boundary.radius + + all_verts = sphere_boundary.all_vertices + l_all = length(all_verts) + count = 0 + count2 = 0 + for (sig,r) in vertices_iterator(mesh,_Cell) + dist = norm(r-center) + count2 += 1 + abs(radius-dist)/radius > 0.5 && continue + count += 1 + if (count>l_all) + l_all *= 2 + resize!(all_verts, l_all) + end + all_verts[count] = r + end + first_r = all_verts[1] + facet_center = sum(view(all_verts,1:count))/count + for kk in 1:count + all_verts[kk] -= facet_center + end + + onb = sphere_boundary.onb + dim = length(sphere_boundary.onb) + onb[dim] = normalize(tree.extended_xs[_Cell]-sphere_boundary.center) + for i in 1:(dim-1) + maximal_angle = 0.0 + max_index = 0 + for kk in 1:count + this_angle = abs(dot(onb[dim],all_verts[kk])) + if this_angle>maximal_angle + maximal_angle = this_angle + max_index = kk + end + end + onb[i] = all_verts[max_index] + for kk in 1:(i-1) + onb[i] -= onb[kk] * dot(onb[kk],onb[i]) + onb[i] -= onb[kk] * dot(onb[kk],onb[i]) + end + onb[i] = normalize(onb[i]) + onb[dim] -= onb[i] * dot(onb[dim],onb[i]) + onb[dim] -= onb[i] * dot(onb[dim],onb[i]) + onb[dim] = normalize(onb[dim]) + end + #= + for i in 1:count + for j in (i+1):count + print("($i,$j,$(dot(all_verts[i]-all_verts[j],onb[dim]))), ") + end + end + println() + for i in 1:count + print("($i,$(dot(all_verts[i],onb[dim]))), ") + end + error()=# + x0 = extended_xs[_Cell] + b_node = x0 + 2*dot(first_r-x0,onb[dim])*onb[dim] + tree.extended_xs[lxs+1] = b_node #reflect(tree.extended_xs[_Cell],tree.domain,plane) + #=println(x0) + println(b_node) + println(length(extended_xs)) + println(lxs) + for i in 1:dim + for j in i:dim + print("($i,$j,$(dot(onb[i],onb[j]))), ") + end + end + println() + for kk in 1:count + all_verts[kk] += facet_center + end + + for i in 1:count + print("$(norm(all_verts[i]-b_node)-norm(all_verts[i]-x0)) - ") + end + error("")=# + + end + end +end + function neighbors_of_cell(_Cell::Int,mesh::Voronoi_MESH,data::IntegrateData, condition = r->true) adj = neighbors_of_cell(_Cell,mesh,adjacents=true) activate_data_cell(data,_Cell,adj) @@ -133,34 +232,38 @@ For each implemented Integrator type this method shall be overwritten. In particular, the passage of calculate and iterate might be modified according to the needs of the respective class. See also Polygon_Integrator and Montecarlo_Integrator for reference. """ -function integrate(Integrator; progress=ThreadsafeProgressMeter(0,true,""), domain=Boundary(), relevant=1:(length(Integrator.Integral)+length(domain)), modified=1:(length(Integrator.Integral))) - _integrate(Integrator; progress=progress,domain=domain, calculate=modified, iterate=relevant) +function integrate(Integrator, domain, relevant, modified, progress) + _integrate(Integrator, domain, modified, relevant, progress) end +integrate(a,b,relevant,d) = integrate(a,b,relevant,d,ThreadsafeProgressMeter(2*length(relevant),false,"")) + + @inline function __integrate_getdata(I_data::Nothing,Integral,domain,Integrator) nn = nodes(mesh(Integral)) IntegrateData(nn,domain,Integrator) end __integrate_getdata(I_data,Integral,domain,Integrator) = I_data +_integrate(Integrator, domain, calculate, iterate,intro::String) = _integrate(Integrator, domain, calculate, iterate,ThreadsafeProgressMeter(2*length(iterate),false,intro)) """ Iterates integrate_cell over all elements of iterate. It thereby passes the information on whether volume, areas, bulk- or surface integrals shall be calculated. """ -function _integrate(Integrator; domain=Boundary(), calculate, iterate, # =1:(length(Integrator.Integral)+length(domain)), iterate=1:(length(Integrator.Integral)), - I_data=nothing, compact=false, intro="$(Integrator_Name(Integrator))-integration over $(length(collect(iterate))) cells:",progress=ThreadsafeProgressMeter(2*length(iterate),false,intro)) +function _integrate(Integrator, domain, calculate, iterate,progress; # =1:(length(Integrator.Integral)+length(domain)), iterate=1:(length(Integrator.Integral)), + I_data=nothing)#, compact=false, intro="$(Integrator_Name(Integrator))-integration over $(length(collect(iterate))) cells:") TODO=collect(iterate) - try #vp_print(0,intro) - position_0 = length(intro)+5 + #position_0 = length(intro)+5 #vp_print(position_0-5," \u1b[0K") - Integral=Integrator.Integral + Integral=Integrator.Integral data = __integrate_getdata(I_data,Integral,domain,Integrator) if length(TODO)==0 - vp_print(position_0,"nothing to integrate") + #vp_print(position_0,"nothing to integrate") return Integrator, data end + try vol=enabled_volumes(Integral) ar=enabled_area(Integral) bulk=enabled_bulk(Integral) @@ -172,12 +275,22 @@ function _integrate(Integrator; domain=Boundary(), calculate, iterate, # =1:(len count=0 bb = typeof(Integrator.Integral)<:ThreadsafeIntegral # println("Hallo") +#println("integrating.... $TODO_count, $TODO") #println(Integrator.Integral.area) for k in 1:TODO_count # initialize and array of length "length(xs)" to locally store verteces of cells #vp_print(position_0,"Cell $(string(TODO[k], base=10, pad=max_string_i)) (in cycle: $(string(k, base=10, pad=max_string_todo)) of $TODO_count)") #@descend integrate_cell(vol,ar,bulk,inter,TODO[k],iterate, calculate, data,Integrator) #error("") - V=integrate_cell(vol,ar,bulk,inter,TODO[k],iterate, calculate, data,Integrator) + #print("$(TODO[k]), ") + V = 0.0 + err = false + try + V=integrate_cell(vol,ar,bulk,inter,TODO[k],iterate, calculate, data,Integrator) + catch + err = true + println(TODO[k]) + end + err && error() #print(Threads.threadid()) if vol vol_sum+=V #Integral.volumes[TODO[k]] @@ -211,15 +324,16 @@ function _integrate(Integrator; domain=Boundary(), calculate, iterate, # =1:(len #vp_line_up(1) #if (!compact) vp_line() end #println("Differenz: $(vol_sum-V1)") - return Integrator,data catch e open("error_log$(Threads.threadid()).txt", "w") do f # Stacktrace speichern Base.showerror(f, e, catch_backtrace()) end + rethrow() #sync(s) #sync(s) end +return Integrator,data end function integrate_cell(vol::Bool,ar::Bool,bulk::Bool,inter::Bool, _Cell, iterate, calculate, data, Integrator::Nothing) @@ -232,7 +346,7 @@ afterwards it calls the true integration function that is provided by the Integr """ function integrate_cell(vol::Bool,ar::Bool,bulk::Bool,inter::Bool, _Cell, iterate, calculate, data, Integrator) I=Integrator.Integral - + #println("Letzte: $(data.extended_xs[81]), $(data.extended_xs[82])") adj = neighbors_of_cell(_Cell,mesh(I),adjacents=true) activate_data_cell(data,_Cell,adj) new_neighbors = neighbors_of_cell(_Cell,mesh(I),extended_xs=data.extended_xs,edgeiterator=data.NFfind, neighbors=adj) diff --git a/src/mcintegrator.jl b/src/mcintegrator.jl index 99089cc..f874589 100644 --- a/src/mcintegrator.jl +++ b/src/mcintegrator.jl @@ -59,8 +59,8 @@ function copy(I::Montecarlo_Integrator) return Montecarlo_Integrator(copy(I.Integral),I.bulk,I.interface,nmc_bulk=I.NMC_bulk,nmc_interface=I.NMC_interface,recycle=I._recycle,cal_area=I.area) end -function integrate(Integrator::Montecarlo_Integrator; progress=ThreadsafeProgressMeter(0,true,""), domain=Boundary(), relevant=1:(length(mesh(Integrator.Integral))+length(domain)), modified=1:(length(mesh(Integrator.Integral)))) - _integrate(Integrator; domain=domain, calculate=relevant, progress=progress, iterate=Base.intersect(relevant,1:(length(mesh(Integrator.Integral))))) +function integrate(Integrator::Montecarlo_Integrator, domain, relevant, modified, progress) + _integrate(Integrator, domain, modified, relevant, progress) end @@ -97,7 +97,7 @@ function integrate(neighbors,_Cell,iterate, calculate, data,Integrator::Montecar vec = Float64[] vecvec = [vec] I=Integrator - xs=data.extended_xs + xs = data.extended_xs if mod(I.cycle[1], I._recycle)==0 new_directions!(I.directions,length(xs[1])) I.cycle[1]=0 @@ -107,6 +107,12 @@ function integrate(neighbors,_Cell,iterate, calculate, data,Integrator::Montecar x = xs[_Cell] d = data.dimension lneigh = length(neighbors) + if _Cell==1 + println(x) + for n in neighbors + println("$n: $(xs[n])") + end + end # Bulk computations: V stores volumes y stores function values in a vector format V = 0.0 diff --git a/src/nodes.jl b/src/nodes.jl index ce1f0a4..823a086 100644 --- a/src/nodes.jl +++ b/src/nodes.jl @@ -251,7 +251,7 @@ SearchTree(nodes::ACN,type=HVKDTree) where {P<:Point,S,ACN<:SubArray{P,S,NodesVi ##################################################################################################################### -## ReflectedNodes +## ReflectedNodes ##################################################################################################################### diff --git a/src/periodicmesh.jl b/src/periodicmesh.jl index 3817cf4..38aa733 100644 --- a/src/periodicmesh.jl +++ b/src/periodicmesh.jl @@ -477,7 +477,7 @@ function periodic_final_integration(Integral,MESH,xs,search,domain,myintegrator, bI = backup_Integrator(II2,true) l1 = public_length(d2) l_2 = (length(mesh(d2))+lboundary) - HighVoronoi.integrate(bI,domain=internal_boundary(domain),relevant=1:l1,modified=1:l_2,progress = ThreadsafeProgressMeter(l1,false,"$(Integrator_Name(bI))-integration over $(l1) cells:")) + HighVoronoi.integrate(bI, internal_boundary(domain), 1:l1, 1:l_2, ThreadsafeProgressMeter(l1,false,"$(Integrator_Name(bI))-integration over $(l1) cells:")) end diff --git a/src/polyintegrator.jl b/src/polyintegrator.jl index 1496a31..9228158 100644 --- a/src/polyintegrator.jl +++ b/src/polyintegrator.jl @@ -30,8 +30,8 @@ function copy(I::Polygon_Integrator) return Polygon_Integrator{typeof(I._function),typeof(I.Integral)}(I._function,I.bulk,Inte, IterativeDimensionChecker(length(Inte.MESH.nodes[1]))) end -@inline function integrate(Integrator::Polygon_Integrator; progress=ThreadsafeProgressMeter(0,true,""),domain=Boundary(), relevant=1:(length(Integrator.Integral)+length(domain)), modified=1:(length(Integrator.Integral))) - _integrate(Integrator; domain=domain, calculate=modified, iterate=relevant,progress=progress) +@inline function integrate(Integrator::Polygon_Integrator, domain, relevant, modified, progress) + _integrate(Integrator, domain, modified, relevant, progress) end diff --git a/src/progress.jl b/src/progress.jl index 4989015..1ce35f4 100644 --- a/src/progress.jl +++ b/src/progress.jl @@ -15,6 +15,11 @@ end function ThreadsafeProgressMeter(n::Int,silence::Bool,intro,::MultiThread) return ThreadsafeProgressMeter(Progress(n,intro),silence,BusyFIFOLock()) end + +function ThreadsafeProgressMeter(intro::String) + return ThreadsafeProgressMeter(Progress(0;desc=intro,showspeed=true),false,BusyFIFOLock()) +end + #import Progress @inline function next!(tpm::TPM) where {TPM<:ThreadsafeProgressMeter} diff --git a/src/queuehashing.jl b/src/queuehashing.jl index 7956758..9e53096 100644 --- a/src/queuehashing.jl +++ b/src/queuehashing.jl @@ -22,7 +22,7 @@ mutable struct QueueHashTable{V<:AbstractVector{HashedQueue}, R} end # Method for inserting into the QueueHashTable -function pushqueue!(ht::Q, key::K, write::Bool=true) where {Q<:QueueHashTable,K} +function pushqueue!(ht::Q, key::K, write::Bool=true,level=0) where {Q<:QueueHashTable,K} value = fnv1a_hash(key, UInt128) value_ = UInt64(value & UInt128(0x7FFFFFFFFFFFFFFF)) index1 = value_ & ht.mylength[1] @@ -30,7 +30,9 @@ function pushqueue!(ht::Q, key::K, write::Bool=true) where {Q<:QueueHashTable,K} i = UInt64(0) ret = false + status(51+level*100) lock(ht.lock) + status(52+level*100) try while true idx = reinterpret(Int64, (index1 + i * index2) & ht.mylength[1] + 1) # try bitcast( ) instead @@ -48,13 +50,15 @@ function pushqueue!(ht::Q, key::K, write::Bool=true) where {Q<:QueueHashTable,K} break # return false elseif data.value == value && data.value2 == index2 ret = true - break # return false + break # return true end i += 1 if i >= ht.mylength[1] if write + status(53+level*100) extend(ht) - ret = pushqueue!(ht, key) + status(54+level*100) + ret = pushqueue!(ht, key,true,level+1) end break end @@ -62,6 +66,7 @@ function pushqueue!(ht::Q, key::K, write::Bool=true) where {Q<:QueueHashTable,K} finally unlock(ht.lock) end + status(55+level*100) return ret end @inline Base.haskey(ht::Q, key::K) where {Q<:QueueHashTable,K} = pushqueue!(ht,key,false) @@ -103,6 +108,7 @@ end # Function to extend the QueueHashTable function extend(ht::QueueHashTable) + #println("ex queue keys") len2 = 2 * (ht.mylength[1] + 1) V2 = similarqueuehash(HashedQueue, reinterpret(Int64, len2)) new_occupied = falses(len2) diff --git a/src/raycast-types.jl b/src/raycast-types.jl index 885b743..b01b4ab 100644 --- a/src/raycast-types.jl +++ b/src/raycast-types.jl @@ -145,6 +145,84 @@ struct RaycastParameter{FLOAT,TREE,R,T,C} copynodes::C end +struct NewRaycastParameter{FLOAT,TREE,R,T,C,LF} + variance_tol::FLOAT + break_tol::FLOAT + b_nodes_tol::FLOAT + plane_tolerance::FLOAT + ray_tol::FLOAT + nntree::TREE + method::R + threading::T + copynodes::C + locking_functions::LF +end +struct NewRaycastParameter_Store_1{FLOAT,TREE,R,T,C,LF} + variance_tol::FLOAT + break_tol::FLOAT + b_nodes_tol::FLOAT + plane_tolerance::FLOAT + ray_tol::FLOAT + nntree::TREE + method::R + threading::T + copynodes::C + locking_functions::LF +end +NewRaycastParameter_Store_1(p::NewRaycastParameter{FLOAT, TREE, R, T, C, LF}) where {FLOAT, TREE, R, T, C, LF} = +NewRaycastParameter_Store_1( + p.variance_tol, + p.break_tol, + p.b_nodes_tol, + p.plane_tolerance, + p.ray_tol, + p.nntree, + p.method, + p.threading, + p.copynodes, + p.locking_functions +) +NewRaycastParameter(p::NewRaycastParameter_Store_1{FLOAT, TREE, R, T, C, LF}) where {FLOAT, TREE, R, T, C, LF} = +NewRaycastParameter( + p.variance_tol, + p.break_tol, + p.b_nodes_tol, + p.plane_tolerance, + p.ray_tol, + p.nntree, + p.method, + p.threading, + p.copynodes, + p.locking_functions +) +NewRaycastParameter(p::RaycastParameter{FLOAT, TREE, R, T, C}) where {FLOAT, TREE, R, T, C} = +NewRaycastParameter( + p.variance_tol, + p.break_tol, + p.b_nodes_tol, + p.plane_tolerance, + p.ray_tol, + p.nntree, + p.method, + p.threading, + p.copynodes, + nothing +) + +JLD2.writeas(::Type{NewRaycastParameter{P, T , TT, X1, X2, X3}}) where {P, T, TT, X1, X2, X3} = NewRaycastParameter_Store_1{P, T , TT, X1, X2, X3} +JLD2.wconvert(::Type{NewRaycastParameter_Store_1{P, T , TT, X1, X2, X3}},domain::NewRaycastParameter{P, T , TT, X1, X2, X3} ) where {P, T, TT, X1, X2, X3} = + NewRaycastParameter_Store_1(domain) +JLD2.rconvert(::Type{NewRaycastParameter{P, T , TT, X1, X2, X3}},store::NewRaycastParameter_Store_1{P, T , TT, X1, X2, X3}) where {P, T, TT, X1, X2, X3} = + NewRaycastParameter(store) + +JLD2.writeas(::Type{RaycastParameter{P, T , TT, X1, X2}}) where {P, T, TT, X1, X2} = NewRaycastParameter_Store_1{P, T , TT, X1, X2, Nothing} +JLD2.wconvert(::Type{NewRaycastParameter_Store_1{P, T , TT, X1, X2, X3}},domain::RaycastParameter{P, T , TT, X1, X2} ) where {P, T, TT, X1, X2, X3} = + NewRaycastParameter_Store_1(NewRaycastParameter(domain)) +JLD2.rconvert(::Type{NewRaycastParameter{P, T , TT, X1, X2, Nothing}},store::RaycastParameter{P, T , TT, X1, X2}) where {P, T, TT, X1, X2} = + NewRaycastParameter(store) + + + @inline variance_tol(::Type{Float64}) = 1.0E-15 @inline break_tol(::Type{Float64}) = 1.0E-5 @inline b_nodes_tol(::Type{Float64}) = 1.0E-7 @@ -165,10 +243,26 @@ end struct Raycast_Original end const RCOriginal = Raycast_Original() +struct Raycast_Original_Safety end +const RCOriginalSafety = Raycast_Original_Safety() struct Raycast_Non_General end const RCNonGeneral = Raycast_Non_General() +struct Raycast_Non_General_Skip end +const RCNonGeneralSkip = Raycast_Non_General_Skip() struct Raycast_Combined end const RCCombined = Raycast_Combined() + +abstract type Raycast_HP end +struct Raycast_Non_General_HP<:Raycast_HP end +const RCNonGeneralHP = Raycast_Non_General_HP() +struct Raycast_Original_HP<:Raycast_HP end +const RCOriginalHP = Raycast_Original_HP() +RaycastHP(::Raycast_Original) = RCOriginalHP +RaycastHP(::Raycast_Non_General) = RCNonGeneralHP +RaycastHP(::T) where T = begin + println("hier") + RCNonGeneralHP +end const RCCopyNodes = statictrue const RCNoCopyNodes = staticfalse @@ -203,36 +297,37 @@ RaycastParameter{FLOAT}(;variance_tol = variance_tol(FLOAT), ray_tol = ray_tol(FLOAT), nntree = VI_KD, method= RCStandard, threading = SingleThread(), - copynodes = RCCopyNodes, kwargs... ) where {FLOAT<:Real} = RaycastParameter{FLOAT,typeof(nntree),typeof(method),typeof(threading),typeof(copynodes)}(FLOAT(variance_tol), + copynodes = RCCopyNodes, locking_functions = nothing, kwargs... ) where {FLOAT<:Real} = NewRaycastParameter{FLOAT,typeof(nntree),typeof(method),typeof(threading),typeof(copynodes),typeof(locking_functions)}( + FLOAT(variance_tol), FLOAT(break_tol), FLOAT(b_nodes_tol), FLOAT(plane_tolerance), - FLOAT(ray_tol), nntree, method, threading, copynodes ) + FLOAT(ray_tol), nntree, method, threading, copynodes, locking_functions ) -RaycastParameter{FLOAT}(RP::RaycastParameter;variance_tol = FLOAT(RP.variance_tol), +RaycastParameter{FLOAT}(RP::NewRaycastParameter;variance_tol = FLOAT(RP.variance_tol), break_tol = FLOAT(RP.break_tol), b_nodes_tol = FLOAT(RP.b_nodes_tol), plane_tolerance = FLOAT(RP.plane_tolerance), ray_tol = FLOAT(RP.ray_tol), nntree = RP.nntree, method=RP.method, threading = RP.threading, - copynodes = RP.copynodes, kwargs... ) where {FLOAT<:Real} = RaycastParameter{FLOAT,typeof(nntree),typeof(method),typeof(threading),typeof(copynodes)}(FLOAT(variance_tol), + copynodes = RP.copynodes, locking_functions = RP.locking_functions, kwargs... ) where {FLOAT<:Real} = NewRaycastParameter{FLOAT,typeof(nntree),typeof(method),typeof(threading),typeof(copynodes),typeof(locking_functions)}(FLOAT(variance_tol), FLOAT(break_tol), FLOAT(b_nodes_tol), FLOAT(plane_tolerance), - FLOAT(ray_tol), nntree, method, threading, copynodes ) + FLOAT(ray_tol), nntree, method, threading, copynodes, locking_functions ) RaycastParameter(::Type{T};kwargs...) where T = RaycastParameter{T}(;kwargs...) -RaycastParameter(::Type{T},RP::RaycastParameter;kwargs...) where T = RaycastParameter{T}(RP;kwargs...) +RaycastParameter(::Type{T},RP::NewRaycastParameter;kwargs...) where T = RaycastParameter{T}(RP;kwargs...) RaycastParameter(::Type{T},NT::NamedTuple;kwargs...) where T = RaycastParameter{T}(;NT...,kwargs...) -RaycastParameter(r1::RaycastParameter{FLOAT},RP::RaycastParameter) where FLOAT = RaycastParameter(FLOAT,r1,variance_tol = FLOAT(RP.variance_tol), +RaycastParameter(r1::NewRaycastParameter{FLOAT},RP::NewRaycastParameter) where FLOAT = RaycastParameter(FLOAT,r1,variance_tol = FLOAT(RP.variance_tol), break_tol = FLOAT(RP.break_tol), b_nodes_tol = FLOAT(RP.b_nodes_tol), plane_tolerance = FLOAT(RP.plane_tolerance), ray_tol = FLOAT(RP.ray_tol), - nntree = RP.nntree, method=RP.method, threading=RP.threading, copynodes = RP.copynodes) -RaycastParameter(r1::RaycastParameter{FLOAT},kwargs::NamedTuple) where FLOAT = RaycastParameter(FLOAT,r1;kwargs...) -Base.merge(r1::RaycastParameter{FLOAT},RP::RaycastParameter) where FLOAT = RaycastParameter(r1,RP) -Base.merge(r1::RaycastParameter{FLOAT},tup::NamedTuple) where FLOAT = RaycastParameter(FLOAT,r1;tup...) + nntree = RP.nntree, method=RP.method, threading=RP.threading, copynodes = RP.copynodes, locking_functions = RP.locking_functions) +RaycastParameter(r1::NewRaycastParameter{FLOAT},kwargs::NamedTuple) where FLOAT = RaycastParameter(FLOAT,r1;kwargs...) +Base.merge(r1::NewRaycastParameter{FLOAT},RP::NewRaycastParameter) where FLOAT = RaycastParameter(r1,RP) +Base.merge(r1::NewRaycastParameter{FLOAT},tup::NamedTuple) where FLOAT = RaycastParameter(FLOAT,r1;tup...) struct MiniRaycast{T,B} tree::T @@ -256,17 +351,17 @@ function getMultiThreadRaycasters(rc::RR,meshes::PP) where {RR,PP} return rcs end -mutable struct RaycastIncircleSkip{T,TTT,TTTTT,TTTTTT,FEI,PA,FLOAT} +mutable struct RaycastIncircleSkip{T,TTT,TTTTT,TTTTTT,FEI,PA,FLOAT,CORRECTOR_FLOAT,HP} tree::T lmesh::Int64 lboundary::Int64 visited::Vector{Int64} ts::Vector{FLOAT} positions::BitVector - vectors::Matrix{Float64} - symmetric::Matrix{FLOAT} - rhs::Vector{Float64} - rhs_cg::Vector{Float64} + vectors::Matrix{CORRECTOR_FLOAT} + symmetric::Matrix{CORRECTOR_FLOAT} + rhs::Vector{CORRECTOR_FLOAT} + rhs_cg::Vector{CORRECTOR_FLOAT} ddd::Vector{FLOAT} domain::Boundary rare_events::Vector{Int64} @@ -278,17 +373,34 @@ mutable struct RaycastIncircleSkip{T,TTT,TTTTT,TTTTTT,FEI,PA,FLOAT} find_general_edgeiterator::TTTTTT FEIStorage_global::FEI parameters::PA + hp_vars::HP end -function RaycastIncircleSkip(xs_::HN,dom,parameters::RaycastParameter{FLOAT,TREE}) where {P,HN<:HVNodes{P},FLOAT,TREE} + +struct HP_corrector_data{CORRECTOR_FLOAT} + vectors::Matrix{CORRECTOR_FLOAT} + symmetric::Matrix{CORRECTOR_FLOAT} + rhs::Vector{CORRECTOR_FLOAT} + rhs_cg::Vector{CORRECTOR_FLOAT} + function HP_corrector_data(dim,::Type{T}) where T + return new{T}(zeros(T,dim,dim),zeros(T,dim,dim),zeros(T,dim),zeros(T,dim)) + end +end +HPCorrector(dim,::RCHP) where {RCHP<:Raycast_HP} = HP_corrector_data(dim,Double64) +HPCorrector(dim,_) = nothing + +CORRECTOR_FLOAT(::RCHP) where {RCHP<:Raycast_HP} = Double64 +CORRECTOR_FLOAT(x) = Float64 +function RaycastIncircleSkip(xs_::HN,dom,parameters::NewRaycastParameter{FLOAT,TREE}) where {P,HN<:HVNodes{P},FLOAT,TREE} xs = copynodes(xs_,parameters.copynodes) lxs=length(xs) + CF = CORRECTOR_FLOAT(parameters.method) dim=size(eltype(xs))[1]#length(xs[1]) z1d_1=zeros(FLOAT,lxs+length(dom)) - z1d_2=zeros(Float64,dim) - z1d_3=zeros(Float64,dim) + z1d_2=zeros(Float64,dim) #CF + z1d_3=zeros(Float64,dim) #CF z1d_4=zeros(FLOAT,dim+1) - z2d_1=zeros(Float64,dim,dim) - z2d_2=zeros(FLOAT,dim,dim) + z2d_1=zeros(Float64,dim,dim) #CF + z2d_2=zeros(Float64,dim,dim) #CF tree = ExtendedTree(xs,dom,parameters.nntree) EI = FastEdgeIterator(zeros(P),1E-8) EI2 = FastEdgeIterator(zeros(P),1E-8) @@ -296,7 +408,7 @@ function RaycastIncircleSkip(xs_::HN,dom,parameters::RaycastParameter{FLOAT,TREE #sizehint!(FEIStorage_global,length(xs)*2^(length(xs[1])-1)) return RaycastIncircleSkip( tree, lxs, length(dom), zeros(Int64,lxs+length(dom)+3), z1d_1, BitVector(zeros(Int8,length(xs))), z2d_1, z2d_2, z1d_2, z1d_3, z1d_4, dom, - zeros(Int64,SRI_max),dim,EI,EI2,xs,General_EdgeIterator(size(eltype(xs))[1]),General_EdgeIterator(size(eltype(xs))[1]),FEIStorage_global,parameters) + zeros(Int64,SRI_max),dim,EI,EI2,xs,General_EdgeIterator(size(eltype(xs))[1]),General_EdgeIterator(size(eltype(xs))[1]),FEIStorage_global,parameters,HPCorrector(dim,parameters.method)) end function copy_RaycastIncircleSkip(original::RaycastIncircleSkip{T, TTT, TTTTT, TTTTTT, FEI, PA, FLOAT}) where {T, TTT, TTTTT, TTTTTT, FEI, PA, FLOAT} diff --git a/src/raycast.jl b/src/raycast.jl index 34e0e80..7f0923f 100644 --- a/src/raycast.jl +++ b/src/raycast.jl @@ -129,6 +129,7 @@ function walkray(full_edge::Sigma, r::Point, xs::Points, searcher, sig, u, edge) end Rest = sig[k] generator, t, r2 = raycast_des(full_edge, r, u, xs, searcher, Rest,edge,sig,Raycast_By_Walkray()) + (generator in sig) && error("hier") exception_raycast(t,r,sig,edge,u,searcher) if t==0.0 return sig, r, false @@ -362,21 +363,59 @@ end ######################################################################################################################################## -function myskips(xs::Points,i::Int64,c::Float64,u::Point,_Cell::Int64,sig::Sigma)#,ts,r) - #b = dot(xs[i], u) <= c -#= if !b - new_t = get_t(r,u,xs[_Cell],xs[i]) - if new_txs[i], sig)/lsig + __x0 = sum(i->xs[i],sig)/lsig + nominator = 2 * u' * (x_new-__x0) + t = 0.0 + for i in sig + x0 = xs[i] + Dx = x_new-x0 + xx = x0+x_new-2*r + t += dot(Dx,xx) / nominator + end + return t/lsig + # r^2 - 2 r x_new + x_new^2 - r^2 + 2 r x0 -x0^2 = 2 r (x0 - x_new) + x_new^2 - x0^2 +end + +function get_t_hp(r,u,x0::P,x_new) where P + #r = SVector{size(P)[1],Double64}(r_) + #u = SVector{size(P)[1],Double64}(u_) + #x0 = SVector{size(P)[1],Double64}(x0_) + #x_new = SVector{size(P)[1],Double64}(x_new_) + Dx = x_new-x0 + xx = x0+x_new-2*r + + return dot(Dx,xx) / (2 * u' * Dx) +# return Float64(dot(Dx,xx) / (2 * u' * Dx)) + # r^2 - 2 r x_new + x_new^2 - r^2 + 2 r x0 -x0^2 = 2 r (x0 - x_new) + x_new^2 - x0^2 end struct Raycast_By_Descend @@ -384,7 +423,7 @@ end struct Raycast_By_Walkray end -function correct_cast(r,r2,u,edge,generator,origin,searcher,cast_type::Raycast_By_Descend) +@inline function correct_cast(r,r2,u,edge,generator,origin,searcher,cast_type::Raycast_By_Descend) return r2 end @@ -394,6 +433,14 @@ function correct_cast(r,r2,u,edge,generator,origin,searcher,cast_type::Raycast_B return r3 end +@inline function correct_cast_hp(edge,xs,searcher,r,u,i,::Raycast_By_Descend,full_mode) + return r +end + +@inline function correct_cast_hp(edge,xs,searcher,r,u,i,::Raycast_By_Walkray,full_mode) + return full_mode ? vertex_calculation_hp(edge,xs,searcher,r,u,i) : r +end + #=function verify_edge(sig,r,u,edge,searcher,origin,xs) ortho = sum(s->abs(dot(u,xs[s]-xs[edge[1]])),sig) # should be almost zero check = sum(s->max(0.0,dot(u,xs[s]-xs[edge[1]])),origin) # should be zero @@ -467,31 +514,46 @@ function raycast_des2(sig::Sigma, old_r, u, xs, searcher::RaycastIncircleSkip, o end -function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,::Raycast_Non_General) +function get_scale(sig, xs,r) + lsig = length(sig) + mid = sum(i->xs[i],sig)/lsig + dist = sum(i->dot(xs[i]-mid,xs[i]-mid),sig)/lsig + return sqrt(dist/dot(xs[sig[1]]-r,xs[sig[1]]-r)) +end + +function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,::RCType,debug=false) where {RCType<:Union{Raycast_Non_General,Raycast_Non_General_Skip}} + shortversion = RCType == Raycast_Non_General_Skip max_int = typemax(Int64) x0 = xs[edge[1]] c1 = maximum(dot(xs[g], u) for g in sig) c2 = abs(c1) c = c1 + c2*searcher.plane_tolerance skip = _i->myskips(xs,_i,c,u,edge[1],sig) - vvv = r + u * dot(u , (x0-r)) + offset = shortversion ? sqrt(norm(r-x0)) : 0.0 + vvv = r + u * (offset +dot(u , (x0-r))) + #vvv = r + u * dot(u , (x0-r)) i, t = _nn(searcher.tree, vvv, skip) t == Inf && return 0, Inf, r x = xs[i] t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) - vvv = r+t*u - i, t = _nn(searcher.tree, vvv, skip) + _vvv = r+t*u + t = get_t(_vvv,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + vvv = _vvv+t*u + i, t = shortversion ? (0, Inf64) : _nn(searcher.tree, vvv, skip) if i!=0 x = xs[i] t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + #__r = r+t*u + #t = t2 + get_t(__r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) end _r = r+t*u measure = maximum(norm(xs[s]-_r) for s in sig) old_measure = measure upper_t = t+2*measure - idss = _inrange(searcher.tree,_r,(1+searcher.b_nodes_tol*100)*measure) + scale = get_scale(sig,xs,_r) + idss = _inrange(searcher.tree,_r,(1+max(1E-12,searcher.b_nodes_tol*100*scale))*measure) lidss = length(idss) # println("1: ",map(k->kupper_t ts[k] = 0.0 + #idss[k] = max_int elseif idss[k]max_dist @@ -544,29 +607,391 @@ function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old , end append!(sig,filter!(i->it +end -function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,::Raycast_Original) +function get__r(vvv,edge,searcher,skip,t,xs,u,r,x0,full_mode,full_error,du,debug,first_t) + _i0 = edge[1] + buffer = MVector{1,Float64}(first_t) + running = true + _modified = false + _r = vvv + ii = 1 + t2__buffer = first_t + while running + ii += 1 + ii>500 && (debug=true) + ii>505 && error("") + # print(ii) + i, _ = _nn(searcher.tree, vvv, i_->skip(i_))# || ts_skip(i_,buffer,xs,edge,x0,r,u,du)) + debug && print("d$i - ") + i==_i0 && break + if i!=0 + _modified = true + x = xs[i] + t2__= get_t_hp(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + t2__buffer = t2__ + t2__ += get_t_hp(r+t2__*u,u,x0,x) + if t2__>=t #-full_error + running = false + break + end + t = t2__ + #__r = r+t*u + #t = t2 + get_t(__r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + else + break + end + __r = r+t*u + #full_mode && prepare_vertex_calculation(edge,xs,searcher) + if full_mode && _modified + _r = vertex_calculation_hp(edge,xs,searcher.hp_vars,__r,u,i) + end + #running = false + end + return _r,t2__buffer +end + +function get__r_2(vvv,edge,searcher,skip,t,xs,u,r,x0,full_mode,_,_,debug,first_t) + i, _ = _nn(searcher.tree, vvv, skip) + _modified = false + if i!=0 + _modified = true + debug && print("h ") + x = xs[i] + t = get_t_hp(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + #t += get_t_hp(r+t*u,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + #__r = r+t*u + #t = t2 + get_t(__r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + end + __r = r+t*u + #full_mode && prepare_vertex_calculation(edge,xs,searcher) + _r = full_mode && _modified ? vertex_calculation_hp(edge,xs,searcher.hp_vars,__r,u,i) : vvv + return _r,t +end + +function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,::Raycast_Non_General_HP,debug=false,du=1E-14) + max_int = typemax(Int64) x0 = xs[edge[1]] + + full_mode = typeof(cast_type)==Raycast_By_Walkray + c1 = maximum(dot(xs[g], u) for g in sig) - c2 = abs(c1) + c2 = abs(c1) c = c1 + c2*searcher.plane_tolerance skip = _i->myskips(xs,_i,c,u,edge[1],sig) vvv = r + u * dot(u , (x0-r)) i, t = _nn(searcher.tree, vvv, skip) + #t==inf && error("") t == Inf && return 0, Inf, r - i2 = i + x = xs[i] + t,full_error = get_t_hp_(r,u,x0,x,du) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + #__t2 = t + first_t = t + _vvv = r+t*u + scale = get_scale(sig,xs,_vvv) + relative_error = full_error / norm(_vvv) + + full_mode &= (relative_error>1E-10 || full_error>1E-8/max(scale,1E-4)) + t += get_t_hp_(_vvv,u,x0,x,du)[1] + _vvv_tu = _vvv+t*u + full_mode && prepare_vertex_calculation(sig,xs,searcher.hp_vars) + vvv = full_mode ? vertex_calculation_hp(edge,xs,searcher.hp_vars,_vvv_tu,u,i) : _vvv_tu + if full_mode + #print("!") + #=if (norm(_vvv_tu-vvv)/(norm(_vvv_tu)+norm(_vvv_tu))>1E-6) + print("!") + elseif (norm(_vvv_tu-vvv)/(norm(_vvv_tu)+norm(_vvv_tu))>1E-8) + print("-") + elseif (norm(_vvv_tu-vvv)/(norm(_vvv_tu)+norm(_vvv_tu))>1E-10) + #print("#$(abs(dot(u,x0-x)/__t2)) - ") + u2 = normalize(vvv-r) + print("#$(__t2), $(norm(r)), $error1 - ") + #print("#$(__t2), $(norm(r)), $(abs((dot(u,x0-x))/__t2)), $error1 - ") + elseif (norm(_vvv_tu-vvv)/(norm(_vvv_tu)+norm(_vvv_tu))<1E-10) && (error1>1E-8)# || __t2>1000) #(norm(vvv)>200 || __t2>200) + print("*") + # print("#$(__t2) - ") + #print("*$(abs(dot(u,x0-x)/__t2)) - ") + end=# + end + #vvv_2 = vertex_calculation_hp(edge,xs,searcher,vvv,u,i) + debug && println(vvv,t) + _r,t = get__r(vvv,edge,searcher,skip,t,xs,u,r,x0,full_mode,full_error,du,debug,first_t) + #_r,t = get__r_2(vvv,edge,searcher,skip,t,xs,u,r,x0,full_mode,full_error,du,debug,first_t) + debug && println(_r,t) + # #= + + # =# + #= + =# + + #=debug && println(full_mode) + if full_mode && _modified && (norm(_r-__r)/(norm(_r)+norm(__r))>1E-6) + println() + println(norm(_r-__r)/(norm(_r)+norm(__r)),", ",norm(_r),", ",norm(__r),", ",dot(_r-__r,u),", ",norm(_r-__r)) + for j in vcat(edge,[i]) + print("$(norm(__r-xs[j])), ") + end + println() + for j in vcat(edge,[i]) + print("$(norm(_r-xs[j])), ") + end + println() + for j in edge + print("$(norm(xs[edge[1]]-xs[j])), ") + end + print("$(norm(xs[edge[1]]-xs[i])), ") + println() + println(_inrange(searcher.tree,__r,norm(__r-xs[edge[1]])*1.00000000001)) + println(_inrange(searcher.tree,_r,norm(_r-xs[edge[1]])*1.00000000001)) + println(dot(u,__r-_r)) + println(typeof(searcher.rhs_cg)) + println(edge) + u2 = normalize(_r-r) + u1 = normalize(__r-r) + println(norm(u2-u1)) + error(norm(_r-__r),", $i, $(norm(vvv-__r)), vv1=$(vertex_variance(vcat(edge,[i]),__r,xs)), vv2=$(vertex_variance(vcat(edge,[i]),_r,xs))") + end=# + measure = maximum(norm(xs[s]-_r) for s in edge) + old_measure = measure + #debug && println(vvv-_r) + upper_t = t*1.0000000001 #+2*measure + idss = _inrange(searcher.tree,_r,(1+max(1E-12,searcher.b_nodes_tol*100*scale))*measure) + lidss = length(idss) + debug && println(idss) + #lidss==0 && error("") +# println("1: ",map(k->k ii∈origin ? 0.0 : get_t_hp(r,u,x0,xs[ii]) ,ts, idss) + if lidss>100 + #print("*") + #=println("*$scale, $measure, $(norm(_r)), $((1+max(1E-12,searcher.b_nodes_tol*scale))*measure)") + for s in edge + print("$s: $(norm(xs[s]-_r)), ") + end + println() + for i in 1:lidss + if ts[i]<=t + print("$i: $(ts[i]) - ") + end + end + #println(idss) + #println(ts) + #println(map(s->norm(xs[s]-_r),idss)) + println(t) + error("")=# + end + k=0 + while kupper_t + idss[k]=max_int + ts[k] = 0.0 + elseif ts[k]upper_t + ts[k] = 0.0 + #idss[k] = max_int + elseif idss[k]max_dist + max_dist = ts[k] + generator = idss[k] + end + end + end + #println("3: ",map(k->k -> -> -> -> ",dot(u,xs[generator]-xs[edge[1]])) + t = get_t_hp(r,u,x0,xs[generator])# (sum(abs2, r - xs[generator]) - sum(abs2, r - x0)) / (2 * u' * (xs[generator]-x0)) + #println(generator,", ",r2) + r2 = correct_cast_hp(edge,xs,searcher.hp_vars,r+t*u,u,generator,cast_type,full_mode) + #=_r = r+t*u + __r = r2 + if full_mode && (norm(_r-__r)/(norm(_r)+norm(__r))>1E-6) + println() + for j in vcat(edge,[generator]) + print("$(norm(__r-xs[j])), ") + end + println() + for j in vcat(edge,[generator]) + print("$(norm(_r-xs[j])), ") + end + println() + for j in edge + print("$(norm(xs[edge[1]]-xs[j])), ") + end + print("$(norm(xs[edge[1]]-xs[generator])), ") + println() + println(_inrange(searcher.tree,__r,norm(__r-xs[edge[1]])*1.00000000001)) + println(_inrange(searcher.tree,_r,norm(_r-xs[edge[1]])*1.00000000001)) + println(dot(u,__r-_r)) + println(typeof(searcher.rhs_cg)) + println(edge) + error(norm(_r-__r),", $i, $(norm(vvv-__r)), vv1=$(vertex_variance(vcat(edge,[i]),__r,xs)), vv2=$(vertex_variance(vcat(edge,[i]),_r,xs))") + end=# + #r2 = correct_cast(r,r+t*u,u,edge,generator,origin,searcher,cast_type) + #if typeof(cast_type)==Raycast_By_Walkray + measure2 = maximum(norm(xs[s]-r2) for s in sig) + measure2 = max(measure2, norm(xs[generator]-r2)) * (1+scale*searcher.b_nodes_tol+1E-16) + k=0 + while kmeasure2 + idss[k] = max_int + elseif idss[k]i1E-6 && error(test,", ",dot(r2-r,u),", ",norm(r2-r)) + return generator, t, r2 +end + + + +function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,::Union{Raycast_Original,Raycast_Original_HP},maxiter=2^length(r)) + x0 = xs[edge[1]] + c1 = maximum(dot(xs[g], u) for g in origin) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + skip = _i->myskips(xs,_i,c,u,edge[1],origin) + offset = sqrt(norm(r-x0)) #10.0#norm(r)>100.0 ? 25.0 : 5.0 + vvv = r + u * (offset +dot(u , (x0-r))) + i, t = _nn(searcher.tree, vvv, skip) + #println(i," ",dot(u,xs[i]),", ",c) + t == Inf && return 0, Inf, r + i2 = i + #if i2 in origin + # HVNearestNeighbors.global_i2[1] = i2 + # i, t = _nn(searcher.tree, vvv, skip) + #end + #i2==195 && print("!") + sk = skip(i2) + #HVNearestNeighbors.global_i2[1] = 0 + (i2 in origin) && error("????????? $i2, $(dot(xs[i2],u)) $c $(sk) $(typeof(searcher.tree.tree))") + count = 0 while i!=0 + count==maxiter && (break) + count += 1 x = xs[i] t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) vvv = r+t*u i, t = _nn(searcher.tree, vvv) - (i==i2 || (i in sig)) && (i=0) - i!=0 && (i2 = i) + if (i==i2 || (i in origin)) + i=0 + end + if i!=0 + i2 = i + end end x2 = xs[i2] t = get_t(r,u,x0,x2) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) @@ -575,21 +1000,119 @@ function raycast_des2(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old , append!(sig,i2) sort!(sig) - return i2, t, _r + return i2, (count==maxiter && i!=0) ? Inf : t, _r end @inline raycast_des(sig, r, u, xs, searcher, old ,edge,origin,cast_type) = raycast_des(sig, r, u, xs, searcher, old ,edge,origin,cast_type,searcher.parameters.method) @inline raycast_des(sig, r, u, xs, searcher, old ,edge,origin,cast_type,method) = raycast_des2(sig, r, u, xs, searcher, old ,edge,origin,cast_type,method) + + + + + ######################################################################################################################################## ######################################################################################################################################## ######################################################################################################################################## -## Raycast several nodes case +## Raycast Hull Search ######################################################################################################################################## ######################################################################################################################################## ######################################################################################################################################## +#= +function raycast_hull(sig::Sigma, r, u, xs, searcher::RaycastIncircleSkip, old ,edge,origin,cast_type,) + max_int = typemax(Int64) + println("---------------------------------------------------------------") + x0 = xs[edge[1]] + c1 = maximum(dot(xs[g], u) for g in sig) + c2 = abs(c1) + c = c1 + c2*searcher.plane_tolerance + skip = _i->myskips(xs,_i,c,u,edge[1],sig) + vvv = r + u * dot(u , (x0-r)) + i, t = _nn(searcher.tree, vvv, skip) + println(i,",",t) + t == Inf && return 0, Inf, r + + x = xs[i] + t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + vvv = r+t*u + i, t = _nn(searcher.tree, vvv, skip) + println(i,",",t) + if i!=0 + x = xs[i] + t = get_t(r,u,x0,x) #(sum(abs2, r - x) - sum(abs2, r - x0)) / (2 * u' * (x-x0)) + end + _r = r+t*u + measure = maximum(norm(xs[s]-_r) for s in sig) + old_measure = measure + + upper_t = t+2*measure + idss = _inrange(searcher.tree,_r,(1+searcher.b_nodes_tol*100)*measure) + lidss = length(idss) +# println("1: ",map(k->k i∈origin ? 0.0 : get_t(r,u,x0,xs[i]) ,ts, idss) + k=0 + while kupper_t + ts[k] = 0.0 + elseif idss[k]max_dist + max_dist = ts[k] + generator = idss[k] + end + end + end + #println("3: ",map(k->k -> -> -> -> ",dot(u,xs[generator]-xs[edge[1]])) + t = get_t(r,u,x0,xs[generator])# (sum(abs2, r - xs[generator]) - sum(abs2, r - x0)) / (2 * u' * (xs[generator]-x0)) + #println(generator,", ",r2) + r2 = correct_cast(r,r+t*u,u,edge,generator,origin,searcher,cast_type) + #if typeof(cast_type)==Raycast_By_Walkray + measure2 = maximum(norm(xs[s]-r2) for s in sig) + measure2 = max(measure2, norm(xs[generator]-r2)) * (1+0.1*searcher.b_nodes_tol) + k=0 + while kmeasure2 + idss[k] = max_int + end + end + append!(sig,filter!(i->ifalse, data) # sortres=false end + +function search_vertex_plane(tree::UnstructuredTree, point::AbstractVector{T}, idx,dist,data) where {T <: Number}#<:Function} + HVNearestNeighbors._knn_plane(tree.tree, point, idx, dist, data) # sortres=false +end diff --git a/src/serialdomain.jl b/src/serialdomain.jl index a2cdea3..22b8760 100644 --- a/src/serialdomain.jl +++ b/src/serialdomain.jl @@ -5,10 +5,10 @@ struct Serial_Domain{P<:Point,M<:AbstractMesh{P}, AM<:AbstractMesh{P}, SI<:HVInt references::SerialVector_Vector{Int64} reference_shifts::SerialVector_Vector{BitVector} internal_boundary::Boundary - _mesh::M - _internal_mesh::AM - _integral::SI - _internal_integral::TT + _mesh::M # public mesh + _internal_mesh::AM # internal mesh + _integral::SI # public integral + _internal_integral::TT # internal integral reflections::SerialVector_Vector{BitVector} _lref::MVector{1,Int64} end @@ -118,6 +118,8 @@ end return (mesh = MeshView(vd.internal_mesh,sv), integral = IntegralView(vd.internal_integral,sv)) end +VDDomain(vd::Serial_Domain) = vd + @inline standardize_mesh(i_mesh::SM) where {SM<:SerialMesh} = MeshView(i_mesh,SortedView(i_mesh.meshes)) @inline standardize_integral(integral::SI,i_mesh) where {SI<:SerialIntegral} = IntegralView(integral,SortedView(i_mesh.meshes)) @inline function standardize(domain::Serial_Domain) @@ -131,6 +133,7 @@ function prepend!(domain::SD,new_xs::ReflectedNodes;kwargs...) where SD<:Serial_ l1 = length(references(domain)) _internal_indeces(domain.mesh,new_xs.references) # transform `new_xs.references` to internal references first domain.internal_mesh = append(domain.internal_mesh,new_xs.data,false) + # recall that append! for SerialVector is different from normal append! append!(domain.references,new_xs.references,domain.internal_mesh.dimensions[end]) append!(domain.reference_shifts,new_xs.reference_shifts,domain.internal_mesh.dimensions[end]) append!(domain.reflections,BitVector[],domain.internal_mesh.dimensions[end]) @@ -144,6 +147,22 @@ function prepend!(domain::SD,new_xs::ReflectedNodes;kwargs...) where SD<:Serial_ #println(sv2 / collect(1:50)) end +function prepend_center!(domain::SD,x) where SD<:Serial_Domain + lnxs = 1 + l1 = length(references(domain)) + domain.internal_mesh = append(domain.internal_mesh,[x],false) + lm = length(domain.internal_mesh) + # recall that append! for SerialVector is different from normal append! + append!(domain.references,[lm],domain.internal_mesh.dimensions[end]) + append!(domain.reference_shifts,[falses(length(boundary(domain)))],domain.internal_mesh.dimensions[end]) + append!(domain.reflections,BitVector[],domain.internal_mesh.dimensions[end]) + i_mesh = domain.internal_mesh + domain.internal_integral = append(domain.internal_integral, i_mesh, false) + domain.lref = domain.lref + lnxs + standardize(domain) + #domain.references.vectors[2].data[1] = lm+1 +end + function append!(domain::SD,new_xs::HV;kwargs...) where {SD<:Serial_Domain, HV<:HVNodes} lnxs = length(new_xs) diff --git a/src/sphere.jl b/src/sphere.jl new file mode 100644 index 0000000..5f1ad02 --- /dev/null +++ b/src/sphere.jl @@ -0,0 +1,165 @@ +abstract type AbstractSphericalDomain{P} <: AbstractDomain{P} end + +struct SphereBoundary{P, M<:AbstractMesh{P}} + mesh::M + center_id::Int64 + center::P + onb::Vector{P} + all_vertices::Vector{P} + radius::Float64 + boundary::Boundary +end +SphereBoundary(a::M,b,c,rad,boundary) where {P,M<:AbstractMesh{P}} = SphereBoundary(a,b,c,[zeros(P) for i in 1:size(P)[1]], [c],rad,boundary) +Base.length(sb::SphereBoundary) = 1 + length(sb.boundary) + +struct SphereReferenceWrapper{T} <: AbstractVector{Int64} + data::T + l::Int64 + SphereReferenceWrapper(data::T) where T = new{T}(data, length(data) - 1) +end +Base.size(d::SphereReferenceWrapper) = (length(d),) +Base.length(d::SphereReferenceWrapper) = d.l +Base.getindex(d::SphereReferenceWrapper, ind::Int64) = d.data[ind] - 1 + + + +struct SimpleSphericalDomain{P,D<:AbstractDomain{P},M,F,TA} <: AbstractSphericalDomain{P} + domain::D + mesh::M + center::P + internal_center_id::Int64 + maps::F + applied_maps::SerialVector_Vector{Int64} + central_sigma::Vector{Int64} + scale::Float64 + radius::Float64 + total_area::TA +end + +function SimpleSphericalDomain(mesh,boundary,center, maps,scale,total_area) + dom = Serial_Domain(mesh,boundary) + ici = length(mesh) + 1 + xs = nodes(mesh) + lxs = length(xs) + app_maps = SerialVector_Vector(Int64[],mesh.dimensions[1]) + prepend_center!(dom,center) + append!(app_maps,Int64[],dom.internal_mesh.dimensions[end]) + c_sig = collect(1:lxs) + + zero_bool_vec = falses(length(boundary)) + if typeof(maps)!=Nothing + for i in 1:length(maps) + new_xs = [maps[i](xs[k]) for k in 1:lxs] + offset = 1 + (i-1)*lxs + references = collect((offset+1):(lxs+offset)) + ref_shifts = [zero_bool_vec for k in 1:lxs] + prepend!(dom,ReflectedNodes(new_xs,references,ref_shifts)) + append!(app_maps,[i for k in 1:lxs],dom.internal_mesh.dimensions[end]) + end + end + SimpleSphericalDomain(dom,mesh,center,ici,maps,app_maps,c_sig,scale,norm(xs[1]-center),total_area) +end + +references(d::SimpleSphericalDomain) = references(d.domain) #SphereReferenceWrapper(references(d.domain)) +mesh(d::SimpleSphericalDomain) = mesh(d.domain) +integral(d::SimpleSphericalDomain) = PublicSphereIntegral(d) #integral(d.domain) +boundary(d::SimpleSphericalDomain) = boundary(d.domain) +shifts(d::SimpleSphericalDomain) = shifts(d.domain) +reference_shifts(d::SimpleSphericalDomain) = reference_shifts(d.domain)#view(reference_shifts(d.domain),1:(reference_length(reference_shifts(domain))-1)) +internal_boundary(d::SimpleSphericalDomain) = internal_boundary(d.domain) + + +struct PublicSphericalDomain{P,Integral<:HVIntegral{P},D} <: AbstractDomain{P} + integral::Integral + domain::D +end +function VDDomain(d::SimpleSphericalDomain{P}) where {P} + return d + integral = PublicSphereIntegral(d) + return PublicSphericalDomain(integral,d) +end +@inline mesh(d::PublicSphericalDomain) = mesh(d.integral) +references(d::PublicSphericalDomain) = SphereReferenceWrapper(references(d.domain)) +integral(d::PublicSphericalDomain) = d.integral +boundary(d::PublicSphericalDomain) = boundary(d.domain) +shifts(d::PublicSphericalDomain) = shifts(d.domain) +reference_shifts(d::PublicSphericalDomain) = view(reference_shifts(d.domain),1:(length(reference_shifts(d.domain))-1)) +internal_boundary(d::PublicSphericalDomain) = internal_boundary(d.domain) + +function internal_boundary(d2::SimpleSphericalDomain,myinte) + m2 = mesh(myinte.Integral) + return SphereBoundary(m2,external_index(m2,d2.internal_center_id),d2.center,d2.radius,boundary(d2.domain)) +end + +@inline function integrate_view(dom::SimpleSphericalDomain) + vd = dom.domain + lim = length(vd.internal_mesh) + sv = CombinedView(SwitchView(length(references(vd))+1,lim),SortedView(vd.internal_mesh.meshes)) + boundary_id = internal_index(vd.internal_mesh,lim+1) + return (mesh = MeshView(SphericalMeshView(vd.internal_mesh,dom.central_sigma,dom.center,dom.scale,dom.internal_center_id,boundary_id),sv), integral = IntegralView(SphericalIntegralView(vd.internal_integral,dom.central_sigma,dom.center,dom.scale,dom.internal_center_id,boundary_id),sv)) +end + + +function get_center(xs) + x0 = xs[1] + + # Initialize A as a static matrix and b as a static vector + A = zeros(MMatrix{size(eltype(xs))[1],size(eltype(xs))[1],Float64} ) + b = zeros(MVector{size(eltype(xs))[1],Float64}) + d = size(eltype(xs))[1] + for i in 1:d + xi = xs[i + 1] + A[i, :] = 2.0 .* (xi .- x0) + b[i] = dot(xi, xi) - dot(x0, x0) + end + + # Solve the linear system A p0 = b + # Using a direct solver since d is typically small + p0_vector = A \ b + + # Convert the result back to an SVector + return SVector{d, Float64}(p0_vector...) +end + +function VoronoiSphere(xs::Points,b=Boundary(); total_area=nothing, transformations = nothing, center = get_center(xs), systematic_error=0.0001, improving=(max_iterations=0, tolerance=1.0,), search_settings::NamedTuple=NamedTuple(), integrator=VI_GEOMETRY,integrand=nothing,mc_accurate=(10000,5,20),silence=false,printevents=false,integrate=true,kwargs...) + oldstd = stdout + result = nothing + check_boundary(xs,b) + + vertex_storage = DatabaseVertexStorage() + + try + myintegrator = replace_integrator(IntegratorType(integrator)) + redirect_stdout(silence ? devnull : oldstd) + + !silence && println(Crayon(foreground=:red,underline=true), "Initialize bulk mesh with $(length(xs)) points",Crayon(reset=true)) + + search=RaycastParameter(eltype(eltype(xs));search_settings...) + mmm_start = cast_mesh(vertex_storage,copy(xs)) + d2 = SimpleSphericalDomain(mmm_start,b, center, transformations, systematic_error,total_area) + lxs = length(xs) + + mmm, integral = integrate_view(d2.domain) + voronoi(mmm,searcher=Raycast(nodes(mmm);domain=b,options=search),Iter = 1:lxs, intro="",printsearcher=printevents, silence=silence) + + m2, i2 = integrate_view(d2) + n = nodes(m2) + #println(n) + #println(length(n)) + #println(center) + #error() + + lboundary = 1 + length(b) + integrate_geo(integrate,d2,myintegrator,integrand,mc_accurate,collect(1:public_length(d2)),collect(1:(length(mesh(d2.domain))+lboundary)),silence) + m,i3 = integrate_view(d2.domain) + println(" ",sum(view(i3.volumes,1:40))) + println(sum(x->i3.area[x][end],1:40)) + result = VoronoiGeometry(myintegrator,d2,integrand,search,mc_accurate,nothing)#NoFile()) + PublicSphereIntegral(d2) + catch err + redirect_stdout(oldstd) + rethrow() + end + redirect_stdout(oldstd) + return result +end diff --git a/src/sphericalmeshview.jl b/src/sphericalmeshview.jl new file mode 100644 index 0000000..c810506 --- /dev/null +++ b/src/sphericalmeshview.jl @@ -0,0 +1,275 @@ +########################################################################################################### + +## SphericalMeshView... + +########################################################################################################### +# Assuming HVNodes is defined elsewhere +# struct HVNodes{P} +# ... +# end + +struct StretchNodesView{P, HVN<:HVNodes{P}} <: HVNodes{P} + data::HVN + center::P + scale::Float64 + center_index::Int64 +end +function StretchNodesView(a,b,c) + ci = 0 + for i in 1:length(a) + ci = i + a[i]==b && break + end + return StretchNodesView(a,b,c,ci) +end + +# Define the size of StretchNodesView to match the underlying data +Base.size(xs::StretchNodesView) = size(xs.data) + +# Define the getindex method to allow array-like indexing +Base.getindex(A::StretchNodesView, inds) = inds != A.center_index ? A.center .+ A.scale .* normalize(A.data[inds] - A.center) : A.center + +# (Optional) Define eltype for better type inference and compatibility +Base.eltype(::Type{StretchNodesView{P, HVN}}) where {P, HVN<:HVNodes{P}} = eltype(HVN) + + +# Define the custom iterator type +mutable struct ConeIterator{I, P} + iterator::I + scale::Float64 + radius::Float64 + center::P + center_id::Int64 + boundary::Int64 + mode::Bool + sigma::Vector{Int64} + r::P + stretch::Float64 +end +function ConeIterator(a,b,c,d,e,f) + center = d + rad = 0.0 + count = 0 + for (sig,r) in a + count += 1 + rad += norm(r-center) + end + rad /= count + + ConeIterator(a,b,c,d,e,f,true,Int64[],d,c/rad) +end + +# Implement the iterate function +function Base.iterate(it::ConeIterator, state=1) + rescale(::Nothing, _ , __) = nothing + function rescale(data,scale,center) + (sig,r) = data[1] + return (sig,center+scale*normalize(r-center)), data[2] + end + decide_first(::Nothing,it,st) = begin + it.mode = true + return nothing + end + function decide_first(data,it,st) + resize!(it.sigma,length(data[1][1])) + it.sigma .= data[1][1] + it.r = data[1][2] + return (it.sigma,it.center + it.scale*(it.r-it.center)), st + end + if it.mode + it.mode = false + return decide_first(iterate(it.iterator, state ), it, state) + else + it.mode = true + it.sigma[findfirst(x->x==it.center_id,it.sigma)] = it.boundary + return (it.sigma,it.center + it.stretch*(it.r-it.center)), state + 1 + end +end +#= +function Base.iterate(it::ConeIterator, state=0) + rescale(::Nothing, _ , __) = nothing + function rescale(data,scale,center) + (sig,r) = data[1] + return (sig,center+scale*normalize(r-center)), data[2] + end + if state == 0 + # First call: yield `first` and set state to 1 + return (it.first, 1) + else + # Delegate to the underlying iterator, adjusting the state + return rescale(iterate(it.iterator, state ),it.scale,it.center) + end +end +=# + +# (Optional) Implement the length method if the underlying iterator is iterable +Base.length(it::ConeIterator{T, I}) where {T, I} = 1 + length(it.iterator) + +# (Optional) Implement eltype for better type inference +Base.IteratorEltype(::Type{ConeIterator{T, I}}) where {T, I} = eltype(I) + + + + +struct SphericalMeshView{P<:Point, VDB <: VertexDB{P}, AM<:AbstractMesh{P,VDB}} <: AbstractMesh{P,VDB} + data::AM + sigma::Vector{Int64} + center::P + radius::Float64 + scale::Float64 + center_id::Int64 + boundary_id::Int64 + # Internal constructor + function SphericalMeshView{P, VDB, AM}(data::AM, sigma, center, scale, center_id, boundary_id) where {P<:Point, VDB <: VertexDB{P}, AM<:AbstractMesh{P,VDB}} + xs = nodes(data) + lxs = length(xs) + lsig = length(sigma) + #=for i in 1:lxs + print("$(norm(xs[i]-center)), ") + end + println() + radius = norm(xs[1]-center) + old_rad = radius + println(center) + for i in 1:lxs + i==lsig+1 && continue + for (s,r) in vertices_iterator(data,i) + sig = copy(s) + for (s2,r2) in vertices_iterator(data,i) + if (abs(dot(normalize(r2)-normalize(r),xs[i]))>1E-10) + for j in sig + print("$(norm(xs[j]-xs[i])), ") + end + for j in s2 + print("$(norm(xs[j]-xs[i])), ") + end + error() + end + end + v = xs[i]-center + new_r = dot(r-center,v)/norm(v) + #print("$new_r -- ")#, $r -- ") + print("$(norm(r-center)) -- ")#, $r -- ") + radius = min(radius,new_r) + #break + end + #radius!=old_rad && break + end + println() + =# + new{P, VDB, AM}( + data, # data + sigma, center, norm(xs[1]-center), scale, center_id, boundary_id + ) + end +end +# External constructor +function SphericalMeshView(data::AM, sigma, center, a, b, c) where {P,VDB,AM<:AbstractMesh{P,VDB}} + SphericalMeshView{P, ptype(data), AM}(data, sigma, center, a, b, c) +end +#@inline copy_sig(mesh::LM,sig) where {LM<:LockMesh} = _copy_indeces(mesh,sig,mesh.sigma) + +@inline Base.getproperty(cd::SphericalMeshView, prop::Symbol) = dyncast_get(cd,Val(prop)) +@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:length_sigma}) = :(getfield(cd,:int_data)[1]) +@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:internal_length_sigma}) = :(getfield(cd,:int_data)[2]) +@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:_Cell}) = :(getfield(cd,:int_data)[3]) +@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:boundary_Vertices}) = :(getfield(cd,:data).boundary_Vertices) +@inline @generated dyncast_get(cd::SphericalMeshView, d::Val{S}) where S = :( getfield(cd, S)) + +@inline Base.setproperty!(cd::SphericalMeshView, prop::Symbol, val) = dyncast_set(cd,Val(prop),val) +@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:length_sigma},val) = :(getfield(cd,:int_data)[1]=val) +@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:internal_length_sigma},val) = :(getfield(cd,:int_data)[2]=val) +@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:_Cell},val) = :(getfield(cd,:int_data)[3]=val) +@inline @generated dyncast_set(cd::SphericalMeshView, d::Val{S},val) where S = :( setfield(cd, S,val)) + +@inline Base.length(mv::SphericalMeshView) = length(mv.data) +@inline internal_length(mv::SphericalMeshView) = internal_length(mv.data) + +@inline internal_index(m::MV,index::Int64) where MV<:SphericalMeshView = internal_index(m.data, index) +@inline external_index(m::MV,index::Int64) where MV<:SphericalMeshView = external_index(m.data,index) +@inline external_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = external_index(m.data,inds) +@inline internal_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = internal_index(m.data,inds) + +@inline internal_sig(mesh::MV,sig::AVI,static::StaticTrue) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = sort!(internal_index(mesh,sig)) +@inline function internal_sig(mesh::MV,sig::AVI,static::StaticFalse) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} + sig .= internal_sig(mesh,sig,statictrue) + return sig +end +@inline external_sig(mesh::MV,sig::AVI,static::StaticTrue) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = sort!(external_index(mesh,sig)) + + +################################################################################################################################################## +################################################################################################################################################## +################################################################################################################################################## + +## basically the only methods that require modification: + +@inline nodes(mv::SphericalMeshView)= StretchNodesView(nodes(mv.data), mv.center, mv.scale) + +@inline function vertices_iterator(mv::MV, index::Int64, internal::StaticTrue) where MV<:SphericalMeshView + return ConeIterator(vertices_iterator(mv.data,index,statictrue), mv.scale, mv.radius, mv.center, mv.center_id, mv.boundary_id) +end + + +@inline number_of_vertices(m::MV, index::Int64, internal::StaticTrue) where MV<:SphericalMeshView = 2*number_of_vertices(m.data,index,statictrue) + 1 # account for the one extra vertex + +################################################################################################################################################## +################################################################################################################################################## +################################################################################################################################################## + +@inline all_vertices_iterator(m::MV, index::Int64, internal::StaticTrue) where MV<:SphericalMeshView = all_vertices_iterator(m.data,index,statictrue) + +@inline push!(mesh::MV, p::Pair{Vector{Int64},T},index) where {T<:Point,MV<:SphericalMeshView{T}} = push!(mesh.data,p,index) +@inline push_ref!(mesh::MV, ref,index) where {T<:Point,MV<:SphericalMeshView{T}} = push_ref!(mesh.data,ref,index) +@inline haskey(mesh::MV,sig::AbstractVector{Int64},index::Int) where MV<:SphericalMeshView = haskey(mesh.data,sig,index) +@inline delete_reference(mesh::MV,s,ref) where MV<:SphericalMeshView = delete_reference(mesh.data,s,ref) + +@inline cleanupfilter!(mesh::MV,i) where MV<:SphericalMeshView = cleanupfilter!(mesh.data,i) +@inline mark_delete_vertex!(mesh::MV,sig,i,ii) where MV<:SphericalMeshView = mark_delete_vertex!(mesh.data,sig,i,ii) + +@inline get_vertex(m::MV,i) where MV<:SphericalMeshView = get_vertex(m.data,i) + + + + + +struct SphericalIntegralView{P<:Point, HVI<:HVIntegral{P}, AM<:AbstractMesh{P}} <: HVIntegral{P} + data::HVI + mesh::AM +end +function SphericalIntegralView(d::HVI, sigma, center, a, b, c) where {P<:Point, HVI<:HVIntegral{P}} + mv = SphericalMeshView(mesh(d),sigma,center, a ,b ,c) + #new{P, HVI, V, typeof(mv)} + SphericalIntegralView( + d, # data + mv + ) +end +#function SphericalIntegralView(d::HVIntegral, v::HV) where HV<:HVView +# SphericalIntegralView{PointType(d), typeof(d), typeof(v)}(d, v) +#end +mesh(iv::SphericalIntegralView) = iv.mesh #MeshView(mesh(iv.data), iv.view) + +@inline Base.getproperty(cd::SphericalIntegralView, prop::Symbol) = dyncast_get(cd,Val(prop)) +@inline @generated dyncast_get(cd::SphericalIntegralView, ::Val{:neighbors}) = :(cd.data.neighbors) +@inline @generated dyncast_get(cd::SphericalIntegralView, ::Val{:volumes}) = :(cd.data.volumes) +@inline @generated dyncast_get(cd::SphericalIntegralView, ::Val{:area}) = :(cd.data.area) +@inline @generated dyncast_get(cd::SphericalIntegralView, ::Val{:interface_integral}) = :(cd.data.interface_integral) +@inline @generated dyncast_get(cd::SphericalIntegralView, ::Val{:bulk_integral}) = :(cd.data.bulk_integral) +@inline @generated dyncast_get(cd::SphericalIntegralView, d::Val{S}) where S = :( getfield(cd, S)) + + +@inline _has_cell_data(I::SphericalIntegralView,_Cell) = _has_cell_data(I.data,_Cell) +@inline cell_data_writable(I::SphericalIntegralView,_Cell,vec,vecvec,::StaticFalse;get_integrals=statictrue) = cell_data_writable(I.data,_Cell,vec,vecvec,staticfalse,get_integrals=get_integrals) +@inline get_neighbors(I::SphericalIntegralView,_Cell,::StaticFalse) = get_neighbors(I.data,_Cell,staticfalse) +@inline set_neighbors(I::SphericalIntegralView,_Cell,new_neighbors,proto_bulk,proto_interface,::StaticFalse) = set_neighbors(I.data,_Cell,new_neighbors,proto_bulk,proto_interface,staticfalse) + +@inline enable(iv::IV;kwargs...) where IV<:SphericalIntegralView = enable(iv.data;kwargs...) +@inline enabled_volumes(Integral::SphericalIntegralView) = enabled_volumes(Integral.data) +@inline enabled_area(Integral::SphericalIntegralView) = enabled_area(Integral.data) +@inline enabled_bulk(Integral::SphericalIntegralView) = enabled_bulk(Integral.data) +@inline enabled_interface(Integral::SphericalIntegralView) = enabled_interface(Integral.data) + +@inline get_area(iv::SphericalIntegralView,c,n,::StaticTrue) = get_area(iv.data,c,n,statictrue) +@inline get_integral(iv::SphericalIntegralView,c,n,::StaticTrue) = get_integral(iv.data,c,n,statictrue) + diff --git a/src/sphericalpublicview.jl b/src/sphericalpublicview.jl new file mode 100644 index 0000000..92e74e9 --- /dev/null +++ b/src/sphericalpublicview.jl @@ -0,0 +1,358 @@ +struct SkipVector{P, AV<:AbstractVector{P}} <: AbstractVector{P} + data::AV + skip::Int64 +end +@inline Base.size(v::SkipVector) = (length(v),) +@inline Base.length(v::SkipVector) = length(v.data) - 1 + +# Implement the `getindex` function +@inline Base.getindex(v::SkipVector, i::Int) = i < v.skip ? v.data[i] : v.data[i + 1] + + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## MultiplyVector + +######################################################################################################################################################################## +######################################################################################################################################################################## + +struct MultiplyVector{T, V<:AbstractVector{T}} <: AbstractVector{T} + data::V + factor::Float64 + + function MultiplyVector(data::V, factor::Float64) where {T, V<:AbstractVector{T}} + new{T,V}(data, factor) + end +end + +# Implement size and length for MultiplyVector +function Base.size(mv::MultiplyVector) + return (length(mv),) +end + +function Base.length(mv::MultiplyVector) + return length(mv.data) +end + +# Define eltype for MultiplyVector +Base.eltype(::Type{MultiplyVector{Float64, V}}) where {V} = Float64 +Base.eltype(::Type{MultiplyVector{V2, V}}) where {T, V2<:AbstractVector{T}, V} = MultiplyVector{T,V2} + +# Implement getindex for MultiplyVector +Base.getindex(v::MultiplyVector{Float64}, i::Int64) = v.factor * v.data[i] +Base.getindex(v::MultiplyVector{T}, i::Int64) where {T} = MultiplyVector(v.data[i], v.factor) + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## TrafoData + +######################################################################################################################################################################## +######################################################################################################################################################################## + + +struct TrafoData + skip1::Int64 + skip2::Int64 + skips::Int64 + + function TrafoData(vec::AbstractVector, taboo1, taboo2) + indices1 = findall(x -> x == taboo1, vec) + indices2 = findall(x -> x == taboo2, vec) + + skip1 = isempty(indices1) ? length(vec) + 1 : minimum(indices1) + skip2 = isempty(indices2) ? length(vec) + 2 : minimum(indices2) + + if skip1 > skip2 + skip1, skip2 = skip2, skip1 + end + + skips = (isempty(indices1) ? 0 : 1) + (isempty(indices2) ? 0 : 1) + + new(skip1, skip2, skips) + end +end + +@inline transform_index(t::TrafoData, ind::Int64) = ind < t.skip1 ? ind : (ind < t.skip2-1 ? ind + 1 : ind + 2) + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## TrafoVector + +######################################################################################################################################################################## +######################################################################################################################################################################## + + + +struct TrafoVector{NV} <: AbstractVector{TrafoData} + data::NV + taboo1::Int64 + taboo2::Int64 + + function TrafoVector(data::AbstractVector{<:AbstractVector{Int64}}, taboo1::Int64, taboo2::Int64) + new{typeof(data)}(data, taboo1, taboo2) + end +end + +# Implement size and length for TrafoVector +function Base.size(tv::TrafoVector) + return (length(tv.data),) +end + +function Base.length(tv::TrafoVector) + return length(tv.data) +end + +# Define eltype for TrafoVector +Base.eltype(::Type{TrafoVector{NV}}) where NV = TrafoData + +# Implement getindex for TrafoVector +@inline Base.getindex(tv::TrafoVector, ind::Int64) = TrafoData(tv.data[ind], tv.taboo1, tv.taboo2) + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## TransformedVector + +######################################################################################################################################################################## +######################################################################################################################################################################## + + +struct TransformedVector{T,AV<:AbstractVector{T}} <: AbstractVector{T} + data::AV + trafo::TrafoData + + function TransformedVector(data::AV2, trafo::TrafoData) where {T,AV2<:AbstractVector{T}} + new{T,AV2}(data, trafo) + end +end + +# Implement size and length for TransformedVector +function Base.size(tv::TransformedVector) + return (length(tv),) +end + +function Base.length(tv::TransformedVector) + return length(tv.data) - tv.trafo.skips +end + +# Define eltype for TransformedVector +Base.eltype(::Type{TransformedVector{AV, T}}) where {AV, T} = T + +# Implement getindex for TransformedVector +function Base.getindex(tv::TransformedVector, ind::Int64) + return tv.data[transform_index(tv.trafo, ind)] +end + + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## TransformedVectorVector + +######################################################################################################################################################################## +######################################################################################################################################################################## + +struct TransformedVectorVector{T, DD<:AbstractVector{T} ,V1<:AbstractVector{DD}, TV} <: AbstractVector{TransformedVector{T, DD}} + data::V1 + trafos::TV + + function TransformedVectorVector(data::V1, trafos::TV) where {T, DD<:AbstractVector{T} ,V1<:AbstractVector{DD}, TV} + new{T,DD,V1,TV}(data, trafos) + end +end + +# Implement size and length for TransformedVectorVector +function Base.size(v::TransformedVectorVector) + return (length(v),) +end + +function Base.length(v::TransformedVectorVector) + return length(v.data) +end + +# Define eltype for TransformedVectorVector +Base.eltype(::Type{TransformedVectorVector{DD, V1, TV}}) where {DD, V1, TV} = TransformedVector{eltype(DD), DD} + +# Implement getindex for TransformedVectorVector +function Base.getindex(v::TransformedVectorVector, ind::Int64) + return TransformedVector(v.data[ind], v.trafos[ind]) +end + + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## PublicVolVector + +######################################################################################################################################################################## +######################################################################################################################################################################## + + +struct PublicVolVector{D<:AbstractVector} <: AbstractVector{Float64} + data::D + skip::Int64 + function PublicVolVector(data::D, skip::Int64) where D<:AbstractVector + new{D}(data, skip) + end +end + +# Implement size and length for PublicVolVector +@inline Base.size(pvv::PublicVolVector) = (length(pvv),) + +@inline Base.length(pvv::PublicVolVector) = length(pvv.data) + +# Define eltype for PublicVolVector +Base.eltype(::Type{PublicVolVector{D}}) where D = Float64 + +# Implement getindex for PublicVolVector +@inline Base.getindex(pvv::PublicVolVector, ind::Int64) = ind <= pvv.skip ? 0.0 : pvv.data[ind][end] + +struct PublicAVVector{AV<:AbstractVector, D<:AbstractVector{AV}} <: AbstractVector{AV} + data::D + skip::Int64 + + function PublicAVVector(data::D, skip::Int64) where {AV<:AbstractVector, D<:AbstractVector{AV}} + new{AV,D}(data, skip) + end +end + +# Implement size and length for PublicAVVector +function Base.size(pavv::PublicAVVector) + return (length(pavv),) +end + +function Base.length(pavv::PublicAVVector) + return length(pavv.data) +end + +# Define eltype for PublicAVVector +Base.eltype(::Type{PublicAVVector{D, AV}}) where {D, AV} = AV + +# Implement getindex for PublicAVVector +function Base.getindex(pavv::PublicAVVector, ind::Int64) + if ind <= pavv.skip + throw(ErrorException("Access to indices <= skip is not allowed.")) + end + return pavv.data[ind] +end + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## PublicSphereMesh + +######################################################################################################################################################################## +######################################################################################################################################################################## + +struct PublicSphereMeshIterator{P,I} + sigma::Vector{Int64} + center_id::Int64 + center::P + radius::Float64 + iterator::I +end + +function Base.iterate(it::PublicSphereMeshIterator,state=1) + decide(::Nothing, _, __) = nothing + function decide(data,it,st) + sig = data[1][1] + r2 = data[1][2] + r = it.center + it.radius * normalize(r2-it.center) + resize(it.sigma,length(sig)-1) + count = 0 + count2 = 0 + for i in 1:length(sig) + count += 1 + sig[count]==it.center_id && continue + count2 +=1 + it.sigma[count2] = sig[count] + end + return (it.sigma,r), st+1 + end + return decide(iterate(it.iterator),it,state) +end + +struct PublicSphereNodes{P,N<:HVNodes{P}} <: HVNodes{P} + data::N + center_id::Int64 + center::P + radius +end +@inline Base.eltype(::PublicSphereNodes{P}) where P = P +@inline Base.getindex(nodes::PublicSphereNodes,ind::Int64) = nodes.center_id==ind ? nodes.center : nodes.center + nodes.radius*normalize(nodes.data[ind]-nodes.center) +@inline Base.size(n::PublicSphereNodes) = size(n.data) + +struct PublicSphereMesh{P,VDB,M<:AbstractMesh{P,VDB}} <: AbstractMesh{P,VDB} + mesh::M + center_id::Int64 + center::P + radius::Float64 +end + +@inline nodes(m::PSM) where {PSM<:PublicSphereMesh} = PublicSphereNodes(nodes(m.mesh),m.center_id,m.center,m.radius) +@inline number_of_vertices(m::PSM, i) where {PSM<:PublicSphereMesh} = number_of_vertices(m.mesh, i) +@inline Base.length(m::PSM) where {PSM<:PublicSphereMesh} = length(m.mesh) +@inline vertices_iterator(m::PS,ind::Int64) where {PS<:PublicSphereMesh} = PublicSphereMeshIterator(Int64[],m.center_id,m.center,m.radius,vertices_iterator(m.mesh,ind)) +@inline external_index(m::M,ind::Int64) where {M<:PublicSphereMesh} = external_index(m.mesh,ind) + +######################################################################################################################################################################## +######################################################################################################################################################################## + +## PublicSphereIntegral + +######################################################################################################################################################################## +######################################################################################################################################################################## + +struct PublicSphereIntegral{P<:Point, HVI<:HVIntegral{P}, AM<:AbstractMesh{P},SVN,SVVol,SVar,SVBI,SVII} <: HVIntegral{P} + data::HVI + mesh::AM + neighbors::SVN + volumes::SVVol + area::SVar + interface_integral::SVII + bulk_integral::SVBI + center_id::Int64 + boundary_id::Int64 +end + +function PublicSphereIntegral(domain) + internal_integral = integral(domain.domain) # the public integral of internal domain + internal_mesh = HighVoronoi.mesh(internal_integral) + + areas = internal_integral.area + + area_factor(::Nothing, _ ) = 1.0 + area_factor(f::Float64, total_ar) = f/total_ar + + lref = length(references(domain.domain)) + total_area = sum(x->areas[x+lref][end],1:(length(internal_mesh)-lref)) + factor = area_factor(domain.total_area, total_area) + dim = 1.0*(length(domain.center)-1) + factor_interface = ( factor^((dim-1)/dim) ) * dim / domain.radius # Ar = Vol * dim / radius + + public_center_id = external_index(internal_mesh,domain.internal_center_id) + boundary_id = length(internal_mesh) + length(internal_boundary(domain)) + 1 + taboo1 = domain.internal_center_id + taboo2 = internal_index(internal_mesh, boundary_id) # official artificial boundary node + + mesh = PublicSphereMesh(internal_mesh,public_center_id,domain.center,domain.radius) + + neighs = internal_integral.neighbors + trafovector = TrafoVector(neighs,taboo1,taboo2) + _neighbors = TransformedVectorVector(neighs,trafovector) + i_ints = internal_integral.interface_integral + tfvv = TransformedVectorVector(i_ints,trafovector) + _interface_integral = MultiplyVector(tfvv,factor_interface) + _area = MultiplyVector(TransformedVectorVector(areas,trafovector),factor_interface) + _bulk_integral = MultiplyVector(i_ints,factor) + _volumes = MultiplyVector(PublicVolVector(areas,public_center_id),factor) + + return PublicSphereIntegral(internal_integral,mesh,_neighbors,_volumes,_area,_interface_integral,_bulk_integral,public_center_id,boundary_id) +end + +mesh(i::PublicSphereIntegral) = i.mesh + + diff --git a/src/statistics.jl b/src/statistics.jl index 3fbebc5..d76a60f 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -41,7 +41,7 @@ function VoronoiStatistics(dim,samples;periodic=nothing,points=1,my_generator=no I,_=voronoi(xs,searcher=Raycast(xs),intro="Run number: $S ",compact=true, Iter=1:number) vp_line_up() I2=Integrator(I.Integral,geodata ? VI_POLYGON : VI_GEOMETRY,integrand=nothing) - _integrate( I2, intro="Run number: $S ", calculate = 1:length(xs), iterate=1:number, compact=true) + _integrate( I2, Boundary(), 1:length(xs), 1:number, "Run number: $S ") mesh = I.Integral.MESH for i in 1:number verteces[(S-1)*number+i] = length(I2.Integral.MESH.All_Vertices[i])+length(I2.Integral.MESH.Buffer_Vertices[i]) diff --git a/src/sysvoronoi.jl b/src/sysvoronoi.jl index 9c0a162..4aef5c9 100644 --- a/src/sysvoronoi.jl +++ b/src/sysvoronoi.jl @@ -338,18 +338,21 @@ function systematic_work_queue(xs,_Cell,mesh,edgecount,searcher,queue,vi_lock,no end end end=# - + function systematic_explore_vertex_multithread(xs::Points,sig,r,_Cell,edgecount,mesh,queue,boundary,searcher,edgeIterator,globallock) k=0 dim = size(eltype(xs))[1] - typeof(edgeIterator)<:FastEdgeIterator && println("$(Threads.threadid())z") + #csig = copy(sig) + #typeof(edgeIterator)<:FastEdgeIterator && println("$(Threads.threadid())z") for (edge,skip) in edgeIterator # print("1") b = pushedge!(edgecount,edge,_Cell) (edge[1]!=_Cell || b ) && continue - full_edge, u = get_full_edge(sig,r,edge,edgeIterator,xs) + #full_edge, u = [0],r + #csig!=sig && error("$sig vs. $csig") + full_edge, u = get_full_edge(sig,r,edge,edgeIterator,xs) sig2, r2, success = walkray(full_edge, r, xs, searcher, sig, u, edge ) # provide missing node "j" of new vertex and its coordinate "r" -# + if sig2 == sig pushray!(mesh,full_edge,r,u,_Cell) continue @@ -439,6 +442,15 @@ function get_full_edge(sig,r,edge,::General_EdgeIterator,xs) u = u_default(sig, xs, i) return Vector{Int64}(edge), u end +#=function get_full_edge_basis(sig,r,edge,::General_EdgeIterator,xs) + i = 0 + while true + i+=1 + !(sig[i] in edge) && break + end + u, base = u_with_base(sig, xs, i) + return Vector{Int64}(edge), u, base +end=# @inline get_full_edge_indexing(sig,r,edge,::General_EdgeIterator,xs) = edge @@ -479,7 +491,10 @@ function fraud_vertex(dim,sig,r,lsig2,searcher,xs) if !verify_vertex(sig,r,xs,searcher) return true end - max_dist = max(map(s->norm(r-xs[s]),sig)) + #println() + #println(sig) + #println(r) + max_dist = maximum(map(s->norm(r-xs[s]),sig)) distance = 0.0 max_distance = 0.0 lsig = length(sig) diff --git a/src/tools.jl b/src/tools.jl index ee6633e..1447f0b 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -781,4 +781,35 @@ function u_qr(sig, xs::HN, i) where {P, HN<:AbstractVector{P}} return eltype(xs)(u) end +function u_qr_onb(onb,x0::P) where {P} + dimension = size(P)[1] + X = MMatrix{dimension, dimension, Float64}(undef) + for j in 1:dimension + X[:, j] .= onb[j] + end + X = SMatrix(X) + Q, R = qr(X) + u = -Q[:,end] * sign(R[end,end]) + return P(u) +end + +function u_with_base(sig, xs::HN, i) where {P, HN<:AbstractVector{P}} + n = length(sig) + dimension = size(P)[1] + X = MMatrix{dimension, dimension, Float64}(undef) + for j in 1:i-1 + X[:, j] = xs[sig[j]] + end + for j in i:n-1 + X[:, j] = xs[sig[j+1]] + end + origin = X[:, end] + X[:, end] = xs[sig[i]] + X .-= origin + X = SMatrix(X) + Q, R = qr(X) + u = -Q[:,end] * sign(R[end,end]) + return eltype(xs)(u), [MVector(Q[:,i]) for i in 1:size(P)[1]] +end + diff --git a/src/voronoi_mesh.jl b/src/voronoi_mesh.jl index 7711a08..0550601 100644 --- a/src/voronoi_mesh.jl +++ b/src/voronoi_mesh.jl @@ -144,9 +144,10 @@ ClassicMesh(xs::Points) = Voronoi_MESH(xs,nothing) @inline set_offset(vdb::VM,i) where {VM<:Voronoi_MESH} = set_offset(vdb.vertices,i) -@inline vertices_iterator(m::VM, i::Int64, ::StaticTrue) where VM<:Voronoi_MESH = vertices_iterator(m.vertices,i,statictrue) +@inline vertices_iterator(m::VM, i::Int64, ::StaticTrue) where {P,VDB,VM<:Voronoi_MESH{P,VDB}} = VertexIterator(m,vertices_iterator(m.vertices,i,statictrue),0) +@inline vertices_iterator(m::VM, i::Int64, ::StaticTrue) where {P,VDB<:VertexDBReference{P}, VM<:Voronoi_MESH{P,VDB}} = vertices_iterator(m.vertices,i,statictrue) #@inline vertices_iterator(m::VM, i::Int64, internal::StaticFalse) where VM<:Voronoi_MESH = VertexIterator(m,vertices_iterator(m.vertices,i,staticfalse)) -@inline all_vertices_iterator(m::VM,i::Int64,static::StaticTrue) where VM<:Voronoi_MESH = all_vertices_iterator(m.vertices,i,static) +@inline all_vertices_iterator(m::VM,i::Int64,static::StaticTrue) where VM<:Voronoi_MESH = all_vertices_iterator(m.vertices,i,static) #VertexIterator(m,all_vertices_iterator(m.vertices,i,static),0) @inline number_of_vertices(m::VM,i::Int64,static::StaticFalse) where VM<:Voronoi_MESH = number_of_vertices(m.vertices,internal_index(m,i),statictrue) @inline number_of_vertices(m::VM,i::Int64,static::StaticTrue) where VM<:Voronoi_MESH = number_of_vertices(m.vertices,i,statictrue) @inline internal_length(m::M) where M<:Voronoi_MESH = m.initial_length diff --git a/src/voronoidata.jl b/src/voronoidata.jl index 8ab9475..865155c 100644 --- a/src/voronoidata.jl +++ b/src/voronoidata.jl @@ -16,6 +16,7 @@ end @inline function transform(bonus::D,neighs,offset,simple=staticfalse) where {D<:DeepNeighborData} lneigh = length(neighs) + #println(typeof(neighs)) vbuffer = simple==false ? _external_indeces(bonus.mesh,neighs,bonus.buffer) : _copy_indeces(bonus.mesh,neighs,bonus.buffer) VoronoiDataArray(vbuffer,offset,bonus.references;lsigma=lneigh,lreferences=offset) return vbuffer @@ -206,10 +207,10 @@ struct Vertices_Vector{M,B} #<: AbstractArray{PublicIterator} offset::Int64 bonus::B function Vertices_Vector(d::D,publicview=false) where {D<:AbstractDomain} - m = mesh(d) ref = references(d) i = integral(d) - b = DeepNeighborData(ref,mesh(i)) + m = mesh(i)#d) + b = DeepNeighborData(ref,m) return new{typeof(m),typeof(b)}(m,publicview*length(ref),b) end end @@ -265,7 +266,7 @@ struct OrientationsVector{P,N<:HVNodes{P},AV,S<:StaticBool} <: AbstractVector{Or onboundary::S lmesh::Int64 function OrientationsVector(d::AD,onboundary,publicview=staticfalse;sorting=nothing) where {P,AD<:AbstractDomain{P}} - m = mesh(d) + m = mesh(integral(d)) n = nodes(m) _nei = DeepVectorNeighbor(d,false) # gets external representation of full neighbor matrix o = publicview==true ? length(references(d)) : 0 @@ -353,10 +354,10 @@ struct BNodesDictDict{P, S, onBoundary<:StaticBool, N,N2,M<:AbstractMesh{P}} <: function BNodesDictDict(d::AD,onboundary=false,publicview=false) where {P,AD<:AbstractDomain{P}} x0 = zeros(P) n = DeepVectorNeighbor(d,false) - _nodes = nodes(mesh(d)) + m = mesh(integral(d)) + _nodes = nodes(m)#mesh(integral(d))) b = boundary(d) buf = MVector{length(b),Bool}(falses(length(b))) - m = mesh(d) lmesh = length(m) myob = StaticBool(onboundary) o = publicview==false ? 0 : length(references(d)) @@ -541,27 +542,32 @@ The call of `VoronoiData(VG)` provides the following options: - `sorted=true`: During the reduction of the internal pseudo periodic mesh to the fully periodic output, the neighbors (jointly with their respective properties) get sorted by their numbers. This is only possible if `getarea`,`getneighbors` and `getinterfaceintegral` are `true`. Otherwise it will be ignored """ function VoronoiData(VG::PGeometry{P};reduce_to_periodic=true,view_only=true,copyall=!view_only,getboundary=copyall,getbulk_integral=copyall,getreferences=copyall,getreference_shifts=copyall,getinterface_integral=copyall,getvolume=copyall,getarea=copyall,getneighbors=copyall,getnodes=copyall,getboundary_vertices=copyall,getorientations=copyall,getvertices=copyall,getboundary_nodes=copyall,onboundary=false,sorted=false) where P - domain = VG.domain + domain = VDDomain(VG.domain) + this_integral = integral(domain) + #println(typeof(this_integral)) + #println(length(this_integral.volumes)) + #println(length(references(domain))) offset = reduce_to_periodic * length(references(domain)) deepversion(::StaticTrue,a) = convert_to_vector(a) deepversion(::StaticFalse,a) = a - ___nodes = nodes(mesh(domain)) + _mesh = mesh(this_integral) + ___nodes = nodes(_mesh) _nodes = deepversion(StaticBool(getnodes),view(___nodes,(offset+1):length(___nodes))) _vertices = deepversion(StaticBool(getvertices),Vertices_Vector(domain,reduce_to_periodic)) - _bv = deepversion(StaticBool(getboundary_vertices),Public_BV_Iterator(mesh(domain))) # referred to by boundary_verteces + _bv = deepversion(StaticBool(getboundary_vertices),Public_BV_Iterator(_mesh)) # referred to by boundary_verteces _bn = deepversion(StaticBool(getboundary_nodes),BNodesDictDict(domain,onboundary,reduce_to_periodic)) #boundary_nodes_on_boundary = onboundary __neighbors = deepversion(StaticBool(getneighbors),DeepVectorNeighbor(domain,reduce_to_periodic)) _neighbors, bonus = VoronoiData_neighbors(StaticBool(sorted),StaticBool(getneighbors),__neighbors) ori = deepversion(StaticBool(getorientations),OrientationsVector(domain,onboundary,StaticBool(reduce_to_periodic),sorting=bonus)) - _volume = deepversion(StaticBool(getvolume),DeepVectorFloat64(integral(domain).volumes,offset)) - _area = deepversion(StaticBool(getarea),DeepVector(integral(domain).area,offset,staticfalse,Val(:deep),bonus)) - _bi = deepversion(StaticBool(getbulk_integral),DeepVectorFloat64Vector(integral(domain).bulk_integral,offset)) - _ii = deepversion(StaticBool(getinterface_integral),DeepVector(integral(domain).interface_integral,offset,staticfalse,Val(:deep),bonus)) - _boun = reduce_to_periodic ? boundary(VG.domain) : internal_boundary(VG.domain) + _volume = deepversion(StaticBool(getvolume),DeepVectorFloat64(this_integral.volumes,offset)) + _area = deepversion(StaticBool(getarea),DeepVector(this_integral.area,offset,staticfalse,Val(:deep),bonus)) + _bi = deepversion(StaticBool(getbulk_integral),DeepVectorFloat64Vector(this_integral.bulk_integral,offset)) + _ii = deepversion(StaticBool(getinterface_integral),DeepVector(this_integral.interface_integral,offset,staticfalse,Val(:deep),bonus)) + _boun = reduce_to_periodic ? boundary(domain) : internal_boundary(domain) boun = getboundary ? deepcopy(_boun) : _boun - _references = deepversion(StaticBool(getreferences),references(VG.domain)) - _reference_shifts = deepversion(StaticBool(getreference_shifts),ShiftVector{P}(reference_shifts(VG.domain),shifts(VG.domain))) + _references = deepversion(StaticBool(getreferences),references(domain)) + _reference_shifts = deepversion(StaticBool(getreference_shifts),ShiftVector{P}(reference_shifts(domain),shifts(domain))) return VoronoiData(_nodes,_vertices,_bv,_bn,onboundary,_neighbors,ori,_volume,_area,_bi,_ii,length(references(domain)),_references,_reference_shifts,VG,boun) end function VoronoiData_neighbors(::StaticTrue,getneighbors::StaticBool,__neighbors) diff --git a/src/voronoidomain.jl b/src/voronoidomain.jl index 9ce5ead..015c633 100644 --- a/src/voronoidomain.jl +++ b/src/voronoidomain.jl @@ -37,6 +37,8 @@ function Voronoi_Domain(mesh,boundary,internaly_precise=true) return Voronoi_Domain(mesh,boundary,shifts,Vector{BitVector}(),Vector{Int64}(),copy(boundary),internaly_precise,bv) end +VDDomain(vd::Voronoi_Domain) = vd + @inline function integrate_view(vd::Voronoi_Domain) sv = SwitchView(length(vd.references)+1,length(mesh(vd))) return (mesh = MeshView(vd._mesh.data,sv), integral = IntegralView(vd._integral,sv)) From 380e95cb8f31cbb51ee7f56dba2c85da6f5d6fa5 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:49:51 +0100 Subject: [PATCH 09/19] Update runtests.jl --- test/runtests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7e0afc4..b3df448 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,4 @@ using Test -using Revise using HighVoronoi using SpecialFunctions From 337c9c1ade84158652c2dc866bf440ef79b19c9e Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 09:54:53 +0100 Subject: [PATCH 10/19] Update HighVoronoi.jl --- src/HighVoronoi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HighVoronoi.jl b/src/HighVoronoi.jl index 42b0a81..b387270 100644 --- a/src/HighVoronoi.jl +++ b/src/HighVoronoi.jl @@ -14,7 +14,7 @@ using GLPK using Plots using Distances using ProgressMeter -using Cthulhu +#using Cthulhu using Base.Threads using Base.Threads: Atomic, atomic_cas! using DoubleFloats From 0c274c23bd6d1a9109efb30eedb829724e9a09c8 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 10:27:02 +0100 Subject: [PATCH 11/19] tree modifecations --- src/NearestNeighborModified/kd_tree.jl | 348 ++++++++++++++++++++++-- src/NearestNeighborModified/tree_ops.jl | 49 +++- 2 files changed, 356 insertions(+), 41 deletions(-) diff --git a/src/NearestNeighborModified/kd_tree.jl b/src/NearestNeighborModified/kd_tree.jl index 5187eea..1860390 100644 --- a/src/NearestNeighborModified/kd_tree.jl +++ b/src/NearestNeighborModified/kd_tree.jl @@ -335,6 +335,153 @@ function knn_kernel_flex!(tree::HVKDTree{V}, return true end +############################################################################################################################### +############################################################################################################################### + +## PEAK + +############################################################################################################################### +############################################################################################################################### + +mutable struct PeakData{T,T2} + c0::Float64 + direction::T + #maxes::T2 + #mins::T2 + vertex::T2 + index::Int64 + function PeakData(a,b,_max,_min) + m1 = MVector(zeros(typeof(b))) + #m2 = MVector(zeros(typeof(b))) + #m3 = MVector(zeros(typeof(b))) + m1 .= _min + l1 = length(m1) + for i in 1:l1 + if b[i]>0 + m1[i] = _max[i] + end + end + #m2 .= d + return new{typeof(b),typeof(m1)}(a,b,m1,0)#,m2,m3,0) + end +end + +function peak_direction(tree::H, + direction, + c0,return_on_find=false) where {A,B,C,H<:HVKDTree{A,B,C}} + d = PeakData(c0,direction,tree.hyper_rec.maxes, tree.hyper_rec.mins) + peak_kernel!(tree, 1, d,return_on_find) + return d.index +end + +function peak_kernel!(tree::HVKDTree{V}, index, data::D,return_on_find) where {V, D} + sd = 0 + v_dim = 0.0 + if index<=length(tree.nodes) + node = tree.nodes[index] + sd = node.split_dim + v_dim = data.vertex[sd] + data.vertex[sd] = data.direction[sd]>0 ? node.hi : node.lo + c_ = dot(data.vertex,data.direction) + data.vertex[sd] = v_dim + c_ < data.c0 && (return false) + else + c_ = dot(data.vertex,data.direction) + c_ < data.c0 && (return false) + end + + if isleaf(tree.tree_data.n_internal_nodes, index) + for z in get_leaf_range(tree.tree_data, index) + @inbounds tiz = tree.indices[z] + idx = tree.reordered ? z : tiz + p = tree.data[idx] + c_ = dot(p,data.direction) + #(tiz==10 || tiz==34) && print("TREFFER: c=$(data.c0), p*u=$(c_), u=$(data.direction), idx=$idx p=$p") + if c_ > data.c0 + #println("+") + data.c0 = c_ + data.index = idx + return_on_find && (return true) + end + #print("-") + end + return false + end + + right = getright(index) # indices left and right + left = getleft(index) + data.vertex[sd] = data.direction[sd]>0 ? node.hi : node.split_val + peak_kernel!(tree,right,data,return_on_find) && return_on_find && (return true) + data.vertex[sd] = data.direction[sd]<0 ? node.lo : node.split_val + peak_kernel!(tree,left,data,return_on_find) && return_on_find && (return true) + + data.vertex[sd] = v_dim + return false +end + + +function search_node_direction(tree::H, + direction, + c0,idx) where {A,B,C,H<:HVKDTree{A,B,C}} + d = PeakData(c0,direction,tree.hyper_rec.maxes, tree.hyper_rec.mins) + search_kernel!(tree, 1, d,idx,1) + return d.index +end + +function search_kernel!(tree::HVKDTree{V}, index, data::D,idx,level) where {V, D} + sd = 0 + v_dim = 0.0 + if index<=length(tree.nodes) + node = tree.nodes[index] + sd = node.split_dim + v_dim = data.vertex[sd] + + data.vertex[sd] = data.direction[sd]>0 ? node.hi : node.lo + c_ = dot(data.vertex,data.direction) + data.vertex[sd] = v_dim + print("$level: $(data.vertex)") + println(" -> $c_") + #c_ < data.c0 && (return false) + else + c_ = dot(data.vertex,data.direction) + print("$level: $(data.vertex)") + println(" -> $c_") + #c_ < data.c0 && (return false) + end + + if isleaf(tree.tree_data.n_internal_nodes, index) + for z in get_leaf_range(tree.tree_data, index) + @inbounds tiz = tree.indices[z] + idx_ = tree.reordered ? z : tiz + p = tree.data[idx_] + if tiz==idx + println("$p") + return true + end + end + return false + end + + right = getright(index) # indices left and right + left = getleft(index) + data.vertex[sd] = data.direction[sd]>0 ? node.hi : node.split_val + search_kernel!(tree,right,data,idx,level+1) && (return true) + data.vertex[sd] = data.direction[sd]<0 ? node.lo : node.split_val + search_kernel!(tree,left,data,idx,level+1) && (return true) + + data.vertex[sd] = v_dim + return false +end + +############################################################################################################################### +############################################################################################################################### + +## KNN + +############################################################################################################################### +############################################################################################################################### + + function _knn(tree::HVKDTree, point::AbstractVector, best_idxs::AbstractVector{Int}, @@ -354,42 +501,185 @@ function knn_kernel!(tree::HVKDTree{V}, best_dists::AbstractVector, min_dist, skip::F) where {V, F} -# At a leaf node. Go through all points in node and add those in range -if isleaf(tree.tree_data.n_internal_nodes, index) -add_points_knn!(best_dists, best_idxs, tree, index, point, false, skip) -return + # At a leaf node. Go through all points in node and add those in range + if isleaf(tree.tree_data.n_internal_nodes, index) + add_points_knn!(best_dists, best_idxs, tree, index, point, false, skip) + return + end + + node = tree.nodes[index] + p_dim = point[node.split_dim] + split_val = node.split_val + lo = node.lo + hi = node.hi + split_diff = p_dim - split_val + M = tree.metric + # Point is to the right of the split value + if split_diff > 0 + close = getright(index) + far = getleft(index) + ddiff = max(zero(eltype(V)), p_dim - hi) + else + close = getleft(index) + far = getright(index) + ddiff = max(zero(eltype(V)), lo - p_dim) + end + # Always call closer sub tree + knn_kernel!(tree, close, point, best_idxs, best_dists, min_dist, skip) + + split_diff_pow = eval_pow(M, split_diff) + ddiff_pow = eval_pow(M, ddiff) + diff_tot = eval_diff(M, split_diff_pow, ddiff_pow) + new_min = eval_reduce(M, min_dist, diff_tot) + if new_min < best_dists[1] + knn_kernel!(tree, far, point, best_idxs, best_dists, new_min, skip) + end + return end -node = tree.nodes[index] -p_dim = point[node.split_dim] -split_val = node.split_val -lo = node.lo -hi = node.hi -split_diff = p_dim - split_val -M = tree.metric -# Point is to the right of the split value -if split_diff > 0 -close = getright(index) -far = getleft(index) -ddiff = max(zero(eltype(V)), p_dim - hi) -else -close = getleft(index) -far = getright(index) -ddiff = max(zero(eltype(V)), lo - p_dim) + +############################################################################################################################### +############################################################################################################################### + +## knn_plane + +############################################################################################################################### +############################################################################################################################### + +#=function _knn_plane(tree::HVKDTree, + point::AbstractVector, + best_idxs::AbstractVector{Int}, + best_dists::AbstractVector, + data) where {F} +init_min = get_min_distance(tree.hyper_rec, point) +knn_kernel_plane!(tree, 1, point, best_idxs, best_dists, init_min, data) +@simd for i in eachindex(best_dists) +@inbounds best_dists[i] = eval_end(tree.metric, best_dists[i]) +end end -# Always call closer sub tree -knn_kernel!(tree, close, point, best_idxs, best_dists, min_dist, skip) -split_diff_pow = eval_pow(M, split_diff) -ddiff_pow = eval_pow(M, ddiff) -diff_tot = eval_diff(M, split_diff_pow, ddiff_pow) -new_min = eval_reduce(M, min_dist, diff_tot) -if new_min < best_dists[1] -knn_kernel!(tree, far, point, best_idxs, best_dists, new_min, skip) +function knn_kernel_plane!(tree::HVKDTree{V}, + index::Int, + point::AbstractVector, + best_idxs::AbstractVector{Int}, + best_dists::AbstractVector, + min_dist, + full_data) where {V} + # At a leaf node. Go through all points in node and add those in range + sd = 0 + v_dim = 0.0 + HighVoronoi.blocked(full_data,index) && (return true) + if index<=length(tree.nodes) + node = tree.nodes[index] + sd = node.split_dim + v_dim = full_data.vertex[sd] + full_data.vertex[sd] = full_data.direction[sd]>0 ? node.hi : node.lo + c_ = dot(full_data.vertex,full_data.direction) + full_data.vertex[sd] = v_dim + c_ < full_data.c0 && (return true) + else + c_ = dot(full_data.vertex,full_data.direction) + c_ < full_data.c0 && (return true) + end + + if isleaf(tree.tree_data.n_internal_nodes, index) + return add_points_knn_plane!(best_dists, best_idxs, tree, index, point, false, full_data) + end + + node = tree.nodes[index] + p_dim = point[node.split_dim] + split_val = node.split_val + lo = node.lo + hi = node.hi + split_diff = p_dim - split_val + M = tree.metric + # Point is to the right of the split value + if split_diff > 0 + close = getright(index) + far = getleft(index) + ddiff = max(zero(eltype(V)), p_dim - hi) + else + close = getleft(index) + far = getright(index) + ddiff = max(zero(eltype(V)), lo - p_dim) + end + # Always call closer sub tree + plane1 = plane2 = knn_kernel_plane!(tree, close, point, best_idxs, best_dists, min_dist, full_data) + plane1 && HighVoronoi.block(full_data,close) + + split_diff_pow = eval_pow(M, split_diff) + ddiff_pow = eval_pow(M, ddiff) + diff_tot = eval_diff(M, split_diff_pow, ddiff_pow) + new_min = eval_reduce(M, min_dist, diff_tot) + if new_min < best_dists[1] + plane2 = knn_kernel_plane!(tree, far, point, best_idxs, best_dists, new_min, full_data) + plane2 && HighVoronoi.block(full_data,far) + else + plane2 = HighVoronoi.blocked(full_data,far) + end + return plane1 && plane2 end -return + +=# + +############################################################################################################################### +############################################################################################################################### + +## REDUCTION + +############################################################################################################################### +############################################################################################################################### + +reduction!(tree::HVKDT) where {HVKDT<:HVKDTree} = reduction_kernel!(tree, 1) + +function reduction_kernel!(tree::HVKDTree{V,M,T}, + index::Int) where {V,M,T} + # At a leaf node. Go through all points in node and add those in range + low = MVector(Inf*ones(V)) + hi = MVector(-Inf*ones(V)) + if isleaf(tree.tree_data.n_internal_nodes, index) + for z in get_leaf_range(tree.tree_data, index) + @inbounds tiz = tree.indices[z] + idx = tree.reordered ? z : tiz + p = tree.data[idx] + for i in 1:size(V)[1] + val = p[i] + if val>hi[i] + hi[i] = val + end + if val < low[i] + low[i] = val + end + end + end + else + close = getright(index) + far = getleft(index) + low1,hi1 = reduction_kernel!(tree, close) + low2,hi2 = reduction_kernel!(tree, far) + for i in 1:size(V)[1] + low[i] = min(low1[i],low2[i]) + hi[i] = max(hi1[i],hi2[i]) + end + end + if index<=length(tree.nodes) + old_node = tree.nodes[index] + sd = old_node.split_dim + new_node = KDNode{T}(low[sd],hi[sd],old_node.split_val,sd) + tree.nodes[index] = new_node + end + return low,hi end + +############################################################################################################################### +############################################################################################################################### + +## INRANGE + +############################################################################################################################### +############################################################################################################################### + function _inrange(tree::HVKDTree, point::AbstractVector, radius::Number, diff --git a/src/NearestNeighborModified/tree_ops.jl b/src/NearestNeighborModified/tree_ops.jl index 34555c3..9757b2e 100644 --- a/src/NearestNeighborModified/tree_ops.jl +++ b/src/NearestNeighborModified/tree_ops.jl @@ -127,26 +127,51 @@ end return #!result # return true iff new_r has not changed end +const global_i2 = MVector{1,Int64}([0]) + # Checks the distance function and add those points that are among the k best. # Uses a heap for fast insertion. @inline function add_points_knn!(best_dists::AbstractVector, best_idxs::AbstractVector{Int}, tree::HVNNTree, index::Int, point::AbstractVector, do_end::Bool, skip::F) where {F} -for z in get_leaf_range(tree.tree_data, index) -idx = tree.reordered ? z : tree.indices[z] -dist_d = myevaluate(tree.metric, tree.data[idx], point, do_end) -if dist_d <= best_dists[1] -if skip(tree.indices[z]) -continue + for z in get_leaf_range(tree.tree_data, index) + idx = tree.reordered ? z : tree.indices[z] + dist_d = myevaluate(tree.metric, tree.data[idx], point, do_end) + if dist_d <= best_dists[1] + tiz = tree.indices[z] + #if tiz==global_i2[1] + # print(tiz,",",skip(tiz),": ") + #end + if skip(tree.indices[z]) + continue + end + + best_dists[1] = dist_d + best_idxs[1] = idx + percolate_down!(best_dists, best_idxs, dist_d, idx) + end + end end -best_dists[1] = dist_d -best_idxs[1] = idx -percolate_down!(best_dists, best_idxs, dist_d, idx) -end -end +#=@inline function add_points_knn_plane!(best_dists::AbstractVector, best_idxs::AbstractVector{Int}, + tree::HVNNTree, index::Int, point::AbstractVector, + do_end::Bool, data) + final_skip = true + for z in get_leaf_range(tree.tree_data, index) + idx = tree.reordered ? z : tree.indices[z] + dist_d = myevaluate(tree.metric, tree.data[idx], point, do_end) + tiz = tree.indices[z] + skip = dot(data.direction,data.xs[tiz]) < data.c0 + final_skip &= skip + if dist_d <= best_dists[1] && !skip + best_dists[1] = dist_d + best_idxs[1] = idx + percolate_down!(best_dists, best_idxs, dist_d, idx) + end + end + return final_skip end - +=# # Add those points in the leaf node that are within range. From 1ed2fb270cf4d2182c887938352333a9426773dd Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 11:04:51 +0100 Subject: [PATCH 12/19] Update make.jl --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 5937cd2..9dcc165 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,7 +19,7 @@ makedocs(; ), pages=[ "Home" => "index.md", - "Convex Hulls" => "convexhull.md", + "Convex Hulls" => "man/convexhull.md", "Manual Voronoi" => Any[ "Examples Voronoi Generation" => "man/short.md", "man/workflowmesh.md", From de60d50368da40dd0a61e1f1e7e44ba5cdd0e581 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 11:14:12 +0100 Subject: [PATCH 13/19] Update CI.yml --- .github/workflows/CI.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0abc879..917361c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -14,11 +14,14 @@ jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} + # Für nightly-Tests Fehler nur als Warnung behandeln + continue-on-error: ${{ matrix.version == 'nightly' }} strategy: fail-fast: false matrix: version: - '1.7' + - '1.9' - 'nightly' os: - ubuntu-latest From 0b6707ccefeb58dfda624a974ea71a53e323be7d Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 11:35:15 +0100 Subject: [PATCH 14/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 851a900..6317e48 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ Crayons = "4.1" GLPK = "1.1" IterTools = "1.4" IterativeSolvers = "0.9.2" -JLD2 = "0.4.31" +JLD2 = "0.4.31-0.5.11" NearestNeighbors = "0.4.13" Polyhedra = "0.7.6" SpecialFunctions = "1,2" From 0f9f9545d33e7b5e36d843bb2d1d69c878810a7b Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 11:41:31 +0100 Subject: [PATCH 15/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6317e48..3df748a 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ Crayons = "4.1" GLPK = "1.1" IterTools = "1.4" IterativeSolvers = "0.9.2" -JLD2 = "0.4.31-0.5.11" +JLD2 = ">=0.4.31, <0.5.12" NearestNeighbors = "0.4.13" Polyhedra = "0.7.6" SpecialFunctions = "1,2" From d04d51a8db24dd4dd9eb8dd007a352ed8ef7cdd3 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 13:32:51 +0100 Subject: [PATCH 16/19] remove unused code --- src/NearestNeighborModified/kd_tree.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NearestNeighborModified/kd_tree.jl b/src/NearestNeighborModified/kd_tree.jl index 1860390..cc78bc0 100644 --- a/src/NearestNeighborModified/kd_tree.jl +++ b/src/NearestNeighborModified/kd_tree.jl @@ -420,7 +420,7 @@ function peak_kernel!(tree::HVKDTree{V}, index, data::D,return_on_find) where {V end -function search_node_direction(tree::H, +#=function search_node_direction(tree::H, direction, c0,idx) where {A,B,C,H<:HVKDTree{A,B,C}} d = PeakData(c0,direction,tree.hyper_rec.maxes, tree.hyper_rec.mins) @@ -471,7 +471,7 @@ function search_kernel!(tree::HVKDTree{V}, index, data::D,idx,level) where {V, D data.vertex[sd] = v_dim return false -end +end=# ############################################################################################################################### ############################################################################################################################### From 075da829ddab303e71fa9f75e6fb515323051b64 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 13:33:56 +0100 Subject: [PATCH 17/19] remove unused code --- src/chull.jl | 20 ++++++++++++-------- src/extended.jl | 4 ++-- src/sphere.jl | 5 +++-- src/sphericalmeshview.jl | 38 +++++++++++++++++--------------------- src/sphericalpublicview.jl | 8 ++++---- 5 files changed, 38 insertions(+), 37 deletions(-) diff --git a/src/chull.jl b/src/chull.jl index bc9f875..4299bff 100644 --- a/src/chull.jl +++ b/src/chull.jl @@ -44,7 +44,9 @@ function Base.getindex(pv::PV, i::Int) where {T<:Int,PV<:PrependedVector{T}} return @inbounds pv.rest[i - 1] + 1 else throw(BoundsError(pv, i)) + return pv.first end + return pv.first end function Base.setindex!(pv::PV, value, i::Int) where {T<:Int,PV<:PrependedVector{T}} @@ -71,6 +73,7 @@ function Base.iterate(pv::PrependedVector, state) end # Optional: print für eine anschauliche Ausgabe +#= Base.show(io::IO, pv::PV) where {PV<:PrependedVector} = begin print(io, "[$(pv.first)") for x in pv.rest @@ -78,7 +81,7 @@ Base.show(io::IO, pv::PV) where {PV<:PrependedVector} = begin end print(io, "]") end - +=# ################################################################################################################## ################################################################################################################## @@ -305,7 +308,6 @@ function __chull(db,indices,queue_position,prog,searcher,edgecount,node_threads, start_time = time_ns() atomic_xchg!(PRINT_STATUS,false) this_sig = [0] - ps_data = PlaneSearchData(P,length(xs.rest),xs.rest) try while !global_end[] status(1) @@ -363,7 +365,7 @@ function __chull(db,indices,queue_position,prog,searcher,edgecount,node_threads, status(5) #@descend explore_chull_vertex(ei,facet_data,xs,onb,sig,r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) #error("") - nv = explore_chull_vertex(ei,facet_data,xs,onb,sig,r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) + nv = explore_chull_vertex(ei,facet_data,xs,onb,sig,r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull)#,ps_data) end status_2(6) #facet_data[1]==[18, 174, 243, 333, 377, 418] && error("KLAPPT!") @@ -384,6 +386,7 @@ end #println("ending: $(Threads.threadid())") end +#= function verify(sig,r,u,xs) ret = true r += u * dot(u,xs[sig[2]]-r) @@ -424,6 +427,7 @@ function verify(sig,r,u,xs) end return ret end +=# @inline function queue_convex_edges_on_find(raw_data,searcher,xs::PV_P,edgecount) where {P,PV_P<:PrependedVector{P}} sig,r = preparate_convex_data!(xs,raw_data) @@ -540,7 +544,7 @@ end end=# -function explore_chull_vertex(edgeIterator,facet_data,xs,onb,sig,original_r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull,ps_data) +function explore_chull_vertex(edgeIterator,facet_data,xs,onb,sig,original_r,edgecount,searcher,indices,db,RADIUS,final_edge,confirmed_hull)#,ps_data) dim = size(eltype(xs))[1] debug = false#facet_data[1]==[477, 1457, 1540, 1605, 1660, 1980] #[33, 104, 210, 285, 392, 458] #println(facet_data) @@ -574,7 +578,7 @@ function explore_chull_vertex(edgeIterator,facet_data,xs,onb,sig,original_r,edge #@descend raycast_des2(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginalSafety,20,debug) #error("") #generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginal,20,debug,ps_data) - generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),searcher.parameters.method,20,debug,ps_data) + generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),searcher.parameters.method,20,debug)#,ps_data) #generator_1, t, r2 = raycast_des3(full_edge, r, u, xs.rest, searcher,0,full_edge,facet_data[1],Raycast_By_Descend(),RCOriginalSafety,20,debug,ps_data) generator_1==0 && continue @@ -733,7 +737,7 @@ function explore_chull_vertex(edgeIterator,facet_data,xs,onb,sig,original_r,edge #error("$(facet_data[1]), e2=$edge2, fe=$full_edge, ir=$ir, m=$m") end end=# - debug && println("Final: ",final_edge,", ",full_edge,", ir=",sort!(copy(inrange(searcher.tree.tree,r3,norm(r3-xs.rest[final_edge[1]])*1.0000001))),", shift=",sort!(copy(inrange(searcher.tree.tree,r3-1*u,norm(r3-1*u-xs.rest[final_edge[1]])*1.0000001)))," ",norm(r3)," ",norm(u)," ",u, verify(final_edge,r3,u,xs.rest)) + #debug && println("Final: ",final_edge,", ",full_edge,", ir=",sort!(copy(inrange(searcher.tree.tree,r3,norm(r3-xs.rest[final_edge[1]])*1.0000001))),", shift=",sort!(copy(inrange(searcher.tree.tree,r3-1*u,norm(r3-1*u-xs.rest[final_edge[1]])*1.0000001)))," ",norm(r3)," ",norm(u)," ",u, verify(final_edge,r3,u,xs.rest)) vvvv = false HighVoronoi.samecount(final_edge,edge2,dim)==false && continue #error("what the fuck??? $final_edge, $(facet_data[1])") #HighVoronoi.samecount(facet_data[1],edge2,dim-1)==false && error("hier schon Mist: $edge2, $(facet_data[1])") @@ -820,7 +824,7 @@ function rotate_outwards(onb,searcher,sig,r_::P,xs,x0,orientation,final_sig,fina full_sig = copy(sig) deleteat!(full_sig,i_) - if debug + #=if debug println("+",sig,inrange(searcher.tree.tree,r,norm(r-x0)*1.00000001),", ",full_sig) for iii in 1:dim dd = dot(onb[dim],onb[iii]) @@ -835,7 +839,7 @@ function rotate_outwards(onb,searcher,sig,r_::P,xs,x0,orientation,final_sig,fina end end println() - end + end=# #if sig==[248, 294, 379, 397, 418, 21, 406] # raycast_hull(full_sig, r, u, xs, searcher,0,full_sig,last_sig,Raycast_By_Walkray()) # println(last_sig) diff --git a/src/extended.jl b/src/extended.jl index 7322ee6..3f3e4af 100644 --- a/src/extended.jl +++ b/src/extended.jl @@ -175,7 +175,7 @@ function search_vertex2(tree::ExtendedTree,point::MV,idx,dist) where {S,FLOAT<:R #error("") end -mutable struct PlaneSearchData{P,P2,HVN} +#=mutable struct PlaneSearchData{P,P2,HVN} blocked::BitVector c0::Float64 direction::P @@ -255,7 +255,7 @@ function search_vertex_plane(tree::ExtendedTree,data) #where {S,FLOAT<:Real} # , #error("") end - +=# function _nn(tree::ExtendedTree,x::Point,skip=(x->false))::Tuple{Int64,Float64} diff --git a/src/sphere.jl b/src/sphere.jl index 5f1ad02..ac3ff2e 100644 --- a/src/sphere.jl +++ b/src/sphere.jl @@ -67,9 +67,9 @@ boundary(d::SimpleSphericalDomain) = boundary(d.domain) shifts(d::SimpleSphericalDomain) = shifts(d.domain) reference_shifts(d::SimpleSphericalDomain) = reference_shifts(d.domain)#view(reference_shifts(d.domain),1:(reference_length(reference_shifts(domain))-1)) internal_boundary(d::SimpleSphericalDomain) = internal_boundary(d.domain) +VDDomain(d::SimpleSphericalDomain{P}) where {P} = d - -struct PublicSphericalDomain{P,Integral<:HVIntegral{P},D} <: AbstractDomain{P} +#=struct PublicSphericalDomain{P,Integral<:HVIntegral{P},D} <: AbstractDomain{P} integral::Integral domain::D end @@ -85,6 +85,7 @@ boundary(d::PublicSphericalDomain) = boundary(d.domain) shifts(d::PublicSphericalDomain) = shifts(d.domain) reference_shifts(d::PublicSphericalDomain) = view(reference_shifts(d.domain),1:(length(reference_shifts(d.domain))-1)) internal_boundary(d::PublicSphericalDomain) = internal_boundary(d.domain) +=# function internal_boundary(d2::SimpleSphericalDomain,myinte) m2 = mesh(myinte.Integral) diff --git a/src/sphericalmeshview.jl b/src/sphericalmeshview.jl index c810506..ef374b7 100644 --- a/src/sphericalmeshview.jl +++ b/src/sphericalmeshview.jl @@ -170,32 +170,28 @@ end #@inline copy_sig(mesh::LM,sig) where {LM<:LockMesh} = _copy_indeces(mesh,sig,mesh.sigma) @inline Base.getproperty(cd::SphericalMeshView, prop::Symbol) = dyncast_get(cd,Val(prop)) -@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:length_sigma}) = :(getfield(cd,:int_data)[1]) -@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:internal_length_sigma}) = :(getfield(cd,:int_data)[2]) -@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:_Cell}) = :(getfield(cd,:int_data)[3]) +#@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:length_sigma}) = :(getfield(cd,:int_data)[1]) +#@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:internal_length_sigma}) = :(getfield(cd,:int_data)[2]) +#@inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:_Cell}) = :(getfield(cd,:int_data)[3]) @inline @generated dyncast_get(cd::SphericalMeshView, ::Val{:boundary_Vertices}) = :(getfield(cd,:data).boundary_Vertices) @inline @generated dyncast_get(cd::SphericalMeshView, d::Val{S}) where S = :( getfield(cd, S)) -@inline Base.setproperty!(cd::SphericalMeshView, prop::Symbol, val) = dyncast_set(cd,Val(prop),val) -@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:length_sigma},val) = :(getfield(cd,:int_data)[1]=val) -@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:internal_length_sigma},val) = :(getfield(cd,:int_data)[2]=val) -@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:_Cell},val) = :(getfield(cd,:int_data)[3]=val) -@inline @generated dyncast_set(cd::SphericalMeshView, d::Val{S},val) where S = :( setfield(cd, S,val)) +#@inline Base.setproperty!(cd::SphericalMeshView, prop::Symbol, val) = dyncast_set(cd,Val(prop),val) +#@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:length_sigma},val) = :(getfield(cd,:int_data)[1]=val) +#@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:internal_length_sigma},val) = :(getfield(cd,:int_data)[2]=val) +#@inline @generated dyncast_set(cd::SphericalMeshView, ::Val{:_Cell},val) = :(getfield(cd,:int_data)[3]=val) +#@inline @generated dyncast_set(cd::SphericalMeshView, d::Val{S},val) where S = :( setfield(cd, S,val)) @inline Base.length(mv::SphericalMeshView) = length(mv.data) -@inline internal_length(mv::SphericalMeshView) = internal_length(mv.data) +#@inline internal_length(mv::SphericalMeshView) = internal_length(mv.data) -@inline internal_index(m::MV,index::Int64) where MV<:SphericalMeshView = internal_index(m.data, index) -@inline external_index(m::MV,index::Int64) where MV<:SphericalMeshView = external_index(m.data,index) -@inline external_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = external_index(m.data,inds) -@inline internal_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = internal_index(m.data,inds) +#@inline internal_index(m::MV,index::Int64) where MV<:SphericalMeshView = internal_index(m.data, index) +#@inline external_index(m::MV,index::Int64) where MV<:SphericalMeshView = external_index(m.data,index) +@inline external_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:Union{Int64,AbstractVector{Int64}}} = external_index(m.data,inds) +@inline internal_index(m::MV,inds::AVI) where {MV<:SphericalMeshView,AVI<:Union{Int64,AbstractVector{Int64}}} = internal_index(m.data,inds) -@inline internal_sig(mesh::MV,sig::AVI,static::StaticTrue) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = sort!(internal_index(mesh,sig)) -@inline function internal_sig(mesh::MV,sig::AVI,static::StaticFalse) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} - sig .= internal_sig(mesh,sig,statictrue) - return sig -end -@inline external_sig(mesh::MV,sig::AVI,static::StaticTrue) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = sort!(external_index(mesh,sig)) +@inline internal_sig(mesh::MV,sig::AVI,static::SS) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64},SS<:Union{StaticTrue,StaticFalse}} = internal_sig(mesh.data,sig,static) +@inline external_sig(mesh::MV,sig::AVI,static::StaticTrue) where {MV<:SphericalMeshView,AVI<:AbstractVector{Int64}} = external_sig(mesh.data,sig,static) ################################################################################################################################################## @@ -217,7 +213,7 @@ end ################################################################################################################################################## ################################################################################################################################################## -@inline all_vertices_iterator(m::MV, index::Int64, internal::StaticTrue) where MV<:SphericalMeshView = all_vertices_iterator(m.data,index,statictrue) +#=@inline all_vertices_iterator(m::MV, index::Int64, internal::StaticTrue) where MV<:SphericalMeshView = all_vertices_iterator(m.data,index,statictrue) @inline push!(mesh::MV, p::Pair{Vector{Int64},T},index) where {T<:Point,MV<:SphericalMeshView{T}} = push!(mesh.data,p,index) @inline push_ref!(mesh::MV, ref,index) where {T<:Point,MV<:SphericalMeshView{T}} = push_ref!(mesh.data,ref,index) @@ -228,7 +224,7 @@ end @inline mark_delete_vertex!(mesh::MV,sig,i,ii) where MV<:SphericalMeshView = mark_delete_vertex!(mesh.data,sig,i,ii) @inline get_vertex(m::MV,i) where MV<:SphericalMeshView = get_vertex(m.data,i) - +=# diff --git a/src/sphericalpublicview.jl b/src/sphericalpublicview.jl index 92e74e9..c50da5d 100644 --- a/src/sphericalpublicview.jl +++ b/src/sphericalpublicview.jl @@ -210,7 +210,7 @@ Base.eltype(::Type{PublicVolVector{D}}) where D = Float64 # Implement getindex for PublicVolVector @inline Base.getindex(pvv::PublicVolVector, ind::Int64) = ind <= pvv.skip ? 0.0 : pvv.data[ind][end] -struct PublicAVVector{AV<:AbstractVector, D<:AbstractVector{AV}} <: AbstractVector{AV} +#=struct PublicAVVector{AV<:AbstractVector, D<:AbstractVector{AV}} <: AbstractVector{AV} data::D skip::Int64 @@ -237,7 +237,7 @@ function Base.getindex(pavv::PublicAVVector, ind::Int64) throw(ErrorException("Access to indices <= skip is not allowed.")) end return pavv.data[ind] -end +end=# ######################################################################################################################################################################## ######################################################################################################################################################################## @@ -261,7 +261,7 @@ function Base.iterate(it::PublicSphereMeshIterator,state=1) sig = data[1][1] r2 = data[1][2] r = it.center + it.radius * normalize(r2-it.center) - resize(it.sigma,length(sig)-1) + resize!(it.sigma,length(sig)-1) count = 0 count2 = 0 for i in 1:length(sig) @@ -272,7 +272,7 @@ function Base.iterate(it::PublicSphereMeshIterator,state=1) end return (it.sigma,r), st+1 end - return decide(iterate(it.iterator),it,state) + return decide(iterate(it.iterator,state),it,state) end struct PublicSphereNodes{P,N<:HVNodes{P}} <: HVNodes{P} From a27002a6a34b2c74f1f906b914504b8bcda642a1 Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 13:34:42 +0100 Subject: [PATCH 18/19] add some testing --- test/runtests.jl | 1 + test/sphere.jl | 10 +++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index b3df448..7e0afc4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using Test +using Revise using HighVoronoi using SpecialFunctions diff --git a/test/sphere.jl b/test/sphere.jl index 824c0ac..ce4fb56 100644 --- a/test/sphere.jl +++ b/test/sphere.jl @@ -8,9 +8,17 @@ A[1,k] = -A[1,k] end end - vs = VoronoiSphere(VoronoiNodes(A),transformations=(x->-x,), integrate=true, integrator=VI_FAST_POLYGON) + vs = VoronoiSphere(VoronoiNodes(A),transformations=(x->-x,), integrate=true, integrand=x->[x[1]], integrator=VI_FAST_POLYGON) #error() vd = VoronoiData(vs) + println(vd.area[1]) + println(vd.interface_integral[1]) + println(vd.bulk_integral[1]) + m = HighVoronoi.mesh(HighVoronoi.integral(vs.domain)) + for (sig,r) in HighVoronoi.vertices_iterator(m,403) + print("$sig, ") + end + println() @test sum(vd.volume)>6.0 end From cccff1eed9c8ee2eea1bc961ff4861d955cdbc0d Mon Sep 17 00:00:00 2001 From: Martin Heida <122048469+martinheida@users.noreply.github.com> Date: Fri, 17 Jan 2025 13:40:19 +0100 Subject: [PATCH 19/19] Update runtests.jl --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7e0afc4..928f1e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ using Test -using Revise +#using Revise using HighVoronoi using SpecialFunctions