diff --git a/src/baryref.jl b/src/baryref.jl
index 4f7749e..052aab1 100644
--- a/src/baryref.jl
+++ b/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)
diff --git a/src/fileio/gid.jl b/src/fileio/gid.jl
index 29c4f4e..8b441e4 100644
--- a/src/fileio/gid.jl
+++ b/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)
diff --git a/src/fileio/gmsh.jl b/src/fileio/gmsh.jl
index 055596d..877c46d 100644
--- a/src/fileio/gmsh.jl
+++ b/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
diff --git a/src/fileio/gmsh3d.jl b/src/fileio/gmsh3d.jl
index 3a2778d..2a79fd1 100644
--- a/src/fileio/gmsh3d.jl
+++ b/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
diff --git a/src/fileio/readmesh.jl b/src/fileio/readmesh.jl
index 496f7ee..c1954e6 100644
--- a/src/fileio/readmesh.jl
+++ b/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)])
diff --git a/src/mesh.jl b/src/mesh.jl
index 6d89e2e..f801417 100644
--- a/src/mesh.jl
+++ b/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,7 +353,7 @@ 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)
@@ -360,7 +361,7 @@ function boundary(mesh)
         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
diff --git a/src/meshes/quadrilateralmesh.jl b/src/meshes/quadrilateralmesh.jl
index bd61c21..a94668e 100644
--- a/src/meshes/quadrilateralmesh.jl
+++ b/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
 
diff --git a/src/primitives/linemeshes/linemeshes.jl b/src/primitives/linemeshes/linemeshes.jl
index da90cf7..fa9e97f 100644
--- a/src/primitives/linemeshes/linemeshes.jl
+++ b/src/primitives/linemeshes/linemeshes.jl
@@ -5,7 +5,7 @@
 function meshsegment(L::T, delta::T, udim=2) where T<:Real
 
     PT = SVector{udim, T}
-    CT = SVector{2,Int}
+    CT = SimplexGraph{2}
 
     num_segments = ceil(Int, L/delta)
     actual_delta = L/num_segments
@@ -20,7 +20,7 @@ function meshsegment(L::T, delta::T, udim=2) where T<:Real
 
     faces = Array{CT}(undef,num_segments)
     for i in 1 : length(faces)
-        faces[i] = SVector(i, i+1)
+        faces[i] = SimplexGraph(i, i+1)
     end
 
     Mesh(vertices, faces)
