From 1be28cb97be8dbf484a0f212763cb678d3920719 Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 10:39:05 +0200 Subject: [PATCH 1/6] adding general `GraphPolytope` --- src/ReferenceFEs/GraphPolytopes.jl | 650 ++++++++++++++++++ src/ReferenceFEs/Polytopes.jl | 1 - src/ReferenceFEs/ReferenceFEs.jl | 6 + test/ReferenceFEsTests/GraphPolytopesTests.jl | 117 ++++ test/ReferenceFEsTests/runtests.jl | 2 + 5 files changed, 775 insertions(+), 1 deletion(-) create mode 100644 src/ReferenceFEs/GraphPolytopes.jl create mode 100644 test/ReferenceFEsTests/GraphPolytopesTests.jl diff --git a/src/ReferenceFEs/GraphPolytopes.jl b/src/ReferenceFEs/GraphPolytopes.jl new file mode 100644 index 000000000..43d925978 --- /dev/null +++ b/src/ReferenceFEs/GraphPolytopes.jl @@ -0,0 +1,650 @@ +""" + struct GraphPolytope{D,Dp,Tp} <: Polytope{D} + + A graph polytope is a general polytope defined by a set of vertices and a rototation + system (a planar oriented graph). This polytopal representation can represent + any polytope in 2 and 3 dimensions. + + In 2 dimensions ([`Polygon`](@ref)), the representation of the polygon is a closed polyline. + + In 3 dimensions ([`Polyhedron`](@ref)), the rotation system generates the connectivities, each facet is a closed cycle of the graph. + This construction allows complex geometrical operations, e.g., intersecting polytopes by halfspaces. + See also, + + > K. Sugihara, "A robust and consistent algorithm for intersecting convex polyhedra", Comput. Graph. Forum 13 (3) (1994) 45–54, doi: [10.1111/1467-8659.1330045](https://doi.org/10.1111/1467-8659.1330045) + + > D. Powell, T. Abel, "An exact general remeshing scheme applied to physically conservative voxelization", J. Comput. Phys. 297 (Sept. 2015) 340–356, doi: [10.1016/j.jcp.2015.05.022](https://doi.org/10.1016/j.jcp.2015.05.022. + + > S. Badia, P. A. Martorell, F. Verdugo. "Geometrical discretisations for unfitted finite elements on explicit boundary representations", J.Comput. Phys. 460 (2022): 111162. doi: [10.1016/j.jcp.2022.111162](https://doi.org/10.1016/j.jcp.2022.111162) +""" +struct GraphPolytope{D,Dp,Tp,Td} <: Polytope{D} + vertices::Vector{Point{Dp,Tp}} + edge_vertex_graph::Vector{Vector{Int32}} + n_m_to_nface_to_mfaces::Matrix{Vector{Vector{Int}}} + facedims::Vector{Int32} + dimranges::Vector{UnitRange{Int}} + dface_nfaces::Vector{Vector{Int}} + isopen::Bool + metadata::Td + function GraphPolytope{D}( + vertices::Vector{Point{Dp,Tp}}, + edge_vertex_graph::Vector{Vector{Int32}}, + n_m_to_nface_to_mfaces::Matrix{Vector{Vector{Int}}}, + facedims::Vector{Int32}, + dimranges::Vector{UnitRange{Int}}, + dface_nfaces::Vector{Vector{Int}}, + isopen::Bool, + metadata::Td) where {D,Dp,Tp,Td} + new{D,Dp,Tp,Td}( + vertices, + edge_vertex_graph, + n_m_to_nface_to_mfaces, + facedims, + dimranges, + dface_nfaces, + isopen, + metadata) + end +end + +""" + Polygon = GraphPolytope{2} + + A polygon is a [`GraphPolytope`](@ref) in 2 dimensions. +""" +const Polygon = GraphPolytope{2} + +""" + Polyhedron = GraphPolytope{3} + + A polyhedron is a [`GraphPolytope`](@ref) in 3 dimensions. +""" +const Polyhedron = GraphPolytope{3} + +# Constructors + +function GraphPolytope{D}( + vertices::Vector{<:Point}, + graph::Vector{Vector{Int32}}, + isopen::Bool, + data) where D + + n = D+1 + n_m_to_nface_to_mfaces = Matrix{Vector{Vector{Int}}}(undef,n,n ) + dimranges = Vector{UnitRange{Int}}(undef,0) + dface_nfaces = Vector{Vector{Int}}(undef,0) + facedims = Vector{Int32}(undef,0) + GraphPolytope{D}( + vertices, + graph, + n_m_to_nface_to_mfaces, + facedims, + dimranges, + dface_nfaces, + isopen, + data) +end + +function GraphPolytope{D}( + vertices::AbstractVector{<:Point}, + graph::Vector{<:Vector}, + isopen::Bool, + data) where D + + GraphPolytope{D}(collect(vertices),graph,isopen,data) +end + +function GraphPolytope{D}( + vertices::AbstractVector{<:Point}, + graph::Vector{<:Vector} + ;isopen=false::Bool, + metadata=nothing) where D + + GraphPolytope{D}(vertices,graph,isopen,metadata) +end + +function GraphPolytope{D}( + vertices::AbstractVector{<:Point}, + graph::Vector{<:Vector}, + data) where D + + isopen = false + GraphPolytope{D}(vertices,graph,isopen,data) +end + +function get_polytope_data(p::Polytope;metadata=nothing) + get_polytope_data(p,metadata) +end + +function get_polytope_data(p::Polytope,::Nothing) + nothing +end + +function Polygon(p::Polytope{2},vertices::AbstractVector{<:Point};kwargs...) + if p == TRI + e_v_graph = [[2,3],[3,1],[1,2]] + perm = [1,2,3] + elseif p == QUAD + e_v_graph = [[2, 3],[4, 1],[1, 4],[3, 2]] + perm = [1,2,4,3] + else + @unreachable + end + vertices = map(Reindex(vertices),perm) + e_v_graph = map(Reindex(e_v_graph),perm) + e_v_graph = map(i->replace(i, Dict(perm .=> 1:length(perm))...),e_v_graph) + e_v_graph = map(i->Int32.(i),e_v_graph) + data = get_polytope_data(p;kwargs...) + Polygon(vertices,e_v_graph,data) +end + +function Polygon(vertices::AbstractVector{<:Point};kwargs...) + graph = map(1:length(vertices)) do i + inext = i == length(vertices) ? 1 : i+1 + iprev = i == 1 ? length(vertices) : i-1 + Int32[iprev,inext] + end + Polygon(vertices,graph;kwargs...) +end + +function Polyhedron(p::Polytope{3},vertices::AbstractVector{<:Point};kwargs...) + if p == TET + e_v_graph = [[2,4,3],[3,4,1],[1,4,2],[1,2,3]] + elseif p == HEX + e_v_graph = [ + [5, 2, 3], + [6, 4, 1], + [7, 1, 4], + [8, 3, 2], + [1, 7, 6], + [2, 5, 8], + [3, 8, 5], + [4, 6, 7] ] + else + @unreachable + end + e_v_graph = map(i->Int32.(i),e_v_graph) + data = get_polytope_data(p;kwargs...) + Polyhedron(vertices,e_v_graph,data) +end + +function GraphPolytope{D}(p::Polytope;kwargs...) where D + GraphPolytope{D}(p,get_vertex_coordinates(p);kwargs...) +end + +# Interface + +num_dims(::T) where T<:GraphPolytope = num_dims(T) + +num_dims(::Type{<:GraphPolytope{D}}) where D = D + +num_cell_dims(a::GraphPolytope) = num_dims(a) + +point_eltype(::T) where T<:GraphPolytope = point_eltype(T) + +point_eltype(::Type{<:GraphPolytope{D,Dp,T}}) where {D,Dp,T} = T + +num_point_dims(::Type{<:GraphPolytope{D,Dp}}) where {D,Dp} = Dp + +num_vertices(a::GraphPolytope) = length(a.vertices) + +get_vertex_coordinates(a::GraphPolytope) = a.vertices + +Base.getindex(a::GraphPolytope,i::Integer) = a.vertices[i] + +""" + get_graph(p::GraphPolytope) -> Vector{Vector{Int32}} + + It returns the edge-vertex graph of the polytope `p`. +""" +@inline get_graph(a::GraphPolytope) = a.edge_vertex_graph + + +""" + get_data(p::GraphPolytope) + + It return the metadata stored in the polytope `p`. +""" +get_data(a::GraphPolytope) = a.metadata + + +""" + isopen(p::GraphPolytope) -> Bool + + In return whether the polytope is watter tight or not. +""" +Base.isopen(a::GraphPolytope) = a.isopen + +""" + isactive(p::GraphPolytope,vertex::Integer) -> Bool + + It returns whether a vertex is connected to any other vertex. +""" +function isactive(p::Polyhedron,vertex::Integer) + !isempty( get_graph(p)[vertex] ) +end + +""" + check_polytope_graph(p::GraphPolytope) -> Bool + + It checks whether the graph is well-constructed. The graph must be oriented + and planar. +""" +function check_polytope_graph(p::GraphPolytope) + check_polytope_graph(get_graph(p)) +end + +function check_polytope_graph(graph::AbstractVector{<:AbstractVector}) + for v in 1:length(graph) + !isempty(graph[v]) || continue + for vneig in graph[v] + vneig > 0 || continue + v ∈ graph[vneig] || return false + end + end + true +end + +is_simplex(::GraphPolytope) = false + +is_n_cube(::GraphPolytope) = false + +function simplexify(p::GraphPolytope{D}) where D + @assert !isopen(p) + X,T = simplexify_interior(p) + @check X == get_vertex_coordinates(p) + T, simplex_polytope(Val{D}()) +end + +simplex_polytope(::Val{0}) = VERTEX + +simplex_polytope(::Val{1}) = SEGMENT + +simplex_polytope(::Val{2}) = TRI + +simplex_polytope(::Val{3}) = TET + +function Polytope{0}(p::GraphPolytope,faceid::Integer) + VERTEX +end + +function Polytope{1}(p::GraphPolytope,faceid::Integer) + SEGMENT +end + +function Polytope{D}(p::GraphPolytope{D},faceid::Integer) where D + p +end + +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 + +Base.:(==)(a::GraphPolytope,b::GraphPolytope) = false + +function Base.:(==)(a::GraphPolytope{D},b::GraphPolytope{D}) where D + a === b || + ( get_vertex_coordinates(a) == get_vertex_coordinates(b) && + get_graph(a) == get_graph(b) ) +end + +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 + +function get_facet_orientations(p::GraphPolytope) + ones(Int,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_facet_normal(p::Polygon) + D = 2 + f_to_v = get_faces(p,D-1,0) + coords = get_vertex_coordinates(p) + @notimplementedif num_components(eltype(coords)) != 2 + map(f_to_v) do v + e = coords[v[2]]-coords[v[1]] + n = VectorValue( e[2], -e[1] ) + 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 + e = coords[v[2]]-coords[v[1]] + e / norm(e) + end +end + +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 + lens = map(i->num_faces(p,i),0:D) + sums = cumsum(lens) + resize!(p.dimranges,D+1) + 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{Int}}(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 + +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::GraphPolytope,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::GraphPolytope,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 = Gridap.Geometry.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::GraphPolytope,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 = Gridap.Geometry.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_face_to_vertices(p::Polygon,d::Integer) + if d == 1 + generate_facet_to_vertices(p) + else + @unreachable + end +end + +function generate_facet_to_vertices(poly::Polyhedron) + 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 > 0 || 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 > 0 || break + push!(k,vcurrent) + end + if length(k) >=D + push!(T,k) + end + end + end + T +end + +function generate_edge_to_vertices(poly::GraphPolytope) + 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 + +function generate_facet_to_vertices(poly::Polygon) + graph = get_graph(poly) + T = Vector{Int32}[] + for v in 1:length(graph) + vnext = v == length(graph) ? 1 : v+1 + @assert vnext ∈ graph[v] + push!(T,[v,vnext]) + end + T +end + +function Base.copy(poly::Polyhedron) + vertices = get_vertex_coordinates(poly) + graph = get_graph(poly) + open = isopen(poly) + data = get_data(poly) + data = !isnothing(data) ? copy(data) : nothing + Polyhedron(vertices,graph,open,data) +end + +function simplexify_interior(p::Polygon) + @assert !isopen(p) + e_to_v = generate_facet_to_vertices(p) + T = Vector{Int32}[] + if length(e_to_v) > 0 + v0 = e_to_v[1][1] + for verts in e_to_v + if v0 ∉ verts + push!(T,[v0,verts[1],verts[2]]) + end + end + end + get_vertex_coordinates(p),T +end + +function simplexify_interior(poly::Polyhedron) + !isopen(poly) || return simplexify_surface(poly) + vstart = fill(UNSET,num_vertices(poly)) + stack = Int32[] + for v in 1:num_vertices(poly) + isactive(poly,v) || continue + vstart[v] == UNSET || continue + vstart[v] = v + empty!(stack) + push!(stack,v) + while !isempty(stack) + vcurrent = pop!(stack) + for vneig in get_graph(poly)[vcurrent] + if vstart[vneig] == UNSET + vstart[vneig] = v + push!(stack,vneig) + end + end + end + end + istouch = map( i -> falses(length(i)), get_graph(poly) ) + vertex_coordinates = get_vertex_coordinates(poly) + T = Vector{Int32}[] + X = vertex_coordinates + 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] + while vnext != v + inext = findfirst( isequal(vcurrent), get_graph(poly)[vnext] ) + !isnothing(inext) || break + inext = ( inext % length( get_graph(poly)[vnext] ) ) + 1 + istouch[vnext][inext] = true + vcurrent = vnext + vnext = get_graph(poly)[vnext][inext] + @assert vcurrent ≠ vnext + vcurrent ≠ vnext || break + if v ∉ (vstart[v],vcurrent,vnext) + k = [vstart[v],v,vcurrent,vnext] + push!(T,k) + end + end + end + end + X,T +end + +function simplexify_surface(poly::Polyhedron) + 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 > 0 || continue + 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 > 0 || break + if v ∉ (vcurrent,vnext) + k = [v,vcurrent,vnext] + push!(T,k) + end + end + end + end + get_vertex_coordinates(poly),T +end diff --git a/src/ReferenceFEs/Polytopes.jl b/src/ReferenceFEs/Polytopes.jl index 318ef406b..e45d584a1 100644 --- a/src/ReferenceFEs/Polytopes.jl +++ b/src/ReferenceFEs/Polytopes.jl @@ -726,4 +726,3 @@ function test_polytope(p::Polytope{D};optional::Bool=false) where D @test isa(is_n_cube(p),Bool) end end - diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl index f7cd472ca..92b874c48 100644 --- a/src/ReferenceFEs/ReferenceFEs.jl +++ b/src/ReferenceFEs/ReferenceFEs.jl @@ -10,6 +10,7 @@ using DocStringExtensions using LinearAlgebra using Combinatorics using FillArrays +using ..Gridap using Gridap.Helpers using Gridap.Arrays @@ -45,6 +46,9 @@ import Base: == export Polytope export ExtrusionPolytope +export GraphPolytope +export Polygon +export Polyhedron export get_extrusion export get_faces export get_dimranges @@ -213,6 +217,8 @@ include("Polytopes.jl") include("ExtrusionPolytopes.jl") +include("GraphPolytopes.jl") + include("Dofs.jl") include("LagrangianDofBases.jl") diff --git a/test/ReferenceFEsTests/GraphPolytopesTests.jl b/test/ReferenceFEsTests/GraphPolytopesTests.jl new file mode 100644 index 000000000..6c9a716dd --- /dev/null +++ b/test/ReferenceFEsTests/GraphPolytopesTests.jl @@ -0,0 +1,117 @@ +module GraphPolytopesTests + +using Test +using Gridap +using Gridap.ReferenceFEs + +using Gridap.ReferenceFEs: check_polytope_graph +using Gridap.ReferenceFEs: simplexify_surface +using Gridap.ReferenceFEs: simplexify_interior + +p = Polygon(TRI) +test_polytope(p,optional=true) +@test check_polytope_graph(p) +@test get_faces(p,0,0) == [[1],[2],[3]] +@test get_faces(p,1,0) == [[1,2],[2,3],[3,1]] +@test get_faces(p,2,0) == [[1,2,3]] +@test get_faces(p,0,1) == [[1,3],[1,2],[2,3]] +@test get_faces(p,1,1) == [[1],[2],[3]] +@test get_faces(p,2,1) == [[1,2,3]] +@test get_faces(p,0,2) == [[1],[1],[1]] +@test get_faces(p,1,2) == [[1],[1],[1]] +@test get_faces(p,2,2) == [[1]] +@test get_facedims(p) == [0,0,0,1,1,1,2] +@test Polytope{2}(p,1) === p +@test Polytope{1}(p,1) == SEGMENT +@test Polytope{0}(p,1) == VERTEX +@test simplexify(p) == ([[1,2,3]],TRI) + +p = Polygon(QUAD) +test_polytope(p,optional=true) +@test check_polytope_graph(p) +@test get_faces(p,0,0) == [[1],[2],[3],[4]] +@test get_faces(p,1,0) == [[1,2],[2,3],[3,4],[4,1]] +@test get_faces(p,2,0) == [[1,2,3,4]] +@test get_faces(p,0,1) == [[1,4],[1,2],[2,3],[3,4]] +@test get_faces(p,1,1) == [[1],[2],[3],[4]] +@test get_faces(p,2,1) == [[1,2,3,4]] +@test get_faces(p,0,2) == [[1],[1],[1],[1]] +@test get_faces(p,1,2) == [[1],[1],[1],[1]] +@test get_faces(p,2,2) == [[1]] +@test get_facedims(p) == [0,0,0,0,1,1,1,1,2] +@test Polytope{2}(p,1) === p +@test Polytope{1}(p,1) == SEGMENT +@test Polytope{0}(p,1) == VERTEX +@test simplexify(p) == ([[1,2,3],[1,3,4]],TRI) + +p = Polyhedron(TET) +test_polytope(p,optional=true) +@test check_polytope_graph(p) +@test get_faces(p,0,0) == [[1],[2],[3],[4]] +@test get_faces(p,1,0) == [[1,2],[1,4],[1,3],[2,3],[2,4],[3,4]] +@test get_faces(p,2,0) == [[1,2,3],[1,4,2],[1,3,4],[2,4,3]] +@test get_faces(p,3,0) == [[1,2,3,4]] +@test get_faces(p,0,1) == [[1,2,3],[1,4,5],[3,4,6],[2,5,6]] +@test get_faces(p,1,1) == [[1],[2],[3],[4],[5],[6]] +@test get_faces(p,2,1) == [[1,4,3],[2,5,1],[3,6,2],[5,6,4]] +@test get_faces(p,3,1) == [[1,2,3,4,5,6]] +@test get_faces(p,0,2) == [[1,2,3],[1,2,4],[1,3,4],[2,3,4]] +@test get_faces(p,1,2) == [[1,2],[2,3],[1,3],[1,4],[2,4],[3,4]] +@test get_faces(p,2,2) == [[1],[2],[3],[4]] +@test get_faces(p,3,2) == [[1,2,3,4]] +@test get_faces(p,0,3) == [[1],[1],[1],[1]] +@test get_faces(p,1,3) == [[1],[1],[1],[1],[1],[1]] +@test get_faces(p,2,3) == [[1],[1],[1],[1]] +@test get_faces(p,3,3) == [[1]] +@test get_facedims(p) == [0,0,0,0,1,1,1,1,1,1,2,2,2,2,3] +@test Polytope{3}(p,1) === p +@test isa(Polytope{2}(p,1),Polygon) +@test Polytope{1}(p,1) == SEGMENT +@test Polytope{0}(p,1) == VERTEX +@test simplexify(p) == ([[1,2,4,3]],TET) +x,t = simplexify_surface(p) +@test t == get_faces(p,2,0) +x,t = simplexify_interior(p) +@test t == [[1,2,4,3]] + +p = Polyhedron(HEX) +test_polytope(p,optional=true) +@test check_polytope_graph(p) +@test get_faces(p,0,0) == [[1],[2],[3],[4],[5],[6],[7],[8]] +@test get_faces(p,1,0) == [ + [1,5],[1,2],[1,3],[2,6],[2,4],[3,7],[3,4],[4,8],[5,7],[5,6],[6,8],[7,8]] +@test get_faces(p,2,0) == [ + [1,5,7,3],[1,2,6,5],[1,3,4,2],[2,4,8,6],[3,7,8,4],[5,6,8,7]] +@test get_faces(p,3,0) == [[1,2,3,4,5,6,7,8]] +@test get_faces(p,0,1) == [ + [1,2,3],[2,4,5],[3,6,7],[5,7,8],[1,9,10],[4,10,11],[6,9,12],[8,11,12]] +@test get_faces(p,1,1) == [[1],[2],[3],[4],[5],[6],[7],[8],[9],[10],[11],[12]] +@test get_faces(p,2,1) == [ + [1,9,6,3],[2,4,10,1],[3,7,5,2],[5,8,11,4],[6,12,8,7],[10,11,12,9]] +@test get_faces(p,3,1) == [[1,2,3,4,5,6,7,8,9,10,11,12]] +@test get_faces(p,0,2) == [ + [1,2,3],[2,3,4],[1,3,5],[3,4,5],[1,2,6],[2,4,6],[1,5,6],[4,5,6]] +@test get_faces(p,1,2) == [ + [1,2],[2,3],[1,3],[2,4],[3,4],[1,5],[3,5],[4,5],[1,6],[2,6],[4,6],[5,6]] +@test get_faces(p,2,2) == [[1],[2],[3],[4],[5],[6]] +@test get_faces(p,3,2) == [[1,2,3,4,5,6]] +@test get_faces(p,0,3) == [[1],[1],[1],[1],[1],[1],[1],[1]] +@test get_faces(p,1,3) == [[1],[1],[1],[1],[1],[1],[1],[1],[1],[1],[1],[1]] +@test get_faces(p,2,3) == [[1],[1],[1],[1],[1],[1]] +@test get_faces(p,3,3) == [[1]] +@test get_facedims(p) == [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3] +@test Polytope{3}(p,1) === p +@test isa(Polytope{2}(p,1),Polygon) +@test Polytope{1}(p,1) == SEGMENT +@test Polytope{0}(p,1) == VERTEX +@test simplexify(p) == ( + [[1,2,4,8],[1,2,8,6],[1,3,7,8],[1,3,8,4],[1,5,6,8],[1,5,8,7]],TET) +x,t = simplexify_surface(p) +@test t == [ + [1,5,7],[1,7,3],[1,2,6],[1,6,5],[1,3,4],[1,4,2], + [2,4,8],[2,8,6],[3,7,8],[3,8,4],[5,6,8],[5,8,7]] +x,t = simplexify_interior(p) +@test t == [ + [1,2,4,8],[1,2,8,6],[1,3,7,8],[1,3,8,4],[1,5,6,8],[1,5,8,7]] + +end # module diff --git a/test/ReferenceFEsTests/runtests.jl b/test/ReferenceFEsTests/runtests.jl index 45faf6f44..bb3c95a1e 100644 --- a/test/ReferenceFEsTests/runtests.jl +++ b/test/ReferenceFEsTests/runtests.jl @@ -6,6 +6,8 @@ using Test @testset "ExtrusionPolytopes" begin include("ExtrusionPolytopesTests.jl") end +@testset "GraphPolytopes" begin include("GraphPolytopesTests.jl") end + @testset "Dofs" begin include("DofsTests.jl") end @testset "MockDofs" begin include("MockDofsTests.jl") end From 688b9363dfe527bb01902603393ebdfce2fa5554 Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 10:42:46 +0200 Subject: [PATCH 2/6] update news --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index 34d86df6e..9e59496c0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added + +- Define `GraphPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#100X](https://github.com/gridap/Gridap.jl/pull/100X). + ## [0.18.2] - 2024-05-02 ### Fixed From e0c627afe26a1b2b752ce0cb41ae56e2bd2fea33 Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 11:25:13 +0200 Subject: [PATCH 3/6] minor changes at GraphPolytope interface --- src/ReferenceFEs/GraphPolytopes.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ReferenceFEs/GraphPolytopes.jl b/src/ReferenceFEs/GraphPolytopes.jl index 43d925978..c01e4b835 100644 --- a/src/ReferenceFEs/GraphPolytopes.jl +++ b/src/ReferenceFEs/GraphPolytopes.jl @@ -112,11 +112,11 @@ function GraphPolytope{D}( GraphPolytope{D}(vertices,graph,isopen,data) end -function get_polytope_data(p::Polytope;metadata=nothing) - get_polytope_data(p,metadata) +function generate_polytope_data(p::Polytope;metadata=nothing) + generate_polytope_data(p,metadata) end -function get_polytope_data(p::Polytope,::Nothing) +function generate_polytope_data(p::Polytope,::Nothing) nothing end @@ -134,7 +134,7 @@ function Polygon(p::Polytope{2},vertices::AbstractVector{<:Point};kwargs...) e_v_graph = map(Reindex(e_v_graph),perm) e_v_graph = map(i->replace(i, Dict(perm .=> 1:length(perm))...),e_v_graph) e_v_graph = map(i->Int32.(i),e_v_graph) - data = get_polytope_data(p;kwargs...) + data = generate_polytope_data(p;kwargs...) Polygon(vertices,e_v_graph,data) end @@ -164,7 +164,7 @@ function Polyhedron(p::Polytope{3},vertices::AbstractVector{<:Point};kwargs...) @unreachable end e_v_graph = map(i->Int32.(i),e_v_graph) - data = get_polytope_data(p;kwargs...) + data = generate_polytope_data(p;kwargs...) Polyhedron(vertices,e_v_graph,data) end @@ -201,11 +201,11 @@ Base.getindex(a::GraphPolytope,i::Integer) = a.vertices[i] """ - get_data(p::GraphPolytope) + get_metadata(p::GraphPolytope) It return the metadata stored in the polytope `p`. """ -get_data(a::GraphPolytope) = a.metadata +get_metadata(a::GraphPolytope) = a.metadata """ @@ -551,7 +551,7 @@ function Base.copy(poly::Polyhedron) vertices = get_vertex_coordinates(poly) graph = get_graph(poly) open = isopen(poly) - data = get_data(poly) + data = get_metadata(poly) data = !isnothing(data) ? copy(data) : nothing Polyhedron(vertices,graph,open,data) end From a1ad7fee649d682f0e30adb159635c8f1f6f7053 Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 11:41:21 +0200 Subject: [PATCH 4/6] add docstrings --- src/ReferenceFEs/GraphPolytopes.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/ReferenceFEs/GraphPolytopes.jl b/src/ReferenceFEs/GraphPolytopes.jl index c01e4b835..5c48666aa 100644 --- a/src/ReferenceFEs/GraphPolytopes.jl +++ b/src/ReferenceFEs/GraphPolytopes.jl @@ -571,6 +571,13 @@ function simplexify_interior(p::Polygon) get_vertex_coordinates(p),T end +""" + simplexify_interior(p::Polyhedron) + + `simplex_interior` computes a simplex partition of the volume inside + the Polyhedron `p`. + It returns a vector of coordinates and an array of connectivitties. +""" function simplexify_interior(poly::Polyhedron) !isopen(poly) || return simplexify_surface(poly) vstart = fill(UNSET,num_vertices(poly)) @@ -621,6 +628,13 @@ function simplexify_interior(poly::Polyhedron) X,T end +""" + simplexify_surface(p::Polyhedron) + + `simplex_surface` computes a simplex partition of the surface bounding + the Polyhedron `p`. + It returns a vector of coordinates and an array of connectivitties. +""" function simplexify_surface(poly::Polyhedron) istouch = map( i -> falses(length(i)), get_graph(poly) ) T = Vector{Int32}[] From 56ec3923df3cbe47993f2803f891fa2f8f422ce5 Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 11:59:58 +0200 Subject: [PATCH 5/6] edit news --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 9e59496c0..8943fd3c4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Define `GraphPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#100X](https://github.com/gridap/Gridap.jl/pull/100X). +- Define `GraphPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#1006](https://github.com/gridap/Gridap.jl/pull/1006). ## [0.18.2] - 2024-05-02 From 8b703464a0dbeacce7f9bf831d9e37be6b5d787d Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Wed, 15 May 2024 13:55:34 +0200 Subject: [PATCH 6/6] renaming GraphPolytope -> GeneralPolytope --- NEWS.md | 2 +- ...{GraphPolytopes.jl => GeneralPolytopes.jl} | 132 +++++++++--------- src/ReferenceFEs/ReferenceFEs.jl | 4 +- ...topesTests.jl => GeneralPolytopesTests.jl} | 2 +- test/ReferenceFEsTests/runtests.jl | 2 +- 5 files changed, 74 insertions(+), 68 deletions(-) rename src/ReferenceFEs/{GraphPolytopes.jl => GeneralPolytopes.jl} (80%) rename test/ReferenceFEsTests/{GraphPolytopesTests.jl => GeneralPolytopesTests.jl} (99%) diff --git a/NEWS.md b/NEWS.md index 8943fd3c4..0298fb970 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Define `GraphPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#1006](https://github.com/gridap/Gridap.jl/pull/1006). +- Define `GeneralPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#1006](https://github.com/gridap/Gridap.jl/pull/1006). ## [0.18.2] - 2024-05-02 diff --git a/src/ReferenceFEs/GraphPolytopes.jl b/src/ReferenceFEs/GeneralPolytopes.jl similarity index 80% rename from src/ReferenceFEs/GraphPolytopes.jl rename to src/ReferenceFEs/GeneralPolytopes.jl index 5c48666aa..959d0ef30 100644 --- a/src/ReferenceFEs/GraphPolytopes.jl +++ b/src/ReferenceFEs/GeneralPolytopes.jl @@ -1,7 +1,7 @@ """ - struct GraphPolytope{D,Dp,Tp} <: Polytope{D} + struct GeneralPolytope{D,Dp,Tp} <: Polytope{D} - A graph polytope is a general polytope defined by a set of vertices and a rototation + The `GeneralPolytope` is definded defined by a set of vertices and a rototation system (a planar oriented graph). This polytopal representation can represent any polytope in 2 and 3 dimensions. @@ -17,7 +17,7 @@ > S. Badia, P. A. Martorell, F. Verdugo. "Geometrical discretisations for unfitted finite elements on explicit boundary representations", J.Comput. Phys. 460 (2022): 111162. doi: [10.1016/j.jcp.2022.111162](https://doi.org/10.1016/j.jcp.2022.111162) """ -struct GraphPolytope{D,Dp,Tp,Td} <: Polytope{D} +struct GeneralPolytope{D,Dp,Tp,Td} <: Polytope{D} vertices::Vector{Point{Dp,Tp}} edge_vertex_graph::Vector{Vector{Int32}} n_m_to_nface_to_mfaces::Matrix{Vector{Vector{Int}}} @@ -26,7 +26,7 @@ struct GraphPolytope{D,Dp,Tp,Td} <: Polytope{D} dface_nfaces::Vector{Vector{Int}} isopen::Bool metadata::Td - function GraphPolytope{D}( + function GeneralPolytope{D}( vertices::Vector{Point{Dp,Tp}}, edge_vertex_graph::Vector{Vector{Int32}}, n_m_to_nface_to_mfaces::Matrix{Vector{Vector{Int}}}, @@ -48,22 +48,22 @@ struct GraphPolytope{D,Dp,Tp,Td} <: Polytope{D} end """ - Polygon = GraphPolytope{2} + Polygon = GeneralPolytope{2} - A polygon is a [`GraphPolytope`](@ref) in 2 dimensions. + A polygon is a [`GeneralPolytope`](@ref) in 2 dimensions. """ -const Polygon = GraphPolytope{2} +const Polygon = GeneralPolytope{2} """ - Polyhedron = GraphPolytope{3} + Polyhedron = GeneralPolytope{3} - A polyhedron is a [`GraphPolytope`](@ref) in 3 dimensions. + A polyhedron is a [`GeneralPolytope`](@ref) in 3 dimensions. """ -const Polyhedron = GraphPolytope{3} +const Polyhedron = GeneralPolytope{3} # Constructors -function GraphPolytope{D}( +function GeneralPolytope{D}( vertices::Vector{<:Point}, graph::Vector{Vector{Int32}}, isopen::Bool, @@ -74,7 +74,7 @@ function GraphPolytope{D}( dimranges = Vector{UnitRange{Int}}(undef,0) dface_nfaces = Vector{Vector{Int}}(undef,0) facedims = Vector{Int32}(undef,0) - GraphPolytope{D}( + GeneralPolytope{D}( vertices, graph, n_m_to_nface_to_mfaces, @@ -85,31 +85,37 @@ function GraphPolytope{D}( data) end -function GraphPolytope{D}( +function GeneralPolytope{D}( vertices::AbstractVector{<:Point}, graph::Vector{<:Vector}, isopen::Bool, data) where D - GraphPolytope{D}(collect(vertices),graph,isopen,data) + GeneralPolytope{D}(collect(vertices),graph,isopen,data) end -function GraphPolytope{D}( +""" + GeneralPolytope{D}(vertices,graph;kwargs...) + + Constructor of a [`GeneralPolytope`](@ref) that generates a polytope of + D dimensions with the given `vertices` and `graph` of connectivities. +""" +function GeneralPolytope{D}( vertices::AbstractVector{<:Point}, graph::Vector{<:Vector} ;isopen=false::Bool, metadata=nothing) where D - GraphPolytope{D}(vertices,graph,isopen,metadata) + GeneralPolytope{D}(vertices,graph,isopen,metadata) end -function GraphPolytope{D}( +function GeneralPolytope{D}( vertices::AbstractVector{<:Point}, graph::Vector{<:Vector}, data) where D isopen = false - GraphPolytope{D}(vertices,graph,isopen,data) + GeneralPolytope{D}(vertices,graph,isopen,data) end function generate_polytope_data(p::Polytope;metadata=nothing) @@ -168,55 +174,55 @@ function Polyhedron(p::Polytope{3},vertices::AbstractVector{<:Point};kwargs...) Polyhedron(vertices,e_v_graph,data) end -function GraphPolytope{D}(p::Polytope;kwargs...) where D - GraphPolytope{D}(p,get_vertex_coordinates(p);kwargs...) +function GeneralPolytope{D}(p::Polytope;kwargs...) where D + GeneralPolytope{D}(p,get_vertex_coordinates(p);kwargs...) end # Interface -num_dims(::T) where T<:GraphPolytope = num_dims(T) +num_dims(::T) where T<:GeneralPolytope = num_dims(T) -num_dims(::Type{<:GraphPolytope{D}}) where D = D +num_dims(::Type{<:GeneralPolytope{D}}) where D = D -num_cell_dims(a::GraphPolytope) = num_dims(a) +num_cell_dims(a::GeneralPolytope) = num_dims(a) -point_eltype(::T) where T<:GraphPolytope = point_eltype(T) +point_eltype(::T) where T<:GeneralPolytope = point_eltype(T) -point_eltype(::Type{<:GraphPolytope{D,Dp,T}}) where {D,Dp,T} = T +point_eltype(::Type{<:GeneralPolytope{D,Dp,T}}) where {D,Dp,T} = T -num_point_dims(::Type{<:GraphPolytope{D,Dp}}) where {D,Dp} = Dp +num_point_dims(::Type{<:GeneralPolytope{D,Dp}}) where {D,Dp} = Dp -num_vertices(a::GraphPolytope) = length(a.vertices) +num_vertices(a::GeneralPolytope) = length(a.vertices) -get_vertex_coordinates(a::GraphPolytope) = a.vertices +get_vertex_coordinates(a::GeneralPolytope) = a.vertices -Base.getindex(a::GraphPolytope,i::Integer) = a.vertices[i] +Base.getindex(a::GeneralPolytope,i::Integer) = a.vertices[i] """ - get_graph(p::GraphPolytope) -> Vector{Vector{Int32}} + get_graph(p::GeneralPolytope) -> Vector{Vector{Int32}} It returns the edge-vertex graph of the polytope `p`. """ -@inline get_graph(a::GraphPolytope) = a.edge_vertex_graph +@inline get_graph(a::GeneralPolytope) = a.edge_vertex_graph """ - get_metadata(p::GraphPolytope) + get_metadata(p::GeneralPolytope) It return the metadata stored in the polytope `p`. """ -get_metadata(a::GraphPolytope) = a.metadata +get_metadata(a::GeneralPolytope) = a.metadata """ - isopen(p::GraphPolytope) -> Bool + isopen(p::GeneralPolytope) -> Bool In return whether the polytope is watter tight or not. """ -Base.isopen(a::GraphPolytope) = a.isopen +Base.isopen(a::GeneralPolytope) = a.isopen """ - isactive(p::GraphPolytope,vertex::Integer) -> Bool + isactive(p::GeneralPolytope,vertex::Integer) -> Bool It returns whether a vertex is connected to any other vertex. """ @@ -225,12 +231,12 @@ function isactive(p::Polyhedron,vertex::Integer) end """ - check_polytope_graph(p::GraphPolytope) -> Bool + check_polytope_graph(p::GeneralPolytope) -> Bool It checks whether the graph is well-constructed. The graph must be oriented and planar. """ -function check_polytope_graph(p::GraphPolytope) +function check_polytope_graph(p::GeneralPolytope) check_polytope_graph(get_graph(p)) end @@ -245,11 +251,11 @@ function check_polytope_graph(graph::AbstractVector{<:AbstractVector}) true end -is_simplex(::GraphPolytope) = false +is_simplex(::GeneralPolytope) = false -is_n_cube(::GraphPolytope) = false +is_n_cube(::GeneralPolytope) = false -function simplexify(p::GraphPolytope{D}) where D +function simplexify(p::GeneralPolytope{D}) where D @assert !isopen(p) X,T = simplexify_interior(p) @check X == get_vertex_coordinates(p) @@ -264,15 +270,15 @@ simplex_polytope(::Val{2}) = TRI simplex_polytope(::Val{3}) = TET -function Polytope{0}(p::GraphPolytope,faceid::Integer) +function Polytope{0}(p::GeneralPolytope,faceid::Integer) VERTEX end -function Polytope{1}(p::GraphPolytope,faceid::Integer) +function Polytope{1}(p::GeneralPolytope,faceid::Integer) SEGMENT end -function Polytope{D}(p::GraphPolytope{D},faceid::Integer) where D +function Polytope{D}(p::GeneralPolytope{D},faceid::Integer) where D p end @@ -283,15 +289,15 @@ function Polytope{2}(p::Polyhedron,faceid::Integer) Polygon(vertices) end -Base.:(==)(a::GraphPolytope,b::GraphPolytope) = false +Base.:(==)(a::GeneralPolytope,b::GeneralPolytope) = false -function Base.:(==)(a::GraphPolytope{D},b::GraphPolytope{D}) where D +function Base.:(==)(a::GeneralPolytope{D},b::GeneralPolytope{D}) where D a === b || ( get_vertex_coordinates(a) == get_vertex_coordinates(b) && get_graph(a) == get_graph(b) ) end -function num_faces(p::GraphPolytope{D},d::Integer) where D +function num_faces(p::GeneralPolytope{D},d::Integer) where D if d == 0 num_vertices(p) elseif d == D @@ -301,12 +307,12 @@ function num_faces(p::GraphPolytope{D},d::Integer) where D end end -function get_faces(p::GraphPolytope,n::Integer,m::Integer) +function get_faces(p::GeneralPolytope,n::Integer,m::Integer) setup_faces!(p,n,m) p.n_m_to_nface_to_mfaces[n+1,m+1] end -function get_facet_orientations(p::GraphPolytope) +function get_facet_orientations(p::GeneralPolytope) ones(Int,num_facets(p)) end @@ -336,7 +342,7 @@ function get_facet_normal(p::Polygon) end end -function get_edge_tangent(p::GraphPolytope) +function get_edge_tangent(p::GeneralPolytope) e_to_v = get_faces(p,1,0) coords = get_vertex_coordinates(p) map(e_to_v) do v @@ -345,27 +351,27 @@ function get_edge_tangent(p::GraphPolytope) end end -function get_dimranges(p::GraphPolytope) +function get_dimranges(p::GeneralPolytope) setup_dimranges!(p) p.dimranges end -function get_dimrange(p::GraphPolytope,d::Integer) +function get_dimrange(p::GeneralPolytope,d::Integer) setup_dimranges!(p) p.dimranges[d+1] end -function get_faces(p::GraphPolytope) +function get_faces(p::GeneralPolytope) setup_face_to_nfaces!(p) p.dface_nfaces end -function get_facedims(p::GraphPolytope) +function get_facedims(p::GeneralPolytope) setup_facedims!(p) p.facedims end -function setup_dimranges!(p::GraphPolytope{D}) where D +function setup_dimranges!(p::GeneralPolytope{D}) where D if length(p.dimranges) < D+1 lens = map(i->num_faces(p,i),0:D) sums = cumsum(lens) @@ -376,7 +382,7 @@ function setup_dimranges!(p::GraphPolytope{D}) where D end end -function setup_face_to_nfaces!(p::GraphPolytope) +function setup_face_to_nfaces!(p::GeneralPolytope) if length(p.dface_nfaces) < num_vertices(p) facedims = get_facedims(p) dface_nfaces = Vector{Vector{Int}}(undef,length(facedims)) @@ -394,7 +400,7 @@ function setup_face_to_nfaces!(p::GraphPolytope) nothing end -function setup_facedims!(p::GraphPolytope) +function setup_facedims!(p::GeneralPolytope) if length(p.facedims) < num_vertices(p) dimranges = get_dimranges(p) n_faces = dimranges[end][end] @@ -403,7 +409,7 @@ function setup_facedims!(p::GraphPolytope) end end -function setup_faces!(p::GraphPolytope{D},dimfrom,dimto) where D +function setup_faces!(p::GeneralPolytope{D},dimfrom,dimto) where D if isassigned(p.n_m_to_nface_to_mfaces,dimfrom+1,dimto+1) return nothing end @@ -421,14 +427,14 @@ function setup_faces!(p::GraphPolytope{D},dimfrom,dimto) where D nothing end -function setup_face_to_vertices!(p::GraphPolytope,d) +function setup_face_to_vertices!(p::GeneralPolytope,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 +function setup_cell_to_faces!(p::GeneralPolytope{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) ] @@ -436,7 +442,7 @@ function setup_cell_to_faces!(p::GraphPolytope{D},d) where D end end -function setup_nface_to_nface!(p::GraphPolytope,n) +function setup_nface_to_nface!(p::GeneralPolytope,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) @@ -465,7 +471,7 @@ function setup_nface_to_mface!(p::Polyhedron,n,m) end end -function setup_nface_to_mface_dual!(p::GraphPolytope,dimto,dimfrom) +function setup_nface_to_mface_dual!(p::GeneralPolytope,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) @@ -523,7 +529,7 @@ function generate_facet_to_vertices(poly::Polyhedron) T end -function generate_edge_to_vertices(poly::GraphPolytope) +function generate_edge_to_vertices(poly::GeneralPolytope) graph = get_graph(poly) T = Vector{Int32}[] for v in 1:length(graph) diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl index 92b874c48..dc3d56205 100644 --- a/src/ReferenceFEs/ReferenceFEs.jl +++ b/src/ReferenceFEs/ReferenceFEs.jl @@ -46,7 +46,7 @@ import Base: == export Polytope export ExtrusionPolytope -export GraphPolytope +export GeneralPolytope export Polygon export Polyhedron export get_extrusion @@ -217,7 +217,7 @@ include("Polytopes.jl") include("ExtrusionPolytopes.jl") -include("GraphPolytopes.jl") +include("GeneralPolytopes.jl") include("Dofs.jl") diff --git a/test/ReferenceFEsTests/GraphPolytopesTests.jl b/test/ReferenceFEsTests/GeneralPolytopesTests.jl similarity index 99% rename from test/ReferenceFEsTests/GraphPolytopesTests.jl rename to test/ReferenceFEsTests/GeneralPolytopesTests.jl index 6c9a716dd..7c9b87b73 100644 --- a/test/ReferenceFEsTests/GraphPolytopesTests.jl +++ b/test/ReferenceFEsTests/GeneralPolytopesTests.jl @@ -1,4 +1,4 @@ -module GraphPolytopesTests +module GeneralPolytopesTests using Test using Gridap diff --git a/test/ReferenceFEsTests/runtests.jl b/test/ReferenceFEsTests/runtests.jl index bb3c95a1e..47cd088da 100644 --- a/test/ReferenceFEsTests/runtests.jl +++ b/test/ReferenceFEsTests/runtests.jl @@ -6,7 +6,7 @@ using Test @testset "ExtrusionPolytopes" begin include("ExtrusionPolytopesTests.jl") end -@testset "GraphPolytopes" begin include("GraphPolytopesTests.jl") end +@testset "GeneralPolytopes" begin include("GeneralPolytopesTests.jl") end @testset "Dofs" begin include("DofsTests.jl") end