[WIP] typed cells (as opposed to SVector{n,Int})

src/baryref.jl

@@ -94,8 +94,8 @@ function barycentric_refinement(mesh::Mesh{U,2}; sort=:spacefillingcurve) where
         c = NV + E
         a = Edges.faces[E][1]
         b = Edges.faces[E][2]
-        edges[2(E-1) + 1] = SVector(a,c)
-        edges[2(E-1) + 2] = SVector(c,b)
+        edges[2(E-1) + 1] = SimplexGraph(a,c)
+        edges[2(E-1) + 2] = SimplexGraph(c,b)
     end
 
     # return the new mesh after refienments
@@ -154,13 +154,15 @@ function barycentric_refinement(mesh::AbstractMesh{U,3}; sort=:spacefillingcurve
     end
 
     Nodes = skeleton(mesh, 0)
-    node_ctrs = [vertices(Nodes)[node][1] for node in cells(Nodes)]
+    # node_ctrs = [vertices(Nodes)[node][1] for node in cells(Nodes)]
+    node_ctrs = [cartesian(center(chart(Nodes, node))) for node in cells(Nodes)]
     Nodes.faces = Nodes.faces[sort_sfc(node_ctrs)]
 
+    fcs = [SimplexGraph{3}(fc) for fc in fcs]
     fine = Mesh(verts, fcs)
     D = connectivity(Nodes, fine)
     rows, vals = rowvals(D), nonzeros(D)
-    sorted_fcs = Vector{indextype(mesh)}()
+    sorted_fcs = Vector{celltype(mesh)}()
     for (i,Node) in enumerate(cells(Nodes))
         for k in nzrange(D,i)
             j = rows[k]
@@ -336,7 +338,8 @@ function bisecting_refinement(mesh::Mesh{U,3}) where U
 
     # add four faces in each coarse face
     nf = 4NF
-    fcs = zeros(indextype(mesh), nf)
+    # fcs = zeros(celltype(mesh), nf)
+    fcs = Vector{celltype(mesh)}(undef, nf)
 
     for F in 1 : numcells(faces)
 
@@ -356,10 +359,10 @@ function bisecting_refinement(mesh::Mesh{U,3}) where U
         r2 = mesh.faces[F][2]
         r3 = mesh.faces[F][3]
 
-        fcs[4(F-1)+1] = index(r1,es[3],es[2])
-        fcs[4(F-1)+2] = index(r2,es[1],es[3])
-        fcs[4(F-1)+3] = index(r3,es[2],es[1])
-        fcs[4(F-1)+4] = index(es[1],es[2],es[3])
+        fcs[4(F-1)+1] = SimplexGraph(r1,es[3],es[2])
+        fcs[4(F-1)+2] = SimplexGraph(r2,es[1],es[3])
+        fcs[4(F-1)+3] = SimplexGraph(r3,es[2],es[1])
+        fcs[4(F-1)+4] = SimplexGraph(es[1],es[2],es[3])
     end
 
     Mesh(verts, fcs)

src/fileio/gid.jl

@@ -33,13 +33,13 @@ function load_gid_mesh(file,T)
     readline(file); # Elements
 
     # Read the triangles
-    triangles = Vector{Pt{3,Int}}(undef,lineCount-vertexCount);
+    triangles = Vector{SimplexGraph{3}}(undef,lineCount-vertexCount);
     triangleCount = 0;
     while !eof(file)
             line = split(readline(file));
             size(line, 1) == 4 || break; # End Elements
             triangleCount += 1;
-            triangles[triangleCount] = Pt{3,Int}(map(Meta.parse, line[2:4]));
+            triangles[triangleCount] = SimplexGraph{3}(map(Meta.parse, line[2:4]));
     end
     triangles = triangles[1:triangleCount]
     # Mesh(vertices[1:vertexCount], triangles[1:triangleCount])
@@ -48,7 +48,7 @@ function load_gid_mesh(file,T)
     Q = Pt{3,T}
     ctrs = Q[]
     for tr in triangles
-            ctr = sum(vertices[tr])/3
+            ctr = sum(vertices[tr.indices])/3
             push!(ctrs,ctr)
     end
     sorted = sort_sfc(ctrs)

src/fileio/gmsh.jl

@@ -87,8 +87,10 @@ function load_gmsh_mesh(meshfile;
 
     # TODO: once we support hexahedrons, we need to catch that one too
     if element == :quadrangle
+        elements = [QuadrilateralGraph(el) for el in elements]
         return QuadMesh(vertices, elements)
     else
+        elements = [SimplexGraph(el) for el in elements]
         return Mesh(vertices, elements)
     end
 end
@@ -374,7 +376,7 @@ function read_gmsh_mesh(io; physical=nothing, dimension=2, sort=true, T=Float64)
         d[2] == type || continue
         entity_tag == d[4] || entity_tag == 0 || continue
         # f[i +=1] = SVector(d[end-2], d[end-1], d[end-0])
-        f[i+=1] = d[end-dimension:end]
+        f[i+=1] = I(d[end-dimension:end])
     end
     resize!(f,i)
 
@@ -396,6 +398,7 @@ function read_gmsh_mesh(io; physical=nothing, dimension=2, sort=true, T=Float64)
 
 
 
+    f = [SimplexGraph(fc) for fc in f]
     return Mesh(v, f)
 
 end

src/fileio/gmsh3d.jl

@@ -92,6 +92,8 @@ function read_gmsh3d_mesh(io; physical=nothing, T=Float64)
             f[i] = @SVector[fc[1],fc[2],fc[4],fc[3]]
         end
     end
+
+    f = [SimplexGraph{4}(fc) for fc in f]
     msh = Mesh(v, f)
     @assert isoriented(msh)
     return msh

src/fileio/readmesh.jl

@@ -52,8 +52,9 @@ function readmesh(filename;T=Float64)
 
         # TODO: remove explicit reference to dimension
         #I = SVector{3,Int}
-        C = SVector{dim1,Int}
-        faces = zeros(C, num_faces)
+        # C = SVector{dim1,Int}
+        C = SimplexGraph{dim1}
+        faces = Vector{C}(undef, num_faces)
         for i in 1 : num_faces
             l = readline(f)
             faces[i] = C([parse(Int,s) for s in split(l)])

src/mesh.jl

@@ -13,26 +13,20 @@ export vertexarray, cellarray
 abstract type AbstractMesh{U,D1,T} end
 
 
+"""
+    struct Mesh
+
+Basic mesh structure representing flat-faceted simplicial meshes.
+"""
 mutable struct Mesh{U,D1,T} <: AbstractMesh{U,D1,T}
     vertices::Vector{SVector{U,T}}
-    faces::Vector{SVector{D1,Int}}
-    # """
-    # maps a face on its index in enumeration
-    # """
-    # dict::Dict{SVector{D1,Int},Int}
+    faces::Vector{SimplexGraph{D1}}
 end
 
-# function Mesh(vertices, faces)
-#     # dict = Dict((f,i) for (i,f) in enumerate(faces))
-#     # dict = Dict{Int,Int}()
-#     T = eltype(eltype(vertices))
-#     U = length(eltype(vertices))
-#     D1 = length(eltype(faces))
-#     Mesh{U,D1,T}(vertices, faces)
-# end
 
-function indices(m::Mesh{U,D1}, cell) where {U,D1}
-    return m.faces[cell]
+
+function indices(m::Mesh{U,D1}, i) where {U,D1}
+    return m.faces[i].indices
 end
 
 vertexarray(m::Mesh) = [ v[i] for v in m.vertices, i in 1:universedimension(m) ]
@@ -45,7 +39,9 @@ cellarray(m::Mesh) = [ k[i] for k in m.faces, i in 1:dimension(m)+1 ]
 Returns an empty mesh with `coordtype` equal to `type`, of dimension `mdim`
 and embedded in a universe of dimension `udim`
 """
-mesh(T, mdim, udim=mdim+1) = Mesh(Pt{udim,T}[], SVector{mdim+1,Int}[])
+function mesh(T, mdim, udim=mdim+1)
+    Mesh(Pt{udim,T}[], CompScienceMeshes.SimplexGraph{mdim+1}[])
+end
 
 
 
@@ -107,7 +103,7 @@ universedimension(m::AbstractMesh) = length(vertextype(m))
 
 Returns an indexable iterable to the vertices of the mesh
 """
-vertices(m::Mesh) = m.vertices
+function vertices(m::Mesh) m.vertices end
 
 
 #     vertices(mesh, i)
@@ -117,10 +113,10 @@ vertices(m::Mesh) = m.vertices
 # form, only statically sized arrays are allowed to discourage
 # memory allocation. The returned vector in that case will also
 # be statically typed and of the same size as `I`.
-vertices(m::Mesh, i::Number) = m.vertices[i]
+function vertices(m::Mesh, i::Number) m.vertices[i] end
 cellvertices(m::AbstractMesh, cell) = vertices(m, indices(m, cell))
-@generated function vertices(m::Mesh, I::SVector)
-    N = length(I)
+@generated function vertices(m::Mesh, I::NTuple)
+    N = fieldcount(I)
     xp = :(())
     for i in 1:N
         push!(xp.args, :(m.vertices[I[$i]]))
@@ -195,7 +191,7 @@ end
 
 Change the orientation of a cell by interchanging the first to indices.
 """
-@generated function flip(cell)
+@generated function flip(cell::SVector)
     # generate `T(cell[2],cell[1],cell[3],...)`, with `T = typeof(cell)`
     N = length(cell)
     if N <= 1
@@ -208,6 +204,11 @@ Change the orientation of a cell by interchanging the first to indices.
     xp
 end
 
+function flip(cell::SimplexGraph)
+    t = flip(cell.indices)
+    return SimplexGraph(t)
+end
+
 """
     flipmesh!(mesh)
 
@@ -305,9 +306,9 @@ function fliporientation(m::Mesh)
     fliporientation!(n)
 end
 
-@generated function fliporientation(I::SVector{N,T}) where {N,T}
+@generated function fliporientation(I::SimplexGraph{N}) where {N}
     @assert N >= 2
-    xp = :(SVector{N,T}(I[2],I[1]))
+    xp = :(SimplexGraph(I[2],I[1]))
     for i in 3:N
         push!(xp.args, :(I[$i]))
     end
@@ -352,15 +353,15 @@ function boundary(mesh)
     rows = rowvals(conn)
     vals = nonzeros(conn)
 
-    I = indextype(edges)
+    I = celltype(edges)
     bnd_edges = Vector{I}(undef, length(edges))
     i = 1
     for (e,edge) in enumerate(edges)
         nzr = nzrange(conn,e)
         length(nzr) != 1 && continue
         relop = vals[nzr[1]]
         inds = indices(edges, edge)
-        bnd_edges[i] = (relop > 0) ? inds : flip(inds)
+        bnd_edges[i] = (relop > 0) ? I(inds) : flip(I(inds))
         i += 1
     end
 
@@ -382,7 +383,7 @@ function interior(mesh::Mesh, edges=skeleton(mesh,1))
     @assert size(C) == (numcells(mesh), numcells(edges))
 
     nn = vec(sum(abs.(C), dims=1))
-    T = CompScienceMeshes.indextype(edges)
+    T = CompScienceMeshes.celltype(edges)
     interior_edges = Vector{T}()
     for (i,edge) in pairs(cells(edges))
         nn[i] > 1 && push!(interior_edges, edge)
@@ -531,7 +532,8 @@ function skeleton(mesh, dim::Int; sort=:spacefillingcurve)
 
     # sort the simplices on a SFC
     simplices = cells(sk)
-    ctrs = [sum(cellvertices(sk,c))/(dim+1) for c in sk]
+    # ctrs = [sum(cellvertices(sk,c))/(dim+1) for c in sk]
+    ctrs = [cartesian(center(chart(sk, cell))) for cell in simplices]
     if length(sk) > 0
         simplices = simplices[sort_sfc(ctrs)]
     end
@@ -570,18 +572,17 @@ end
 
 
 function skeleton_fast(mesh, dimtype::Type)
-    C = celltype(mesh)
-    I = indextype(mesh, dimtype)
-    E = celltype(mesh, dimtype)
-    faces = Vector{I}()
-    for cell in mesh
-        idcs = indices(mesh, cell)
-        for face in skeleton(C(idcs), dimtype)
-            push!(faces, face)
+    # C = celltype(mesh)
+    # I = indextype(mesh, dimtype)
+    F = celltype(mesh, dimtype)
+    fcs = Vector{F}()
+    for cell in cells(mesh)
+        for face in faces(cell, dimtype)
+            push!(fcs, face)
         end
     end
-    faces = unique!(x -> E(x), faces)
-    Mesh(vertices(mesh), faces)
+    fcs = unique!(fcs)
+    Mesh(vertices(mesh), fcs)
 end
 
 function skeleton_fast(mesh, dim::Int)
@@ -780,8 +781,8 @@ function cellpairs(mesh, edges; dropjunctionpair=false)
 
     k = 1
     Cells = cells(mesh)
-    for e in edges
-        edge = indices(edges, e)
+    for edge in cells(edges)
+        # edge = indices(edges, e)
 
         # neighborhood of startvertex
         v = edge[1]
@@ -864,7 +865,16 @@ end
 
 Return a chart describing the supplied cell of `mesh`.
 """
-chart(mesh::Mesh, cell) = simplex(vertices(mesh, indices(mesh,cell)))
+function chart(mesh::Mesh, c)
+    cell = cells(mesh)[c]
+    chart(mesh, cell)
+end
+
+function chart(mesh::Mesh, cell::SimplexGraph)
+    # cell = cells(mesh)[c]
+    verts = vertices(mesh, cell.indices.data)
+    simplex(verts)
+end
 
 parent(mesh::AbstractMesh) = nothing
 # parent(mesh::Mesh) = nothing

src/meshes/quadrilateralmesh.jl

@@ -2,10 +2,14 @@
 
 struct QuadMesh{T} <: AbstractMesh{3,4,T}
     vertices::Vector{SVector{3,T}}
-    faces::Vector{SVector{4,Int}}
+    faces::Vector{QuadrilateralGraph}
 end
 
-function indices(m::QuadMesh, cellptr) m.faces[cellptr] end
+function QuadMesh(V::Vector{SVector{3,T}}, F::Vector{SVector{4,Int}}) where {T}
+    return QuadMesh{T}(V, [QuadrilateralGraph(f) for f in F])
+end
+
+function indices(m::QuadMesh, cellptr) m.faces[cellptr].indices end
 function vertextype(m::QuadMesh) eltype(m.vertices) end
 function dimension(m::QuadMesh) 2 end
 
@@ -23,7 +27,7 @@ function celltype(m::QuadMesh, ::Type{Val{2}}) QuadrilateralGraph end
 function celltype(m::QuadMesh) QuadrilateralGraph end
 
 function chart(m::QuadMesh, cellptr)
-    verts = m.vertices[m.faces[cellptr]]
+    verts = m.vertices[m.faces[cellptr].indices]
     Quadrilateral(verts...)
 end