Add Mesh2Dual functionality

src/Metis.jl

@@ -7,7 +7,7 @@ using METIS_jll: libmetis
 
 # Metis C API
 include("metis_h.jl")
-const options = fill(Cint(-1), METIS_NOPTIONS)
+const options = fill(idx_t(-1), METIS_NOPTIONS)
 options[METIS_OPTION_NUMBERING] = 1
 
 # Julia interface
@@ -82,6 +82,65 @@ function graph(G::LightGraphs.AbstractSimpleGraph)
     return Graph(idx_t(N), xadj, adjncy)
 end
 
+"""
+    graph(elements::AbstractMatrix{<:Integer}, ncommon=1)
+
+Construct the 1-based CSR representation of the mesh described by the
+connectivity matrix `elements`.
+
+`elements[i,j]` contains the `i`-th node of the `j`-th element. Thus the size of
+the first dimension `elements` is equal to the number of nodes per element, and
+the size of the second dimension is the number of elements.
+
+The nodes are assumed to be numbered consecutively, starting from 1.
+"""
+function graph(elements::AbstractMatrix{<:Integer}, ncommon=1)
+    # Assume the number of nodes in the mesh is equal to the highest node
+    # number.
+    nnodes = maximum(elements)
+
+    # Number of nodes per element.
+    n_npe = size(elements, 1)
+
+    # Number of elements in the mesh.
+    ne = size(elements, 2)
+
+    # Get the element connectivity data in the form METIS needs.
+    eptr = convert(Array{Metis.idx_t}, collect(0:ne) .* n_npe)
+    eind = convert(Array{Metis.idx_t}, elements[:])
+
+    # With 1-based indexing, need to add 1 to eptr.
+    eptr .+= 1
+
+    # Create some pointers and stuff that will be needed for the METIS library call.
+    ne = Metis.idx_t(ne)
+    nn = Metis.idx_t(nnodes)
+    ncommon = Metis.idx_t(ncommon)
+    num_flag = Metis.options[Metis.METIS_OPTION_NUMBERING]
+    xadj = [Ptr{Metis.idx_t}()]
+    adjncy = [Ptr{Metis.idx_t}()]
+
+    # Do it!
+    Metis.METIS_MeshToDual(ne, nn, eptr, eind, ncommon, num_flag, xadj, adjncy)
+
+    # Create a Julian copy of the xadj array. The length of xadj is one more
+    # than the number of graph verticies (here, the number of mesh elements).
+    xadj_out = copy(unsafe_wrap(Vector{Metis.idx_t}, xadj[1], ne+1))
+
+    # Free the METIS memory.
+    Metis.METIS_Free(xadj[1])
+
+    # Create a Julian copy of the adjncy array.  The length of adjncy is twice
+    # the number of graph edges. Turns out that's the last entry in xadj minus
+    # 1 if we're using one-based indices.
+    adjncy_out = copy(unsafe_wrap(Vector{Metis.idx_t}, adjncy[1], xadj_out[end]-1))
+
+    # Free the METIS memory.
+    Metis.METIS_Free(adjncy[1])
+
+    return Graph(ne, xadj_out, adjncy_out)
+end
+
 """
     perm, iperm = Metis.permutation(G)
 

src/metis_h.jl

@@ -122,14 +122,14 @@ function METIS_PartGraphKway(nvtxs, ncon, xadj, adjncy, vwgt, vsize, adjwgt, npa
     return
 end
 
-# function METIS_MeshToDual(ne, nn, eptr, eind, ncommon, numflag, r_xadj, r_adjncy)
-#     r = ccall((:METIS_MeshToDual, libmetis), Cint,
-#               (Ref{idx_t}, Ref{idx_t}, Ptr{idx_t}, Ptr{idx_t}, Ref{idx_t}, Ref{idx_t},
-#                Ptr{idx_t}, Ptr{idx_t}),
-#               ne, nn, eptr, eind, ncommon, numflag, r_xadj, r_adjncy)
-#     r == METIS_OK || throw(MetisError(r))
-#     return
-# end
+function METIS_MeshToDual(ne, nn, eptr, eind, ncommon, numflag, r_xadj, r_adjncy)
+    r = ccall((:METIS_MeshToDual, libmetis), Cint,
+              (Ref{idx_t}, Ref{idx_t}, Ptr{idx_t}, Ptr{idx_t}, Ref{idx_t}, Ref{idx_t},
+               Ptr{idx_t}, Ptr{idx_t}),
+              ne, nn, eptr, eind, ncommon, numflag, r_xadj, r_adjncy)
+    r == METIS_OK || throw(MetisError(r))
+    return
+end
 
 # function METIS_MeshToNodal(ne, nn, eptr, eind, numflag, r_xadj, r_adjncy)
 #     r = ccall((:METIS_MeshToNodal, libmetis), Cint,
@@ -172,11 +172,11 @@ function METIS_NodeND(nvtxs, xadj, adjncy, vwgt, options, perm, iperm)
     return
 end
 
-# function METIS_Free(ptr)
-#     r = ccall((:METIS_Free, libmetis), Cint, (Ptr{idx_t},), ptr)
-#     r == METIS_OK || throw(MetisError(r))
-#     return
-# end
+function METIS_Free(ptr)
+    r = ccall((:METIS_Free, libmetis), Cint, (Ptr{idx_t},), ptr)
+    r == METIS_OK || throw(MetisError(r))
+    return
+end
 
 # function METIS_SetDefaultOptions(options)
 #     r = ccall((:METIS_SetDefaultOptions, libmetis), Cint, (Ptr{idx_t},), options)

test/runtests.jl

@@ -33,3 +33,23 @@ end
         @test extrema(parts) == (1, 3)
     end
 end
+
+@testset "Metis.MeshToDual" begin
+    # Connectivity matrix of 2D mesh of 4 triangles.
+    cells = [1 2 2 3;
+             2 5 3 6;
+             4 4 5 5]
+    # Tell Metis to consider a cell connected if it shares two vertices.
+    ncommon = 2
+
+    # Get a Graph object. The number of graph verticies is the number of mesh
+    # elements.
+    G = Metis.graph(cells, ncommon)
+
+    # Partition the graph.
+    parts = Metis.partition(G, 2)
+    # First two elements should be in the same partition.
+    @test parts[1] == parts[2]
+    # Last two elements should be in the same partition.
+    @test parts[3] == parts[4]
+end