Skip to content

Commit

Permalink
[polytope] add polytope connectivities with deferred computations
Browse files Browse the repository at this point in the history
  • Loading branch information
@pmartorell
pmartorell committed May 13, 2024
1 parent 50f1302 commit 1c642bb
Show file tree
Hide file tree
Showing 3 changed files with 296 additions and 16 deletions.
295 changes: 285 additions & 10 deletions src/Polyhedra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,10 @@ end
struct Polyhedron{Dp,Tp,Td} <: GraphPolytope{3}
vertices::Vector{Point{Dp,Tp}}
edge_vertex_graph::Vector{Vector{Int32}}
n_m_to_nface_to_mfaces::Matrix{Vector{Vector{Int32}}}
facedims::Vector{Int32}
dimranges::Vector{UnitRange{Int}}
dface_nfaces::Vector{Vector{Int32}}
isopen::Bool
data::Td
end
Expand All @@ -30,6 +34,12 @@ end

# Constructors

function Polygon(p::Polytope{2},vertices::AbstractVector{<:Point})
@notimplemented
data = polyhedron_data(p)
Polygon(vertices,data)
end

function Polyhedron(stl::GridTopology)
𝓖 = compute_graph(stl)
X = get_vertex_coordinates(stl)
Expand All @@ -39,6 +49,29 @@ function Polyhedron(stl::GridTopology)
p
end

function Polyhedron(
vertices::Vector{<:Point},
graph::Vector{Vector{Int32}},
isopen::Bool,
data)

D = 3
n = D+1
n_m_to_nface_to_mfaces = Matrix{Vector{Vector{Int32}}}(undef,n,n )
dimranges = Vector{UnitRange{Int}}(undef,0)
dface_nfaces = Vector{Vector{Int32}}(undef,0)
facedims = Vector{Int32}(undef,0)
Polyhedron(
vertices,
graph,
n_m_to_nface_to_mfaces,
facedims,
dimranges,
dface_nfaces,
isopen,
data)
end