@@ -32,8 +32,9 @@ end
 """
 function meshcircle(radius::T, delta::T, udim=2) where T<:Real
 
-PT = SVector{udim,T}
-CT = SVector{2,Int}
+    PT = SVector{udim,T}
+    # CT = SVector{2,Int}
+    CT = SimplexGraph{2}
 
     circumf = 2 * pi *radius
     num_segments = ceil(Int, circumf / delta)
@@ -49,11 +50,11 @@ CT = SVector{2,Int}
         vertices[i] = PT(a)
     end
 
-    faces = Array{CT}(undef,num_segments)
+    faces = Vector{CT}(undef,num_segments)
     for i in 1 : length(faces)-1
-        faces[i] = SVector{2,Int}(i, i+1)
+        faces[i] = SimplexGraph(i, i+1)
     end
-    faces[end] = SVector{2,Int}(num_segments, 1)
+    faces[end] = SimplexGraph(num_segments, 1)
 
     return Mesh(vertices, faces)
 end
@@ -107,5 +108,6 @@ function meshcurve(curve, delta::T; udim=2, tstart=T(0.0), tend=T(2π), order=1)
         faces[end] = SVector{2,Int}(num_segments, num_segments + 1)
     end
 
+    faces = SimplexGraph.(faces)
     return Mesh(vertices, faces)
 end
\ No newline at end of file
diff --git a/src/primitives/surfacemeshes/mesh_cuboid.jl b/src/primitives/surfacemeshes/mesh_cuboid.jl
index 608984a..611ce9f 100644
--- a/src/primitives/surfacemeshes/mesh_cuboid.jl
+++ b/src/primitives/surfacemeshes/mesh_cuboid.jl
@@ -512,6 +512,7 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F
                 end
             end
         end
+    faces = [SimplexGraph(f) for f in faces]
     Γ = Mesh(nodes, faces)
     return flipmesh!(Γ)
 end
diff --git a/src/primitives/surfacemeshes/mesh_icosphere.jl b/src/primitives/surfacemeshes/mesh_icosphere.jl
index e0b6e66..645e349 100644
--- a/src/primitives/surfacemeshes/mesh_icosphere.jl
+++ b/src/primitives/surfacemeshes/mesh_icosphere.jl
@@ -50,6 +50,7 @@ function meshicosphere(nu::Int = 1, r::F = 1.0, nr_verts = nothing) where F
             vertices[i] = r.*vertices[i]./sqrt.(sum(x->x.^2, vertices[i]))
         end
     end 
+    faces = [SimplexGraph(f) for f in faces]
     return Mesh(vertices, faces)
 end
 
diff --git a/src/primitives/surfacemeshes/mesh_rectangle.jl b/src/primitives/surfacemeshes/mesh_rectangle.jl
index 17d838f..1d117fc 100644
--- a/src/primitives/surfacemeshes/mesh_rectangle.jl
+++ b/src/primitives/surfacemeshes/mesh_rectangle.jl
@@ -88,7 +88,7 @@ function mesh_rectangle_impl(a::F, b::F, h::F, udim) where F
     m = Int(round(b/h))  # number of elements along b
     
     nodes = zeros(SVector{udim, F}, (m + 1)*(n + 1))
-    faces = Vector{SVector{3, Int64}}(undef, 2*m*n)
+    faces = Vector{SimplexGraph{3}}(undef, 2*m*n)
     
     for ix in 0 : n - 1
         for iy in 1 : m
@@ -99,14 +99,14 @@ function mesh_rectangle_impl(a::F, b::F, h::F, udim) where F
             end
             
             nodes[(ix)*(m + 1) + iy] = node
-            face = SVector(
+            face = SimplexGraph(
                 (ix)*(m + 1) + (iy),
                 (ix)*(m + 1) + (iy + 1),
                 (ix + 1)*(m + 1) + (iy)
             )
             faces[(ix)*2*m + (2*iy - 1)] = face
             
-            face = SVector(
+            face = SimplexGraph(
                 (ix)*(m + 1) + (iy + 1), 
                 (ix + 1)*(m + 1) + (iy + 1), 
                 (ix + 1)*(m + 1) + (iy)
@@ -244,6 +244,7 @@ function meshrectangle_unit_tri_2graded(h::T, udim) where {T}
         end
     end
 
+    faces = [SimplexGraph(f) for f in faces]
     return Mesh(vertices, faces)
 end
 
diff --git a/src/primitives/volumemeshes/mesh_cuboid.jl b/src/primitives/volumemeshes/mesh_cuboid.jl
index 8066283..26e7ab2 100644
--- a/src/primitives/volumemeshes/mesh_cuboid.jl
+++ b/src/primitives/volumemeshes/mesh_cuboid.jl
@@ -170,6 +170,7 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F
         end
     end
 
+    faces = [SimplexGraph(f) for f in faces]
     return Mesh(vertices, faces)
 end
 
diff --git a/src/subdMesh.jl b/src/subdMesh.jl
index d2d3a62..4daf6f2 100644
--- a/src/subdMesh.jl
+++ b/src/subdMesh.jl
@@ -71,7 +71,7 @@ function GSubdMesh(mesh::Mesh)
         Svertices[V] = SV
     end
     for (E, Edge) in enumerate(cells(edges))
-        EndVerticesIndices = Edge
+        EndVerticesIndices = Edge.indices
         FacesIndex = []
         RelatedVerticesIndices = []
         orientation = []
@@ -82,7 +82,7 @@ function GSubdMesh(mesh::Mesh)
                 append!(orientation,connectMat[F,E])
                     ## Loop over vertices in the faces find out related vertices
                 for iv = 1 : length(Face)
-                    if Face[iv] != Edge[1] && Face[iv] != Edge[2]
+                    if Face[iv] != Edge.indices[1] && Face[iv] != Edge.indices[2]
                             append!(RelatedVerticesIndices, Face[iv])
                     end
                 end
diff --git a/src/submesh.jl b/src/submesh.jl
index e00999a..9b2e1b5 100644
--- a/src/submesh.jl
+++ b/src/submesh.jl
@@ -45,7 +45,7 @@ mutable struct SubMesh{U,D1,T} <: AbstractMesh{U,D1,T}
     supermesh::AbstractMesh{U,D1,T}
     sub2sup::Vector{Int}
     sup2sub::Vector{Int}
-    cells::Vector{SVector{D1,Int}}
+    # cells::Vector{SVector{D1,Int}}
 end
 
 # function celltype(m::SubMesh{U,D1}) where {U,D1} SimplexGraph{D1} end
@@ -64,7 +64,8 @@ function SubMesh(supermesh, sub2sup)
     # Cells = cells(supermesh)[sub2sup]
     Cells = [indices(supermesh, i) for i in supermesh]
     Cells = Cells[sub2sup]
-    SubMesh(supermesh, sub2sup, sup2sub, Cells)
+    # SubMesh(supermesh, sub2sup, sup2sub, Cells)
+    SubMesh(supermesh, sub2sup, sup2sub)
 end
 
 function vertextype(m::SubMesh) vertextype(m.supermesh) end
@@ -84,9 +85,13 @@ vertices(m::SubMesh, cell) = vertices(m.supermesh, cell)
 numvertices(m::SubMesh) = numvertices(m.supermesh)
 
 indices(m::SubMesh, p) = indices(m.supermesh, m.sub2sup[p])
-cells(m::SubMesh) = m.cells #m.supermesh.faces[m.sub2sup]
 numcells(m::SubMesh) = length(m.sub2sup)
 
+function cells(m::SubMesh)
+    return cells(parent(m))[m.sub2sup]
+    # m.supermesh.faces[m.sub2sup]
+end
+
 issubmesh(sub, sup) = (sub == sup)
 issubmesh(sub::SubMesh, sup) = (sub.supermesh == sup)
 issubmesh(sub::SubMesh, sup::SubMesh) = (sub == sup)
diff --git a/src/tpredicates.jl b/src/tpredicates.jl
index ef68351..4213efe 100644
--- a/src/tpredicates.jl
+++ b/src/tpredicates.jl
@@ -1,9 +1,12 @@
 import Base.in
 function in(mesh::CompScienceMeshes.AbstractMesh)
-    cells_mesh = sort.(cells(mesh))
+    # cells_mesh = sort.(cells(mesh))
+    cells_mesh = cells(mesh)
     function f(m,p)
-        cell = indices(m,p)
-        sort(cell) in cells_mesh
+        # cell = indices(m,p)
+        cell = cells(m)[p]
+        # sort(cell) in cells_mesh
+        cell in cells_mesh
     end
     return f
 end
diff --git a/src/utils/localgraphs.jl b/src/utils/localgraphs.jl
index 6a42a0c..5c4ee8f 100644
--- a/src/utils/localgraphs.jl
+++ b/src/utils/localgraphs.jl
@@ -1,6 +1,6 @@
 abstract type AbstractGraph end
 function Base.getindex(g::AbstractGraph, i::Int) g.indices[i] end
-function Base.hash(g::AbstractGraph, a::UInt) hash(sort(g.indices), a) end
+function Base.hash(g::AbstractGraph, a::UInt) hash(sort(SVector(g.indices)).data, a) end
 function Base.length(g::AbstractGraph) length(g.indices) end
 
 struct QuadrilateralGraph <: AbstractGraph
@@ -14,53 +14,88 @@ struct SimplexGraph{N} <: AbstractGraph
     indices::SVector{N,Int}
 end
 
+# function SimplexGraph{1}(i::Int)
+#     return SimplexGraph{1}(SVector(i))
+# end
+
+function SimplexGraph(indices::Vararg{Int,N}) where {N}
+    return SimplexGraph{N}(SVector{N}(indices))
+end
+
+function SimplexGraph(i::Int)
+    return SimplexGraph{1}(SVector(i))
+end
+
 function dimension(h::SimplexGraph{N}) where {N} N-1 end
 function dimtype(g::SimplexGraph{N}) where {N} Val{N-1} end
 
 function Base.isequal(g::G, h::G) where {G <: AbstractGraph}
-    p = indexin(g.indices, h.indices)
+    p = indexin(SVector(g.indices), SVector(h.indices))
     any(p .== nothing) && return false
     return true
 end
 
-function skeleton(g::SimplexGraph{N}, ::Type{Val{M}}) where {N,M}
-    collect(idcs for idcs in combinations(g.indices, M+1))
+# function skeleton(g::SimplexGraph{N}, ::Type{Val{M}}) where {N,M}
+#     collect(idcs for idcs in combinations(g.indices, M+1))
+# end
+
+function faces(g::SimplexGraph{N}, ::Type{Val{M}}) where {N,M}
+    # F = SimplexGraph{M+1}
+    F = SimplexGraph
+    collect(F(idcs...) for idcs in combinations(g.indices, M+1))
 end
 
 # function skeleton(g::SimplexGraph{N}, ::Type{Val{N-1}}) where {N}
 #     SVector(g.indices)
 # end
 
-function skeleton(g::SimplexGraph{3}, ::Type{Val{1}})
-    return SVector(
-        (g[2],g[3]),
-        (g[3],g[1]),
-        (g[1],g[2]),
+# function skeleton(g::SimplexGraph{3}, ::Type{Val{1}})
+#     return SVector(
+#         (g[2],g[3]),
+#         (g[3],g[1]),
+#         (g[1],g[2]),
+#     )
+# end
+
+function faces(g::SimplexGraph{3}, ::Type{Val{1}})
+    return (
+        SimplexGraph(g[2],g[3]),
+        SimplexGraph(g[3],g[1]),
+        SimplexGraph(g[1],g[2]),
     )
 end
 
 
 @testitem "skeleton SimplexGraph" begin
     g = CompScienceMeshes.SimplexGraph{4}((1,2,3,4))
-    edges = skeleton(g,Val{1})
+    edges = faces(g,Val{1})
     @test length(edges) == 6
     # @test eltype(edges) == CompScienceMeshes.SVector{2,Int}
 end
 
-function skeleton(g::QuadrilateralGraph, ::Type{Val{1}})
+# function skeleton(g::QuadrilateralGraph, ::Type{Val{1}})
+#     return SVector(
+#         (g[1],g[2]),
+#         (g[2],g[3]),
+#         (g[3],g[4]),
+#         (g[4],g[1]),)
+# end
+
+function faces(g::QuadrilateralGraph, ::Type{Val{1}})
     return SVector(
-        (g[1],g[2]),
-        (g[2],g[3]),
-        (g[3],g[4]),
-        (g[4],g[1]),)
+        SimplexGraph(g[1],g[2]),
+        SimplexGraph(g[2],g[3]),
+        SimplexGraph(g[3],g[4]),
+        SimplexGraph(g[4],g[1]),)
 end
 
 function relorientation(g::G, h::AbstractGraph) where {G <: SimplexGraph}
     dimension(g) > dimension(h) && return relorientation(h, g)
-    faces = collect(skeleton(h, dimtype(g)))
-    i = something(findfirst(face -> isequal(G(face), g), faces), 0)
+    fcs = collect(faces(h, dimtype(g)))
+    i = something(findfirst(face -> isequal(face, g), fcs), 0)
     i == 0 && return 0
-    p = Permutations.Permutation(indexin(g.indices, SVector(faces[i])))
+    # p = Permutations.Permutation(indexin(g.indices, SVector(fcs[i])))
+    p = Permutations.Permutation(indexin(g.indices, fcs[i].indices))
     σ = Permutations.sign(p)
     return σ * i
 end
diff --git a/src/weld.jl b/src/weld.jl
index fd319ed..84f1886 100644
--- a/src/weld.jl
+++ b/src/weld.jl
@@ -76,8 +76,10 @@ function weld(Γ₁, Γ₂; boundary=false, glueop=unique)
 
     V = [V1; V2]
 
-    F1 = [indices(Γ₁,c) for c in Γ₁]
-    F2 = [indices(Γ₂,c) for c in Γ₂]
+    # F1 = [indices(Γ₁,c) for c in Γ₁]
+    # F2 = [indices(Γ₂,c) for c in Γ₂]
+    F1 = collect(cells(Γ₁))
+    F2 = collect(cells(Γ₂))
     F = [F1; F2]
 
     nc1 = numcells(Γ₁)
@@ -188,8 +190,8 @@ function weld(G1::SComplex2D, G2::SComplex2D; seam)
         vcat(G1.faces, Faces2))
 end
 
-@generated function map_ids(c::SVector{N,T}, idmap) where {N,T}
-    xp = :(SVector{$N,$T}())
+@generated function map_ids(c::SimplexGraph{N}, idmap) where {N,T}
+    xp = :(SimplexGraph())
     for i in 1:N
         push!(xp.args, :(idmap[c[$i]]))
     end
diff --git a/test/test_flipped.jl b/test/test_flipped.jl
index e54db31..c771ab6 100644
--- a/test/test_flipped.jl
+++ b/test/test_flipped.jl
@@ -8,9 +8,10 @@ for U in [Float32, Float64]
         point(U,1,1,0),
         point(U,0,1,0)]
 
+    C = CompScienceMeshes.SimplexGraph{3}
     global F = [
-        index(1,2,3),
-        index(1,3,4)]
+        C(index(1,2,3)),
+        C(index(1,3,4))]
 
     global m = Mesh(V,F)
 
@@ -19,11 +20,11 @@ for U in [Float32, Float64]
     global g = -f
     @test g === m
 
-    global T = [index(2,1,3), index(3,1,4)]
+    global T = [C(index(2,1,3)), C(index(3,1,4))]
     # for (i,_c) in enumerate(cells(f))
-    for i in f
-        inds = CompScienceMeshes.indices(f,i)
-        @test inds == T[i]
+    for (i,cell) in pairs(cells(f))
+        # inds = CompScienceMeshes.indices(f,i)
+        @test cell == T[i]
     end
 end
 # @test f.faces[1] == index(2,1,3)
diff --git a/test/test_geometry.jl b/test/test_geometry.jl
index d48ca48..a941246 100644
--- a/test/test_geometry.jl
+++ b/test/test_geometry.jl
@@ -56,10 +56,14 @@ for T in [Float32, Float64]
     #for i in 1:numcells(m)
     for cl in cells(m)
         #p = chart(m, cells(m,i))
-        v = vertices(m, cl)
-        p = simplex(v)
-        _c = (p.vertices[1] + p.vertices[2]) / 2
-        @test dot(_c, p.normals[1]) > 0
+        # v = vertices(m, cl)
+        # p = simplex(v)
+        # _c = (p.vertices[1] + p.vertices[2]) / 2
+        p = center(chart(m, cl))
+        _c = cartesian(p)
+        # _c = cartesian(center(chart(m, cl)))
+        # @test dot(_c, p.normals[1]) > 0
+        @test dot(cartesian(p), normal(p)) > 0
     end
 end
 
@@ -73,8 +77,8 @@ for T in [Float32, Float64]
     p3 = point(T,-1,0,0)
     p4 = point(T,0,-1,0)
 
-    i1 = index(1,2,3)
-    i2 = index(1,4,2)
+    i1 = CompScienceMeshes.SimplexGraph(1,2,3)
+    i2 = CompScienceMeshes.SimplexGraph(1,4,2)
 
     local m = Mesh([p1,p2,p3,p4], [i1,i2])
     e = skeleton(m,1)
@@ -133,7 +137,7 @@ for T in [Float32, Float64]
                 if val == 1
                     @test e1_cells[row] == E1_cells[col]
                 else
-                    @test e1_cells[row] == reverse(E1_cells[col])
+                    @test e1_cells[row].indices == reverse(E1_cells[col].indices)
                 end
             end
         end
diff --git a/test/test_gmsh_multiple_dims.jl b/test/test_gmsh_multiple_dims.jl
index a6538c7..7735dc5 100644
--- a/test/test_gmsh_multiple_dims.jl
+++ b/test/test_gmsh_multiple_dims.jl
@@ -12,7 +12,7 @@ sea_vertices    = skeleton(sea, 0)
 
 # interior_vertices = submesh(v -> !(v in border_vertices || v in coast_vertices), sea_vertices)
 interior_vertices = submesh(sea_vertices) do mesh, i
-    v = CompScienceMeshes.indices(mesh, i)
+    v = CompScienceMeshes.cells(mesh)[i]
     v in cells(border_vertices) && return false
     v in cells(coast_vertices) && return false
     return true
diff --git a/test/test_jctweld.jl b/test/test_jctweld.jl
index 9cca060..4087f1b 100644
--- a/test/test_jctweld.jl
+++ b/test/test_jctweld.jl
@@ -8,7 +8,7 @@ vrts = [
 ]
 
 idcs = [
-    index(1,2,3)
+    CompScienceMeshes.SimplexGraph(1,2,3)
 ]
 
 axis = 0.5π * [1,0,0]
@@ -23,9 +23,9 @@ G = weld(G1, G2, G3)
 @test numvertices(G) == 5
 
 edges = skeleton(G,1)
-pairs = cellpairs(G, edges, dropjunctionpair=true)
+prs = cellpairs(G, edges, dropjunctionpair=true)
 
-@test size(pairs) == (2,8)
+@test size(prs) == (2,8)
 
 mesh = meshrectangle(1.0, 1.0, 0.5)
 edges = skeleton(mesh,1)
diff --git a/test/test_mapper.jl b/test/test_mapper.jl
index ac58e47..bf16d08 100644
--- a/test/test_mapper.jl
+++ b/test/test_mapper.jl
@@ -1,7 +1,10 @@
 using CompScienceMeshes
+
+using SparseArrays
 using Test
 
 #include(Pkg.dir("CompScienceMeshes","examples","matlab_patches.jl"))
+# T = Float64
 for T in [Float32, Float64]
     sphere = readmesh(joinpath(dirname(@__FILE__),"assets","sphere2.in"),T=T)
     index_in_sphere = mapper(sphere)
@@ -15,6 +18,7 @@ for T in [Float32, Float64]
     # build the restriction operator
     R = spzeros(Int, numcells(southern_hemisphere), numcells(sphere))
     for (i,_c) in enumerate(cells(southern_hemisphere))
+        # @show _c
         j = index_in_sphere[_c]
         R[i,j] = 1
     end
diff --git a/test/test_mesh.jl b/test/test_mesh.jl
index 545910d..749375d 100644
--- a/test/test_mesh.jl
+++ b/test/test_mesh.jl
@@ -69,14 +69,14 @@ T = Float64
     ]
 
     F = [
-        index(1,2,3)
+        CompScienceMeshes.SimplexGraph{3}(index(1,2,3))
     ]
 
     m = Mesh(V,F)
     n = fliporientation(m)
 
-    @test m.faces[1] == index(1,2,3)
-    @test n.faces[1] == index(2,1,3)
+    @test m.faces[1].indices == index(1,2,3)
+    @test n.faces[1].indices == index(2,1,3)
 
     # test the construction of transposed connectivity matricees
     Σᵀ = connectivity(edges, faces, identity)
diff --git a/test/test_overlap.jl b/test/test_overlap.jl
index 6f784ae..437fca0 100644
--- a/test/test_overlap.jl
+++ b/test/test_overlap.jl
@@ -44,10 +44,11 @@ vt = skeleton(l1,0)
 bd = boundary(l1)
 
 overlaps = overlap_gpredicate(bd)
-pred1 = c -> overlaps(simplex(vertices(vt,c)))
-@test pred1(SVector(1))
-@test !pred1(SVector(2))
-@test pred1(SVector(3))
+# pred1 = c -> overlaps(simplex(vertices(vt,c)))
+pred1 = c -> overlaps(chart(vt,c))
+@test pred1(CompScienceMeshes.SimplexGraph(1))
+@test !pred1(CompScienceMeshes.SimplexGraph(2))
+@test pred1(CompScienceMeshes.SimplexGraph(3))
 
 
 # test a case where the segments are:
diff --git a/test/test_patches.jl b/test/test_patches.jl
index 87ede45..403c156 100644
--- a/test/test_patches.jl
+++ b/test/test_patches.jl
@@ -3,9 +3,10 @@ using CompScienceMeshes
 
 for T in [Float32, Float64]
     local mesh = meshrectangle(T(1.0), T(1.0), T(1.0))
-    local faces = skeleton(mesh, 2)
-    local verts = vertices(mesh, faces.faces[1])
-    p = simplex(verts)
+    # local faces = skeleton(mesh, 2)
+    # local verts = vertices(mesh, faces.faces[1])
+    p = chart(mesh, cells(mesh)[1])
+    # p = simplex(verts)
 
     @test p.vertices == [
         point(T, 0.0, 0.0, 0.0),