rigetti/DiscreteExteriorCalculus.jl

Discrete differential geometry on simplicial complexes

DiscreteExteriorCalculus.jl

DiscreteExteriorCalculus.jl is a package implementing Discrete Exterior Calculus. Data structures for cell complexes, primal, and dual meshes are provided along with implementations of the exterior derivative, hodge star, codifferential, and Laplace-de Rham operators.

Installation

Clone the repository from GitHub and install Julia 1.1. No build is required beyond the default Julia compilation.

Example usage: modes of the Laplace-de Rham operator on a rectangle with Dirichlet boundary conditions

See test/test_laplacian_rectangle.jl for a more complete version of this example. Also see DiscretePDEs.jl for more examples.

Import packages.

using DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus
using LinearAlgebra: eigen
using AdmittanceModels: sparse_nullbasis
using Base.Iterators: product
using PlotlyJS: plot, heatmap

Create a rectangular grid that is r1×r2, then subdivide each rectangle into two triangles. Collect these into a cell complex.

r1, r2 = .5, .4
num = 40
points, tcomp = DEC.triangulated_lattice([r1, 0], [0, r2], num, num)

Orient the cell complex positively and compute the dual mesh.

orient!(tcomp.complex)
m = Metric(2)
mesh = Mesh(tcomp, circumcenter(m))

Compute the Laplace-de Rham operator for 0-forms and a sparse nullbasis for the constraint that the 0-form goes to 0 on the boundary.

laplacian = differential_operator(m, mesh, "Δ", 1, true)
comp = tcomp.complex
_, exterior = boundary_components_connected(comp)
constraint = zero_constraint(comp, exterior.cells[1], 1)
nullbasis = sparse_nullbasis(constraint)

Compute the eigenvalues and eigenvectors of the restricted Laplace-de Rham operator.

vals, vects = eigen(collect(transpose(nullbasis) * laplacian * nullbasis))
inds = sortperm(vals)
vals, vects = vals[inds], vects[:, inds]

Lift the eigenvectors back to the original space and reshape for plotting.

comp_points = [c.points[1] for c in comp.cells[1]]
ordering = [findfirst(isequal(p), comp_points) for p in points]
vs = [collect(transpose(reshape((nullbasis * vects[:,i])[ordering],
    num+1, num+1))) for i in 1:size(vects, 2)]
for i in 1:length(vs)
    j = argmax(abs.(vs[i]))
    vs[i] /= vs[i][j]
end

Plot the first eigenvector.

plot(heatmap(z=vs[1]))

Plot the second eigenvector.

plot(heatmap(z=vs[2]))

Plot the third eigenvector.

plot(heatmap(z=vs[3]))

Plot the fourth eigenvector.

plot(heatmap(z=vs[4]))

And so on.