function Polyhedron(
vertices::AbstractVector{<:Point},
graph::Vector{<:Vector},
Expand All @@ -48,7 +81,6 @@ function Polyhedron(
Polyhedron(collect(vertices),graph,isopen,data)
end


function Polyhedron(
vertices::AbstractVector{<:Point},
graph::Vector{<:Vector}
Expand Down Expand Up @@ -1302,23 +1334,266 @@ simplex_polytope(::Val{2}) = TRI

simplex_polytope(::Val{3}) = TET


function get_faces(::GraphPolytope,args...)
@notimplemented
function Polytope{0}(p::GraphPolytope,faceid::Integer)
VERTEX
end

function get_dimranges(::GraphPolytope,args...)
@notimplemented
function Polytope{1}(p::GraphPolytope,faceid::Integer)
SEGMENT
end

function get_facedimms(::GraphPolytope)
@notimplemented
function Polytope{D}(p::GraphPolytope{D},faceid::Integer) where D
p
end

function Polytope{N}(p::GraphPolytope,faceid::Integer) where N
@notimplemented
function Polytope{2}(p::Polyhedron,faceid::Integer)
f_to_v = get_faces(p,2,0)
coords = get_vertex_coordinates(p)
vertices = coords[f_to_v[faceid]]
Polygon(vertices)
end

import Base: ==

(==)(a::GraphPolytope,b::GraphPolytope) = false

function num_faces(p::GraphPolytope{D},d::Integer) where D
if d == 0
num_vertices(p)
elseif d == D
1
else
length( get_faces(p,d,0) )
end
end

function get_faces(p::GraphPolytope,n::Integer,m::Integer)
setup_faces!(p,n,m)
p.n_m_to_nface_to_mfaces[n+1,m+1]
end

import Gridap.ReferenceFEs: get_facet_orientations
import Gridap.ReferenceFEs: get_facet_normal
import Gridap.ReferenceFEs: get_edge_tangent
function get_facet_orientations(p::GraphPolytope)
Ones(num_facets(p))
end

function get_facet_normal(p::Polyhedron)
D = 3
f_to_v = get_faces(p,D-1,0)
coords = get_vertex_coordinates(p)
map(f_to_v) do v
v1 = coords[v[2]]-coords[v[1]]
v2 = coords[v[3]]-coords[v[1]]
v1 /= norm(v1)
v2 /= norm(v2)
n = v1 × v2
n /= norm(n)
end
end

function get_edge_tangent(p::GraphPolytope)
e_to_v = get_faces(p,1,0)
coords = get_vertex_coordinates(p)
map(e_to_v) do v
v = coords[v[2]]-coords[v[1]]
v / norm(v)
end
end

import Gridap.ReferenceFEs: get_dimranges
import Gridap.ReferenceFEs: get_dimrange


function get_dimranges(p::GraphPolytope)
setup_dimranges!(p)
p.dimranges
end

function get_dimrange(p::GraphPolytope,d::Integer)
setup_dimranges!(p)
p.dimranges[d+1]
end

function get_faces(p::GraphPolytope)
setup_face_to_nfaces!(p)
p.dface_nfaces
end

function get_facedims(p::GraphPolytope)
setup_facedims!(p)
p.facedims
end

function setup_dimranges!(p::GraphPolytope{D}) where D
if length(p.dimranges) < D+1
resize!(p.dimranges,D+1)
lens = map(i->num_faces(p,i),0:D)
sums = cumsum(lens)
for (i,(l,s)) in enumerate(zip(lens,sums))
p.dimranges[i] = s-l+1 : s
end
end
end

function setup_face_to_nfaces!(p::GraphPolytope)
if length(p.dface_nfaces) < num_vertices(p)
facedims = get_facedims(p)
dface_nfaces = Vector{Vector{Int32}}(undef,length(facedims))
ofsets = get_offsets(p)
for (f,d) in enumerate(facedims)
df = f - ofsets[d+1]
nfs = Int[]
for n in 0:d
union!(nfs,get_faces(p,d,n)[df] .+ ofsets[n+1])
end
dface_nfaces[f] = nfs
end
copy!(p.dface_nfaces,dface_nfaces)
end
nothing
end

using Gridap.ReferenceFEs: _nface_to_nfacedim
function setup_facedims!(p::GraphPolytope)
if length(p.facedims) < num_vertices(p)
dimranges = get_dimranges(p)
n_faces = dimranges[end][end]
facedims = _nface_to_nfacedim(n_faces,dimranges)
copy!(p.facedims,facedims)
end
end

function setup_faces!(p::GraphPolytope{D},dimfrom,dimto) where D
if isassigned(p.n_m_to_nface_to_mfaces,dimfrom+1,dimto+1)
return nothing
end
if dimfrom == dimto
setup_nface_to_nface!(p,dimfrom)
elseif dimfrom == D
setup_cell_to_faces!(p,dimto)
elseif dimto == 0
setup_face_to_vertices!(p,dimfrom)
elseif dimfrom > dimto
setup_nface_to_mface!(p,dimfrom,dimto)
else
setup_nface_to_mface_dual!(p,dimto,dimfrom)
end
nothing
end

function setup_face_to_vertices!(p::Polyhedron,d)
if !isassigned(p.n_m_to_nface_to_mfaces,d+1,1)
df_to_v = generate_face_to_vertices(p,d)
p.n_m_to_nface_to_mfaces[d+1,1] = df_to_v
end
end

function setup_cell_to_faces!(p::GraphPolytope{D},d) where D
if !isassigned(p.n_m_to_nface_to_mfaces,D+1,d+1)
num_f = num_faces(p,d)
c_to_f = [ collect(1:num_f) ]
p.n_m_to_nface_to_mfaces[D+1,d+1] = c_to_f
end
end

function setup_nface_to_nface!(p::Polyhedron,n)
if !isassigned(p.n_m_to_nface_to_mfaces,n+1,n+1)
num_nf = num_faces(p,n)
nf_to_nf = map(i->Int[i],1:num_nf)
p.n_m_to_nface_to_mfaces[n+1,n+1] = nf_to_nf
end
end

function setup_nface_to_mface!(p::Polyhedron,n,m)
if !isassigned(p.n_m_to_nface_to_mfaces,n+1,m+1)
@notimplementedif n != 2 && m != 1
nf_to_v = get_faces(p,n,0)
v_to_mf = get_faces(p,0,m)
nf_to_ftype = map( length, nf_to_v )
ftype_to_lmf_to_lv = map(1:maximum(nf_to_ftype)) do ftype
map(1:ftype) do i
inext = i == ftype ? 1 : i+1
Int[i,inext]
end
end
nf_to_mf = find_cell_to_faces(
Table(nf_to_v),
ftype_to_lmf_to_lv,
nf_to_ftype,
Table(v_to_mf))
p.n_m_to_nface_to_mfaces[n+1,m+1] = Vector(nf_to_mf)
end
end

function setup_nface_to_mface_dual!(p::Polyhedron,dimto,dimfrom)
if !isassigned(p.n_m_to_nface_to_mfaces,dimfrom+1,dimto+1)
@assert dimfrom < dimto
nf_to_mf = get_faces(p,dimto,dimfrom)
n_mf = num_faces(p,dimfrom)
mf_to_nf = generate_cells_around(Table(nf_to_mf),n_mf)
p.n_m_to_nface_to_mfaces[dimfrom+1,dimto+1] = Vector(mf_to_nf)
end
end


function generate_face_to_vertices(p::Polyhedron,d::Integer)
if d == 1
generate_edge_to_vertices(p)
elseif d == 2
generate_facet_to_vertices(p)
else
@unreachable
end
end

function generate_facet_to_vertices(poly::Polyhedron{3})
D = 3
istouch = map( i -> falses(length(i)), get_graph(poly) )
T = Vector{Int32}[]
for v in 1:num_vertices(poly)
isactive(poly,v) || continue
for i in 1:length(get_graph(poly)[v])
!istouch[v][i] || continue
istouch[v][i] = true
vcurrent = v
vnext = get_graph(poly)[v][i]
vnext ∉ (OPEN,UNSET) || continue
k = [v]
while vnext != v
inext = findfirst( isequal(vcurrent), get_graph(poly)[vnext] )
inext = ( inext % length( get_graph(poly)[vnext] ) ) + 1
istouch[vnext][inext] = true
vcurrent = vnext
vnext = get_graph(poly)[vnext][inext]
vnext ∉ (OPEN,UNSET) || break
push!(k,vcurrent)
end
if length(k) >=D
push!(T,k)
end
end
end
T
end

function generate_edge_to_vertices(poly::Polyhedron{3})
graph = get_graph(poly)
T = Vector{Int32}[]
for v in 1:length(graph)
for vneig in graph[v]
if vneig > v
push!(T,[v,vneig])
end
end
end
T
end



# get_edge_tangent
# get_facet_normal
# get_facet_orientations
# get_vertex_permutations # not implemented
2 changes: 2 additions & 0 deletions src/STLCutters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ using PartitionedArrays

using Gridap.Geometry: get_face_to_parent_face
using Gridap.Geometry: get_cell_to_parent_cell
using Gridap.Geometry: find_cell_to_faces
using Gridap.Geometry: generate_cells_around
using GridapEmbedded.Distributed: consistent_bgcell_to_inoutcut!
using GridapEmbedded.Distributed: consistent_bgfacet_to_inoutcut!
using GridapEmbedded.Distributed: DistributedEmbeddedDiscretization
Expand Down
15 changes: 9 additions & 6 deletions test/PolyhedraTests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,13 @@ using STLCutters: merge_nodes
using STLCutters: simplexify_cell_boundary
using STLCutters: compute_grid

p = Polyhedron(TRI)
@test check_graph(p)
#@test volume(p) ≈ 1/2
#@test surface(p) ≈ 2+√2
# p = Polyhedron(TRI)
# @test check_graph(p)
# #@test volume(p) ≈ 1/2
# #@test surface(p) ≈ 2+√2

p = Polyhedron(QUAD)
@test check_graph(p)
# p = Polyhedron(QUAD)
# @test check_graph(p)
#@test volume(p) ≈ 1
#@test surface(p) ≈ 4

Expand All @@ -39,6 +39,9 @@ p = Polyhedron(TET)
@test surface(p) ≈ 3/2 + (√3)/2

p = Polyhedron(HEX)
STLCutters.generate_facets(p)
STLCutters.generate_edges(p)

@test check_graph(p)
@test volume(p) ≈ 1
@test surface(p) ≈ 6
Expand Down

0 comments on commit 1c642bb

Please sign in to comment.